Capturing soil-water and groundwater interactions with an iterative feedback coupling scheme: New HYDRUS package for MODFLOW

Authors deal with the scale-mismatching problem when coupling the soil-water and groundwater models. A range of numerical cases were employed to address three concerns arose using the iterative feedback coupling method. This work successfully present its advantages in reducing computational cost, coupling errors, and maintaining the numerical stabilities of the sub-models at disparate scales. The method presented here will be promising in the application of large scale problems. This paper is of significant contribution to scientific progress regarding the coupling of soil-water and groundwater systems. I am interested in this paper and recommend some minor revisions before its acceptance for publication.


Introduction
Numerical modeling of the soil-water 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 soil-water and saturated groundwater flows with similar properties are usually integrated into a whole modeling system.Although physically consistent and numerically rigorous, methods involving the 3-D 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 soil-water and groundwater interactions at regional scales.
The original RE, also known as the mixed-form 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 mixed-form RE, either h or θ , or a switching of both, is assigned as the primary variable.The h-form RE is widely employed for unsaturatedsaturated flow simulation, especially in heterogeneous soils, such as the HYDRUS package (Šimůnek et al., 2016).Signif-Published by Copernicus Publications on behalf of the European Geosciences Union.
J. Zeng et al.: Capturing soil-water and groundwater interactions icant 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 Celia-format h-form 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 quasi-3-D 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 quasi-3-D methods are categorized into (1) the fully coupling scheme, which simultaneously builds the nodal hydraulic connections of sub-models at both sides and implicitly solves the assembled matrices; (2) the one-way coupling scheme, which delivers the soilwater model solutions onto the groundwater model without feedback mechanism; and (3) the feedback (or two-way) 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 non-convergence for the iterative solvers (Edwards, 1996).Owing to high nonlinearity in the soil-water sub-models, the assembled matrices can only be solved with unified small time steps, which adds to the computational expense.The one-way 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 soil-water and groundwater sub-models can be built with governing equations and numerical schemes at different scales.For flow processes with multi-scale components, such as boundary geometries, parameter heterogeneities, and hydrologic stresses, the scale-separation strategy can be implemented easily.Although the feedback coupling method is numerically more rigorous than a one-way 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 non-iterative 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 time-step 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 time-stepping 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 scale-mismatching 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 large-scale parameter, while for soil-water models (Niswonger et al., 2006;Šimůnek et al., 2009;Thoms et al., 2006), the small-scale phreatic water release is influenced by the water-table depth and the unsaturated soil moisture profile (Dettmann and Bechtold, 2016;Nachabe, 2002).Delivering small-scale solutions of the soilwater models onto the large-scale 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 sub-models causes significant coupling errors and instability.A multi-scale 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 2-D groundwater model usually require additional efforts to build water budget equations in each subdivision represented by the soil columns.A moving-boundary strategy helps to avoid the saturated lateral flow in the groundwater body.
In this work, the hand θ -form of the 1-D RE are switched at the equation level to obtain a new HYDRUS package.To handle three of the aforementioned concerns, a multi-scale 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 dis-Hydrol.Earth Syst.Sci., 23, 637-655, 2019 www.hydrol-earth-syst-sci.net/23/637/2019/ parate scales.The saturated lateral fluxes between the soil columns are fully removed from the interfacial water balance equation, making it a moving-interface coupling framework.
In this paper, the governing equations at different scales, the multi-scale 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.

