Appendix E — Numerical scheme soil water boundary conditions

Top boundary condition

Appropriate criteria for the procedure with respect to the top boundary condition are important for accurate simulation of rapidly changing soil water fluxes near the soil surface. This is for instance the case with infiltration/runoff events during intensive rain showers or when the soil occasionally gets flooded in areas with shallow groundwater tables. At moderate weather and soil wetness conditions the soil top boundary condition will be flux-controlled. In either very wet or very dry conditions the prevailing water pressure head at the soil surface starts to govern the boundary condition.

Flux controlled top boundary

In case of a flux controlled top boundary the term \(-K_{i-1/2}^{j+\kappa,\kappa p} \left(\frac{h_{i-1}^{j+1,p}-h_i^{j+1,p}} {0.5\left(\Delta z_{i-1}+\Delta z_i\right)} +1\right)\) is replaced by the flux through the soil surface \(q_\text{top}\) (cm d-1) which yields the following expression:

\[ F_1 = \frac{\Delta z_1}{\Delta t^j} \left(\theta^{j+1}_1 - \theta^j_1\right) + q_{top} + K^{j+\kappa}_{1+0.5} \frac{h_1^{j+1}-h_2^{j+1}} {0.5\left(\Delta z_1+\Delta z_2\right)}+ K^{j+\kappa}_{1+0.5} + \Delta Z_1 S^j_{d,1} + \Delta Z_1 S^{j+1}_{m,1} \tag{E.1}\]

where \(q_{top}\) is calculated from external driving factors as net precipitation (\(q_{prec}\)), irrigation (\(q_{irri}\)), melt of a snow pack (\(q_{melt}\)) runon originating from adjacent fields (\(q_{runon}\)) and inundation from adjacent water courses (\(q_{inun}\)):

\[ q_{top} = - q_{prec} - q_{irri} - q_{melt} - q_{runon} - q_{inun} \tag{E.2}\]

Head controlled top boundary

In case of a head controlled top boundary the term \(-K_{i-1/2}^{j+\kappa,\kappa p} \left(\frac{h_{i-1}^{j+1,p}-h_i^{j+1,p}} {0.5\left(\Delta z_{i-1}+\Delta z_i\right)} +1\right)\) is replaced by \(-K_{1/2}^{j+1} \left(\frac{h_0^{j+1}-h_1^{j+1}} {0.5\Delta z_1} +1\right)\), where \(h^{j+1}_0\) is the pressure head at the soil surface at the new time level. The internodal conductivity \(K_{1/2}^{j+1}\) is always treated implicitly.

Within each iteration round and also within each backtracking sub-cycle it is tested whether the combination of \(q_{top}\) and \(h_1^{j+1}\) would lead to \(h_1^{j+1}-0.5\Delta z_1\left(\frac{q_{top}} {K_{1/2}^{j+1}} +1\right) > 0\). In such a case it is decided that the head boundary condition holds and the water balance of the so-called ponding layer is calculated which includes the surface runoff flux and the ponding depth at time level \(j\)+1. The value of \(h_0^{j+1}\)is set to the ponding depth at time level \(j\)+1. The water balance of the ponding layer reads as:

\[ \begin{align} \frac{h_0^{j+1} - h_0^j}{\Delta t} = &-K_{1/2}^{j+1} \left(\frac{h_0^{j+1}-h_1^{j+1}} {0.5\Delta z_1} +1\right) \\ &+ q_{prec} + q_{irri} + q_{melt} + q_{runon} + q_{inun} - q_{e,pond} - q_{runoff} - q_{tl} \end{align} \tag{E.3}\]

Where \(q_{tl}\) is the runoff into the macropores (Section 6.3.1).

The surface runoff flux $q_{runoff} is defined as a function of the ponding height:

\[ q_{runoff} = \begin{cases} 0 & \text{for} \space h^{j+1}_0 \le h_{0,threshold}\\ \alpha \left(h^{j+1}_0 - h_{0,threshold}\right)^\beta & \text{for} \space h^{j+1}_0 > h_{0,threshold}\\ \end{cases} \tag{E.4}\]

