the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Capturing soilwater and groundwater interactions with an iterative feedback coupling scheme: new HYDRUS package for MODFLOW
Jicai Zeng
Jinzhong Yang
Yuanyuan Zha
Liangsheng Shi
Accurately capturing the complex soilwater and groundwater interactions is vital for describing the coupling between subsurface–surface–atmospheric systems in regionalscale models. The nonlinearity of Richards' equation (RE) for water flow, however, introduces numerical complexity to large unsaturated–saturated modeling systems. An alternative is to use quasi3D methods with a feedback coupling scheme to practically join submodels with different properties, such as governing equations, numerical scales, and dimensionalities. In this work, to reduce the nonlinearity in the coupling system, two different forms of RE are switched according to the soilwater content at each numerical node. A rigorous multiscale water balance analysis is carried out at the phreatic interface to link the soilwater and groundwater models at separated spatial and temporal scales. For problems with dynamic groundwater flow, the nontrivial coupling errors introduced by the saturated lateral fluxes are minimized with a movingboundary approach. It is shown that the developed iterative feedback coupling scheme results in significant error reduction and is numerically efficient for capturing drastic flow interactions at the water table, especially with dynamic local groundwater flow. The coupling scheme is developed into a new HYDRUS package for MODFLOW, which is applicable to regionalscale problems.
 Article
(5044 KB)  Fulltext XML
 BibTeX
 EndNote