Methodology
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 soil-water sub-models 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 soil-water and groundwater models at independent scales (Sect.2.4).As for the second concern, a multi-scale 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 soil-water 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.

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 soil-water capacity (C = ∂θ/∂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 1-D RE is represented by a switchable format, where ψ is the primary variable.For an h-form RE, ψ = h, Ĉ = C, and K = K, while for a θ -form RE, ψ = θ , Ĉ = 1, and K = D; where D (L 2 T −1 ) is the hydraulic diffusivity, D = K/C.

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 h-form 1-D RE for temporal approximation.Both forms of RE are handled with a temporally backward finite difference discretization (Zha et al., 2013b(Zha et al., , 2017)).Each sub-model 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 + 1/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 + 1/2) is the number of node (or element), and z i+1/2 is the length of the element i + 1/2, z i+1/2 = (z i+1 − z i ).When a soil interface exists at node i, for example, the soil moisture contents in elements i − 1/2 and i + 1/2 are discontinuous at node i, thus dissatisfying the θ-form RE.To address this problem, the correction term ε j +1,k i+1/2 , suggested by Zha et al. (2013b), is employed to handle the heterogeneous interface at nodes i and i + 1, where variables (scalar) are the continuously distributed ψ within element i + 1/2, i.e., between the vertices i and i + 1.

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 soil-water capacity in the storage term, the θ -form RE is more robust than the h-form RE, especially when dealing with rapidly changing atmospheric boundary conditions (Zeng et al., 2018).In our work, the hand θ -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 hand θ -form REs are stable and efficient.When Se ≥ Se crit , the soil moisture is closer to saturation, so the h-form RE is chosen as the governing equation; otherwise, when it undergoes dry soil condition, the θ -form RE is preferred.
For element i + 1/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 ψ i+1 = h i+1 as an example, the nodal fluxes calculated by Eq. ( 5) for different forms of RE have to be carefully handled by substituting . The resulting equivalent nodal fluxes q h i+1/2 and q θ i+1/2 are then weighted to obtain an approximation by where ω is the weighting factor, 0 ≤ ω ≤ 1.In our work, ω = 0.5 is applied to implicitly maintain the unknown variables of both h j +1,k+1 i+1 and θ j +1,k+1 i . 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 h-form 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 switching-RE approach is incorporated into a new HYDRUS package.

Iterative feedback coupling scheme
The Dirichlet and Neumann boundaries are iteratively transferred across the phreatic interface.The groundwater head solution serves as the head-specified lower boundary of the soil columns, while the unsaturated solution is converted into the flux-specified upper boundary of the groundwater model.Due to moderate variation of the groundwater flow, the predicted water-table solution is usually adopted in advance as the Dirichlet lower boundary of the fine-scale 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 coarse-scale groundwater flow model.Section A1 provides the method for a moving Dirichlet lower boundary, while Sect.A2 presents the Neumann upper boundary for the 3-D groundwater model.In Sect.A3, the relaxed iterative feedback coupling scheme is used to solve the unsaturatedsaturated sub-models at two sides of the coupling interface.