Where \(\alpha\) and \(\beta\) are coefficients of the surface runoff equation employed by the SWAP model (see Section 4.1). Substitution of surface runoff expression into water balance equation for the ponding layer yields the relation which is solved in the iteration procedure:

\[ \begin{align} h^{j+1}_0 \left(\frac{1}{\Delta t} + \frac{K^{j+1}_{1/2}}{0.5\Delta z_1} \right) + \alpha \left(h^{j+1}_0 - h_{0,threshold}\right)^\beta &= \frac{h^j_0}{\Delta t}+\frac{K^{j+1}_{1/2}}{0.5\Delta z_1} h^{j+1}_1 -K^{j+1}_{1/2} \\ &+ q_{prec} + q_{melt} + q_{runon} + q_{inun} - q_{e,pond} - q_{tl} \end{align} \tag{E.5}\]

which is solved in the iteration procedure.

Bottom boundary condition

The SWAP model provides a number of options to describe the relation between saturated shallow soil layers with deep groundwater (see Chapter 5). Beside handling the flux controlled boundary condition and head controlled boundary condition, the model has additional capabilities to combine these basic types of conditions. Additional options comprise the handling of:

Flux controlled bottom boundary

For \(i\)=\(n\), the \(K_{i+1/2}^{j+\kappa,\kappa p} \left(\frac{h_i^{j+1,p}-h_{i+1}^{j+1,p}} {0.5\left(\Delta z_i+\Delta z_{i+1}\right)} +1\right)\)-term in Equation 2.54 is replaced by -\(q_{bot}\), which yields:

\[ F_n = \frac{\Delta z_n}{\Delta t^j} \left(\theta^{j+1}_n - \theta^j_n\right) - K^{j+\kappa}_{n-1/2} \frac{h_{n-1}^{j+1}-h_n^{j+1}} {0.5\left(\Delta z_{n-1}+\Delta z_n\right)} - q_{bot} + \Delta z_n \left(S^{j+\kappa}_{a,n} + S^j_{d,n} + S^{j+1}_{m,n}\right) \tag{E.6}\]

Beside the flux boundary condition, the SWAP model has options to handle groundwater level dependent bottom fluxes. The flux can be formulated as an exponential function of groundwater level, or as the difference between groundwater level and hydraulic head in deep groundwater outside the flow domain divided by a flow resistance. Such a flux is calculated explicitly at the start of the current time step and is treated as a flux condition in the numerical scheme.

Head controlled bottom boundary

For \(i\)=\(n\), \(\Delta z_{i+1}\) is set to zero and \(h_{i+1}^{j+1,p} = h_{bot}\), which leads to the following expression:

\[ \begin{align} F_n &= \frac{\Delta z_n}{\Delta t^j} \left(\theta^{j+1}_n - \theta^j_n\right) - K^{j+\kappa}_{n-1/2} \left(\frac{h_{n-1}^{j+1}-h_n^{j+1}} {0.5\left(\Delta z_{n-1}+\Delta z_n\right)}+1\right) \\ &+ K^{j+\kappa}_{n+1/2}\left( \frac{h_n^{j+1}-h_{bot}} {0.5\Delta z_n}+1\right) + \Delta z_n \left(S^{j+\kappa}_{a,n} + S^j_{d,n} + S^{j+1}_{m,n}\right) \end{align} \tag{E.7}\]

Predefined groundwater levels

First, the lowest partially unsaturated compartment is searched for and is called \(n^*\). The set of \(n^*\) non-linear equations for \(F(h)\) is then solved for the unsaturated compartments. The bottom boundary condition for this set of equations is defined by:

\[ \begin{align} F_{n^*} &= \frac{\Delta z_{n^*}}{\Delta t^j} \left(\theta^{j+1}_{n^*} - \theta^j_{n^*}\right) - K^{j+\kappa}_{{n^*}-1/2} \left(\frac{h_{{n^*}-1}^{j+1}-h_{n^*}^{j+1}} {0.5\left(\Delta z_{{n^*}-1}+\Delta z_{n^*}\right)}+1\right) \\ &+ K^{j+\kappa}_{{n^*}+1/2}\left( \frac{h_{n^*}^{j+1}-h^{j+1}_{n^*+1}} {0.5\left(\Delta z_{n^*} + \Delta z_{n^* + 1}\right)}+1\right) + \Delta z_n \left(S^{j+\kappa}_{a,n^*} + S^j_{d,n^*} + S^{j+1}_{m,n^*}\right) \end{align} \tag{E.8}\]

The groundwater level is situated between the nodal points \(n^*\) and \(n^{*+1}\). The pressure head of nodal point \(n^*+1\) is approximated \(h_{n^*}^{j+1}\) by:

\[ h_{n^*+1}^{j+1} = - h_{n^*}^{j+1} \frac{gwl-z_{n^*+1}}{z_{n^*}-gwl} \tag{E.9}\]

where \(z_{n^*}\) is the height of the node in compartment \(n^*\) and \(gwl\) is the groundwater level. Substitution into Equation 2.55 yields:

\[ \begin{align} F_{n^*} &= \frac{\Delta z_{n^*}}{\Delta t^j} \left(\theta^{j+1}_{n^*} - \theta^j_{n^*}\right) - K^{j+\kappa}_{{n^*}-1/2} \left(\frac{h_{{n^*}-1}^{j+1}-h_{n^*}^{j+1}} {0.5\left(\Delta z_{{n^*}-1}+\Delta z_{n^*}\right)}+1\right) \\ &+ K^{j+\kappa}_{{n^*}+1/2}\left( \frac{h_{n^*}^{j+1} \left(\frac{z_{n^*}-z_{n^*+1}} {Z_{n^*}-gwl}\right) } {0.5\left(\Delta z_{n^*} + \Delta z_{n^* + 1}\right)}+1\right) + \Delta z_n \left(S^{j+\kappa}_{a,n^*} + S^j_{d,n^*} + S^{j+1}_{m,n^*}\right) \end{align} \tag{E.10}\]

After iteratively solving the set of equations for \(i \le n^*\), the pressure head profile of the compartments \(i > n^*\) can be calculated from the pressure head of the two adjacent upward nodes.

\[ \begin{align} h_{i+1}^{j+1} &= h_{i}^{j+1} \left(1+\frac{K^{j+\kappa}_{i-1/2} \left(\Delta z_i + \Delta z_{i+1}\right)}{K^{j+\kappa}_{i+1/2} \left(\Delta z_{i-1} + \Delta z_i\right)}\right) - h_{i-1}^{j+1} \frac{K^{j+\kappa}_{i-1/2} \left(\Delta z_i + \Delta z_{i+1}\right)}{K^{j+\kappa}_{i+1/2} \left(\Delta z_{i-1} + \Delta z_i\right)} \\ &+ 0.5\left(\Delta z_i + \Delta z_{i+1}\right) \left(1-\frac{K^{j+\kappa}_{i-1/2}}{K^{j+\kappa}_{i+1/2}}\right) \\ &+ \frac{0.5\left(\Delta z_i + \Delta z_{i+1}\right)}{K^{j+\kappa}_{i+1/2}} \left(\frac{\Delta z_i}{\Delta t^j} \left(\theta_{sat,i} - \theta^j_i\right) + \Delta z_i\left(S^{j+\kappa}_{a,i} + S^j_{d,i} + S^{j+1}_{m,i} \right)\right) \end{align} \tag{E.11}\]

Cauchy relation for the bottom boundary

The flux through the bottom boundary is defined by the difference the hydraulic head at the lower boundary and the hydraulic head \(\phi\) (cm) of the regional groundwater specified by the user, divided by a flow resistance \(c\) (d). The hydraulic head at the lower boundary is approximated by the pressure head of the lowest nodal point plus the elevation head of node \(n\).

\[ q_{bot} = \frac{h^{j+1}_{n} + \Delta z_n - \phi} {\frac{0.5\Delta z_n}{K^{j+\kappa}_{n+1/2}}+c} \tag{E.12}\]

