4. Model Equations

In this chapter, we discuss the model equations used by ParFlow for its fully and variably saturated flow, overland flow, and multiphase flow and transport models. First, section Steady-State, Saturated Groundwater Flow describes steady-state, groundwater flow (specified by solver IMPES). Next, section Richards’ Equation describes the Richards’ equation model (specified by solver RICHARDS) for variably saturated flow as implemented in ParFlow. Section Terrain Following Grid describes the terrain following grid formulation. Next, the overland flow equations are presented in section Overland Flow. In section Multi-Phase Flow Equations we describe the multi-phase flow equations (specified by solver IMPES), and in section Transport Equations we describe the transport equations. Finally, section Notation and Units presents some notation and units and section Water Balance presents some basic water balance equations.

4.1. Steady-State, Saturated Groundwater Flow

Many groundwater problems are solved assuming steady-state, fully-saturated groundwater flow. This follows the form often written as:

(4.1)\[\begin{aligned} \nabla \cdot\textbf{q} = Q(x) \end{aligned}\]

where \(Q\) is the spatially-variable source-sink term (to represent wells, etc) and \(\textbf{q}\) is the Darcy flux \([L^{2}T^{-1}]\) which is commonly written as:

(4.2)\[\begin{aligned} \textbf{q}=- \textbf{K} \nabla H \end{aligned}\]

where \(\textbf{K}\) is the saturated, hydraulic conductivity tensor \([LT^{-1}]\) and \(H\) \([L]\) is the head-potential. Inspection of (4.18) and (4.19) show that these equations agree with the above formulation for a single-phase (\(i=1\)), fully-saturated (\(S_i=S=1\)), problem where the mobility, \({\lambda}_i\), is set to the saturated hydraulic conductivity, \(\textbf{K}\), below. This is accomplished by setting the relative permeability and viscosity terms to unity in (4.20) as well as the gravity and density terms in (4.19). This is shown in the example in Annotated Input Scripts, but please note that the resulting solution is in pressure-head, \(h\), not head potential, \(H\), and will still contain a hydrostatic pressure gradient in the \(z\) direction.

4.2. Richards’ Equation

The form of Richards’ equation implemented in ParFlow is given as,

(4.3)\[\begin{aligned} S(p)S_s\frac{\partial p}{\partial t} - \frac{\partial (S(p)\rho(p)\phi)}{\partial t} - \nabla \cdot(\textbf{K}(p)\rho(p)(\nabla p - \rho(p) {\vec g})) = Q, \; {\rm in} \; \Omega, \end{aligned}\]

where \(\Omega\) is the flow domain, \(p\) is the pressure-head of water \([L]\), \(S\) is the water saturation, \(S_s\) is the specific storage coefficient \([L^{-1}]\), \(\phi\) is the porosity of the medium, \(\textbf{K}(p)\) is the hydraulic conductivity tensor \([LT^{-1}]\), and \(Q\) is the water source/sink term \([L^{3}T^{-1}]\) (includes wells and surface fluxes). The hydraulic conductivity can be written as,

(4.4)\[\begin{aligned} K(p) = \frac{{\bar k}k_r(p)}{\mu} \end{aligned}\]

Boundary conditions can be stated as,

(4.5)\[\begin{split}\begin{align} p & = & p_D, \; {\rm on} \; \Gamma^D, \\ -K(p)\nabla p \cdot {\bf n} & = & g_N, \; {\rm on} \; \Gamma^N, \end{align}\end{split}\]

where \(\Gamma^D \cup \Gamma^N = \partial \Omega\), \(\Gamma^D \neq \emptyset\), and \({\bf n}\) is an outward pointing, unit, normal vector to \(\Omega\). This is the mixed form of Richards’ equation. Note here that due to the constant (or passive) air phase pressure assumption, Richards’ equation ignores the air phase except through its effects on the hydraulic conductivity, \(K\). An initial condition,

(4.6)\[\begin{aligned} p = p^0(x), \; t = 0, \end{aligned}\]

completes the specification of the problem.

4.3. Terrain Following Grid

The terrain following grid formulation transforms the ParFlow grid to conform to topography [Max13]. This alters the form of Darcy’s law to include a topographic slope component:

(4.7)\[\begin{aligned} q_x=\textbf{K}(p)\rho(p)(\frac{\partial p}{\partial x}\cos \theta_x + \sin \theta_x) \end{aligned}\]