Multi-scale 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 time-splitting strategy (see Fig. 1) is adopted to separate sub-models at different scales.That is, the soil-water models are established by z = 10 −3 -10 0 m and t = 10 −5 -10 0 days, while for the saturated model, the grid sizes are x = 10 0 -10 3 m, and time-step 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 sub-model 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 multi-scale water balance analysis is conducted within each step of the large-scale saturated flow model.At small spatial and temporal scales, e.g., within a macro-time step T = T J +1 − T J and at a local area of interest (with thickness of M = z s − z 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 one-dimensional expression (Dettmann and Bech-Figure 1. Schematic of the space-and time-splitting strategy for coupling models at two independent scales.For a groundwater model, spatial discretization is expected to be large ( x = 10 0 -10 3 m), while for soil-water models, it is small ( x = 10 −3 -10 0 m).Multiple levels of temporal discretization are common for regional problems.For the groundwater model, the stress periods (SP) and macro-time-step sizes ( T ) appear in months and days (10 0 -10 1 days).For soil-water models, the time-step sizes are about 10 −5 -10 0 days.told, 2016), where w (L) is the amount of unsaturated water in the moving balancing domain (see Fig. 2b), w(t) = z s z t (t) θ(t, z)dz;  2b) during a small-scale 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 ( x,y,z∈ dxdydz 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 top = q top (t)dt.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 bot = T J +1 T J q bot (t)dt.ε l (L) is the cumulative lateral input water into the moving balancing domain, where N is the number of time steps for the small-scale soilwater model within a macro-time step T , and ε l is the nontrivial saturated later flux produced by a stationary-boundary method (Seo et al., 2007;Xu et al., 2012).By taking R top as the specific recharge at z s , the small-scale specific yield S y is derived from Eqs. ( 8) and ( 10) as Supposing that z t is linearly fluctuating in time, i.e., z t = a • t+ b (where a and b are constants), we get the watertable change during a small-scale step (dt) by dz t = a • dt; thus, ε l = o(dt 2 ), which means that linearly refining the local time-step size (dt) in the soil-water model brings about at least quadratic approximation of ε l towards zero.Thus ε l can be neglected from the small-scale mass balance analysis.
In the developed model, the large-scale specific yield, S y in Eq. (A3), represents the water release in the phreatic aquifer, while the small-scale 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 scale-separation strategy is employed in Eq. ( 13).The specific yields at two different scales are explicitly linked in F top .The large-scale properties in the groundwater model (MODFLOW) are thus fully maintained.

Numerical experiments
In this section, a range of 1-D, 2-D, 3-D, and regional numerical test cases are presented.The 1-D tests are benchmarked by the globally refined solutions from the HYDRUS-1D code (Šimůnek et al., 2009).The 2-D and 3-D "truth"  solutions are obtained from the fully 3-D 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 (i3-4160).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 root-mean-square 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 .

Case 1: rapidly changing atmospheric boundaries
The 1-D case is used to investigate the benefit brought by switching Richards' equation in the unsaturated zone.A soil column is initialized with a hydrostatic water-table depth of 800 cm.That is, h(t = 0, z) = 200 − z cm, with z = 0 at the bottom and z = 1000 cm on the top.The lower boundary is set to non-flux 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 non-flux 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 1-D Richards' equation are investigated.Solutions obtained from the HYDRUS-1D model with z = 1 cm and t = 0.05 days are taken as the truth.

Case 2: dynamic groundwater flow
A 2-D case is analyzed with sharp groundwater flow (see Fig. 4).To minimize the unsaturated lateral flow, the soil surface is set with a non-flux boundary.The bottom and lateral boundaries are also non-flux.Two pumping stresses are applied to the cross-sectional field with x × z = 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 0 (x, z) = 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 1-D soil-water models are discretized with a segmental thickness of z = 1 cm.The fully 2-D unsaturated-saturated solutions from the VSF model are taken as the truth.

Case 3: pumping and irrigation
Case 3 is used to investigate the efficiency and applicability of a quasi-3-D coupling model in comparison of the fully 3-D approaches.A phreatic aquifer with x × y × z = 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 non-flux boundaries.The aquifer is horizontally discretized with x = y = 40 m for the coupled saturated model as well as for the VSF model to obtain the truth solution.The top-down thicknesses of the fully 3-D grid are z = 0.10 m (×30) , 0.4 m (×5) , 1 m (×5) , and 2 m (×5) .
For the 1-D 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.

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 x = y = 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).

Reducing the complexity of a feedback coupling system
The numerical difficulty in a coupled unsaturated-saturated flow system originates from the nonlinearity of the soil-water 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 quasi-3-D model could be lowered by ( 1) taking full advantage of the hand θ -form REs to reduce the nonlinearity in the soil-water 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 switching-form RE for lowering the nonlinearity in the soil-water models.To evaluate the benefits brought by a switching-form RE, the numerical stability was first considered, as shown in Fig. 7.The coupled model in our work was tested with h-form and switching-form REs.Compared with the HYDRUS-1D model (also based on an h-form RE), the switching-form method was numerically more robust, i.e., with larger minimal time-step sizes ( t min ) and a lower computational cost, where a minimal time-step 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 h-form methods, including HYDRUS-1D and the coupled h-form method, t min was constrained to 10 −8 days before reaching a painstaking convergence.In Fig. 8, the soil-water content solution by the coupled switching-form method and the HYDRUS-1D 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 switching-form 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 HYDRUS-1D 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 non-sensitive to the numerical efficiency.A wide range of Se crit [0.3, 0.9] was suggested according to substantial trial-and-error 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 time-varying specific yield calculated by the small-scale soilwater models, S y in Eq. ( 12), introduced significant variability to the large-scale groundwater model, thus caused extra iterations.A large-scale S y reduced the nonlinearity of the storage term in the groundwater equation.In case 1, using an 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 large-scale S y , the nonlinearity introduced by the small-scale soil-water models could be quickly smoothed, as shown in Eq. ( 12).