Substitution into Equation 2.55 yields:

\[ \begin{align} F_n &=\frac{\Delta z_n}{\Delta t^j} \left(\theta^{j+1}_n - \theta^j_n\right) - K^{j+\kappa}_{n-1/2} \left(\frac{h_{n-1}^{j+1}-h_n^{j+1}} {0.5\left(\Delta z_{n-1}+\Delta z_n\right)}+1\right) \\ &+ K^{j+\kappa}_{n+1/2}\left( \frac{h_n^{j+1}-\left(\phi - z_n\right)} {0.5\Delta z_n + cK^{j+\kappa}_{n+1/2}}+1\right) + \Delta z_n \left(S^{j+\kappa}_{a,n} + S^j_{d,n} + S^{j+1}_{m,n}\right)\\ \end{align} \tag{E.13}\]

Seepage face

The seepage face option is used to simulate the soil moisture flow in a lysimeter with an open outlet at the bottom. No outflow occurs when the bottom soil layer is still unsaturated. Since the flow resistance of the outlet is negligible small, no positive pressure head values will be build up at the bottom when the soil water percolates at the bottom. Within the iteration cycle for solving the numerical expression of the Richards equation, it is checked whether the flux or the head controlled boundary condition prevails. When \(h_n^{j+1}+1/2\Delta z_n > 0\) the bottom flux q_bot is set to zero, but when \(h_n^{j+1}+ 0.5\Delta z_n\) tends to take values larger than zero, the pressure at the bottom is set to zero (\(h_{bot}\)=0). The numerical implementation is as follows:

\[ F_n = \begin{cases} \begin{array}{ll} \frac{\Delta z_n}{\Delta t^j}&\left(\theta^{j+1}_n - \theta^j_n\right) - K^{j+\kappa}_{n-1/2} \left(\frac{h_{n-1}^{j+1}-h_n^{j+1}} {0.5\left(\Delta z_{n-1}+\Delta z_n\right)}+1\right)\\ & +\Delta z_n \left(S^{j+\kappa}_{a,n} + S^j_{d,n} + S^{j+1}_{m,n}\right) \end{array} & \text{for} \space h^{j+1}_n + 0.5 \Delta z_n < 0\\ \begin{array}{ll} \frac{\Delta z_n}{\Delta t^j}&\left(\theta^{j+1}_n - \theta^j_n\right) - K^{j+\kappa}_{n-1/2} \left(\frac{h_{n-1}^{j+1}-h_n^{j+1}} {0.5\left(\Delta z_{n-1}+\Delta z_n\right)}+1\right) \\ & + K^{j+\kappa}_{n+1/2}\left( \frac{h_n^{j+1}} {0.5\Delta z_n}+1\right) +\Delta z_n \left(S^{j+\kappa}_{a,n} + S^j_{d,n} + S^{j+1}_{m,n}\right) \end{array} & \text{for} \space h^{j+1}_n + 0.5\Delta z_n \ge 0\\ \end{cases} \tag{E.14}\]

Free drainage

The free drainage option is applied for soil profiles with deep groundwater levels. The bottom flux is only provoked by gravity flow and the head pressure gradient equals zero:

\[ q_{bot} = -K^{j+1}_{n+1/2} \left(\left(\frac{\partial h}{\partial z}\right)_{n+1/2}+1\right) \to q_{bot} = -K^{j+1}_{n+1/2} \tag{E.15}\]

Substitution into Equation 2.55 yields:

\[ \begin{array}{ll} F_n =&\frac{\Delta z_n}{\Delta t^j} = \left(\theta^{j+1}_n - \theta^j_n\right) - K^{j+\kappa}_{n-1/2} \left(\frac{h_{n-1}^{j+1}-h_n^{j+1}} {0.5\left(\Delta z_{n-1}+\Delta z_n\right)}+1\right) + K^{j+\kappa}_{n+1/2}\\ &+ \Delta z_n \left(S^{j+\kappa}_{a,n} + S^j_{d,n} + S^{j+1}_{m,n}\right)\\ \end{array} \tag{E.16}\]