where \(\theta_x = \arctan(S_0,x)\) and \(\theta_y = \arctan(S_0,y)\) which are assumed to be the same as the TopoSlope keys assigned for overland flow, described below. The terrain following grid formulation can be very useful for coupled surface-subsurface flow problems where groundwater flow follows the topography. As cells are distributed near the ground surface and can be combined with the variable \(\delta Z\) capability, the number of cells in the problem can be reduced dramatically over the orthogonal formulation. For complete details on this formulation, the stencil used and the function evaluation developed, please see Maxwell [Max13]. NOTE: in the original formulation, \(\theta_x\) and \(\theta_y\) for a cell face is calculated as the average of the two adjacent cell slopes (i.e. assuming a cell centered slope calculation). The TerrainFollowingGrid.SlopeUpwindFormulation key provide options to use the slope of a grid cell directly (i.e. assuming face centered slope calculations) and removing the sine term from (4.7). The Upwind and UpwindSine options for this key will provide consistent results with OverlandKinematic and OverlandDiffusive boundary conditions while the Original option is consistent with the standard OverlandFlow boundary condition.

4.4. Flow Barriers

The the flow barrier multipliers allow for the reduction in flow across a cell face. This slightly alters Darcy’s law to include a flow reduction in each direction, show here in x:

(4.8)\[\begin{aligned} q_x=FB_x\textbf{K}(p)\rho(p)(\frac{\partial p}{\partial x}\cos \theta_x + \sin \theta_x) \end{aligned}\]

where \(FB_x\), \(FB_y\) and \(FB_z\) are a dimensionless multipliers specified by the FBx, FBy and FBz keys. This creates behavior equivalent to the Hydraulic Flow Barrier (HFB) or ITFC (flow and transport parameters at interfaces) conditions in other models.

4.5. Overland Flow

As detailed in Kollet and Maxwell [KM06], ParFlow may simulate fully-coupled surface and subsurface flow via an overland flow boundary condition. While complete details of this approach are given in that paper, a brief summary of the equations solved are presented here. Shallow overland flow is now represented in ParFlow by the kinematic wave equation. In two spatial dimensions, the continuity equation can be written as:

(4.9)\[\begin{aligned} \frac{\partial \psi_s}{\partial t} = \nabla \cdot({\vec v}\psi_s) + q_r(x) \end{aligned}\]

where \({\vec v}\) is the depth averaged velocity vector \([LT^{-1}]\); \(\psi_s\) is the surface ponding depth \([L]\) and \(q_r(x)\) is the a general source/sink (e.g. rainfall) rate \([LT^{-1}]\). If diffusion terms are neglected the momentum equation can be written as:

(4.10)\[\begin{aligned} S_{f,i} = S_{o,i} \end{aligned}\]

which is commonly referred to as the kinematic wave approximation. In Equation (4.10) \(S_{o,i}\) is the bed slope (gravity forcing term) \([-]\), which is equal to the friction slope \(S_{f,i}\) \([L]\); \(i\) stands for the \(x\)- and \(y\)-direction. Manning’s equation is used to establish a flow depth-discharge relationship:

(4.11)\[\begin{aligned} v_x=- \frac{\sqrt{S_{f,x}}}{n}\psi_{s}^{2/3} \end{aligned}\]

and

(4.12)\[\begin{aligned} v_y=- \frac{\sqrt{S_{f,y}}}{n}\psi_{s}^{2/3} \end{aligned}\]

where \(n\) \([TL^{-1/3}]\) is the Manning’s coefficient. Though complete details of the coupled approach are given in Kollet and Maxwell [KM06], brief details of the approach are presented here. The coupled approach takes Equation eq:kinematic and adds a flux for subsurface exchanges, \(q_e(x)\).

(4.13)\[\begin{aligned} \frac{\partial \psi_s}{\partial t} = \nabla \cdot({\vec v}\psi_s) + q_r(x) + q_e(x) \end{aligned}\]

We then assign a continuity of pressure at the top cell of the boundary between the surface and subsurface systems by setting pressure-head, \(p\) in (4.3) equal to the vertically-averaged surface pressure, \(\psi_s\) as follows:

(4.14)\[\begin{aligned} p = \psi_s = \psi \end{aligned}\]

If we substitute this relationship back into Equation (4.13) as follows:

(4.15)\[\begin{aligned} \frac{\partial \parallel\psi,0\parallel}{\partial t} = \nabla \cdot({\vec v}\parallel\psi,0\parallel) + q_r(x) + q_e(x) \end{aligned}\]

