the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Coupling saturated and unsaturated flow: comparing the iterative and the noniterative approach
Natascha Brandhorst
Daniel Erdal
Insa Neuweiler
Fully integrated threedimensional (3D) physically based hydrologic models usually require many computational resources. For many applications, simplified models can be a costeffective alternative. The 3D models of subsurface flow are often simplified by coupling a 2D groundwater model with multiple 1D models for the unsaturated zone. The crucial part of such models is the coupling between the two model compartments. In this work we compare two approaches for the coupling. One is iterative where the 1D unsaturated zone models go down to the impervious bottom of the aquifer, and the other one is noniterative and uses a moving lower boundary for the unsaturated zone. In this context we also propose a new way of treating the specific yield, which plays a crucial role in linking the unsaturated and the groundwater model. Both models are applied to three test cases with increasing complexity and analyzed in terms of accuracy and speed compared to fully integrated model runs. The noniterative approach is faster but does not yield a good accuracy for the model parameters in all applied test cases, whereas the iterative one gives good results in all cases. Which strategy is applied depends on the requirements: computational speed vs. model accuracy.
Threedimensional (3D) hydrologic models for subsurface flow are an important tool to gain a better understanding of the processes in a hydrologic system and to estimate the impact of possible future changes to it. However, these models tend to be very computationally demanding especially when applied to settings on large length scales that require a fine grid resolution (Kollet et al., 2010). Therefore, a lot of approaches to reduce the complexity of such models have been developed. One very common approach is to divide the model domain into two components: the saturated (groundwater) part and the unsaturated part. This is often complemented by the assumption that flow in the unsaturated part occurs predominantly in the vertical direction. Thus, the unsaturated zone, which is usually the more complex one due to the strong nonlinearity of the Richards equation, can be represented by one or multiple onedimensional soil columns that can be solved independently from each other. Groundwater flow is then either modeled as 3D or as 2D depthintegrated flow. The coupled model can then be seen as 2.5D or quasi 3D.
In this work we focus on phreatic aquifers where the crucial part is the coupling of the two model components. Usually, they are linked by the exchange of information on groundwater recharge and the position of the groundwater table. This can be done either in an iterative manner to account for the feedback between the two models (Refsgaard and Storm, 1995; Shen and Phanikumar, 2010; Zhu et al., 2011; Kuznetsov et al., 2012; Xie et al., 2012; Xu et al., 2012; Mao et al., 2019; Zeng et al., 2019) or in a noniterative manner without feedback (Pikul et al., 1974; Yakirevich et al., 1998; Beegum et al., 2018; Erdal et al., 2019). Both methods have been applied successfully. While iterative coupling in general yields a higher accuracy, noniterative coupling requires less computational resources (Zeng et al., 2019).
Differences among the coupling strategies exist also in terms of the spatial extent of the unsaturated zone model. Mostly, the unsaturated zone is modeled to span from the groundwater table to the land surface (Pikul et al., 1974; Refsgaard and Storm, 1995; Yakirevich et al., 1998; Zhu et al., 2011; Kuznetsov et al., 2012; Xie et al., 2012; Zhu et al., 2012; Renard and Tognelli, 2016; Erdal et al., 2019; Zeng et al., 2019). As the groundwater table is moving with time, a remesh is needed after each coupling time step. Others assign a small overlap of the two model compartments which allows for a better calculation of the recharge flux (Xu et al., 2012; Beegum et al., 2018; Mao et al., 2019). A third approach is to consider 1D soil columns that reach down to the bottom of the groundwater domain (Shen and Phanikumar, 2010). Thus, the two compartments overlap in the entire saturated zone. The drawback of this approach is that the size of the unsaturated model grid is increased. The advantage is that lateral fluxes in the groundwater can be passed directly to the unsaturated zone model and no remeshing is needed.
All these different coupling strategies have been tested successfully against real data or fully integrated 3D models. One has to keep in mind that there are limitations to the applicability that are not related to the coupling strategy but the general 2.5D or quasi 3D approach. Due to the 1D representation of the unsaturated zone, these models cannot be applied in the presence of large lateral fluxes in the unsaturated zone which might be the case in steep hillslopes. Furthermore, they neglect lateral water movement in the capillary fringe, which can play a significant role.
One shared problem of these coupled models is due to the fact that the exchange quantities, e.g., recharge and water table position, are usually assumed to be constant over one coupling time step. This might lead to nonconsistent model states when using a coupling scheme without feedback. While the consistency of the two model compartments is not necessary when the compartments are spatially separated, it is fundamental when they overlap. The easiest way to prevent these inconsistencies is to choose a small enough coupling time step, but this might lead to unfeasible compute times. Beegum et al. (2018) overcame this problem by updating the pressure head profile in the unsaturated zone after each time step in such a way that the new groundwater table position and the current recharge flux are both matched. By this procedure, they prevent sudden inflow or outflow of the unsaturated zone due to a stepwise changing groundwater table position. However, this requires the application of an optimization routine which can be rather time consuming. With an iterative coupling scheme, consistency can be achieved. In the MIKESHE model (Refsgaard and Storm, 1995), the water table is adjusted stepwise during the iterations until the mass balance is met with a predefined accuracy. But there still remains an error due to neglection of the temporal changes in the exchange information. Zeng et al. (2019) use linear extrapolation of the groundwater table position to reduce this error. Zhu et al. (2012) solve the two models jointly linked via the pressure head values at the interface. However, this is only possible if the groundwater model is 3D and if the time steps in both compartments are equal. Often, the time step of the groundwater model is chosen to be larger than in the unsaturated model where smaller time steps are needed to achieve convergence for the nonlinear system of equations (Beegum et al., 2018).
Regarding the model consistency, another important factor is the quantification of the specific yield. It relates the amount of water being added or released to the fluctuations of the groundwater table, thus accounting for the neglected unsaturated zone. To get a consistent solution in both model parts in terms of water table position and mass conservation, it is obvious that a correct estimation of the specific yield in the groundwater model is needed.
However, the specific yield is not easy to quantify as it depends on soil properties, change of the groundwater table and the water content profile above the groundwater table (Pikul et al., 1974; Sophocleous, 1985; Crosbie et al., 2005; Hilberts et al., 2005). The common approach to define it as drainage porosity and assign a temporally and often even spatially constant value is a strong simplification in most applications. Beegum et al. (2018) performed a small sensitivity analysis on this parameter and detected a strong influence on the dynamics of the groundwater table position. Some research was dedicated to getting a better estimation of this value in the context of groundwater modeling (Hilberts et al., 2005) and coupled saturated–unsaturated zone modeling (Pikul et al., 1974; Xu et al., 2012; Zeng et al., 2019). Hilberts et al. (2005) presented an analytical expression for the specific yield depending on the groundwater table position and the soil hydraulic parameters, assuming hydrostatic equilibrium in the unsaturated zone. A similar approach was used by Crosbie et al. (2005), who propose a soil and depthdependent formulation of the specific yield. Pikul et al. (1974) were the first to use a dynamic heterogeneous specific yield for linking the 1D unsaturated zone model with a 2D Boussinesq equation for the groundwater. They calculate it as the difference between saturated water content and the minimum water content at a specified depth. This approach was criticized by Vachaud and Vauclin (1975), who could confirm the dependency of the specific yield on space and time with physical experiments, but they did not agree with its ambiguous definition proposed in the work by Pikul et al. (1974). Xu et al. (2012) implemented three different methods to define the specific yield when coupling the SWAP model (Kroes and Van Dam, 2003) with MODFLOW2000 (Harbaugh et al., 2000) but do not give a comparison of these methods or a suggestion. Zeng et al. (2019) distinguish between a smallscale and a largescale specific yield. The smallscale specific yield is calculated from the unsaturated zone model as the change of water content up to a userspecified elevation above the groundwater table; the largescale specific yield is the commonly used parameter of the groundwater model. The difference between those two is used to correct the recharge flux. So the quantification of the specific yield remains an issue which is not yet fairly solved.
In this work we compare two different coupling methods for linking 1D unsaturated flow with 2D depthintegrated groundwater flow. Both are quite simple, straightforward to implement and do not build on any additional software. One of them can keep the consistency of the two model compartments while still being fast and efficient, whereas the other one is essentially designed to be very fast. The latter approach builds on the setup developed by Erdal et al. (2019) and uses a noniterative, fast (but also rather approximate) coupling scheme, in which the unsaturated columns vary individually in length over the simulation time. Hence, each column is resized every time step to fit the distance between groundwater table and surface. In the other method the unsaturated zone model covers the entire depth to the bottom of the groundwater domain and is coupled in an iterative manner to the groundwater model. Thus, the solution in both model compartments is guaranteed to be consistent. A new way of calculating the specific yield during the iteration is introduced and passed as coupling information between the model components. This prevents a prior calibration of this parameter and can account for its temporal and spatial variability.
We perform a sensitivity analysis with the iteratively coupled model to examine whether the calculated specific yield behaves reasonably and to identify those parameters which have a high impact on the groundwater dynamics in the coupled model. The number of parameters that need to be calibrated can thus be reduced, saving computation time. As the equations describing unsaturated flow are nonlinear, a global sensitivity analysis, where the sensitivity is measured over a larger (predefined) range of parameter values, is preferable. As global sensitivity metric, we choose activity scores based on active subspaces introduced by Constantine and Diaz (2017). This is a rather new method that has been given increasing attention since its publication (Palar and Shimoyama, 2017; Fritz et al., 2019; Leon et al., 2019) and has also been used with subsurface models (Erdal and Cirpka, 2019). The idea of the method is to find the most influential directions in parameter space, in which the directions are linear combinations of the individual parameters, socalled active variables.
The remainder of this paper is structured as follows: the coupled model with the two different coupling strategies as well as the activity scores used for the sensitivity analysis are described in detail in the following section. Afterwards, the used test cases are introduced. We then present and discuss the results of the coupling scheme comparison and the sensitivity analysis. In the final part we give a summary of our experiments and conclusions.
2.1 Governing equations
The coupled model consists of two components, namely, groundwater and unsaturated zone. Groundwater flow is modeled by the 2D Boussinesq equation for an unconfined aquifer:
where S_{y} [–] is specific yield, h [m] is hydraulic head, t [s] is time, $\stackrel{\mathrm{\u203e}}{{K}_{\mathrm{s}}}$ [ms^{−1}] is the depthaveraged saturated hydraulic conductivity, z_{0} [m] is the bottom elevation of the aquifer and R [ms^{−1}] is groundwater recharge.
The mixed formulation of the 1D Richards equation is used to describe flow in the unsaturated zone columns:
where S [–] is saturation, S_{s} [m^{−1}] is specific storage, h_{p} [m] is pressure head, ϕ [–] is porosity, K [ms^{−1}] is the unsaturated hydraulic conductivity, z [m] is the geodetic height above the bottom of the aquifer and Q [s^{−1}] is a source or sink term. The specific storage is here defined as ${S}_{\mathrm{s}}=\frac{\mathit{\varphi}}{{V}_{\mathrm{t}}}\frac{\partial {V}_{\mathrm{t}}}{\partial {h}_{\mathrm{p}}}$ with V_{t} [m^{3}] being the total volume and assuming that $\frac{\partial \mathit{\varphi}}{\partial {h}_{\mathrm{p}}}$ is negligible. We use this formulation of the Richards equation as it is consistent with Parflow (Kollet and Maxwell, 2006), which we use as fully integrated reference model, even though there might be more common formulations and the first term describing compressibility of the whole medium is not relevant for our examples.
The relation between pressure head, saturation and unsaturated hydraulic conductivity is represented by the van Genuchten–Mualem model (Van Genuchten, 1980):
where
[–] is the effective saturation, and the model parameter m [–] is usually given by $m=\mathrm{1}\mathrm{1}/n$. The model parameters α [m^{−1}] and n [–] are related to the poresize distribution, and S_{sat} and S_{r} [–] are the saturated and residual saturation, respectively.
2.2 Coupling
2.2.1 General approach
The coupled model is a simplified subsurface model that can model flow in the groundwater and in the vadose zone, but it has no overland flow. It consists of a 2D depthaveraged model for horizontal groundwater flow (Eq. 1) and multiple 1D models for vertical unsaturated flow (Eq. 2). Each 1D model is located at a cell of the groundwater model grid. To reduce the computational burden, not every groundwater grid cell has a 1D model. Instead, several grid cells of the groundwater model can be lumped into zones that are assigned one 1D soil column. A schematic example of the model grid is shown in Fig. 1. The atmospheric forcing of the unsaturated zone model as well as the recharge that is calculated from the unsaturated zone model and passed to the groundwater model are assumed to be constant within each zone. Therefore, the definition of the zones should be geared to land use, depth to groundwater table or other factors that affect the recharge flux as also done by Xu et al. (2012), Renard and Tognelli (2016) and Erdal et al. (2019) and investigated by Zhu et al. (2012) and Zeng et al. (2019).
Both components are discretized using finite volumes and an implicit Euler time integration. The nonlinear systems of equations are solved using the Newton–Raphson scheme and line search. For the groundwater model, a rectangular grid with uniform grid spacing is used. The 1D grid for the unsaturated zone can be nonuniform.
The boundary conditions on the lateral sides of the groundwater model can be either Neumann or Dirichlet boundary conditions. The upper boundary condition of the unsaturated zone models is a Neumann boundary representing the flux across the land surface. In this work, the flux will be either precipitation or the net flux from precipitation and evapotranspiration depending on the test case. The feedback from the groundwater model to the unsaturated zone models as well as the lower boundary condition of the unsaturated zone models depends on the coupling method described in the following sections.
The time step size of the unsaturated zone model is smaller than that of the groundwater model to guarantee convergence in the former but also numerical efficiency for the full model. The coupling time step, which defines the interval when the two model components exchange their information, is equal to the time step size of the groundwater model. All exchanged information such as groundwater recharge is assumed to be constant during one coupling time step.
Both components as well as the coupling routine are implemented in MATLAB.
2.2.2 Noniterative method
The first coupling method aims at providing a fast and iterativefree (but also rather approximate) solution to the numerical problem. A previous version of the model is published in Erdal et al. (2019). In this method, specific yield is a parameter of the groundwater model that the user needs to specify during model setup. Each unsaturated zone column is set up to range from the groundwater table to the land surface, with a Dirichlet boundary h=0 m at the bottom. A coupling time step is then performed as follows:

Run all unsaturated 1D column models for the coupling time step, collect the computed recharge (i.e., fluxes leaving over the bottom boundary for all 1D time steps during the coupling time step) and interpolate the 2D map of groundwater recharge.

Run the 2D groundwater model for the coupling time step with the computed recharge.

For each 1D column, resize all cells by a uniformed ratio computed such that the resized grid fits the new groundwater table to land surface distance at the column location. Thus, the cell sizes Δz are changed while all other properties of the model are kept as they are, including saturation and hydraulic head which means a gain or loss of water in the 1D models.

Compute the amount of water lost (or gained) in each unsaturated model by resizing the grid.

Add (or subtract) a ratio r of this water to the recharge computed in the next coupling time step.

Move to the next coupling time step and repeat from point 1.
The ratio r can in principle range from 0 to 1. At 1 the coupled model is locally mass conservative, as no water is lost or gained during the grid resizing. However, a local massconservation in a noniterative coupling method may not be wanted. Reasons for this are, for example, (1) the water content in the unsaturated zone is also effected by lateral flows in the groundwater, (2) the added or subtracted water is only added in the next time step, (3) the dynamics around the rising or falling groundwater table is only captured in rough terms, and (4) the column resizing is distributed across the full column. This can cause numerical problems like oscillating groundwater tables. In this work we compute r as the ratio between the incoming water volume from the 1D column (RΔt_{c} [m]) and the volume of water added to (or subtracted from) the groundwater (ϕ⋅ΔH_{GW} [m]):
where i is the index of the column, Δt_{c} [s] is the coupling time step and r is limited to a maximum value of 1. The current formulation gives an indication to what extent the change of groundwater elevation can be attributed to the incoming groundwater recharge. If recharge is the main contribution, the water lost or gained in the resizing could be assumed to have entered or left the groundwater, while if the recharge is small but the change in groundwater table large, we find it more likely that the water may still be in the unsaturated zone. We note that this is an ad hoc solution and that it differs to the one presented in Erdal et al. (2019), where r was simply set to 0. However, numerical tests show good stability and the best overall performance for the applied test cases when r is computed as in this work even though a value of r=1 would be needed to keep strict mass conservation.
2.2.3 Iterative method
The second method couples the groundwater and the unsaturated model parts iteratively. Here, the 1D soil columns span the entire soil thickness from the impervious bottom of the aquifer up to the land surface. Thus, the 1D and the 2D model overlap in the saturated part. The bottom boundary condition of the 1D columns is then a noflow boundary.
The groundwater recharge R [ms^{−1}], that is passed from the unsaturated zone models to the groundwater model, is defined – similar as in Crosbie et al. (2005) – as the amount of water being added to or subtracted from the groundwater over time and calculated from each unsaturated flow model as
where H_{UZ} [m] is the position of the groundwater table in the unsaturated zone model, S_{y} [–] is the specific yield, Δt_{c} [s] is the coupling time step and the superscript (ν) is the iteration counter for the coupling time step loop.
The information passed from the groundwater model to the unsaturated zone models is the lateral fluxes Q_{lat} [ms^{−1}] into or out of the columns cells. These data are calculated at the locations where the unsaturated and groundwater models overlap as the amount of water being added to or subtracted from the groundwater over time minus recharge:
with H_{GW} [m] being the position of the groundwater table in the groundwater model. This cumulative flux is spread over the saturated part of the corresponding 1D model by taking into account the variable cell sizes of the 1D grid:
The local flux q_{lat,i} [ms^{−1}] is given as a source or sink term to the ith cell of the 1D model. This procedure is shown exemplarily for one 1D model in Fig. 2.
In this approach, the position of the groundwater table in both models is not used as exchange information but as closure criterion for the iteration. With the lateral fluxes from groundwater flow being accounted for in the unsaturated zone models, the calculated groundwater table position at the next time step should be equal in both models. However, this is only the case if the specific yield is chosen to correctly represent the behavior of the unsaturated zone in the groundwater model. As was already discussed in the introduction, the specific yield is dynamic. Therefore, it needs to be determined at every time step. In the literature, different approaches of how to define and calculate the specific yield can be found. Most of them lack the dependency on the fluctuations of the groundwater table ΔH. The method proposed by Zeng et al. (2019) takes this into account but requires the definition of a threshold elevation above the groundwater table. Here, we propose a method that calculates the specific yield depending on ΔH without the need for further userdefined parameters:
It represents the ratio of the lateral fluxes and the corresponding change of the groundwater table position in the unsaturated zone model. In our case this is equivalent to relating all incoming fluxes (${Q}_{\mathrm{lat}}^{\mathit{\nu}}+{R}^{\mathit{\nu}}$) to the total change of the groundwater table position ($\mathrm{\Delta}{H}_{\mathrm{UZ}}^{\mathit{\nu}}$). So the specific yield is no longer a fixed model parameter but the ratio of the calculated water gain or loss per calculated groundwater level change where the water level change is calculated using the 1D Richards equation, which does not have the specific yield as a parameter. We thereby to a certain extent impose the unsaturated zone model upon the groundwater model. As the altered specific yield may lead to different lateral fluxes and consequently to a different specific yield, several iteration steps are needed to achieve a consistent solution in both model compartments. The convergence of the specific yield is – as in any numerical scheme – bound to an adequate choice of the numerical grid and in this case additionally to a sufficient number of unsaturated zone models which can both be found by a convergence study. Note that if $\mathrm{\Delta}{H}_{\mathrm{UZ}}^{\mathit{\nu}}\mathrm{\Delta}{H}_{\mathrm{UZ}}^{\mathrm{1}}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}$m, then there are no lateral fluxes, which means that there is no iteration required and that the specific yield stays at its current value.
The integration of the coupled model to the next time step, which is illustrated in Fig. 3, can be thus summarized as follows:

Run all unsaturated 1D column models for the coupling time step with ${q}_{\mathrm{lat},i}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}$ms^{−1} ∀i, calculate the computed recharge according to Eq. (7) and interpolate the 2D map of groundwater recharge R^{1}.

Run the 2D groundwater model for the coupling time step with the computed recharge R^{1}.

Calculate the lateral flux source or sink terms ${q}_{\mathrm{lat},i}^{\mathit{\nu}}$ for the 1D models (Eq. 9).

Run all unsaturated 1D column models for the coupling time step with ${q}_{\mathrm{lat},i}^{\mathit{\nu}}$, calculate the computed recharge according to Eq. (7) and interpolate the 2D map of groundwater recharge R^{ν}.

Calculate and interpolate the 2D map of specific yield ${S}_{\mathrm{y}}^{\mathit{\nu}}$ using Eq. (10).

Run the 2D groundwater model for the coupling time step with the updated recharge R^{ν} and specific yield ${S}_{\mathrm{y}}^{\mathit{\nu}}$.

Repeat steps 3 to 6 until ${H}_{\mathrm{UZ}}^{\mathit{\nu}}{H}_{\mathrm{GW}}^{\mathit{\nu}}\le {\mathit{\epsilon}}_{H}$.

Move to the next coupling time step and repeat from point 1.
Note that an assumption is made in this approach. The source or sink terms q_{lat,i} have an effect on the recharge ($R({Q}_{\mathrm{lat}}=\mathrm{0})\ne R({Q}_{\mathrm{lat}}\ne \mathrm{0})$), which due to the nonlinearity of the Richards equation cannot be quantified. This effect is assumed to be negligible, and the fluctuations due to recharge ($\mathrm{\Delta}{H}_{\mathrm{UZ}}^{\mathrm{1}}$) are kept constant during the iterations, whereas the remainder of the total fluctuations ($\mathrm{\Delta}{H}_{\mathrm{UZ}}^{\mathit{\nu}}\mathrm{\Delta}{H}_{\mathrm{UZ}}^{\mathrm{1}}$) are assigned entirely to the lateral fluxes. Changes of the recharge R^{ν} during the iteration are then only caused by changes in the specific yield (see Eq. 7).
Besides, as groundwater recharge and specific yield are calculated from the unsaturated zone models, they are spatially uniform for all cells of the groundwater model that belong to the zone assigned to that unsaturated zone model. This leads to a simplified recharge and specific yield pattern but is necessary to keep the model computationally efficient.
2.3 Activity scores
We perform a global sensitivity analysis to identify the most influential parameters for the groundwater dynamics in the coupled model. As metric, we use activity scores (Constantine and Diaz, 2017). For this, we define the vector of m model parameters x normalized to a range of [0,1] with the probability density function ρ(x) of the parameters and the scalar model output f(x). In our case, the parameters are the soil hydraulic parameters and the upper boundary condition, defined in more detail in Sect. 3.4. The model output f is either the groundwater table position, its daily fluctuation or the calculated specific yield in the iterative model (see Sect. 4.4).
The active subspace of f is defined by the eigenvectors of the matrix
where W is the matrix of eigenvectors, and Λ is the diagonal matrix of the corresponding eigenvalues λ_{i} in decreasing order. If there exists a λ_{n} that fulfills ${\mathit{\lambda}}_{n}\gg {\mathit{\lambda}}_{n+\mathrm{1}}$, Λ and W can be partitioned into
with Λ_{1} containing the first n eigenvalues on its diagonal and W_{1} consisting of the first n columns of W. The column span of W_{1} is the active subspace of dimension n of the function f. The product ${\mathbf{W}}_{\mathrm{1}}^{T}\mathit{x}$ gives the n first active variables which are linear combinations of the normalized input parameters. A variation of an active variable causes a stronger variation in f(x) than a variation of the inactive variables given by ${\mathbf{W}}_{\mathrm{2}}^{T}\mathit{x}$. Then, the individual sensitivity for the ith parameter can be calculated as
where λ_{j} is the jth eigenvalue, and w_{ij} is the value corresponding to the ith parameter in the jth eigenvector. This metric is called the activity score. For a transient problem, the activity score of a parameter can change with time, leading to a transient activity score a_{i}(t).
In this work, the activity score is computed independently for each time step to create a transient score and normalized such that the sum of all scores at one time step is equal to one:
Often, the integration over the parameter space required in Eq. (11) is not feasible or not possible. Therefore, a Monte Carlo method, where M samples are drawn independently from ρ(x), is used to compute the matrix C as proposed in Constantine et al. (2015):
Another issue is related to the determination of ∇f. As this gradient is often unavailable, a common approach is to fit a polynomial model to relate the model output to the input parameters (Jefferson et al., 2015; Erdal and Cirpka, 2019). Following the suggestion by Erdal and Cirpka (2019), a secondorder polynomial depending on the input parameters is fitted to the data f(x) using standard multiple regression, from which the gradient ∇f can be computed easily:
We test and compare the above described coupling methods using three different test cases. The first two test cases are taken from Beegum et al. (2018) where the second one was originally set up by Morway et al. (2013). As these are rather artificial, we added a third test case which is more complex and more realistic. The first test case is a flow model in a bucket setup where all boundary sides are closed. It represents a 1D flow problem consisting of two homogeneous soil layers with dynamic precipitation evapotranspiration (PET) data. It is used to investigate the coupling in the absence of any lateral fluxes where fluctuations of the groundwater table are only due to groundwater recharge. The second test case is conceptualized with a 3D model domain. However, as all parameters are constant over the length of the domain, it is effectively a 2D flow problem. Here again the soil is homogeneous, but it contains lateral fluxes. The third test case is similar to the second but the dimensions and the precipitation are modified such that the groundwater table falls during a certain period and the soil is made heterogeneous. This test case is a true 3D flow problem. The main features of the test cases are summarized in Table 1. For the 1D bucket model, one 1D soil column that is not coupled to groundwater is solved for comparison. The results of the other two test cases are compared to those obtained from a fully integrated 3D model set up with Parflow (Kollet and Maxwell, 2006).
^{*} Only needed in the coupled model. Using the iterative coupling, this value is the initial value ${S}_{\mathrm{y}}^{\mathrm{0}}.$
For the global sensitivity analysis, a fourth test case is created which is mainly a combination of the second and the third test case. They are combined and modified in such a way that the groundwater table rises and falls, the simulation is long enough to finish the model spinup and parameters are homogeneous within each model compartment.
We performed a convergence study on the spatial grid resolution (Δx=Δy and Δz) and coupling time step size (Δt_{c}) for the three applied test cases using the iteratively coupled model. As we want to compare our results to those of Beegum et al. (2018), we chose the grid sizes corresponding to the accuracy we obtained for the grid resolution used in their work (mean absolute error of the groundwater table over time of $\mathrm{MAE}\le \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}$m) for our final calculations. To guarantee the comparability among the fully integrated and the coupled model, we always use the same grid for the reference and for the coupled model using the two coupling schemes although one has to keep in mind that the vertical grid resolution in the noniteratively coupled model will change during the remeshing. The convergence study showed that in the vertical direction a much finer grid is needed with form factors up to 500. Besides, the third (heterogeneous) test case requires a finer horizontal grid to achieve the same accuracy as in the other two test cases. For all three test cases, the influence of the coupling time step Δt_{c} was negligible. Therefore, we decided to use the maximal possible value Δt_{c}=1 d (as we have daily changing forcing data in test case 1) for all test cases. The grid for the fourth test case used for the sensitivity analysis was not determined from a convergence study since its accuracy does not need to be comparable to the other test cases and since the convergence study might show different results for different parameter combinations. Instead, we chose a grid resolution between the second and the third test case, of which this test case is a mixture.
3.1 Test case 1: 1D flow
This model is taken from Beegum et al. (2018) and covers a domain of 1 m × 1 m × 10 m. All lateral boundaries are noflow boundaries. The net flux of daily varying precipitation and evapotranspiration (PET) shown in Fig. 4 is applied as a Neumann boundary at the soil surface. These data were collected at a weather station in Gdańsk, Poland, in 2011. The initial condition assumes a hydrostatic pressure profile in the lower 6.5 m of the soil corresponding to an initial position of the groundwater table at 3.95 m below the surface. In the upper 3.5 m of the soil, unit gradient flow is imposed by setting the pressure head to −0.283 m. The soil consists of two layers: a sandy soil in the top 2.5 m and a loamy sand soil underneath (see Fig. 5). The soil hydraulic parameters for the two layers are given in Table 2. Since the groundwater table stays below the interface to the sandy soil layer during the entire simulation, the groundwater model is assigned the loamy sand parameters.
The total simulation time is 1 year. The 1D model used for comparison (pure Richards model, not coupled to groundwater) uses an adaptive time stepping scheme with a minimum value of 1 s and a maximum value of 1 h. The grid is uniform with 1000 cells of 0.01 m height. The same spatial resolution is used for the 1D models in the coupled model. The groundwater domain is divided into a 2×2 grid. Each groundwater cell is assigned a 1D model. The coupling time step as well as the time step of the groundwater model is set to 1 d for both coupling methods. The time step for the 1D models is adaptive with a minimum value of 1 s and a maximum value of 1 d.
3.2 Test case 2: 2D flow
The dimensions of this model (Morway et al., 2013; Beegum et al., 2018) are 8000 m × 4000 m × 15 m. Since there is no variability along the 8000 m side, flow is effectively 2D in this test case. There is a constant slope of 0.1 % at the 4000 m side. Two noflow boundaries are assigned along the shorter side walls. On the other two sides there are Dirichlet boundaries corresponding to a groundwater table elevation of H=7 m and H=0.9 m at the top and the bottom of the slope, respectively (Fig. 6). The initial groundwater table elevation is interpolated linearly between the two Dirichlet boundaries. The initial condition for the 1D models is generated by applying hydrostatic equilibrium and assigning a minimal initial pressure head of −1.25 m at locations where the pressure head is below this value. Monthly varying rainfall (Fig. 7) is used as Neumann boundary condition for the land surface. The forcing is repeated for five cycles summing up to a total simulation time of 5 years. The soil parameters are listed in Table 2.
The 3D Parflow reference model uses a uniform grid with $\mathrm{80}\times \mathrm{40}\times \mathrm{75}$ cells with grid size $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{100}\phantom{\rule{0.125em}{0ex}}$m and Δz=0.2 m. The spatial resolution in the coupled model is the same with a 80×40 grid for the groundwater model and 75 cells for the unsaturated zone models. The groundwater domain is divided into 40 zones all consisting of slices along the 8000 m side as shown in Fig. 6. With the flow problem being 2D, this means that the entire domain is actually covered by 1D models. The coupling and groundwater time step is Δt_{c}=1 d, and the time step for the 1D models is again adaptive with a minimum value of 1 s and a maximum value of 1 d. Parflow uses an adaptive time scheme as well, with a maximum value of 1 h.
3.3 Test case 3: 3D flow
The first two test cases are rather simple. Therefore, a third test case is designed where the properties that are relevant for the coupling are closer to those of a model for a real hydrological system. The setup is similar to that of the second test case but with some modifications. The domain size is reduced to 800 m × 400 m × 15 m, and the slope is maintained. The lateral boundaries are the same: the flux boundary at the land surface now consists of three cycles of the precipitation shown in Fig. 7, followed by 2 noflux years and then again 2 years using the original precipitation data. Thus, the total simulation time is 7 years. With these changes, the groundwater table is forced to eventually fall and then rise again. Apart from the more dynamic groundwater table, the influence of heterogeneities should be investigated with this test case. Therefore, instead of the homogeneous soil, three different soil types, representing loam, loamy sand and sandy clay, are distributed throughout the domain as depicted in Fig. 8. The distribution is synthetic and the corresponding parameter values are given in Table 2. In the coupled model, the parameters for the groundwater model are averaged arithmetically over the depth of the saturated zone in that grid cell. This averaging takes place at every time step adapting to the new position of the groundwater table. The initial condition is the same as in the 2D flow case.
The horizontal spatial resolution is increased to $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{5}\phantom{\rule{0.125em}{0ex}}$m. The vertical grid consists of 200 cells and is nonuniform with grid sizes ranging between Δz_{min}=0.75 mm at the surface and Δz_{max}=0.3 m at the bottom. The same resolution is used for Parflow and for the coupled model. With the heterogeneous soil, each grid cell of the groundwater model would need a 1D model to cover the entire domain. As this is computationally not feasible, the groundwater domain is again divided into 40 zones but now spread evenly across the domain each covering an area of 20×16 grid cells. The 1D models are placed at the center of each zone and assigned the soil hydraulic parameters at this location. Thus, also the influence of a simplified recharge (and specific yield) pattern, which is constant within each zone and calculated from the corresponding 1D model, can be investigated. The time step sizes are the same as in the 2D flow model.
3.4 Test case 4: sensitivity analysis
As already pointed out, for this test case the 2D and 3D flow test cases were combined. The domain size and the boundary conditions as well as the initial condition are the same as in the 3rd test case. This was done to maintain the rising and falling groundwater table. To get rid of spinup effects, the simulation time was increased, running four times the precipitation data used for the 3D flow case, adding up to a total simulation time of 28 years. The soil structure was simplified to be homogeneous within each model compartment, i.e., the groundwater model and the unsaturated zone models. For a heterogeneous soil structure, a global sensitivity analysis would not be feasible nor very meaningful. The parameters included in the sensitivity analysis are listed in Table 3. Here, K_{GW} and K_{UZ} are the saturated hydraulic conductivity of the groundwater model and the unsaturated zone models, respectively, whereas the saturated part of the 1D models, that overlaps with the groundwater model, is also assigned the K_{GW} value. This means that the hydraulic conductivity of an unsaturated model cell is changed when the cell becomes saturated or unsaturated. Thus, inconsistencies between the two model compartments that may lead to numerical problems are avoided. The parameters are sampled from uniform distribution functions ρ to cover their entire physically plausible range. Table 3 shows the limits of these ranges. Additionally, the influence of the upper boundary condition is investigated by multiplying the precipitation data with a constant factor r_{PET}. The residual saturation ${S}_{\mathrm{r}}={\mathit{\theta}}_{\mathrm{r}}/{\mathit{\theta}}_{\mathrm{s}}$ and the specific storage S_{s} are excluded from the analysis and set to 0.01 and 0.0015, respectively, as a previous smaller sensitivity analysis had shown no impact of these parameters (not shown). The specific yield values are calculated by the model, so the range is defined by these calculated values. Even though this is not a predefined parameter, its influence on the model output, i.e., the groundwater table dynamics, is of great interest and should be comprehensible and consistent over the entire simulation time.
The horizontal spatial resolution is $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{10}\phantom{\rule{0.125em}{0ex}}$m, and the vertical resolution is Δz=0.1 m. As the sensitivity analysis requires a lot of model realizations, the number of unsaturated zone models and corresponding zones is decreased to eight to save computation time. Similar to the second test case, they are placed evenly along the 800 m side of the domain since the soil is homogeneous and flow thus 2D again. The time step sizes are the same as in the previously described test cases.
The coupling methods are compared in terms of accuracy and efficiency. The spatial and temporal distribution of the groundwater table is used for the evaluation of the accuracy. The wall clock time serves as indicator for speed and is given in Table 4. No information about run times are given in Beegum et al. (2018). Therefore, their results are not included in the speed comparison. All MATLAB calculations are done on a local desktop PC with an Intel Core i5 using two cores for the coupled model runs. The Parflow runs used for the other test cases are performed on the JUWELS supercomputer at the Research Centre Jülich which provides 2567 compute nodes with Dual Intel Xeon CPUs. Twelve cores were used for the simulations.
The iterative coupling uses a dynamic formulation of the specific yield. The resulting values are investigated in terms of plausibility. The groundwater table position, the groundwater table fluctuations and the specific yield calculated from the iterative coupling (all taken at the center of the domain) are used as observations f(x) for the global sensitivity analysis. The resulting transient activity scores (Eq. 14) of the parameters listed in Table 3 are analyzed to identify the most influential parameters and to see whether the calculated specific yield behaves reasonably.
4.1 Results 1D flow
Figure 9 shows the water table elevation over time for the single column reference run and the two coupling schemes. The reference elevation z=0 m is located at 5 m below the surface. This is done to facilitate the comparison to the work by Beegum et al. (2018). Both the iterative and the noniterative coupling scheme show a good match with the single column run. The accuracy of the iterative approach is slightly higher; no difference to the single column run can be seen in Fig. 9. A visual comparison indicates that the coupling applied by Beegum et al. (2018) yields a comparable accuracy.
The dynamic adaptation of the specific yield in the iterative coupling scheme is active only in the presence of lateral fluxes, because otherwise due to the formulation of the recharge in Eq. (7) the groundwater table fluctuations in the two model compartments are already forced to be the same in the first iteration step. Therefore, in this test case the value stays constant at its initial value, and further investigation is not possible.
For the 1D flow problem, the iterative coupling takes roughly three times as much wall clock time as the noniterative coupling (see Table 4). This is due to the larger grid of the iterative model which also includes the saturated part.
4.2 Results 2D flow
The temporal evolution of the groundwater table position at the center of the domain is shown in the upper plot of Fig. 10. While both coupling approaches tend to overestimate the elevation of the groundwater table, increasing with time, the noniterative coupling shows a rather large offset compared to the Parflow run. This difference is also evident in the spatial plot of the groundwater table in the lower part of Fig. 10. The overestimation of the iterative coupling scheme is also visible there but is clearly smaller. Besides, there is a shift in the shape of the groundwater table. The iterative coupling scheme shows the largest difference to the Parflow model close to the boundary at x=4000 m which will be discussed further in Sect. 4.5. The largest difference between the noniterative coupling and the Parflow run occurs at around x=1000 m. The results by Beegum et al. (2018) have a similar accuracy and shape as the results of the iterative coupling approach.
When considering the noniterative model, it is notable that initial time steps are an issue, as the model starts almost immediately with a too high groundwater table. However, the difference also grows slowly over time. Both of these issues occur as well when using smaller time step sizes and may be related to the reference model essentially acting as a bucket without any plausible steadystate solution (i.e., steady state for the groundwater model would have groundwater tables far above the top of the domain). As such, an initial error in the model can never be compensated and will also lead to feedback effects, such as an overestimation of the groundwater recharge due to too short unsaturated zone model columns. In this respect, the 2D flow test case is of limited use when evaluating simplified models for realistic setups.
The lateral fluxes in this setup cause a temporally and spatially changing value of the specific yield in the iterative model, which can be seen in the left plot of Fig. 11. However, the differences among the predefined zones (different lines) and also the temporal fluctuations are rather small. All values are smaller than the proposed value of 0.28, although the difference is less than 0.03. The temporal fluctuations of S_{y} suggest a dependence on the fluctuations of the groundwater table which also show a repeating pattern due to the cyclic forcing, whereas its spatial differences indicate the dependence on the water table position. The right plot of Fig. 11 shows the specific yield over the groundwater table position. The least squares regression line reveals a linear relationship between these two variables, where S_{y} is increasing with increasing H_{GW}. Nevertheless, the temporal and the spatial differences of S_{y} are small here, which might be due to the homogeneity of the soil.
For this setup, again the noniterative coupling is faster than the iterative coupling (see Table 4), but the difference is very small, especially compared to the high run time used for the full 3D Parflow model.
4.3 Results 3D flow
In Fig. 12 the position of the groundwater table at the center of the domain is plotted over time (upper part). Both coupling schemes show a good agreement with the fully integrated 3D model. Similar behavior is found at other positions in the domain. The noniterative scheme has a slight delay when the groundwater table is rising more rapidly, e.g., around t=1000 d and t=2250 d. The groundwater table along the width of the domain (y=400 m) at the end of the simulation is plotted in the lower part of Fig. 12. Again, both coupling approaches match the Parflow result very well with small differences between x=300 m and x=370 m. The mismatch between the fully integrated and the coupled model is even smaller than in the less complex 2D flow case, which will be explained in the discussion part (Sect. 4.5).
The heterogeneity of the soil causes more variability of the groundwater table throughout the domain. Figure 13 shows the difference between the groundwater table resulting from the coupled models and the full 3D model at the last time step in the entire domain. It can be seen that in the major part of the domain the differences are small with an average deviation of $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}{H}_{\mathrm{GW}}}=\mathrm{0.12}\phantom{\rule{0.125em}{0ex}}$m for the iterative model and an average deviation of $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}{H}_{\mathrm{GW}}}=\mathrm{0.16}\phantom{\rule{0.125em}{0ex}}$m for the noniterative model. Areas with larger differences appear at similar locations for both coupling schemes showing the largest deviations close to the lower Dirichlet boundary (between x=300 m and x=400 m). Within these areas the differences between the iterative coupling and the Parflow model are smaller compared to the noniterative coupling. The reason for these differences can be seen in Fig. 14. The left plot shows the standard deviation of H_{GW} at the last time step within each zone corresponding to one 1D model obtained with the fully integrated model. The right plot shows the difference to this standard deviation obtained from the iteratively coupled model. It is evident that the variation of H_{GW} within one zone is much larger between x=300 and x=400 m than in the other parts of the domain. The lumping into zones of equal recharge, specific yield, etc. reduces this variability in the coupled model, leading to larger differences compared to the fully integrated model in these areas. A higher accuracy could be achieved by using more 1D models, but this would also increase the run times.
The heterogeneity also leads to a larger variability of the specific yield in the iteratively coupled model as can be seen in the left plot of Fig. 15. There is a significant variation in space and time, with values ranging between approximately 0.01 and 0.33. A clear relationship between the groundwater table and the specific yield with an increasing S_{y} for increasing H_{GW}, as in the 2D flow case, cannot be seen here. Instead, the soil heterogeneities have a stronger impact on the S_{y} values. The right plot of Fig. 15 shows the normalized average specific yield and the normalized average saturated hydraulic conductivity at the water table position ${K}_{i}=\stackrel{\mathrm{\u203e}}{{K}_{\mathrm{s}}\left({H}_{\mathrm{GW}}\right)}$. Both lines show a very similar temporal pattern and, except for the time span between t=1400 d and t=2500 d, their values are even close to each other.
For this case, the iteratively coupled model needs roughly twice as much run time as the noniteratively coupled model (see Table 4). Both coupled models are notably faster than the fully integrated 3D model.
4.4 Results global sensitivity analysis (GSA)
A global sensitivity analysis is performed to identify the most influential parameters of the coupled model regarding the groundwater dynamics, i.e., the groundwater table elevation H_{GW} and its daily fluctuation ΔH_{GW}, both taken at the center of the domain. The resulting normalized transient activity score (NAS) data are plotted in the lower plots of Figs. 16 and 17. The upper plots of Figs. 16 and 17 show the water table elevation and fluctuation of the coupled model averaged over all realizations, respectively. Low values of the saturated hydraulic conductivity of the groundwater model K_{GW} result in a higher curvature of the groundwater table which in some cases reaches the land surface at the center of the domain. All these realizations are removed from the analysis as the model is not able to represent overland flow properly. Therefore, 316 out of initially 500 realizations are used. The polynomial representation (Eq. 16) was crossvalidated with 26 randomly chosen parameters and found to be appropriate (result not shown).
Figure 16 shows that the saturated hydraulic conductivity of the groundwater domain K_{GW} is by far the most influential parameter, except for in the beginning. There, the van Genuchten parameter α is dominant but its influence decreases with time. This is probably due to a spinup effect. When the model reaches the dynamic steady state (last two cycles), the influence of α has basically vanished. Additionally, porosity ϕ and the specific yield S_{y} still have a small influence, especially when the water table starts to rise again after a draining period. Note again that the specific yield in the iteratively coupled model is calculated by the model during the simulation. All other parameters play a minor role.
The activity scores concerning the water table fluctuations look different as can be seen in Fig. 17. Here, the most influential parameters are porosity ϕ, the saturated hydraulic conductivity of the groundwater model K_{GW} and the specific yield S_{y}. The dominance of these parameters is alternating. The saturated hydraulic conductivity K_{GW} is the most important parameter during norain periods. The fluctuations are then mainly caused by lateral fluxes, which depend strongly on this parameter. In the iteratively coupled model, when a high or low peak in ΔH_{GW} is reached, the recharge is the main driving force which is then mostly influenced by porosity. When looking at Eqs. (1) and (7), one sees that S_{y} can be eliminated when lateral fluxes are negligible compared to the recharge fluxes, which explains why there is no influence of S_{y} under these conditions. These are the extreme conditions and in between the specific yield has the largest influence. At these times, there is also an increase in the influence of K_{UZ}. As will be shown later, this is because the calculated specific yield mainly depends on this parameter.
To analyze the sensitivity for different flow conditions, the observations (H_{GW} and ΔH_{GW} at the center of the domain) during the last 7year cycle, when the spinup is finished, are divided into three categories based on their average groundwater table position during that period: dry (3–7 m), medium (7–11 m) and wet (12–15 m). For each category, the normalized activity scores are calculated and averaged over time. These averaged scores are shown in Fig. 18. The scores shown in Figs. 16 and 17 correspond to the medium case. One can see that, for both the water table position and its daily fluctuation, the sensitivity towards the saturated hydraulic conductivity of the groundwater model is increasing with decreasing depth to groundwater table. Under drier conditions, the influence of porosity and the specific yield is stronger, especially regarding the water table fluctuation. Here, those two parameters become the most influential parameters, while for the water table position the dominance of K_{GW} remains. The influence of K_{UZ} on the water table fluctuations is also increasing with increasing depth to water table. Under drier conditions, the unsaturated zone parameters have a stronger impact on the recharge which explains this behavior. Interestingly, there is no increase of sensitivity towards the van Genuchten parameters α and n.
Overall, we see the importance of the specific yield for the groundwater table dynamics. With the iterative calculation of the specific yield, a prior calibration of one of the most influential parameters is not needed. To see whether this calculated parameter behaves reasonably, a sensitivity analysis towards the calculated S_{y} at the center of the domain was conducted. The results are shown in Fig. 19. Note that the S_{y} in the parameter list is now an initial value to which no sensitivity is expected. The average calculated S_{y} value shows some smaller fluctuations, but overall it converges to a value around S_{y}=0.17, which is a plausible value. The activity scores are also on the whole constant in time with the two most influential parameters being porosity ϕ and the saturated hydraulic conductivity of the unsaturated zone model K_{UZ}. The strong dependency on K_{UZ} is consistent with the results of the heterogeneous test case (see Fig. 15). These two parameters are followed by the van Genuchten parameters α and n, whose influence is already much smaller, though. Besides, there is a small influence of the saturated hydraulic conductivity of the groundwater model K_{GW}. This means that the specific yield is mainly depending on the unsaturated zone parameters. This is reasonable as its intention is to represent the missing unsaturated zone in the groundwater model.
4.5 Discussion
The abovedescribed numerical experiments show that both coupling methods presented in this work are well suitable to substitute a full 3D model under certain conditions. They can correctly capture the dynamics of the groundwater table in the presence of lateral groundwater flow and soil heterogeneities. The iterative model outperforms the fast but rather approximate noniterative model in terms of accuracy for all applied test cases, especially in the second one. Its accuracy is comparable to that of other coupling approaches (Beegum et al., 2018), but it is easy to implement and fast, and the model compartments are consistent.
One general problem of the coupled model can be seen in the second test case. Because of the low groundwater table and the large gradient at the lower Dirichlet boundary (x=400 m), water is also leaving the domain laterally via the unsaturated zone. This cannot be captured by the coupled model. The water thus has to first move vertically through the 1D model before it can leave the domain through the groundwater. This leads to a pileup of water and therefore an overestimation of the groundwater table along that boundary. This can also be seen in Beegum et al. (2018), which supports the assumption that it is a general issue. In the 3D flow case (the third test case) the gradient at this boundary is smaller, causing smaller lateral fluxes in the groundwater and thus also in the unsaturated part. Therefore, the results of the coupled model are on average more accurate even though this test case is more complex than the 2D flow case. Here, on the other hand, we see the influence of the simplified recharge and specific yield pattern, which is however unavoidable for the efficiency of the coupled model.
The global sensitivity analysis shows that the groundwater table position H_{GW} is mainly depending on the saturated hydraulic conductivity in the groundwater model K_{GW}. This was to be expected as it is the main parameter in the steadystate formulation of Eq. (1). The initially high and then decreasing influence of the van Genuchten parameter α indicates larger fluxes in the unsaturated zone at the beginning of the simulation. This is most likely due to the rather artificial initial condition in the unsaturated zone which enforces unit gradient flow down to 1.25 m above the groundwater table and no flow below. The temporal fluctuations of the groundwater table are influenced by more parameters, namely, the saturated hydraulic conductivity in the groundwater model K_{GW}, the specific yield S_{y} and porosity ϕ. Which parameter is dominating depends on the current flow conditions as described in Sect. 4.4.
The dynamic formulation of the specific yield in the iterative coupling approach allows the groundwater model and the unsaturated zone model to respond consistently to the given fluxes. So far this has been a challenge when using an overlapping coupling approach (Zeng et al., 2019). The specific yield values obtained in these experiments confirm the assumption that it is depending on the saturation profile, the fluctuations of the groundwater table and the soil properties. However, these values should be interpreted with some caution. In this model, the specific yield also needs to compensate the influence of the lateral fluxes on the recharge, which cannot be quantified properly. The sensitivity analysis for the parameters’ influence on the calculated specific yield shows that it behaves reasonably, though. It depends on the parameters of the unsaturated zone models, especially on porosity ϕ and the saturated hydraulic conductivity K_{UZ}. The aim of the specific yield is to represent the missing unsaturated zone in the groundwater model; therefore, a strong dependency on the unsaturated zone model’s parameters is plausible.
In this work, two different ways of coupling a 2D groundwater model with multiple 1D unsaturated zone models are presented and compared in terms of accuracy and efficiency. The first approach uses a moving boundary formulation and couples the two domains noniteratively. In the second approach they are coupled iteratively, and the 1D models overlap with the groundwater model throughout the entire saturated part. A recalculation of the specific yield value during the iteration is applied to get to the same water table position in both model compartments.
The results and run times of the two coupling methods are compared to those of fully integrated model runs. While the iterative coupling shows a higher accuracy, the noniterative approach is faster. Both coupled models are significantly faster than the fully integrated model. Besides that, the iterative model is more robust, having a decent accuracy for all applied test cases, which is not the case for the noniterative model, and it can keep the consistency of the two model compartments.
A global sensitivity analysis revealed that the groundwater table position and its dynamics depend mainly on three parameters: the saturated hydraulic conductivity of the groundwater model, the specific yield and porosity. The specific yield is, however, hard to quantify, especially because it is usually not constant over time, and fluctuations can be large in the presence of soil heterogeneities. In the iteratively coupled model the specific yield is calculated by the model itself during the iterations, accounting for this variability. It depends mainly on one additional parameter which is the saturated hydraulic conductivity of the unsaturated zone model. This is an advantage of this method over other existing coupling strategies as this parameter is easier to be quantified and constant over time.
In general, both models qualify for substituting a fully integrated model when computational efficiency is more important than highly accurate results. The iterative approach is more robust and accurate and is recommended for most applications. When speed is essential, the usage of the noniterative model is also an option. Erdal et al. (2019) used a version of the heredescribed noniteratively coupled model successfully to reduce the compute time needed for model spinup. Another application where such surrogate models can play an important role is data assimilation. Data assimilation often requires an ensemble of model realizations which comes with a strong computational burden. How well the herepresented 2.5D model can replace a fully integrated model in this context is part of ongoing research.
The MATLAB code to run the coupled model using the two presented coupling strategies and to visualize the output is available at https://doi.org/10.5281/zenodo.4737010 (Brandhorst and Erdal, 2021). This code also includes the input and reference data for the three test cases.
Simulations and code development were performed by DE and NB. IN acquired the funding. All authors contributed to the design of the experiments, the analysis of the results and writing of the paper.
Insa Neuweiler is a member of the editorial board of the journal.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Computing time has been provided by the Jülich Supercomputing Centre (http://www.fzjuelich.de, last access: 3 May 2021).
This research has been supported by the Deutsche Forschungsgemeinschaft in the framework of research unit FOR 2131 (grant no. NE 824/121) and the Collaborative Research Center
CAMPOS (CC 1253 CAMPOS – Catchments as Reactors) (grant no. CRC 1253).
The publication of this article was funded by the openaccess fund of Leibniz Universität Hannover.
This paper was edited by Philippe Ackerer and reviewed by Gerrit H. de Rooij, Orgogozo Laurent, and one anonymous referee.
Beegum, S., Šimůnek, J., Szymkiewicz, A., Sudheer, K., and Nambi, I. M.: Updating the Coupling Algorithm between HYDRUS and MODFLOW in the HYDRUS Package for MODFLOW, Vadose Zone J., 17, 180034–180034, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p
Brandhorst, N. and Erdal, D.: Coupled saturated and unsaturated flow model, Zenodo, https://doi.org/10.5281/zenodo.4737010, 2021. a
Constantine, P. G. and Diaz, P.: Global sensitivity metrics from active subspaces, Reliab. Eng. Syst. Safe., 162, 1–13, 2017. a, b
Constantine, P. G., Emory, M., Larsson, J., and Iaccarino, G.: Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet, J. Computat. Phys., 302, 1–20, 2015. a
Crosbie, R. S., Binning, P., and Kalma, J. D.: A time series approach to inferring groundwater recharge using the water table fluctuation method, Water Resour. Res., 41, W01008, https://doi.org/10.1029/2004WR003077, 2005. a, b, c
Erdal, D. and Cirpka, O. A.: Global sensitivity analysis and adaptive stochastic sampling of a subsurfaceflow model using active subspaces, Hydrol. Earth Syst. Sci., 23, 3787–3805, https://doi.org/10.5194/hess2337872019, 2019. a, b, c
Erdal, D., Baroni, G., SánchezLeón, E., and Cirpka, O. A.: The value of simplified models for spin up of complex models with an application to subsurface hydrology, Comput. Geosci., 126, 62–72, 2019. a, b, c, d, e, f, g
Fritz, M., Lima, E. A., Tinsley Oden, J., and Wohlmuth, B.: On the unsteady Darcy–Forchheimer–Brinkman equation in local and nonlocal tumor growth models, Math. Mod. Meth. Appl. Sci., 29, 1691–1731, 2019. a
Harbaugh, A. W., Banta, E. R., Hill, M. C., and McDonald, M. G.: MODFLOW2000, The U. S. Geological Survey Modular GroundWater ModelUser Guide to Modularization Concepts and the GroundWater Flow Process, Openfile Report, U. S. Geological Survey, 92, 134, 2000. a
Hilberts, A. G., Troch, P. A., and Paniconi, C.: Storagedependent drainable porosity for complex hillslopes, Water Resour. Res., 41, W06001, https://doi.org/10.1029/2004WR003725, 2005. a, b, c
Jefferson, J. L., Gilbert, J. M., Constantine, P. G., and Maxwell, R. M.: Active subspaces for sensitivity analysis and dimension reduction of an integrated hydrologic model, Comput. Geosci., 83, 127–138, 2015. a
Kollet, S. J. and Maxwell, R. M.: Integrated surface–groundwater flow modeling: A freesurface overland flow boundary condition in a parallel groundwater flow model, Adv. Water Resour., 29, 945–958, 2006. a, b
Kollet, S. J., Maxwell, R. M., Woodward, C. S., Smith, S., Vanderborght, J., Vereecken, H., and Simmer, C.: Proof of concept of regional scale hydrologic simulations at hydrologic resolution utilizing massively parallel computer resources, Water Resour. Res., 46, W04201, https://doi.org/10.1029/2009WR008730, 2010. a
Kroes, J. and Van Dam, J.: Reference Manual SWAP; version 3.0. 3, Tech. rep., AlterraReport 773, Alterra, Green World Research, Wageningen, the Netherlands, 2003. a
Kuznetsov, M., Yakirevich, A., Pachepsky, Y., Sorek, S., and Weisbrod, N.: Quasi 3D modeling of water flow in vadose zone and groundwater, J. Hydrol., 450, 140–149, 2012. a, b
Leon, L. S., Miles, P. R., Smith, R. C., and Oates, W. S.: Active subspace analysis and uncertainty quantification for a polydomain ferroelectric phasefield model, J. Intell. Mat. Syst. Str., 30, 2027–2051, 2019. a
Mao, W., Zhu, Y., Dai, H., Ye, M., Yang, J., and Wu, J.: A comprehensive quasi3D model for regionalscale unsaturated–saturated water flow, Hydrol. Earth Syst. Sci., 23, 3481–3502, https://doi.org/10.5194/hess2334812019, 2019. a, b
Morway, E. D., Niswonger, R. G., Langevin, C. D., Bailey, R. T., and Healy, R. W.: Modeling variably saturated subsurface solute transport with MODFLOWUZF and MT3DMS, Groundwater, 51, 237–251, 2013. a, b
Palar, P. S. and Shimoyama, K.: Exploiting active subspaces in global optimization: how complex is your problem?, in: Proceedings of the Genetic and Evolutionary Computation Conference Companion, ACM Press, 1487–1494, 2017. a
Pikul, M. F., Street, R. L., and Remson, I.: A numerical model based on coupled onedimensional Richards and Boussinesq equations, Water Resour. Res., 10, 295–302, 1974. a, b, c, d, e, f
Refsgaard, J. and Storm, B.: MIKE SHE, in: Computer models of watershed hydrology, edited by: Singh, V. P., Water Resources Publications, Colorado, USA, 809–846, 1995. a, b, c
Renard, F. and Tognelli, A.: A new quasi3D unsaturated–saturated hydrogeologic model of the Plateau de Saclay (France), J. Hydrol., 535, 495–508, 2016. a, b
Shen, C. and Phanikumar, M. S.: A processbased, distributed hydrologic model based on a largescale method for surface–subsurface coupling, Adv. Water Resour., 33, 1524–1541, 2010. a, b
Sophocleous, M.: The role of specific yield in groundwater recharge estimations: A numerical study, Groundwater, 23, 52–58, 1985. a
Vachaud, G. and Vauclin, M.: Comments on “A numerical model based on coupled onedimensional Richards and Boussinesq equations” by Mary F. Pikul, Robert L. Street, and Irwin Remson, Water Resour. Res., 11, 506–509, 1975. a
Van Genuchten, M. T.: A closedform equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892–898, 1980. a
Xie, Z., Di, Z., Luo, Z., and Ma, Q.: A quasithreedimensional variably saturated groundwater flow model for climate modeling, J. Hydrometeorol., 13, 27–46, 2012. a, b
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, 170–181, 2012. a, b, c, d, e
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, 1998. a, b
Zeng, J., Yang, J., Zha, Y., and Shi, L.: Capturing soilwater and groundwater interactions with an iterative feedback coupling scheme: new HYDRUS package for MODFLOW, Hydrol. Earth Syst. Sci., 23, 637–655, https://doi.org/10.5194/hess236372019, 2019. a, b, c, d, e, f, g, h, i
Zhu, Y., Zha, Y.y., Tong, J.x., and Yang, J.z.: Method of coupling 1D unsaturated flow with 3D saturated flow on large scale, Water Sci. Eng., 4, 357–373, 2011. a, b
Zhu, Y., Shi, L., Lin, L., Yang, J., and Ye, M.: A fully coupled numerical modeling for regional unsaturated–saturated water flow, J. Hydrol., 475, 188–203, 2012. a, b, c