Multi-scale water balance analysis
The traditional non-iterative 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 sub-model at either side of the phreatic interface, the non-iterative feedback methods usually employed a predicted interfacial boundary without correction, which inevitably introduced coupling errors.In traditional non-iterative methods (Seo et al., 2007;Xu et al., 2012), such shortcomings could be alleviated by refining the macro-time-step 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 soil-water 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 non-iteratively (niter = 0), was easier for approaching the truth (see Fig. 9).Even with refined macro-time-step 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 one-dimensional case with constant upper boundary flux, avoiding interference from lateral fluxes, illustrated the importance of a temporal scalematching analysis for coupling the soil-water 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 scale-mismatching problem.That is, when the coarse-grid groundwater flow solutions are converted into the  vertically distributed fine-scale source and sink terms for the soil columns, an extra downscaling approach is needed to ensure their accuracy.Here we carried out a multi-scale water balance analysis above the phreatic surface.The fine-scale saturated lateral flows were thus excluded from Eq. ( 10).The benefits of the moving-boundary approach can be seen in case 2, which produces significant saturated lateral flux.We carried out a comparative analysis against the traditional stationary-boundary methods (Seo et al., 2007;Xu et al., 2012).The 2-D solution of VSF was taken as the truth.Figure 10 presents the effectiveness of the moving-boundary 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 cross-sectional 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 stationary-boundary 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.

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 non-iteratively (Seo et al., 2007;Shen and Phanikumar, 2010;Xu et al., 2012), is numerically more rigorous than the one-way coupling scheme.The main difference between the above two methods lies in the ability to conserve mass within a single step for back-and-forth 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 back-and-forth 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 "non-iterative".The The HYDRUS-1D solution is taken as the truth.Compared with the stepwise extended method (Seo et al., 2007), the coupling error is significantly reduced by a linear prediction.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 macro-time 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.

Parsimonious decision-making
The feedback coupling schemes, either iteratively or noniteratively, increase the degree of freedom of the users to manage the sub-models 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 scale-disparate 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 1-D soil-water models.Although parallelled computation could further reduce the numerical cost, representing the vadose zone with three sequentially calculated soil-water models achieved acceptable accuracy, as presented in Fig. 13a and b.The computational cost for obtaining the fully 3-D solution with VSF was 15.561 s, which was more than 11 times larger than an iterative feedback coupling method with soil-water models sequentially solved.Problems in more complicated real-world situations can thus be simplified to achieve higher numerical efficiency.

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 two-way 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 sub-models are no longer constrained.In our method, the soil-water 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.

Summary and conclusions
Fully 3-D numerical models are available but are numerically expensive for simulating the regional unsaturated-saturated flow.The quasi-3-D 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 soil-water retention curve, as well as the variability in realistic boundary stresses of the vadose and saturated zones, usually results in the scale-mismatching problems when attempting numerical coupling.In this work, the soil-water and models were coupled with an iterative feedback (two-way) coupling scheme.Three concerns about the multi-scale water balance at the phreatic interface are addressed using a range of numerical cases in multiple dimensionalities.We conclude with the following: 1.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 hand θ-form REs to achieve higher numerical stability and computational efficiency.The one-dimensional switching RE was em-ployed 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 time-step size ( t min ) and conserved mass better.
Lowering the nonlinearity of soil-water models with this switching scheme was promising for coupling complex flow modeling systems at the regional scale.
2. Stringent multi-scale water balance analysis at the water table was conducted to handle scale-mismatching 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 soil-  water and groundwater interactions, the coupling errors were significantly reduced by using an iterative feedback coupling scheme.The multi-scale water balance analysis mathematically maintained numerical stabilities in the sub-models at disparate scales.
3. 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 soil-water 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 soil-water 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, Hydrol.Earth Syst.Sci., 23, 637-655, 2019 www.hydrol-earth-syst-sci.net/23/637/2019/ ∂z| z=z b (positive into the balancing domain and negative outside); dz t = z t (t) − z t (t − dt) 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 sub-domain, the lateral flux l =