Where the \(\parallel\psi,0\parallel\) operator chooses the greater of the two quantities, \(\psi\) and \(0\). We may now solve this term for the flux \(q_e(x)\) which we may set equal to flux boundary condition shown in Equation (4.5). This yields the following equation, which is referred to as the overland flow boundary condition [KM06]:

(4.16)\[\begin{aligned} -K(\psi)\nabla \psi \cdot {\bf n} = \frac{\partial \parallel\psi,0\parallel}{\partial t} - \nabla \cdot({\vec v}\parallel\psi,0\parallel) - q_r(x) \end{aligned}\]

This results a version of the kinematic wave equation that is only active when the pressure at the top cell of the subsurface domain has a ponded depth and is thus greater than zero. This method solves both systems, where active in the domain, over common grids in a fully-integrated, fully-mass conservative manner.

The depth-discharge relationship can also be written as

(4.17)\[\begin{aligned} v_x=- \frac{S_{f,x}}{n\sqrt{\overline{S_{f}}}}\psi_{s}^{2/3} \end{aligned}\]

where \(\overline{S_{f}}\) is the magnitude of the friction slope. This formulation for overland flow is used in the OverlandKinematic and OverlandDiffusive boundary conditions. In OverlandKinematic case the friction slope equals the bed slope following Equation (4.10). For the OverlandDiffusive case the friction slope also includes the pressure gradient. The solution for both of these options is formulated to do the upwinding internally and assumes that the user provides face centered bedslopes (\(S_{o,i}\)). This is different from the original formulation which assumes the user provides grid cenered bedslopes.

4.6. Multi-Phase Flow Equations

The flow equations are a set of mass balance and momentum balance (Darcy’s Law) equations, given respectively by,

(4.18)\[\frac{\partial}{\partial t} ( \phi S_i) ~+~ \nabla\cdot {\vec V}_i ~-~ Q_i~=~ 0 ,\]

(4.19)\[{\vec V}_i~+~ {\lambda}_i\cdot ( \nabla p_i~-~ \rho_i{\vec g}) ~=~ 0 ,\]

for \(i = 0, \ldots , \nu- 1\) \((\nu\in \{1,2,3\})\), where

(4.20)\[\begin{split}\begin{aligned} {\lambda}_i& = & \frac{{\bar k}k_{ri}}{\mu_i} , \\ {\vec g}& = & [ 0, 0, -g ]^T ,\end{aligned}\end{split}\]

Table 5.1 defines the symbols in the above equations, and outlines the symbol dependencies and units.

Table 4.1 Notation and units for flow equations.

symbol	quantity	units
\(\phi({\vec x},t)\)	porosity	[]
\(S_i({\vec x},t)\)	saturation	[]
\({ \vec V}_i({\vec x},t)\)	Darcy velocity vector	[\(L T^{-1}\)]
\(Q_i({\vec x},t)\)	source/sink	[\(T^{-1}\)]
\({\lambda}_i\)	mobility	[\(L^{3} T M^{-1}\)]
\(p_i({\vec x},t)\)	pressure	[\(M L^{-1} T^{-2}\)]
\(\rho_i\)	mass density	[\(M L^{-3}\)]
\({\vec g}\)	gravity vector	[\(L T^{-2}\)]
\({ \bar k}({\vec x},t)\)	intrinsic permeability tensor	[\(L^{2}\)]
\(k_{ri}({\vec x},t)\)	relative permeability	[]
\(\mu_i\)	viscosity	[\(M L^{-1} T^{-1}\)]
\(g\)	gravitational acceleration	[\(L T^{-2}\)]

Here, \(\phi\) describes the fluid capacity of the porous medium, and \(S_i\) describes the content of phase \(i\) in the porous medium, where we have that \(0 \le \phi\le 1\) and \(0 \le S_i\le 1\). The coefficient \({\bar k}\) is considered a scalar here. We also assume that \(\rho_i\) and \(\mu_i\) are constant. Also note that in ParFlow, we assume that the relative permeability is given as \(k_{ri}(S_i)\). The Darcy velocity vector is related to the velocity vector, \({\vec v}_i\), by the following:

(4.21)\[{\vec V}_i= \phi S_i{\vec v}_i.\]

To complete the formulation, we have the following \(\nu\) consititutive relations

(4.22)\[\sum_i S_i= 1 ,\]

(4.23)\[p_{i0} ~=~ p_{i0} ( S_0 ) , ~~~~~~ i = 1 , \ldots , \nu- 1 .\]