Numerical modeling of the soilwater and groundwater interactions has to deal with flow components and governing equations at different scales. This adds significant complexity to model development and calibration. Unsaturated soilwater and saturated groundwater flows with similar properties are usually integrated into a whole modeling system. Although physically consistent and numerically rigorous, methods involving the 3D Richards' equation (RE; Richards, 1931) tend to be computationally expensive and numerically unstable due to the large nonlinearity and the demand for dense discretization (Kumar et al., 2009; Maxwell and Miller, 2005; Panday and Huyakorn, 2004; Thoms et al., 2006; Zha et al., 2013a), especially for problems with multiscale properties. In this work, parsimonious approaches, which appear in different governing equations and coupling schemes, are developed for modeling the soilwater and groundwater interactions at regional scales.
Simplifying the soilwater flow details into upper flux boundaries has been widely used to simulate largescale saturated flow dynamics, such as the MODFLOW package and its variants (Langevin et al., 2017; Leake and Claar, 1999; McDonald and Harbaugh, 1988; Niswonger et al., 2011; Panday et al., 2013; Zeng et al., 2017). At the local scale, in contrast, the unsaturated flow processes are usually approximated with reasonable simplifications and assumptions in RE (Bailey et al., 2013; Liu et al., 2016; Paulus et al., 2013; Šimůnek et al., 2009; van Dam et al., 2008; Yakirevich et al., 1998; Zha et al., 2013b).
The original RE, also known as the mixedform RE, takes the pressure head (h) as the driving force variable, while soil moisture content (θ) is the mass accumulation variable (Krabbenhøft, 2007). To solve the mixedform RE, either h or θ, or a switching of both, is assigned as the primary variable. The hform RE is widely employed for unsaturated–saturated flow simulation, especially in heterogeneous soils, such as the HYDRUS package (Šimůnek et al., 2016). Significant improvement in mass conservation has been achieved with Celia's modification (Celia et al., 1990). Then, efforts were made to combine the advantages in the original and the Celiaformat hform REs by switching their storage terms (Hao et al., 2005; Zadeh, 2011). However, these models still suffer from a high computational cost and low numerical robustness when dealing with rapidly changing atmospheric boundary conditions (Crevoisier et al., 2009; Zha et al., 2017). The θform RE, addressing the above problems, is inherently mass conservative and less nonlinear in the averaged nodal hydraulic diffusivity when the soil is dry (Warrick, 1991; Zha et al., 2013b). However, the θform RE is not applicable to saturated and heterogeneous soils (Crevoisier et al., 2009; Zha et al., 2013b). In this work, to take advantage of both forms of RE, the governing equations, rather than primary variables (Diersch and Perrochet, 1999; Forsyth et al., 1995; Zha et al., 2013a), are switched at each node according to its saturation degree.
For regional problems, the vadose zone is usually conceptualized into paralleled soil columns without lateral connections. The resulting quasi3D coupling scheme (Kuznetsov et al., 2012; Seo et al., 2007; Xu et al., 2012; Zhu et al., 2012) significantly reduces the dimensionality and complexity. According to how the messages are transferred across the coupling interface, the quasi3D methods are categorized into (1) the fully coupling scheme, which simultaneously builds the nodal hydraulic connections of submodels at both sides and implicitly solves the assembled matrices; (2) the oneway coupling scheme, which delivers the soilwater model solutions onto the groundwater model without feedback mechanism; and (3) the feedback (or twoway) coupling scheme, which explicitly exchanges the head and/or flux solutions in vicinity of the interface nodes.
The fully coupling scheme (Gunduz and Aral, 2005; Zhu et al., 2012) is numerically rigorous but tends to increase the computational burden under practical conditions. For example, the potentially conditional diagonal dominance causes nonconvergence for the iterative solvers (Edwards, 1996). Owing to high nonlinearity in the soilwater submodels, the assembled matrices can only be solved with unified small time steps, which adds to the computational expense. The oneway coupling scheme, as adopted by the UZF1 package for MODFLOW (Grygoruk et al., 2014; Niswonger et al., 2006) as well as the free drainage mode in SWAP package for MODFLOW (Xu et al., 2012), assumes that the watertable depth has a minor influence on flow interactions at the phreatic interface and is thus problem specific.
The feedback coupling methods, in contrast, are widely used (Kuznetsov et al., 2012; Seo et al., 2007; Shen and Phanikumar, 2010; Stoppelenburg et al., 2005; Xie et al., 2012; Xu et al., 2012) as a compromise of numerical accuracy and computational cost. In a feedback coupling scheme, the soilwater and groundwater submodels can be built with governing equations and numerical schemes at different scales. For flow processes with multiscale components, such as boundary geometries, parameter heterogeneities, and hydrologic stresses, the scaleseparation strategy can be implemented easily. Although the feedback coupling method is numerically more rigorous than a oneway coupling method, and tends to reduce the inconsistency of head and/or flux interfacial boundaries, some concerns arise.
The first concern is the numerical efficiency of the iterative and noniterative feedback coupling methods. The noniterative approach (Twarakavi et al., 2008; Xu et al., 2012) usually leads to significant error accumulation when dealing with a dynamically fluctuating water table, especially with large timestep sizes. The iterative methods, in contrast (Kuznetsov et al., 2012; Stoppelenburg et al., 2005; Xie et al., 2012), by converging the head and/or flux solutions at the coupling interface, are numerically rigorous but computationally expensive, especially when solving the coupled submodels with a unified timestepping scheme (Kuznetsov et al., 2012). Good balance between cost and effect is needed to maintain practical utility of the iterative feedback coupling scheme.
The second concern lies in the scalemismatching problem. For groundwater models (Harbaugh et al., 2005; Langevin et al., 2017; Lin et al., 2010; McDonald and Harbaugh, 1988), the specific yield at the phreatic surface is usually represented by a simple largescale parameter, while for soilwater models (Niswonger et al., 2006; Šimůnek et al., 2009; Thoms et al., 2006), the smallscale phreatic water release is influenced by the watertable depth and the unsaturated soil moisture profile (Dettmann and Bechtold, 2016; Nachabe, 2002). Delivering smallscale solutions of the soilwater models onto the largescale interfacial boundary of the groundwater model, as well as maintaining the global mass balance, usually introduces significant nonlinearity to the entire coupling system (Stoppelenburg et al., 2005). Conditioned by this, the mismatch of numerical scales in the coupled submodels causes significant coupling errors and instability. A multiscale water balance analysis at the phreatic surface helps to relieve such difficulties.
The third concern is the nontrivial lateral fluxes between the saturated regions controlled by the vertical soil columns, which are usually not considered in previous study (Seo et al., 2007; Xu et al., 2012). Though rigorous water balance analysis is conducted to address such inadequacy (Shen and Phanikumar, 2010), the lateral fluxes solved with a 2D groundwater model usually require additional efforts to build water budget equations in each subdivision represented by the soil columns. A movingboundary strategy helps to avoid the saturated lateral flow in the groundwater body.
In this work, the h and θform of the 1D RE are switched at the equation level to obtain a new HYDRUS package. To handle three of the aforementioned concerns, a multiscale water balance analysis is carried out at the phreatic surface to conserve head and/or flux consistent at the coupling interface. An iterative feedback coupling scheme is developed for linking the unsaturated and saturated flow models at disparate scales. The saturated lateral fluxes between the soil columns are fully removed from the interfacial water balance equation, making it a movinginterface coupling framework. The head solution of MODFLOW 2005 (Harbaugh et al., 2005; Langevin et al., 2017) and flux solution of HYDRUS1D (Šimůnek et al., 2009) are relaxed to meet consistency at the phreatic surface.
In this paper, the governing equations at different scales, the multiscale water balance analysis at the phreatic surface, and the iterative feedback coupling scheme for solving the whole system, are presented in Sect. 2. Synthetic numerical experiments are described in Sect. 3. Numerical performance of the developed model is investigated in Sect. 4. Conclusions are drawn in Sect. 5.
To address the aforementioned first concern, governing equations for subsurface flow are given at different levels of complexity (Sect. 2.1), numerical solutions of these equations are presented (Sect. 2.2), and nonlinearity in the soilwater submodels is reduced by a generalized switching scheme that chooses appropriate forms of RE according to the hydraulic conditions at each numerical node (Sect. 2.3). Then, an iterative feedback coupling scheme is developed to solve the soilwater and groundwater models at independent scales (Sect. 2.4). As for the second concern, a multiscale water balance analysis is conducted to deal with the scalemismatching problem at the phreatic surface (Sect. 2.5). To cope with the third concern, a moving Dirichlet boundary at the groundwater table is assigned to the soilwater submodels (see Appendix Sect. A1); the Neumann upper boundary for the saturated model is provided in Sect. A2, and the relaxed iterative feedback process is presented in Sect. A3.
2.1 Governing equations
The mass conservation equation for unsaturated–saturated flow is given by
where t is time (T); θ (L^{3} L^{−3}) is volumetric moisture content; h (L) is the pressure head; β is 1 for the saturated region but zero for the unsaturated region; C (L^{−1}) is the soilwater capacity ($C=\partial \mathit{\theta}/\partial h$) for the unsaturated region but zero for saturated region; μ_{s} (L^{−1}) is the specific elastic storage; and q (L T^{−1}) is the vector for Darcian flux calculated by the following:
where K (L T^{−1}) is the hydraulic conductivity, K=K(θ), and H (L) is the potentiometric head, $H=h+z$, in which z is the vertical location with a coordinate positive upward. Combining Eqs. (1) and (2) results in the governing equation for groundwater flow,
With the assumption that the horizontal unsaturated flows are negligible, the regional vadose zone is usually represented by an assembly of paralleled soil columns. The generalized 1D RE is represented by a switchable format,
where ψ is the primary variable. For an hform RE, ψ=h, $\widehat{C}=C$, and $\widehat{K}=K$, while for a θform RE, ψ=θ, $\widehat{C}=\mathrm{1}$, and $\widehat{K}=D$; where D (L^{2} T^{−1}) is the hydraulic diffusivity, $D=K/C$.
2.2 Numerical approximation
The governing equation for the saturated zone (Eq. 3) is spatially and temporally approximated in the same form as the MODFLOW 2005 model (Harbaugh et al., 2005; Langevin et al., 2017). Celia's modification (Celia et al., 1990; Šimůnek et al., 2009) is applied to the hform 1D RE for temporal approximation. Both forms of RE are handled with a temporally backward finite difference discretization (Zha et al., 2013b, 2017). Each submodel is solved by a Picard iteration scheme, which is widely used in some popular codes or software packages (van Dam et al., 2008; Šimůnek et al., 2016).
The spatial discretization of Eq. (4), as well as the water balance analysis for each node, is based on the nodal flux in element $i+\mathrm{1}/\mathrm{2}$ (bounded by nodes i and i+1), which is
where the superscripts j and k are the levels of time and inner iteration, the subscript i (or $i+\mathrm{1}/\mathrm{2}$) is the number of node (or element), and $\mathrm{\Delta}{z}_{i+\mathrm{1}/\mathrm{2}}$ is the length of the element $i+\mathrm{1}/\mathrm{2}$, $\mathrm{\Delta}{z}_{i+\mathrm{1}/\mathrm{2}}=({z}_{i+\mathrm{1}}{z}_{i})$. When a soil interface exists at node i, for example, the soil moisture contents in elements $i\mathrm{1}/\mathrm{2}$ and $i+\mathrm{1}/\mathrm{2}$ are discontinuous at node i, thus dissatisfying the θform RE. To address this problem, the correction term ${\mathit{\epsilon}}_{i+\mathrm{1}/\mathrm{2}}^{j+\mathrm{1},k}$, suggested by Zha et al. (2013b), is employed to handle the heterogeneous interface at nodes i and i+1,
where variables (scalar) ${\overrightarrow{\mathit{\psi}}}_{i+\mathrm{1}}^{j+\mathrm{1},k}$ and ${\overleftarrow{\mathit{\psi}}}_{i}^{j+\mathrm{1},k}$ are the continuously distributed ψ within element $i+\mathrm{1}/\mathrm{2}$, i.e., between the vertices i and i+1.
When ψ=h, or when ψ=θ, but no heterogeneity occurs, we get ${\mathit{\psi}}_{i+\mathrm{1}}^{j+\mathrm{1},k}={\overrightarrow{\mathit{\psi}}}_{i+\mathrm{1}}^{j+\mathrm{1},k}$ and ${\mathit{\psi}}_{i}^{j+\mathrm{1},k}={\overleftarrow{\mathit{\psi}}}_{i}^{j+\mathrm{1},k}$, so ${\mathit{\epsilon}}_{i+\mathrm{1}/\mathrm{2}}^{j+\mathrm{1},k}=\mathrm{0}$. When ψ=θ, with soil interfaces at node i or i+1, ${\overrightarrow{\mathit{\psi}}}_{i+\mathrm{1}}^{j+\mathrm{1},k}=\mathit{\theta}\left({h}_{i+\mathrm{1}}^{j+\mathrm{1},k},{\mathit{p}}_{i+\mathrm{1}/\mathrm{2}}\right)$ and ${\overleftarrow{\mathit{\psi}}}_{i}^{j+\mathrm{1},k}=\mathit{\theta}\left({h}_{i}^{j+\mathrm{1},k},{\mathit{p}}_{i+\mathrm{1}/\mathrm{2}}\right)$. It is obvious that ${\mathit{\psi}}_{i}^{j+\mathrm{1},k}\ne {\overleftarrow{\mathit{\psi}}}_{i}^{j+\mathrm{1},k}$ (or ${\mathit{\psi}}_{i+\mathrm{1}}^{j+\mathrm{1},k}\ne {\overrightarrow{\mathit{\psi}}}_{i+\mathrm{1}}^{j+\mathrm{1},k})$, so ${\mathit{\epsilon}}_{i+\mathrm{1}/\mathrm{2}}^{j+\mathrm{1},k}\ne \mathrm{0}$.
Hereinafter, ${\mathit{P}}_{i+\mathrm{1}/\mathrm{2}}$ represents the soil parameters in element $i+\mathrm{1}/\mathrm{2}$. For example, in a van Genuchten model (van Genuchten, 1980), ${\mathit{P}}_{i+\mathrm{1}/\mathrm{2}}=$ (θ_{r}, θ_{s}, n, m, α, k_{s}), where θ_{r} (L^{3} L^{−3}) and θ_{s} (L^{3} L^{−3}) are the residual and saturated soil moisture contents; α (L^{−1}), n, and m are the poresize distribution parameters, $m=\mathrm{1}\mathrm{1}/n$; and k_{s} (L T^{−1}) is the saturated hydraulic conductivity.
2.3 Switching Richards' equation
Due to lower nonlinearity of hydraulic diffusivity (D) for dry soils (Zha et al., 2013b) and the avoidance of the mass balance error by removing the soilwater capacity in the storage term, the θform RE is more robust than the hform RE, especially when dealing with rapidly changing atmospheric boundary conditions (Zeng et al., 2018). In our work, the h and θform REs are switched at each node according to its effective saturation, Se. The resulting hybrid matrix equation set is solved by the Picard iteration. The empirical effective saturation for doing switching varies with soil type and is suggested to be Se^{crit}=0.4–0.9, the state when both the h and θform REs are stable and efficient. When Se ≥ Se^{crit}, the soil moisture is closer to saturation, so the hform RE is chosen as the governing equation; otherwise, when it undergoes dry soil condition, the θform RE is preferred.
For element $i+\mathrm{1}/\mathrm{2}$, when the governing equations for nodes i and i+1 are identical, the spatial approximation of nodal flux is given by Eq. (5). When the governing equations differ at nodes i and i+1, a switched element is produced. When taking ψ_{i}=θ_{i} and ${\mathit{\psi}}_{i+\mathrm{1}}={h}_{i+\mathrm{1}}$ as an example, the nodal fluxes calculated by Eq. (5) for different forms of RE have to be carefully handled by substituting ${\mathit{\theta}}_{i+\mathrm{1}}^{j+\mathrm{1},k+\mathrm{1}}$ with ${\mathit{\theta}}_{i+\mathrm{1}}^{j+\mathrm{1},k}$, while ${h}_{i}^{j+\mathrm{1},k+\mathrm{1}}$ is replaced by ${h}_{i}^{j+\mathrm{1},k}$. When ψ_{i}=h_{i} and ${\mathit{\psi}}_{i+\mathrm{1}}={\mathit{\theta}}_{i+\mathrm{1}}$, in contrast, ${h}_{i+\mathrm{1}}^{j+\mathrm{1},k+\mathrm{1}}$ is replaced by ${h}_{i+\mathrm{1}}^{j+\mathrm{1},k}$, while ${\mathit{\theta}}_{i}^{j+\mathrm{1},k+\mathrm{1}}$ is replaced by ${\mathit{\theta}}_{i}^{j+\mathrm{1},k}$. The resulting equivalent nodal fluxes ${q}_{i+\mathrm{1}/\mathrm{2}}^{h}$ and ${q}_{i+\mathrm{1}/\mathrm{2}}^{\mathit{\theta}}$ are then weighted to obtain an approximation by
where ω is the weighting factor, $\mathrm{0}\le \mathit{\omega}\le \mathrm{1}$. In our work, ω=0.5 is applied to implicitly maintain the unknown variables of both ${h}_{i+\mathrm{1}}^{j+\mathrm{1},k+\mathrm{1}}$ and ${\mathit{\theta}}_{i}^{j+\mathrm{1},k+\mathrm{1}}$. Specifically, when ω=1, the hform RE is used at nodes i and i+1; when ω=0, the θform RE is employed instead. A detailed study on the switching of RE between two ends of the soil moisture condition, as well as the description of the numerical formation, can be found in Zeng et al. (2018).
Note that the equation switching method takes full advantage of the θ and hform REs, which is different from the traditional primary variable switching schemes (Diersch and Perrochet, 1999; Forsyth et al., 1995; Zha et al., 2013a). In our work, the switchingRE approach is incorporated into a new HYDRUS package.
2.4 Iterative feedback coupling scheme
The Dirichlet and Neumann boundaries are iteratively transferred across the phreatic interface. The groundwater head solution serves as the headspecified lower boundary of the soil columns, while the unsaturated solution is converted into the fluxspecified upper boundary of the groundwater model. Due to moderate variation of the groundwater flow, the predicted watertable solution is usually adopted in advance as the Dirichlet lower boundary of the finescale soilwater flow models (Seo et al., 2007; Shen and Phanikumar, 2010; Xu et al., 2012), which in sequence then provides the Neumann upper boundary for successively solving the coarsescale groundwater flow model. Section A1 provides the method for a moving Dirichlet lower boundary, while Sect. A2 presents the Neumann upper boundary for the 3D groundwater model. In Sect. A3, the relaxed iterative feedback coupling scheme is used to solve the unsaturated–saturated submodels at two sides of the coupling interface.
2.5 Multiscale water balance analysis
Coupling models at different scales requires consistency in their spatial and temporal scales at the interface (Downer and Ogden, 2004; Rybak et al., 2015). A space and timesplitting strategy (see Fig. 1) is adopted to separate submodels at different scales. That is, the soilwater models are established by $\mathrm{\Delta}z={\mathrm{10}}^{\mathrm{3}}$–10^{0} m and $\mathrm{\Delta}t={\mathrm{10}}^{\mathrm{5}}$–10^{0} days, while for the saturated model, the grid sizes are Δx=10^{0}–10^{3} m, and timestep sizes are Δt=10^{0}–10^{1} days. Water balance at one side of the interface is conserved by scale matching of boundary conditions provided by the submodel on the other side. For unsaturated flow, Richards' equation requires fine discretization of space and time (Miller et al., 2006; Vogel and Ippisch, 2008), while for saturated flow, coarse spatial and temporal grids produce adequate solutions at a large scale (Mehl and Hill, 2004; Zeng et al., 2017). To approximate the upper boundary flux of the groundwater flow model, a multiscale water balance analysis is conducted within each step of the largescale saturated flow model. At small spatial and temporal scales, e.g., within a macrotime step $\mathrm{\Delta}T={T}^{J+\mathrm{1}}{T}^{J}$ and at a local area of interest (with thickness of $\stackrel{\mathrm{\u203e}}{M}={z}_{\mathrm{s}}{z}_{\mathrm{0}}$, where z_{s} and z_{s} are defined in Sect. A2), the specific storage term in Eq. (1) is vertically integrated into a transient onedimensional expression (Dettmann and Bechtold, 2016),
where w (L) is the amount of unsaturated water in the moving balancing domain (see Fig. 2b), $w\left(t\right)={\int}_{{z}_{\mathrm{t}}\left(t\right)}^{{z}_{\mathrm{s}}}\mathit{\theta}(t,z)\mathrm{d}z$; $\mathrm{\Delta}{z}_{\mathrm{t}}=\sum _{j=\mathrm{1}}^{N}\mathrm{d}{z}_{\mathrm{t}}^{j}={z}_{\mathrm{t}}\left({T}^{J+\mathrm{1}}\right){z}_{\mathrm{t}}\left({T}^{J}\right)$ is the total fluctuation of the phreatic surface during $\mathrm{\Delta}T=\sum _{j=\mathrm{1}}^{N}\mathrm{d}{t}^{j}={T}^{J+\mathrm{1}}{T}^{J}$, and the subscript t in z_{t} here means the groundwater table; and θ_{s} is the saturated soilwater content. Approaching a transient state at time t, the water balance in a moving water balancing domain (see $z\in ({z}_{\mathrm{t}},{z}_{\mathrm{s}}$) in Fig. 2b) during a smallscale time step dt (defined in Fig. 1b) is given by
where q_{top}(t) and q_{bot}(t) (L T^{−1}) are the nodal fluxes into and out of the moving balancing domain at a fixed top boundary (z_{s}) and a moving bottom boundary (z_{b}= min(${z}_{\mathrm{t}}\left(t\right),{z}_{\mathrm{t}}(t\mathrm{d}t)\left)\right)$, i.e., ${q}_{\mathrm{top}}={\left(K\left(h\right)\cdot \partial (h+z)/\partial z\right}_{z={z}_{\mathrm{s}}}$ and ${q}_{\mathrm{bot}}={\left(K\left(h\right)\cdot \partial (h+z)/\partial z\right}_{z={z}_{\mathrm{b}}}$ (positive into the balancing domain and negative outside); d${z}_{\mathrm{t}}={z}_{\mathrm{t}}\left(t\right){z}_{\mathrm{t}}(t\mathrm{d}t)$ is the transient fluctuation of the phreatic surface during dt; and l (T^{−1}) is the saturated lateral flux into the balancing domain at time t (see Fig. 2b). Taking Γ as the lateral boundary of a subdomain, the lateral flux $l=\underset{x,y,z\in \mathrm{\Omega}}{\iiint}\left[\frac{\partial}{\partial x}\left(K\frac{\partial H}{\partial x}\right)+\frac{\partial}{\partial y}\left(K\frac{\partial H}{\partial y}\right)\right]\mathrm{d}x\mathrm{d}y\mathrm{d}z/\underset{x,y,z\in \mathrm{\Omega}}{\iiint}\mathrm{d}x\mathrm{d}y\mathrm{d}z$ is supposed to be constant during ΔT; Ω is the volume of the saturated domain controlled by a soil column, which is horizontally projected into Π. Temporally integrating Eq. (9) from time T^{J} to T^{J+1} produces
where R_{top} (L) is the cumulative water flux at z_{s}, ${R}_{\mathrm{top}}={\int}_{{T}^{J}}^{{T}^{J+\mathrm{1}}}{q}_{\mathrm{top}}\left(t\right)\mathrm{d}t$. Note that R_{top} equals F_{top} in Eq. (A3) (Sect. A2). R_{bot} (L) is the cumulative water flux out of the moving balancing domain, ${R}_{\mathrm{bot}}={\int}_{{T}^{J}}^{{T}^{J+\mathrm{1}}}{q}_{\mathrm{bot}}\left(t\right)\mathrm{d}t$. ε_{l} (L) is the cumulative lateral input water into the moving balancing domain,
where N is the number of time steps for the smallscale soilwater model within a macrotime step ΔT, and ${\mathit{\epsilon}}_{\mathrm{l}}^{\prime}$ is the nontrivial saturated later flux produced by a stationaryboundary method (Seo et al., 2007; Xu et al., 2012). By taking R_{top} as the specific recharge at z_{s}, the smallscale specific yield ${\stackrel{\mathrm{\u0303}}{S}}_{y}$ is derived from Eqs. (8) and (10) as
Supposing that z_{t} is linearly fluctuating in time, i.e., ${z}_{\mathrm{t}}=a\cdot t+$ b (where a and b are constants), we get the watertable change during a smallscale step (dt) by d${z}_{\mathrm{t}}=a\cdot \mathrm{d}t$; thus, ε_{l}=o(dt^{2}), which means that linearly refining the local timestep size (dt) in the soilwater model brings about at least quadratic approximation of ε_{l} towards zero. Thus ε_{l} can be neglected from the smallscale mass balance analysis. In the developed model, the largescale specific yield, ${\stackrel{\mathrm{\u203e}}{S}}_{y}$ in Eq. (A3), represents the water release in the phreatic aquifer, while the smallscale ${\stackrel{\mathrm{\u0303}}{S}}_{y}$ in Eq. (12) denotes the dynamically changing water yield caused by the fluctuation of the water table. The upper boundary flux F_{top} in the phreatic flow equation (Eq. A3) is therefore corrected to
Differing from previous studies (Seo et al., 2007; Shen and Phanikumar, 2010; Xu et al., 2012), a scaleseparation strategy is employed in Eq. (13). The specific yields at two different scales are explicitly linked in F_{top}. The largescale properties in the groundwater model (MODFLOW) are thus fully maintained.
In this section, a range of 1D, 2D, 3D, and regional numerical test cases are presented. The 1D tests are benchmarked by the globally refined solutions from the HYDRUS1D code (Šimůnek et al., 2009). The 2D and 3D “truth” solutions are obtained from the fully 3D variably saturated flow (VSF; Thoms et al., 2006) model. At the regional scale, a synthetic case study suggested by Twarakavi et al. (2008) is reproduced. The codes are run on a personal computer with a 16 GB RAM and 3.6 GHz Intel Core (i34160). A maximal number of feedback iterations is set at 20. Soil parameters for the van Genuchten model (van Genuchten, 1980) are given in Table 1. The rootmeansquare error (RMSE) of the solution ψ at time t is given by
where ψ is the numerical solution of either the pressure head or water content, ψ^{ref} is the corresponding reference solution, and subscript i is the number of nodes, i=1, 2, …, N.
3.1 Case 1: rapidly changing atmospheric boundaries
The 1D case is used to investigate the benefit brought by switching Richards' equation in the unsaturated zone. A soil column is initialized with a hydrostatic watertable depth of 800 cm. That is, $h(t=\mathrm{0},z)=\mathrm{200}z$ cm, with z=0 at the bottom and z=1000 cm on the top. The lower boundary is set to nonflux to avoid the extra computational burden caused by variation of the groundwater model. Two scenarios from literature are reproduced with rapidly changing upper boundaries, as well as extreme flow interactions between the unsaturated and saturated zones.
Miller et al.'s problem (Miller et al., 1998) is reproduced in scenario 1. A dry sandy soil column (see soil no. 1 in Table 1) experiences a large constant flux infiltration at the soil surface of q_{top}=30 cm days^{−1}, which ceases at t=4 days.
In scenario 2, Hills et al.'s problem (Hills et al., 1989) is considered. The soils no. 2 and no. 3 from Table 1 are alternatively layered with a thickness of 20 cm within the first 80 cm of depth. Below 80 cm (z=0–920 cm) is soil no. 2, with a nonflux bottom boundary. The atmospheric upper boundary conditions, rainfall and evaporation, change rapidly with time (see Fig. 3), over 365 days.
The coupled unsaturated model is discretized into a fine grid with Δz=1 cm, while the saturated model is discretized into two layers with thickness of 500 cm. The impact of different numbers of feedback iteration, closure criteria, and different forms of the 1D Richards' equation are investigated. Solutions obtained from the HYDRUS1D model with Δz=1 cm and Δt=0.05 days are taken as the truth.
3.2 Case 2: dynamic groundwater flow
A 2D case is analyzed with sharp groundwater flow (see Fig. 4). To minimize the unsaturated lateral flow, the soil surface is set with a nonflux boundary. The bottom and lateral boundaries are also nonflux. Two pumping stresses are applied to the crosssectional field with $x\times z=\mathrm{5000}$ cm ×1000 cm. Well no. 1 is located at x=2500 cm, with a pumping screen at z=0–200 cm, while well no. 2 is at x=5000 cm, with a pumping screen of z=0–200 cm. Pumping rates for wells no. 1 and no. 2 respectively are 2×10^{4} and 1×10^{4} cm^{2} days^{−1} per unit width. The initial hydrostatic head of the cross section is ${h}_{\mathrm{0}}(x,z)=\mathrm{700}$ cm. Soil no. 4 in Table 1 fills the entire cross section. The total simulation lasts 50 days. For the coupled saturated submodel, as well as the reference model (VSF Thoms et al., 2006), the cross section is discretized horizontally into uniform segments with a width Δx=50 cm, while vertically (bottom up) refined into segments with thicknesses Δz=200 cm_{(×1)}, 100 cm_{(×2)}, 50 cm_{(×2)}, 25 cm_{(×2)}, 12.5 cm_{(×4)}, and 5 cm_{(×200)}, where the subscripts hereinafter (×N) are the numbers of discretized segments. The 1D soilwater models are discretized with a segmental thickness of Δz=1 cm. The fully 2D unsaturated–saturated solutions from the VSF model are taken as the truth.
3.3 Case 3: pumping and irrigation
Case 3 is used to investigate the efficiency and applicability of a quasi3D coupling model in comparison of the fully 3D approaches. A phreatic aquifer with $x\times y\times z=\mathrm{1000}$ m ×1000 m ×20 m is stressed by constant irrigation and pumping wells. The infiltration rate is 3 mm days^{−1} in (x, y) = (0–440 and 560–1000 m), while it is 5 mm days^{−1} in (x, y) = (560–1000 and 0–440 m). Screens for three of the pumping wells are located at (x, y, and z) = (220, 220, and 5–10 m), (500, 500, and 5–10 m), and (780, 780, and 5–10 m). The pumping rates are constant at 30 m^{3} days^{−1}. The initial hydrostatic head of the aquifer is 18 m. Around and below the aquifer are nonflux boundaries. The aquifer is horizontally discretized with $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{40}$ m for the coupled saturated model as well as for the VSF model to obtain the truth solution. The topdown thicknesses of the fully 3D grid are Δz=0.10 m_{(×30)}, 0.4 m_{(×5)}, 1 m_{(×5)}, and 2 m_{(×5)}. For the 1D soil columns, they are Δz=0.1 m_{(×30)} and 0.4 m_{(×5)}, which means that no soil column reaches the bottom. Different numbers of the subzones represented by soil columns, as well as their different geometries, are given in Fig. 5. The soil parameters for a sandy loam (soil no. 5) are given in Table 1. The total simulation lasts 60 days.
3.4 Case 4: synthetic regional case study
A hypothetical test case from literature (Niswonger et al., 2006; Prudic et al., 2004; Twarakavi et al., 2008) for a largescale simulation is reproduced here. The overall alluvial basin is divided into uniform grids with $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{1524}$ m. The coupled saturated model is conceptualized into a single layer. The initial head, as well as the elevations of land surface and bedrock, are presented in Fig. 6a, b, and c. The precipitation, evaporation, and pumping rates for 12 stress periods, each lasting 1∕12 of 365 days, are given in Table 2. The infiltration factors (see Fig. 6d) are used to approximate the spatial variability of precipitation. The initial head in the vadose zone is set with hydrostatic status. Twenty soil columns, coinciding with the subzones in Fig. 6d, are discretized separately with a range of gradually refined segments with a thickness (Δz) from 30.48 to 0.3048 cm (bottom up). A comparative analysis is conducted with the solutions obtained from the original HYDRUS package for MODFLOW (taken as HPM for short; Seo et al., 2007).
4.1 Reducing the complexity of a feedback coupling system
The numerical difficulty in a coupled unsaturated–saturated flow system originates from the nonlinearity of the soilwater models, heterogeneity of the parameters, and the variability of the hydrologic stresses (Krabbenhøft, 2007; Zha et al., 2017). In our work, the overall complexity of an iteratively coupled quasi3D model could be lowered by (1) taking full advantage of the h and θform REs to reduce the nonlinearity in the soilwater models and (2) smoothing the variability in the exchanged interfacial messages.
Two scenarios in case 1 were selected to address the first point. Sudden infiltration into a dry sandy soil and the rapidly altering atmospheric upper boundaries were tested to illustrate the importance of applying a switchingform RE for lowering the nonlinearity in the soilwater models. To evaluate the benefits brought by a switchingform RE, the numerical stability was first considered, as shown in Fig. 7. The coupled model in our work was tested with hform and switchingform REs. Compared with the HYDRUS1D model (also based on an hform RE), the switchingform method was numerically more robust, i.e., with larger minimal timestep sizes (Δt_{min}) and a lower computational cost, where a minimal timestep size 10^{−3} days was stable for convergence. Notably at the beginning of the sudden infiltration into a dry sandy soil, in Fig. 7a, the Δt_{min} for a switching method was 10^{−3} days, even at early infiltration times, while for the hform methods, including HYDRUS1D and the coupled hform method, Δt_{min} was constrained to 10^{−8} days before reaching a painstaking convergence. In Fig. 8, the soilwater content solution by the coupled switchingform method and the HYDRUS1D method (taken as the truth) were compared at depths of 0, 50, and 200 cm. The RMSEs of the soil moisture solutions (θ) at three different depths are respectively 0.0189, 0.0032, and 0.0013. To finish the calculation, the coupled switchingform RE method took 17 s, while it took 41 s for the HYDRUS code. When solving the same problem, the matrix equation set was solved 4903 times with the switching scheme, while it was solved 10 925 times for the HYDRUS1D code. The switched governing equations contribute toward cutting the computational cost by half for problems with rapidly changing upper boundary conditions. Here, the threshold for choosing an appropriate form of RE was nonsensitive to the numerical efficiency. A wide range of Se^{crit}ϵ [0.3, 0.9] was suggested according to substantial trialanderror tests.
Reducing the complexity of a coupling system can also be attained by smoothing the exchanged information in space and time. As suggested by Stoppelenburg et al. (2005), a timevarying specific yield calculated by the smallscale soilwater models, ${\stackrel{\mathrm{\u0303}}{S}}_{y}$ in Eq. (12), introduced significant variability to the largescale groundwater model, thus caused extra iterations. A largescale ${\stackrel{\mathrm{\u203e}}{S}}_{y}$ reduced the nonlinearity of the storage term in the groundwater equation. In case 1, using an ${\stackrel{\mathrm{\u203e}}{S}}_{y}$ of 0.1–0.2 in the groundwater model produced best numerical stability for the sandy soil with dramatically uprising water table. With a largescale ${\stackrel{\mathrm{\u203e}}{S}}_{y}$, the nonlinearity introduced by the smallscale soilwater models could be quickly smoothed, as shown in Eq. (12).
4.2 Multiscale water balance analysis
The traditional noniterative feedback coupling methods cannot maintain sound mass balance near the phreatic surface, especially for problems with drastic flow interactions.
One reason is that, to launch a new step of a submodel at either side of the phreatic interface, the noniterative feedback methods usually employed a predicted interfacial boundary without correction, which inevitably introduced coupling errors. In traditional noniterative methods (Seo et al., 2007; Xu et al., 2012), such shortcomings could be alleviated by refining the macrotimestep size (ΔT). However, the Dirichlet head predicted for the soil columns with a stepwise extension method (see Fig. 2a) was easy to implement but tended to suffer from significant coupling error. In this work, we proposed a linear extrapolation method for the lower boundary head prediction for the soilwater models (see Eq. A2). Here, we used niter to indicate the maximal number of feedback iteration. Compared with a traditional stepwise method, the solution obtained by a linear method, either iteratively (with niter =3) or noniteratively (niter =0), was easier for approaching the truth (see Fig. 9). Even with refined macrotimestep sizes (ΔT from 0.2 to 0.005 days), the stepwise method made a thorough effort to minimize the coupling errors. Notably, three feedback iterations (niter =3) were sufficient for reducing the coupling error significantly. Such a onedimensional case with constant upper boundary flux, avoiding interference from lateral fluxes, illustrated the importance of a temporal scalematching analysis for coupling the soilwater and groundwater models.
The other factor contributing to the coupling errors in the traditional method lies in neglecting the saturated lateral flux between adjacent soil columns (Seo et al., 2007; Stoppelenburg et al., 2005; Xu et al., 2012). In practical applications, the fluxes in and out of the saturated parts of the soil columns differ, which adds to the complexity of the coupling scheme. Although a strict water balance equation is established (Shen and Phanikumar, 2010), the concern centers on the spatial scalemismatching problem. That is, when the coarsegrid groundwater flow solutions are converted into the vertically distributed finescale source and sink terms for the soil columns, an extra downscaling approach is needed to ensure their accuracy. Here we carried out a multiscale water balance analysis above the phreatic surface. The finescale saturated lateral flows were thus excluded from Eq. (10). The benefits of the movingboundary approach can be seen in case 2, which produces significant saturated lateral flux. We carried out a comparative analysis against the traditional stationaryboundary methods (Seo et al., 2007; Xu et al., 2012). The 2D solution of VSF was taken as the truth. Figure 10 presents the effectiveness of the movingboundary method. Five stationary soil columns with three different lengths (L=1000, 500, and 300 cm) were compared with an adaptively moving soil column within the iterative feedback coupling scheme. The crosssectional RMSE of the phreatic surface and the head at the bottom layer (z=0) are presented in Fig. 10a and b. The soil columns with bottom nodes fixed deeply into the aquifer, instead of moving with the phreatic surface, introduced large coupling errors. This was caused by the nontrivial saturated lateral fluxes between the adjacent soil columns. With a traditional stationaryboundary method, such problems can be alleviated by avoiding large saturated lateral fluxes between the soil columns. However, for some spatiotemporally varying local events in a regional aquifer (e.g., pumping or flooding irrigation), such problems increased the burden for subzone partitioning. A movingboundary method, instead, was numerically more efficient in minimizing the size of the matrix equation and reducing the coupling errors.
4.3 Regulating the feedback iterations
In coupling two complicated modeling systems, a common agreement has been reached that; feedback coupling, either iteratively (Markstrom et al., 2008; Mehl and Hill, 2013; Stoppelenburg et al., 2005; Xie et al., 2012) or noniteratively (Seo et al., 2007; Shen and Phanikumar, 2010; Xu et al., 2012), is numerically more rigorous than the oneway coupling scheme. The main difference between the above two methods lies in the ability to conserve mass within a single step for backandforth information exchange. In an iterative method, the head and/or flux boundaries are iteratively exchanged. There is a cost–benefit tradeoff for obtaining higher numerical efficiency.
During the late stages of the recharge in scenario 1 of case 1, the groundwater table rises quickly, which increases the burden on the coupling scheme. In our work, feedback iteration was conducted to eliminate the coupling error within the backandforth boundary exchange. To investigate how the feedback iteration influences the numerical accuracy as well as computational cost, solutions were compared with different closure criteria, instead of different maximal numbers of feedback iterations. For this purpose, scenario 1 in case 1 is tested with a range of closure criteria indicated by closure =0.001, 0.01, 0.1, 1, 5, and 20. Specifically, closure =20 (i.e., ε_{H}=20 cm) is too large to regulate any feedback iteration and is thus replaced by “noniterative”. The ε_{F}, indicating the closure of the Neumann boundary feedback iteration, is usually related to the phreatic Darcian flux. To avoid its impact on the discussion below, we assume that ${\mathit{\epsilon}}_{\mathrm{F}}=+\mathrm{\infty}$, which means no regulation from the flux boundary exchange. Due to less dynamic in the groundwater submodel, the empirical relaxation factors were both set by 1.0 to have a straightforward update of the interfacial boundaries, i.e., z_{t} and F_{top}.
When the wetting front approached the phreatic surface (at t=2.4 days), the number of feedback iterations increased dramatically (see Fig. 11a). This was caused by the intensive rise of the water table within each macrotime step ΔT. The head and/or flux interfacial boundaries were thus not easy to approximate the truth. With several attempts to exchange the head and/or flux boundaries, the head solution was effectively drawn towards the truth, (see Fig. 11b). With closure < 2, i.e., ε_{H} < 2 cm, the coupling errors were significantly reduced (see Fig. 11c). The cost–benefit curve, which was quantified by the number of feedback iteration and the RMSE, was indicative of problems at larger scales and higher dimensionalities.
4.4 Parsimonious decisionmaking
The feedback coupling schemes, either iteratively or noniteratively, increase the degree of freedom of the users to manage the submodels with different governing equations, numerical algorithms, and the heterogeneities in parameters and variabilities in hydrologic stresses. For practical purposes, a significant concern is how to efficiently handle the complicated and scaledisparate systems.
For problems with rapid changes in groundwater flows, as in case 2, the hydraulic gradient at the phreatic surface was large. Using a single soil column for such a complex situation introduced significant coupling errors at the water table, (see Fig. 12a). Although more subzones portioned means higher accuracy for the coupling method, five or more soil columns were adequate enough in approaching the truth. Furthermore, for the saturated nodes deep in the aquifer, such differences in coupling errors were of minor influence (see Fig. 12b).
In case 3, a simple pumped and irrigated region was simulated with different numbers of soil columns. A range of tests with total numbers of 16, 12, 9, 5, and 3 soil columns were carried out to obtain a cost–benefit curve, shown in Fig. 13c. When partitioning the vadose zone into more than 12 soil columns, there was a slight reduction in solution errors (RMSE) and a significant increase in the computational cost caused by solving more 1D soilwater models. Although parallelled computation could further reduce the numerical cost, representing the vadose zone with three sequentially calculated soilwater models achieved acceptable accuracy, as presented in Fig. 13a and b. The computational cost for obtaining the fully 3D solution with VSF was 15.561 s, which was more than 11 times larger than an iterative feedback coupling method with soilwater models sequentially solved. Problems in more complicated realworld situations can thus be simplified to achieve higher numerical efficiency.
4.5 Regional application
Prudic et al.'s problem was originally designed to validate a streamflow routing package (Prudic et al., 2004). Stressed by soil surface infiltration, pumping wells, and the general head boundary, the synthetic case was used to evaluate several unsaturated flow packages for MODFLOW (Twarakavi et al., 2008). Based on their studies, in case 4, we compared the developed iterative feedback coupling method with HPM. The hydraulic conductivity, as well as its heterogeneity, was forced to be consistent within the saturated and unsaturated zones, which is different from the case in Twarakavi et al. (2008). Figure 14a gives the contours for the final phreatic head solutions, indicating a good match of the phreatic surface with the HYDRUS package. Figure 14b–e present the absolute head difference of the method developed here and the HYDRUS package at the end of stress periods 3, 6, 9, and 12. The dark color blocks indicated the largest difference in the head solution. According to Fig. 6d, the saturated grid cells controlled by the soil columns of no. 3, no. 9, no. 10, no. 15, and no. 19 were suffering from the largest deviation, although with the same horizontal partitioning of the unsaturated zone. The strict iteratively twoway coupling contributes to such an accuracy improvement.
For unsaturated–saturated flow situations, the vadose zone flow is important. Figure 15 presents the water content profiles at subzones 1, 3, 5, 7, and 9 as examples. The solution obtained from the developed model matched well with the original HPM solution. For practical purposes, the manually controlled stress periods for the unsaturated submodels are no longer constrained. In our method, the soilwater models run at disparate numerical scales, which makes it possible to handle daily or hourly observed information rather than a stress period lasting 2 or more days in traditional coupled models.
Fully 3D numerical models are available but are numerically expensive for simulating the regional unsaturated–saturated flow. The quasi3D method presented here, in contrast, with horizontally disconnected adjacent unsaturated nodes, significantly reduces the dimensionality and complexity of the problem. Such a simplification brings about computational cost saving and flexibility for better manipulation of the submodels. However, the nonlinearity in the soilwater retention curve, as well as the variability in realistic boundary stresses of the vadose and saturated zones, usually results in the scalemismatching problems when attempting numerical coupling. In this work, the soilwater and groundwater models were coupled with an iterative feedback (twoway) coupling scheme. Three concerns about the multiscale water balance at the phreatic interface are addressed using a range of numerical cases in multiple dimensionalities. We conclude with the following:

A new HYDRUS package for MODFLOW was developed by switching the θ and h forms of Richards' equation (RE) at each numerical node. The switching RE circumvents the disadvantages of the h and θform REs to achieve higher numerical stability and computational efficiency. The onedimensional switching RE was employed to simulate the rapid infiltration into a dry sandy soil, and the swiftly altering atmospheric upper boundaries in a layered soil column. Compared with the hform RE, the switching RE used a 10^{5} times larger minimal timestep size (Δt_{min}) and conserved mass better. Lowering the nonlinearity of soilwater models with this switching scheme was promising for coupling complex flow modeling systems at the regional scale.

Stringent multiscale water balance analysis at the water table was conducted to handle scalemismatching problems and to smooth the information delivered back and forth across the interface. In our work, the errors originating from inadequate phreatic boundary predictions were reduced firstly by a linear extrapolation method and then by an iterative feedback. Compared with the traditional stepwise extension method, the linear extrapolation significantly reduced the coupling errors caused by scale mismatching. For problems with severe soilwater and groundwater interactions, the coupling errors were significantly reduced by using an iterative feedback coupling scheme. The multiscale water balance analysis mathematically maintained numerical stabilities in the submodels at disparate scales.

When a moving phreatic boundary was assigned to the soil columns during the phreatic water balance analysis, it avoided the coupling errors originating from the saturated lateral fluxes. In practical applications for regional problems, the fluxes into and out of the saturated parts of the soil columns differed, which added to the complexity and phreatic water balance error of the coupling scheme. With a moving Dirichlet lower boundary, the saturated regions of the soilwater models were removed. The coupling error was significantly reduced for problems with major groundwater flow. Extra cost saving was achieved by minimizing the matrix sizes of the soilwater models.
Future investigation will focus on regional solute transport modeling based on the developed coupling scheme. Surface flow models, as well as the crop models, which appear to be less nonlinear than the subsurface models, will be coupled in an objectoriented modeling system. The RS and GISbased data class can then be used to handle more complicated largescale problems.
All the data and codes used in this study can be requested by email to the corresponding author Yuanyuan Zha at zhayuan87@gmail.com.
A1 The moving Dirichlet lower boundary
The bottom node of a soil column is adaptively located at the phreatic surface, which makes it an areaaveraged moving Dirichlet boundary:
where z_{t}(T) (L) is the elevation of the water table; Π is the control domain of a soil column; H(T) (L) is potentiometric head solution, as well as the elevation of the phreatic surface, which is obtained by solving the groundwater model; and s is the horizontal area.
To simulate the multiscale flow process within a macrotime step $\mathrm{\Delta}{T}^{J+\mathrm{1}}={T}^{J+\mathrm{1}}{T}^{J}$, the lower boundary head of a soil column is temporally predicted either by stepwise extension of z_{t}(T^{J}) (Seo et al., 2007; Shen and Phanikumar, 2010; Xu et al., 2012) or by linear extrapolation from z_{t}(T^{J+1}) and z_{t}(T^{J}). In Fig. 2a, the stepwise extension method (${z}_{t}^{\prime}\left({T}^{J}\right)$) potentially causes large deviation from the truth. In our study, the linear extrapolation is resorted to for reducing the coupling errors and accelerating the convergence of the feedback iteration. The smallscale lower boundary head at time t (${T}^{J}\phantom{\rule{0.125em}{0ex}}\mathit{<}\phantom{\rule{0.125em}{0ex}}t\le {T}^{J+\mathrm{1}}$) is given by
A2 The Neumann upper boundary
The moving Dirichlet boundary introduces the need for water balance of a moving balancing domain above the water table (see Fig. 2b), which is bounded by a specific elevation above the phreatic surface, z_{s} (L), and the dynamically changing phreatic surface, z_{t}(t) (L).
Assuming that the activated top layer in a threedimensional groundwater model is conceptualized into a phreatic aquifer, the governing equation for this layer is given by
where $\stackrel{\mathrm{\u203e}}{M}$ (L) is the thickness of the phreatic layer, which is numerically defined as the layer below the vadose zone, $\stackrel{\mathrm{\u203e}}{M}={z}_{\mathrm{s}}{z}_{\mathrm{0}}$; z_{0} is the bottom elevation of the top phreatic layer, z_{0 ≪}z_{s}; F_{top} (L T^{−1}) is the groundwater recharge into the activated top layer of the phreatic aquifer, ${F}_{\mathrm{top}}={\left(K\cdot \partial H/\partial z\right)}_{z={z}_{\mathrm{s}}}$; and F_{base} is the leakage into an underlying numerical layer, ${F}_{\mathrm{base}}={\left(K\cdot \partial H/\partial z\right)}_{z={z}_{\mathrm{0}}}$ (positive downward, like F_{top}). The longterm regionalscale parameter indicating the water yield caused by the fluctuation of the water table (Nachabe, 2002), ${\stackrel{\mathrm{\u203e}}{S}}_{y}$, (–), is calculated by
where V_{w} (L^{3}) is the amount of water release caused by fluctuation of the phreatic surface (ΔH (L)), and A (L^{2}) is the area of interest.
A3 The relaxed iterative feedback coupling
The relaxed feedback iteration method (Funaro et al., 1988; Mehl and Hill, 2013) is used to improve the convergence of head and/or flux at the phreatic surface (see Fig. A1). The Dirichlet lower boundary head for the soil columns, z_{t}, as well as the Neumann upper boundary flux for the phreatic surface, F_{top}, are updated within each iterative step (niter):
where superscripts old (or new) indicate the previous (or newly calculated) head and/or flux boundaries at the coupling interface, and λ_{h} and λ_{f} are the empirical relaxation factors for head and/or flux boundaries respectively. Their values are suggested to be within $\mathrm{0}<\mathit{\lambda}\le \mathrm{1}$. The iteration ends when agreements are reached at
where ε_{H} (L) and ε_{F} (L T^{−1}) are residuals for the feedback iteration of the interfacial head and flux.
JZ, YZ, and JY developed the new package for soilwater movement based on a switching Richards' equation; JZ and YZ developed the coupling methods for efficiently joining the submodels. JZ, YZ, JY, and LS made nonnegligible efforts in preparing the paper.
The authors declare that they have no conflict of interest.
This work was funded by the Chinese National Natural Science Foundation
(no. 51479143 and 51609173). The authors would thank Ian White and
Wenzhi Zeng for their laborious revisions and helpful suggestions for the
paper as well as the editorial team for the HESS Journal.
Edited by: Bob Su
Reviewed by: two anonymous
referees
Bailey, R. T., Morway, E. D., Niswonger, R. G., and Gates, T. K.: Modeling variably saturated multispecies reactive groundwater solute transport with MODFLOWUZF and RT3D, Groundwater, 51, 752–761, https://doi.org/10.1111/j.17456584.2012.01009.x, 2013.
Celia, M. A., Bouloutas, E. T., and Zarba, R. L.: A general massconservative numerical solution for the unsaturated flow equation, Water Resour. Res., 26, 1483–1496, https://doi.org/10.1029/WR026i007p01483, 1990.
Crevoisier, D., Chanzy, A., and Voltz, M.: Evaluation of the Ross fast solution of Richards' equation in unfavourable conditions for standard finite element methods, Adv. Water Resour., 32, 936–947, https://doi.org/10.1016/j.advwatres.2009.03.008, 2009.
Dettmann, U. and Bechtold, M.: Onedimensional expression to calculate specific yield for shallow groundwater systems with microrelief, Hydrol. Process., 30, 334–340, https://doi.org/10.1002/hyp.10637, 2016.
Diersch, H. J. G. and Perrochet, P.: On the primary variable switching technique for simulating unsaturatedsaturated flows, Adv. Water Resour., 23, 271–301, https://doi.org/10.1016/S03091708(98)000578, 1999.
Downer, C. W. and Ogden, F. L.: Appropriate vertical discretization of Richards' equation for twodimensional watershedscale modelling, Hydrol. Process., 18, 1–22, https://doi.org/10.1002/hyp.1306, 2004.
Edwards, M. G.: Elimination of Adaptive Grid Interface Errors in the Discrete Cell Centered Pressure Equation, J. Comput. Phys., 126, 356–372, https://doi.org/10.1006/jcph.1996.0143, 1996.
Forsyth, P. A., Wu, Y. S., and Pruess, K.: Robust numerical methods for saturatedunsaturated flow with dry initial conditions in heterogeneous media, Adv. Water Resour., 18, 25–38, https://doi.org/10.1016/03091708(95)00020J, 1995.
Funaro, D., Quarteroni, A., and Zanolli, P.: An iterative procedure with interface relaxation for domain decomposition methods, Siam J. Numer. Anal., 25, 1213–1236, https://doi.org/10.1137/0725069, 1988.
Grygoruk, M., Batelaan, O., MirosławŚwiltek, D., Szatyłowicz, J., and Okruszko, T.: Evapotranspiration of bush encroachments on a temperate mire meadow – A nonlinear function of landscape composition and groundwater flow, Ecol. Eng., 73, 598–609, https://doi.org/10.1016/j.ecoleng.2014.09.041, 2014.
Gunduz, O. and Aral, M. M.: River networks and groundwater flow: A simultaneous solution of a coupled system, J. Hydrol., 301, 216–234, https://doi.org/10.1016/j.jhydrol.2004.06.034, 2005.
Hao, X., Zhang, R., and Kravchenko, A.: A massconservative switching method for simulating saturatedunsaturated flow, J. Hydrol., 301, 254–265, https://doi.org/10.1016/j.jhydrol.2005.01.019, 2005.
Harbaugh, A. W.: MODFLOW2005, the U.S. Geological Survey modular groundwater model: The groundwater flow process, US Department of the Interior, US Geological Survey Reston, VA, USA, 2005.
Hills, R. G., Porro, I., Hudson, D. B., and Wierenga, P. J.: Modeling onedimensional infiltration into very dry soils: 1. Model development and evaluation, Water Resour. Res., 25, 1259–1269, https://doi.org/10.1029/WR025i006p01259, 1989.
Krabbenhøft, K.: An alternative to primary variable switching in saturatedunsaturated flow computations, Adv. Water Resour., 30, 483–492, https://doi.org/10.1016/j.advwatres.2006.04.009, 2007.
Kumar, M., Duffy, C. J., and Salvage, K. M.: A SecondOrder Accurate, Finite VolumeBased, Integrated Hydrologic Modeling (FIHM) Framework for Simulation of Surface and Subsurface Flow, Vadose Zone J., 8, 873–890, https://doi.org/10.2136/vzj2009.0014, 2009.
Kuznetsov, M., Yakirevich, A., Pachepsky, Y. A., Sorek, S., and Weisbrod, N.: Quasi 3D modeling of water flow in vadose zone and groundwater, J. Hydrol., 450–451, 140–149, https://doi.org/10.1016/j.jhydrol.2012.05.025, 2012.
Langevin, C. D., Hughes, J. D., Banta, E. R., Provost, A. M., Niswonger, R. G., and Panday, S.: MODFLOW 6 Groundwater Flow (GWF) Model Beta version 0.9.03, U.S. Geol. Surv. Provisional Softw. Release, https://doi.org/10.5066/F76Q1VQV, 2017.
Leake, S. A. and Claar, D. V: Procedures and computer programs for telescopic mesh refinement using MODFLOW, Citeseer, available at: http://az.water.usgs.gov/MODTMR/tmr.html (last access: 29 January 2019), 1999.
Lin, L., Yang, J.Z., Zhang, B., and Zhu, Y.: A simplified numerical model of 3D groundwater and solute transport at large scale area, J. Hydrodyn. Ser. B, 22, 319–328, https://doi.org/10.1016/S10016058(09)600615, 2010.
Liu, Z., Zha, Y., Yang, W., Kuo, Y., and Yang, J.: LargeScale Modeling of Unsaturated Flow by a Stochastic Perturbation Approach, Vadose Zone J., 15, https://doi.org/10.2136/vzj2015.07.0103, 2016.
Markstrom, S. L., Niswonger, R. G., Regan, R. S., Prudic, D. E., and Barlow, P. M.: GSFLOW – Coupled GroundWater and SurfaceWater Flow Model Based on the Integration of the PrecipitationRunoff Modeling System (PRMS) and the Modular GroundWater Flow Model (MODFLOW2005), Geological Survey (US), available at: https://pubs.usgs.gov/tm/tm6d1/ (last access: 29 January 2019), 2008.
Maxwell, R. M. and Miller, N. L.: Development of a coupled land surface and groundwater model, J. Hydrometeorol., 6, 233–247, https://doi.org/10.1175/JHM422.1, 2005.
McDonald, M. G. and Harbaugh, A. W.: A modular threedimensional finitedifference groundwater flow model, available at: http://pubs.er.usgs.gov/publication/twri06A1 (last access: 29 January 2019), 1988.
Mehl, S. and Hill, M. C.: Threedimensional local grid refinement for blockcentered finitedifference groundwater models using iteratively coupled shared nodes: a new method of interpolation and analysis of errors, Adv. Water Resour., 27, 899–912, https://doi.org/10.1016/j.advwatres.2004.06.004, 2004.
Mehl, S. and Hill, M. C.: MODFLOW–LGR – Documentation of ghost node local grid refinement (LGR2) for multiple areas and the boundary flow and head (BFH2) package, 2013th ed., available at: https://pubs.usgs.gov/tm/6a44 (last access: 29 January 2019), 2013.
Miller, C. T., Williams, G. A., Kelley, C. T., and Tocci, M. D.: Robust solution of Richards' equation for nonuniform porous media, Water Resour. Res., 34, 2599–2610, https://doi.org/10.1029/98WR01673, 1998.
Miller, C. T., Abhishek, C., and Farthing, M. W.: A spatially and temporally adaptive solution of Richards' equation, Adv. Water Resour., 29, 525–545, https://doi.org/10.1016/j.advwatres.2005.06.008, 2006.
Nachabe, M. H.: Analytical expressions for transient specific yield and shallow water table drainage, Water Resour. Res., 38, 111–117, https://doi.org/10.1029/2001WR001071, 2002.
Niswonger, R. G., Prudic, D. E., and Regan, S. R.: Documentation of the UnsaturatedZone Flow (UZF1) Package for Modeling Unsaturated Flow Between the Land Surface and the Water Table with MODFLOW2005, US Department of the Interior, US Geological Survey, available at: https://pubs.er.usgs.gov/publication/tm6A19 (last access: 29 January 2019), 2006.
Niswonger, R. G., Panday, S., and Ibaraki, M.: MODFLOWNWT, a Newton formulation for MODFLOW2005, US Geol. Surv. Tech. Methods, 6, A37, available from: https://pubs.usgs.gov/tm/tm6a37 (last access: 29 January 2019), 2011.
Panday, S. and Huyakorn, P. S.: A fully coupled physicallybased spatiallydistributed model for evaluating surface/subsurface flow, Adv. Water Resour., 27, 361–382, https://doi.org/10.1016/j.advwatres.2004.02.016, 2004.
Panday, S., Langevin, C. D., Niswonger, R. G., Ibaraki, M., and Hughes, J. D.: MODFLOW–USG version 1: An unstructured grid version of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finitedifference formulation, US Geological Survey, https://doi.org/10.3133/tm6A45, 2013.
Paulus, R., Dewals, B. J., Erpicum, S., Pirotton, M., and Archambeau, P.: Innovative modelling of 3D unsaturated flow in porous media by coupling independent models for vertical and lateral flows, J. Comput. Appl. Math., 246, 38–51, https://doi.org/10.1016/j.cam.2012.07.032, 2013.
Prudic, D. E., Konikow, L. F., and Banta, E. R.: A New StreamflowRouting (SFR1) Package to Simulate StreamAquifer Interaction with MODFLOW2000, available at: http://pubs.er.usgs.gov/publication/ofr20041042 (last access: 29 January 2019), 2004.
Richards, L. A.: Capillary conduction of liquids through porous mediums, J. Appl. Phys., 1, 318–333, https://doi.org/10.1063/1.1745010, 1931.
Rybak, I., Magiera, J., Helmig, R., and Rohde, C.: Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems, Comput. Geosci., 19, 299–309, https://doi.org/10.1007/s1059601594698, 2015.
Seo, H. S., Šimůnek, J., and Poeter, E. P.: Documentation of the HYDRUS Package for MODFLOW2000, the U.S. Geological Survey Modular GroundWater Model, GWMI 200701, Int. Ground Water Modeling Center, Colorado School of Mines, Golden, CO, 96 pp., 2007.
Shen, C. and Phanikumar, M. S.: A processbased, distributed hydrologic model based on a largescale method for surfacesubsurface coupling, Adv. Water Resour., 33, 1524–1541, https://doi.org/10.1016/j.advwatres.2010.09.002, 2010.
Šimůnek, J., Sejna, M., Saito, H., Sakai, M., and van Genuchten, M. T.: The HYDRUS1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variablysaturated media, Version 4.08, 2009.
Šimůnek, J., Van Genuchten, M. T., and Šejna, M.: Recent Developments and Applications of the HYDRUS Computer Software Packages, Vadose Zone J., 15, https://doi.org/10.2136/vzj2016.04.0033, 2016.
Stoppelenburg, F. J., Kovar, K., Pastoors, M. J. H., and Tiktak, A.: Modelling the interactions between transient saturated and unsaturated groundwater flow. Offline coupling of LGM and SWAP, RIVM Rep., 500026001, 70, 2005.
Thoms, R. B., Johnson, R. L., and Healy, R. W.: User's guide to the variably saturated flow (VSF) process to MODFLOW, available at: http://pubs.er.usgs.gov/publication/tm6A18 (last access: 29 January 2019), 2006.
Twarakavi, N. K. C., Šimůnek, J., and Seo, S.: Evaluating Interactions between Groundwater and Vadose Zone Using the HYDRUSBased Flow Package for MODFLOW, Vadose Zone J., 7, 757–768, https://doi.org/10.2136/vzj2007.0082, 2008.
van Dam, J. C., Groenendijk, P., Hendriks, R. F. A., and Kroes, J. G.: Advances of Modeling Water Flow in Variably Saturated Soils with SWAP, Vadose Zone J., 7, 640–653, https://doi.org/10.2136/vzj2007.0060, 2008.
van Genuchten, M. T.: A Closedform Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980.
Vogel, H.J. and Ippisch, O.: Estimation of a Critical Spatial Discretization Limit for Solving Richards' Equation at Large Scales, Vadose Zone J., 7, 112–114, https://doi.org/10.2136/vzj2006.0182, 2008.
Warrick, A. W.: Numerical approximations of darcian flow through unsaturated soil, Water Resour. Res., 27, 1215–1222, https://doi.org/10.1029/91WR00093, 1991.
Xie, Z., Di, Z., Luo, Z., and Ma, Q.: A QuasiThreeDimensional Variably Saturated Groundwater Flow Model for Climate Modeling, J. Hydrometeorol., 13, 27–46, https://doi.org/10.1175/JHMD1005019.1, 2012.
Xu, X., Huang, G., Zhan, H., Qu, Z., and Huang, Q.: Integration of SWAP and MODFLOW2000 for modeling groundwater dynamics in shallow water table areas, J. Hydrol., 412–413, 170–181, https://doi.org/10.1016/j.jhydrol.2011.07.002, 2012.
Yakirevich, A., Borisov, V., and Sorek, S.: A quasi threedimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi twodimensional case, Adv. Water Resour., 21, 679–689, https://doi.org/10.1016/S03091708(97)000316, 1998.
Zadeh, K. S.: A massconservative switching algorithm for modeling fluid flow in variably saturated porous media, J. Comput. Phys, 230, 664–679, https://doi.org/10.1016/j.jcp.2010.10.011, 2011.
Zeng, J., Zha, Y., Zhang, Y., Shi, L., Zhu, Y., and Yang, J.: On the submodel errors of a generalized oneway coupling scheme for linking models at different scales, Adv. Water Resour., 109, 69–83, https://doi.org/10.1016/j.advwatres.2017.09.005, 2017.
Zeng, J., Zha, Y., and Yang, J.: Switching the Richards' equation for modeling soil water movement under unfavorable conditions, J. Hydrol., 563, 942–949, https://doi.org/10.1016/j.jhydrol.2018.06.069, 2018.
Zha, Y., Shi, L., Ye, M., and Yang, J.: A generalized Ross method for two and threedimensional variably saturated flow, Adv. Water Resour., 54, 67–77, https://doi.org/10.1016/j.advwatres.2013.01.002, 2013a.
Zha, Y., Yang, J., Shi, L., and Song, X.: Simulating OneDimensional Unsaturated Flow in Heterogeneous Soils with Water ContentBased Richards Equation, Vadose Zone J., 12, https://doi.org/10.2136/vzj2012.0109, 2013b.
Zha, Y., Yang, J., Yin, L., Zhang, Y., Zeng, W., and Shi, L.: A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil, J. Hydrol., 551, 56–69, https://doi.org/10.1016/j.jhydrol.2017.05.053, 2017.
Zhu, Y., Shi, L., Lin, L., Yang, J., and Ye, M.: A fully coupled numerical modeling for regional unsaturatedsaturated water flow, J. Hydrol., 475, 188–203, https://doi.org/10.1016/j.jhydrol.2012.09.048, 2012.