Figure 2 .
Figure 2. The Dirichlet-Neumann coupling of the soil-water and groundwater flow models at different scales.(a) Linear or stepwise prediction of Dirichlet lower boundary for the soil-water flow model.(b) Water balance analysis based on a balancing domain with moving lower boundary.Blue dashed line is the linearly extrapolated groundwater table as an alternative for prediction of Dirichlet lower boundary.J (or j ), T (or t), and T (or dt) are the time level, time, and time-step size at coarse (or fine) scale.At any of the transient states (t), the balancing domain is bounded by a user-specified top elevation (z s ) and the moving phreatic surface (z t ).At a transient time t (or T J ), the total mass volume in the moving balancing domain is indicated by w(t) (or w(T J )).The saturated lateral flux of the moving domain is indicated by l(t), while the unsaturated lateral flux is neglected as the assumption of quasi-3-D models.The water flux into and out of the balancing domain is indicated by q top and q bot .

Figure 4 .
Figure 4. Schematic of the cross section for test case 2. Two pumping wells with screens of z = 0-200 cm are located at x = 2500 and 5000 cm.The pumping rates per unit width at well no. 1 and no. 2 are 2 × 10 4 and 1 × 10 4 cm 2 days −1 , respectively.

Figure 6 .
Figure 6.Input of the synthetic regional problem including (a) land surface elevation, (b) initial head, (c) bedrock elevation of the aquifer, and (d) the subzones and boundaries.

Figure 7 .
Figure 7.The time-step sizes through the simulation of (a) sudden infiltration into a dry sandy soil column and (b) rapidly changing atmospheric upper boundary conditions with a layered soil column.

Figure 8 .
Figure 8.Comparison of soil moisture content at z = 0, 50, and 200 cm for the layered soil column with rapidly changing upper boundary conditions (scenario 2, case 1).Taking the HYDRUS-1D solution as the truth, RMSEs of solution of the developed model are provided at different soil depths.

Figure 9 .
Figure 9. Water table changing with time for different macro-time-step sizes ( T = 0.005, 0.05, 0.1, and 0.2 days), in scenario 1, case 1.The HYDRUS-1D solution is taken as the truth.Compared with the stepwise extended method(Seo et al., 2007), the coupling error is significantly reduced by a linear prediction.

Figure 10 .
Figure 10.Comparison of RMSE of (a) the phreatic surface and (b) the head solution (at z = 0) between the moving-boundary and the stationary-boundary methods.Three different lengths of the stationary soil columns, L = 1000, 500, and 300 cm, are considered.

Figure 11 .
Figure 11.(a) The number of feedback iterations and (b) phreatic surface solutions changing with different closure criteria.In the legend, "closure = 0.001" means that ε H = 0.001 cm is used to regulate the feedback iteration.The HYDRUS-1D solution is taken as truth.Tested in scenario 1, case 1.

Figure 12 .
Figure 12.Comparison of (a) water table and (b) head solution (at z = 0) that are changing by the number of soil columns.Solutions obtained with a moving-boundary method in case 2.

Figure 13 .
Figure 13.(a) Comparison of contours of the phreatic surface solution obtained with the fully 3-D and quasi-3-D methods.(b) Comparison of the phreatic surface at A-A' cross section; (c) computational cost and RMSE changing by the number of total soil columns.

Figure 14 .
Figure 14.(a) Comparison of elevation of the water table calculated by the HYDRUS package for MODFLOW(Seo et al., 2007) and the developed method (t = 365 days).(b) The absolute head difference of the phreatic head solution by the method developed here and HYDRUS package at the end of stress periods 3, 6, 9, and 12 (case 4).

Figure 15 .
Figure15.Comparison of water content profiles obtained from the HYDRUS package for MODFLOW(Seo et al., 2007) and the developed iterative feedback coupling method.Subzones 1, 3, 5, 7, and 9 are shown as an example.(t = 365 days in case 4).

Table 1 .
Soil parameters used in the test cases.

Table 2 .
The precipitation, evaporation, and pumping rates in 12 stress periods.