where, \(p_{ij} = p_i - p_j\) is the capillary pressure between phase \(i\) and phase \(j\). We now have the \(3 \nu\) equations, (4.18), (4.19), (4.22), and (4.23), in the \(3 \nu\) unknowns, \(S_i, {\vec V}_i\), and \(p_i\).

For technical reasons, we want to rewrite the above equations. First, we define the total mobility, \({\lambda}_T\), and the total velocity, \({\vec V}_T\), by the relations

(4.24)\[\begin{aligned} {\lambda}_T~=~ \sum_{i} {\lambda}_i, \end{aligned}\]

(4.25)\[\begin{aligned} {\vec V}_T~=~ \sum_{i} {\vec V}_i. \end{aligned}\]

After doing a bunch of algebra, we get the following equation for \(p_0\):

(4.26)\[-~ \sum_{i} \left \{ \nabla\cdot {\lambda}_i \left ( \nabla( p_0 ~+~ p_{i0} ) ~-~ \rho_i{\vec g}\right ) ~+~ Q_i \right \} ~=~ 0 .\]

After doing some more algebra, we get the following \(\nu- 1\) equations for \(S_i\):

(4.27)\[\frac{\partial}{\partial t} ( \phi S_i) ~+~ \nabla\cdot \left ( \frac{{\lambda}_i}{{\lambda}_T} {\vec V}_T~+~ \sum_{j \neq i} \frac{{\lambda}_i{\lambda}_j}{{\lambda}_T} ( \rho_i - \rho_j ) {\vec g} \right ) ~+~ \sum_{j \neq i} \nabla\cdot \frac{{\lambda}_i{\lambda}_j}{{\lambda}_T} \nabla p_{ji} ~-~ Q_i ~=~ 0 .\]

The capillary pressures \(p_{ji}\) in (4.27) are rewritten in terms of the constitutive relations in (4.23) so that we have

(4.28)\[p_{ji} ~=~ p_{j0} ~-~ p_{i0} ,\]

where by definition, \(p_{ii} = 0\). Note that equations (4.27) are analytically the same equations as in (4.18). The reason we rewrite them in this latter form is because of the numerical scheme we are using. We now have the \(3 \nu\) equations, (4.26), (4.27), (4.25), (4.19), and (4.23), in the \(3 \nu\) unknowns, \(S_i, {\vec V}_i\), and \(p_i\).

4.7. Transport Equations

The transport equations in ParFlow are currently defined as follows:

(4.29)\[\begin{split}\begin{aligned} \left ( \frac{\partial}{\partial t} (\phi c_{i,j}) ~+~ \lambda_j~ \phi c_{i,j}\right ) & + \nabla\cdot \left ( c_{i,j}{\vec V}_i\right ) \nonumber \\ & = \\ -\left ( \frac{\partial}{\partial t} ((1 - \phi) \rho_{s}F_{i,j}) ~+~ \lambda_j~ (1 - \phi) \rho_{s}F_{i,j}\right ) & + \sum_{k}^{n_{I}} \gamma^{I;i}_{k}\chi_{\Omega^{I}_{k}} \left ( c_{i,j}- {\bar c}^{k}_{ij}\right ) ~-~ \sum_{k}^{n_{E}} \gamma^{E;i}_{k}\chi_{\Omega^{E}_{k}} c_{i,j}\nonumber\end{aligned}\end{split}\]

where \(i = 0, \ldots , \nu- 1\) \((\nu\in \{1,2,3\})\) is the number of phases, \(j = 0, \ldots , n_c- 1\) is the number of contaminants, and where \(c_{i,j}\) is the concentration of contaminant \(j\) in phase \(i\). Recall also, that \(\chi_A\) is the characteristic function of set \(A\), i.e. \(\chi_A(x) = 1\) if \(x \in A\) and \(\chi_A(x) = 0\) if \(x \not\in A\). Table 5.2 defines the symbols in the above equation, and outlines the symbol dependencies and units. The equation is basically a statement of mass conservation in a convective flow (no diffusion) with adsorption and degradation effects incorporated along with the addition of injection and extraction wells.

Table 4.2 Notation and units for transport equation.

symbol	quantity	units
\(\phi({\vec x})\)	porosity	[]
\(c_{i,j}({\vec x},t)\)	concentration fraction	[]
\({\vec V}_i({\vec x},t)\)	Darcy velocity vector	[\(L T^{-1}\)]
\(\lambda_j\)	degradation rate	[\(T^{-1}\)]
\(\rho_{s}({\vec x})\)	density of the solid mass	[\(M L^{-3}\)]]
\(F_{i,j}({\vec x}, t)\)	mass concentration	[\(L^{3} M^{-1}\)]
\(n_{I}\)	number of injection wells	[]
\(\gamma^{I;i}_{k}(t)\)	injection rate	[\(T^{-1}\)]
\(\Omega^{I}_{k}({\vec x})\)	injection well region	[]
\({\bar c}^{k}_{ij}()\)	injected concentration fraction	[]
\(n_{E}\)	number of extraction wells	[]
\(\gamma^{E;i}_{k}(t)\)	extraction rate	[\(T^{-1}\)]
\(\Omega^{E}_{k}({\vec x})\)	extraction well region	[]

These equations will soon have to be generalized to include a diffusion term. At the present time, as an adsorption model, we take the mass concentration term (\(F_{i,j}\)) to be instantaneous in time and a linear function of contaminant concentration :

(4.30)\[F_{i,j}= K_{d;j}c_{i,j},\]

where \(K_{d;j}\) is the distribution coefficient of the component ([\(L^{3} M^{-1}\)]). If (4.30) is substituted into (4.29) the following equation results (which is the current model used in ParFlow) :

(4.31)\[\begin{split}\begin{aligned} (\phi+ (1 - \phi) \rho_{s}K_{d;j}) \frac{\partial}{\partial t} c_{i,j} & ~+~ \nabla\cdot \left ( c_{i,j}{\vec V}_i\right ) \nonumber \\ & ~=~ \nonumber \\ -~(\phi+ (1 - \phi) \rho_{s}K_{d;j}) \lambda_jc_{i,j} & ~+~ \sum_{k}^{n_{I}} \gamma^{I;i}_{k}\chi_{\Omega^{I}_{k}} \left ( c_{i,j}- {\bar c}^{k}_{ij}\right ) ~-~ \sum_{k}^{n_{E}} \gamma^{E;i}_{k}\chi_{\Omega^{E}_{k}} c_{i,j}\end{aligned}\end{split}\]

4.8. Notation and Units

In this section, we discuss other common formulations of the flow and transport equations, and how they relate to the equations solved by ParFlow.

We can rewrite equation (4.19) as

(4.32)\[{\vec V}_i~+~ {\bar K}_i\cdot ( \nabla h_i~-~ \frac{\rho_i}{\gamma} {\vec g}) ~=~ 0 ,\]

where

(4.33)\[\begin{split}\begin{aligned} {\bar K}_i& = & \gamma{\lambda}_i, \\ h_i& = & ( p_i~-~ \bar{p}) / \gamma.\end{aligned}\end{split}\]

Table 5.3 defines the symbols and their units.

Table 4.3 Notation and units for reformulated flow equations.

symbol	quantity	units
\({\vec V}_i\)	Darcy velocity vector	[\(L T^{-1}\)]
\({\bar K}_i\)	hydraulic conductivity tensor	[\(L T^{-1}\)]
\(h_i\)	pressure head	[\(L\)]
\(\gamma\)	constant scale factor	[\(M L^{-2} T^{-2}\)]
\({\vec g}\)	gravity vector	[\(L T^{-2}\)]

We can then rewrite equations (4.26) and (4.27) as

(4.34)\[-~ \sum_{i} \left \{ \nabla\cdot {\bar K}_i \left ( \nabla( h_0 ~+~ h_{i0} ) ~-~ \frac{\rho_i}{\gamma} {\vec g}\right ) ~+~ Q_i \right \} ~=~ 0 ,\]

(4.35)\[\frac{\partial}{\partial t} ( \phi S_i) ~+~ \nabla\cdot \left ( \frac{{\bar K}_i}{{\bar K}_T} {\vec V}_T~+~ \sum_{j \neq i} \frac{{\bar K}_i{\bar K}_j}{{\bar K}_T} \left ( \frac{\rho_i}{\gamma} - \frac{\rho_j}{\gamma} \right ) {\vec g} \right ) ~+~ \sum_{j \neq i} \nabla\cdot \frac{{\bar K}_i{\bar K}_j}{{\bar K}_T} \nabla h_{ji} ~-~ Q_i ~=~ 0 .\]

Note that \({\bar K}_i\) is supposed to be a tensor, but we treat it as a scalar here. Also, note that by carefully defining the input to ParFlow, we can use the units of equations (4.34) and (4.35). To be more precise, let us denote ParFlow input symbols by appending the symbols in table 5.1 with \((I)\), and let \(\gamma= \rho_0 g\) (this is a typical definition). Then, we want:

(4.36)\[\begin{split}\begin{aligned} {\bar k}(I) & = & \gamma{\bar k}/ \mu_0 ; \\ \mu_i(I) & = & \mu_i/ \mu_0 ; \\ p_i(I) & = & h_i; \\ \rho_i(I) & = & \rho_i/ \rho_0 ; \\ g (I) & = & 1 . \end{aligned}\end{split}\]

By doing this, \({\bar k}(I)\) represents hydraulic conductivity of the base phase \({\bar K}_0\) (e.g. water) under saturated conditions (i.e. \(k_{r0} = 1\)).

4.9. Water Balance

ParFlow can calculate a water balance for the Richards’ equation, overland flow and clm capabilities. For a schematic of the water balance in ParFlow please see Maxwell [Max10]. This water balance is computes using pftools commands as described in Manipulating Data: PFTools. There are two water balance storage components, subsurface and surface, and two flux calculations, overland flow and evapotranspiration. The storage components have units [\(L^3\)] while the fluxes may be instantaneous and have units [\(L^3T^{-1}\)] or cumulative over an output interval with units [\(L^3\)]. Examples of water balance calculations and errors are given in the scripts water_balance_x.tcl and water_balance_y.tcl. The size of water balance errors depend on solver settings and tolerances but are typically very small, \(<10^{-10}\) [-]. The water balance takes the form:

(4.37)\[\begin{aligned} \frac{\Delta [Vol_{subsurface} + Vol_{surface}]}{\Delta t} = Q_{overland} + Q_{evapotranspiration} + Q_{source sink} \end{aligned}\]

where \(Vol_{subsurface}\) is the subsurface storage [\(L^3\)]; \(Vol_{surface}\) is the surface storage [\(L^3\)]; \(Q_{overland}\) is the overland flux [\(L^3 T^{-1}\)]; \(Q_{evapotranspiration}\) is the evapotranspiration flux passed from clm or other LSM, etc, [\(L^3 T^{-1}\)]; and \(Q_{source sink}\) are any other source/sink fluxes specified in the simulation [\(L^3 T^{-1}\)]. The surface and subsurface storage routines are calculated using the ParFlow toolset commands pfsurfacestorage and pfsubsurfacestorage respectively. Overland flow out of the domain is calculated by pfsurfacerunoff. Details for the use of these commands are given in PFTCL Commands and Common examples using ParFlow TCL commands (PFTCL). \(Q_{evapotranspiration}\) must be written out by ParFlow as a variable (as shown in Code Parameters) and only contains the external fluxes passed from a module such as clm or WRF. Note that these volume and flux quantities are calculated spatially over the domain and are returned as array values, just like any other quantity in ParFlow. The tools command pfsum will sum these arrays into a single value for the enrite domain. All other fluxes must be determined by the user.

The subsurface storage is calculated over all active cells in the domain, \(\Omega\), and contains both compressible and incompressible parts based on Equation (4.3). This is computed on a cell-by-cell basis (with the result being an array of balances over the domain) as follows:

(4.38)\[\begin{aligned} Vol_{subsurface} = \sum_\Omega [ S(\psi)S_s \psi \Delta x \Delta y \Delta z + S(\psi)(\psi)\phi \Delta x \Delta y \Delta z] \end{aligned}\]

The surface storage is calculated over the upper surface boundary cells in the domain, \(\Gamma\), as computed by the mask and contains based on Equation [eq:kinematic]. This is again computed on a cell-by-cell basis (with the result being an array of balances over the domain) as follows:

(4.39)\[\begin{aligned} Vol_{surface} = \sum_\Gamma \psi \Delta x \Delta y \end{aligned}\]

For the overland flow outflow from the domain, any cell at the top boundary that has a slope that points out of the domain and is ponded will remove water from the domain. This is calculated, for example in the y-direction, as the multiple of Equation [eq:manningsy] and the area:

(4.40)\[\begin{aligned} Q_{overland}=vA= -\frac{\sqrt{S_{f,y}}}{n}\psi_{s}^{2/3}\psi \Delta x=- \frac{\sqrt{S_{f,y}}}{n}\psi_{s}^{5/3}\Delta x \end{aligned}\]