A methodology for simulating perched conditions in multilayer aquifer systems with 2D variably saturated flow

Simulatingflow through multilayer aquifer systems is essential for groundwater management applications such evaluating natural recharge, assessing managed aquifer recharge (MAR), the characterization of surface water–groundwater interactions, and groundwater pollution risk analysis. However, dealing with flow through variably saturated porous media with perched water tables is a difficult task. Usual codes for groundwater flow modeling such as MODFLOW, TOUGH2, and FEFLOW require tedious procedures based on simplified assumptions and approximations. This paper presents a simple approach to simulate two‐dimensional flow through perched aquifers and aquitards in deep vadose zones by combining the Dirichlet and Neumann boundary conditions that may apply to any code simulating unsaturated flow. The proposed approach is illustrated with the unsaturated flow code VS2DTI. The main strengths of the proposed methodology are (a) variably saturated flow in the multilayer system is solved using Richards’ equation, taking into account the lateral flow that sustains observed water levels in perched aquifers and variations in recharge; (b) a Dirichlet‐type boundary condition at the top perched water table is substituted by a Neumann‐type boundary condition, allowing for the representation of any disturbance (e.g., infiltration from ponds, pumping from wells, etc.), regardless of duration and intensity; and (c) the impact of the disturbance is evaluated by comparing the responses of the undisturbed and the disturbed systems. The versatility of this methodology is applied to a MAR case study of a deep aquifer in a sedimentary basin where aquitards limit its feasibility.


INTRODUCTION
Current modeling in many worldwide groundwater systems faces added difficulties caused by intensive pumping and reduced recharge due to climate change (de Graaf, Gleeson et al., 2019;Famiglietti, 2014). Consequently, the existence of deeper vadose zones, together with highly heterogeneous sediments and fractured rocks, brings about an increase in the number of aquitards, leading to the development of saturated perched water zones (Oostrom et al., 2013). Depending on the hydraulic properties of the soil layers, the prevailing hydrologic conditions, and the effect of anthropogenic modifications such as pumping, perched water zones may be temporary or permanent. Furthermore, water in these zones may be stagnant or flow laterally, affecting groundwater recharge and contaminant transport (Oostrom et al., 2013). Assessment of the impact of perched water on subsurface flow is one of the most challenging aspects of vadose zone modeling because it contradicts the common assumption of vertical-only flow and may lead to the lateral dispersion of infiltrating water (Niswonger & Fogg, 2008;Robinson et al., 2005). Flow in multilayer aquifer systems is complex and highly nonlinear, so suitable numerical groundwater models are required for proper simulation. One straightforward procedure for quantifying perching conditions consists of solving the unsaturated flow equation for a stratified system in which the spatial variation of hydraulic properties determines the distribution of perched water in the profile (Hinds et al., 1999). Although most popular groundwater flow models such as MODFLOW (Harbaugh, 2005) allow to solving coupled unsaturated-saturated water flow schemes, they do not provide integrated schemes to deal with flow through heterogeneous multilayer aquifer systems alternating unsaturated and saturated layers. MODFLOW 6 provides options to simulate problems involving perched aquifers by changing how the vertical conductance between adjacent cells is calculated and using the weight of the overlying water column to drive flow (Langevin et al., 2017). Bedekar et al. (2012) tested different approaches and modifications of the base MODFLOW code to prevent dry cells from becoming inactive and to achieve a stable solution focused on formulations of the unconfined, partially saturated groundwater flow equation. They managed to solve a synthetic two-dimensional (2D) perched aquifer problem where recharge could be applied on top of an aquitard present in the central part of the domain. The best results were obtained with a modified version of MODFLOW-2000 implementing upstream weighting, the Newton-Raphson method with linear smoothing, backtracking, under-relaxation, and an ORTHOMIN solver. However, these approaches do not consider soil water content gradients for unsaturated flow and may be prone to convergence issues.
On the other hand, FEFLOW (Diersch & Kolditz, 2002) can handle variably saturated flow modeling approaches, and a perched water table problem can be solved using Richards' equation. This way, the perched water height, as well as the water pressure head contours following infiltration, can be determined (Diersch, 2014). However, when several aquitards and perched water tables need to be modeled, the convergence to the numerical solution is prone to numerical instability and computationally expensive (Su et al., 2020;Zha et al., 2019). To avoid such issues, Peleg and Gvirtzman (2010) carried out two independent saturated flow models to assess water flow through two perched karstic aquifers, excluding the flow modeling of the unsaturated zone between both aquifers.
Vertical flow of perched water through a low-permeability perching layer is a function of perched water height and thickness and hydraulic conductivity of the low-permeability

Core Ideas
• A new approach to model multilayer systems with perched aquifers using Richards' equation is presented. • A Neumann boundary condition substitutes the Dirichlet-type one in the perched water table. • The impact of disturbances is evaluated through the variation of state variables (e.g., water pressure head). • A case study to assess the efficiency of managed aquifer recharge is discussed.
layer (Oostrom et al., 2013). Several researchers have made specific approximations to perched flow based on Darcian formulations. Robinson et al. (2005) included a permeability reduction factor between two layers to simulate perched water. However, this method implies extensive knowledge of the three-dimensional structure of the perching horizon. Flint and Ellett (2004) and Flint et al. (2012) used an unsaturated water flow model, TOUGH2, to determine the suitability of the multilayered alluvial system of the San Gorgonio Pass area in southern California for managed aquifer recharge (MAR). They calibrated the model using intensive field-monitored data and setting a constant water table head at the perching layer to estimate the saturated hydraulic conductivity. However, the estimated value did not match the water level observations above the perching layer, and manual adjustment was required. Oostrom et al. (2013), in addition to a one-dimensional (1D) vadose zone flow numerical approach, applied a simple Darcian scoping model to estimate perched-water height as a function of recharge rate and saturated hydraulic conductivity of the perching layer for steady-state conditions. The results helped to define the conditions which were consistent with the observed perched water. Greco et al. (2018) modeled water flow through the slope of Cervinara in southern Italy, coupling the solution of 2D Richards' equation applied to the soil cover with the 1D Darcy equation simulating a perched aquifer developing in the epikarst. Nimmo et al. (2017) presented a more complex model incorporating explicitly distinct formulations for diffuse and preferential flow through alternating sedimentary and fractured-rock materials, respectively. For the diffuse case, they deduced a Darcian-based equation to estimate the effective hydraulic conductivity of the limiting perching layer in terms of the hydraulic properties of the layer above (effective specific yield, uniform thickness, and an exponential decline of head vs. time) and the thickness of the perching layer. The model estimates vertical travel times, the occurrence of perched water at layer interfaces, and the thickness and lateral extent of perched water.
Variations of the water table in an unconfined perched aquifer are usually related to lateral flows and recharge, depending on the prevailing climatic conditions and the depth of the aquitard. However, none of the approaches mentioned above consider using any well-suited BC to explicitly represent the observed water level variations in the perched aquifer while performing a predictive analysis of the system. The objective of this study is to present a quick, straightforward methodological framework for the 2D simulation of the impact of a groundwater flow disturbance, such as that of MAR, in multilayer aquifers containing shallow and deep perched water tables. The methodology was tested with the current implementation of the VS2DTI software (Hsieh et al., 1999). This code is widely applied for 2D unsaturated flow modeling (Healy, 2008;Heilweil et al., 2015). The approach provides an exact solution of the variably saturated flow equation by defining a suitable BC to represent the top perched water table. The impact of the disturbance is evaluated by quantifying the different responses of the disturbed and undisturbed systems. The undisturbed system is defined as a modeling scenario where inputs, outputs, and values of state variables are known and constant variations in time are achieved. Therefore, the undisturbed system does not necessarily represent a steady state. The disturbed system represents the addition of new input or output to the previously defined undisturbed situation.

Theoretical background and methodological approach
The 2D circulation of flow in a porous medium with variable saturation can be described by the robust mixed form of Richards' equation (Celia et al., 1990), allowing a straightforward treatment of saturated and unsaturated conditions. Following the formulation of Szymkiewicz et al. (2015), based on Richards (1931) and Lappala et al. (1987), in its extended form, it can be expressed as where θ is the volumetric water content (dimensionless); is the time [T]; q is a source-sink term [L 3 T −1 ]; e = (θ − θ r )∕(θ s − θ r ), is the degree of saturation (dimensionless); θ s is the saturated water content (dimensionless); θ r is the residual water content (dimensionless); s is the spe- where z is the position [L] of the point with respect to a reference height. The model solves for as the principal dependent variable. The determination of θ(ℎ) and r (ℎ) can be carried out using different models (Brooks & Corey, 1964;Haverkamp et al., 1977;Rossi & Nimmo, 1994;van Genuchten, 1980), the van Genuchten model being the most frequently used.
The infiltrated water circulating through a multilayer aquifer system composed of perched aquifers flows alternatively through unsaturated and saturated zones (Figure 1). A simple approach to simulate this flow consists of defining hydraulic parameters and initial moisture conditions for each layer, imposing recharge and ponding infiltration at the surface, considering zero flux through the lateral boundaries, and letting the model solve unsaturated and saturated flow ( Figure 2a). Nevertheless, this would lead to a misrepresentation of the system. The problem is that the lateral flow sustaining the water table of the top perched aquifer would be missed, and desaturation of the perched water table would follow due to vertical flow through the aquitard (Figure 2b). Another option to model the water table of the top perched aquifer would be imposing a predefined head at the lateral boundaries of the model domain to simulate a constant water table or even a hydraulic gradient between both sides to F I G U R E 2 Different modeling approaches to solve water flow through the top perched aquifer. (a-b) Despite considering an external recharge, an instant after the initial moisture content condition of the system, desaturation of the perched aquifer due to the downward vertical flow takes place. (c) The water table in the top perched aquifer is simulated with lateral total head boundary conditions, causing a bulging effect on the water table in the central part of the domain. (d) The water table in the top perched aquifer is simulated with the inner Dirichlet-Neumann (H-Q) boundary condition proposed in this study (methodological Steps 1-5 in Section 2.1) represent the horizontal flow. However, this would cause nonuniform piezometric surfaces in the central part of the domain with either an upward or downward bulging effect depending on the existence or not of an input recharge, respectively ( Figure 2c).
The methodology presented here to evaluate the impact of a disturbance on a multilayer aquifer with perched water tables is based on the contrast between the responses of the undisturbed system and the disturbed system. The undisturbed system refers to where inputs, outputs, and values of state variables are known and constant variations in time are achieved. The disturbed system represents the addition of new input or output to the previously defined undisturbed situation. The water table in the top perched aquifer is represented with suitable BCs (Figure 2d). The approach provides an exact solution of Equation 1 and requires only a simple implementation in a 2D numerical model.
The following steps describe the methodological approach to represent the undisturbed system in a variably saturated media with perched water tables and then simulate the impact of a new disturbance using a numerical unsaturated flow model: 1. In the model that represents the undisturbed system, a specified total head Dirichlet BC is imposed at the nodes that describe the piezometric surface of the top unconfined perched aquifer ( Figure 3): The bottom BC of the model representing the piezometric level of the deep regional aquifer is set to h = 0. No-flow lateral BCs are defined at both sides of the model. Arbitrary, realistic initial moisture conditions are specified according to the soil properties of the different layers (e.g., field capacity in unsaturated layers and saturation in perched layers). The undisturbed system is represented in the model as a predefined scenario where mean inputs (e.g., rainfall recharge, irrigation return flow) and outputs (e.g., regional aquifer discharge, pumping, river discharge) are known and considered. These inputs and outputs may show variations in time and, therefore, the undisturbed system does not necessarily represent a steady state.

2.
A suitable representation of the undisturbed system is achieved when the response becomes stable in time. "Stability" is assumed when both constant variations in the water balance (i.e., inputs -outputs ≈ constant) and constant variations in total heads and degrees of saturation at every node are reached. 3. The magnitude and direction of water flow (in/out), Q pg [L 3 T −1 ], associated with the Dirichlet BC specified in Step 1 are identified through the output of the global water balance once stability, as defined in Step 2, is reached. This Q pg is transformed into specific flow q p [L 2 T -1 ]: where L x is the horizontal length of the domain ( Figure 3).
4. The specified head BC, BC = H p , is replaced by the Neumann-type BC of specified flux, BC = Q p , at all the nodes in the top perched water table ( Figure 2). This step is necessary to represent the impact of a subsequent disturbance to the system; otherwise, the BC = H p would drain all the infiltrating water at the height of the water table of the top perched aquifer. Q p is calculated by multiplying the size of the grid cell at each node (it can be variable in space) by the specific flow obtained in Step 3: where ΔL X is the length of a given cell in the X direction (Figure 3) 5. For the suitable representation of the undisturbed system, the obtained numerical solution in Steps 1 and 4 has to be equivalent. This equivalence is evaluated in terms of the state variables in Equation 1 (H, h, θ, and w ) and the water balance. 6. The new disturbance to be studied (e.g., MAR action) is introduced in the verified model in Step 5 with a suitable BC (e.g., Neumann-type or Dirichlet-type for MAR or Neumann-type pumping). 7. Finally, the disturbance's impact is evaluated through the difference between the model's results and those of the undisturbed system.
It should be stressed the relevance of Step 4, where the BC H p is replaced by Q p , because if the new disturbance is imposed (Step 6), keeping BC = H p (Step 1), no impact on the system would be detected since the condition H p would keep the response of the undisturbed system. Once the system is disturbed (i.e., pond infiltration or pumping), the BC in the piezometric surface of the top perched aquifer has to allow for a dynamic response.

Case study
The proposed methodological approach is applied to simulate a plausible MAR scenario of the Medina del Campo groundwater body, located in the Duero River basin in centralwestern Spain (Figure 4). The area represents a typical geological and climatic sedimentary basin in central Spain. The stratigraphic layout comprises Quaternary deposits near the surface, primarily alluvial, that show interbedded perching layers of fine silty clay materials, followed by Upper Miocene arkosic coarse grain silty sediments with interbedded layers of marsh carbonates that constitute the deeper regional aquifer.
The climate is continental Mediterranean with an annual average potential evapotranspiration of 750 mm and annual average precipitation of about 400 mm. The fluvial network of the region is composed of small and ephemeral streams, characterized by irregular inter-and intra-annual flow regimes with frequent seasonal null flows and significant autumn and spring floods (Llorente & Bejarano, 2018).
Agriculture has a tremendous economic value in the region, and it has been the primary water consumer since the 1960s. Groundwater extractions have exceeded recharge by far, leading to the progressive decline of groundwater levels in the deep regional aquifer and most likely to degrade many of the typical steppe wetlands (Llorente & Bejarano, 2018). Consequently, several corrective measures have been considered to meet the Water Framework Directive 2000/60/EC requirements. These measures include controlling and limiting groundwater pumping and spreading MAR from streambeds. Assessing MAR feasibility requires simulating groundwater flow through the multilayer system.
The simulation considers three different water sources at the top BC: natural recharge as a fraction of average precipitation, seasonal irrigation return flows, and the volume of water added to the streambed for the MAR action. Wet, average, and dry years were defined according to the quartile distribution of annual rainfall data in the period 1941-2012: wet years are those that have annual precipitation (P) in the fourth quartile (P > 476.31 mm), average years are those where annual precipitation lies within the second and third quartiles (374.11 mm < P < 476.31 mm), and dry years are those that have annual precipitation in the first quartile (P < 331.76 mm). The Duero River Basin Authority estimated both annual natural recharge and irrigation return flow for the period 1941-2012 (http://www.mirame.chduero.es).  These values and MAR were weighted for wet and dry years according to the mean precipitation of average years (Table 1). Two 10-yr climatic scenarios were simulated: (a) an average decade (AD) where all years are average-type and have the same annual natural recharge apportioned from October to April (212 d) and the same irrigation return flow apportioned from May to September (153 d); and (b) a modal decade (MD), in which the most frequent distribution of wet, average, and dry years observed in the reference period 1941-2012 was considered (Table 1). In this last scenario, the natural recharge flow rate varies depending on the type of year, whereas the irrigation return flow rate is assumed to be constant independently of the type of year.

Model structure and implementation in VS2DTI
The methodology was tested with the graphical software VS2DTI Version 1.3 for simulating fluid flow and solute transport in variably saturated porous media (Hsieh et al., 1999). The numerical model used to simulate flow is the USGS's computer model VS2DT (Lappala et al., 1987). This code uses the finite difference method to solve Richards's equation for fluid flow (Equation 1). The VS2DTI package contains a graphical preprocessor to facilitate the construction of a simulation and a post-processor to display simulation results graphically (Healy, 2008). Evaluation of axisymmetric infiltration is one of the most typical applications of VS2DTI (Heilweil et al., 2015;Izbicki, 2002;Szymkiewicz et al., 2015). However, these applications focused either on homogeneous (Heilweil et al., 2015;Szymkiewicz et al., 2015) or heterogeneous (Izbicki, 2002) unconfined aquifers without perched water tables. Healy (2008) compiled other applications of VS2DTI, including the analysis of unconfined aquifer tests, estimation of ground and surface water exchange, and the analysis of water movement in soils in response to evapotranspiration. Nevertheless, the current implementation of the VS2DTI software package allows describing flow through a multilayer aquifer system with perched water tables and adequately implementing the BCs and disturbances described in Section 2.1. Similarly, the methodology is extensible to other unsaturated flow models solving Richard's equation.

2.3.1
Model domain and discretization The disturbed system, a MAR spreading operation through a streambed, can be simulated as a symmetric problem with respect to the streambed axis, taking a 2D section perpendicular to the streambed. The half of the longitudinal section of the streambed equals to 4 m. The average depth to the water table in the top perched aquifer is 2 m, and the average depth to the groundwater level of the lower regional aquifer is 55.5 m ( Figure 5). For the simulation of the undisturbed system (numerical domain not shown), the numerical domain extends 1,000 m in the X direction and 55.5 m in the Z direction, with a constant cell size of 100 m in the X direction to check the suitable representation of the undisturbed system following Step 5 of the methodological approach in Section 2.1 (i.e., the coincidence T A B L E 2 Depth intervals and soil hydraulic parameters of the soil layers defined within the model domain. Parameters are defined in the text

Material
Depth Alluvial soil a 0-1 2. Note. K sX , saturated hydraulic conductivity in the X direction; K sZ /K sX , anisotropy ratios of the vertical (K sZ ) to horizontal (K sX ) hydraulic conductivity; S s , specific storage coefficient; ϕ (θ s ), porosity (ϕ) and saturated water content (θ s ); θ r , residual water content; α, inverse of the air-entry pressure (van Genuchten parameter); n, empirical parameter on the pore size distribution (van Genuchten parameter); θ 0 , initial moisture content.
of the obtained results in Step 1, with BC = H p , and Step 4, in which BC = Q p ). This simple grid provides a good representation of the system as all inputs and outputs are evenly distributed along the boundaries, so flow is mainly vertical and, all the same, a 1D model would also be valid. The 2D representation was chosen for coherency, and the cell size in the X direction (100 m) equals that of the largest cells in the model of the case study. For the disturbed system (i.e., MAR case study [numerical domain illustrated in Figure 5]), the model domain is extended to 4,000 m in the X direction, as the MAR inflow through the streambed originates horizontal flow components. Therefore, the spatial domain of the 2D model comprises, in the X direction, half of the streambed and only one of the margins, up to 4,000 m from the axis of the river, and 55.5 m in the Z direction ( Figure 5). As the largest gradient moves downward with the wetting front, the mesh should be refined around it (Zha et al., 2019). The cell size in the X direction is refined towards the streambed: ∆x = 1 m, from 0 to 10 m; ∆x = 2.5 m, from 10 to 35 m; ∆x = 5 m, from 35 to 70 m; ∆x = 10 m, from 70 to 200 m; ∆x = 25 m, from 200 to 400 m; ∆x = 50 m, from 400 to 1,000 m; and ∆x = 100 m, from 1,000 to 4,000 m ( Figure 5).
Discretization in the Z direction for both the undisturbed and disturbed systems is the same and varies with depth to ensure accuracy and convergence ( Figure 5). Under unfavorable conditions such as infiltration with large hydraulic gradients, the mesh size should be equal or less than a few centimeters (Zha et al., 2019). Here, the mesh is refined in the topsoil in the contact between the perched aquifer and the upper unsaturated area (∆z = 1 cm), and, to avoid numerical instability issues, between formations with large permeability contrast. The thickest elements are the inner elements of the most permeable formations (∆z = 1 m). Cell sizes in the X and Z directions were adjusted through trial and error, considering numerical stability and computational cost issues.

Model parameters and zonation
The model parameters are shown in Table 2. The Mualemvan Genuchten soil hydraulic functions (Mualem, 1976;van Genuchten, 1980) implemented in VS2DTI were used: where α is the inverse of the air-entry pressure [L −1 ]; n is the empirical parameter on the pore size distribution (dimensionless); and = 1 − 1∕ . In VS2DTI, θ s is taken equal to the porosity (ϕ, dimensionless) Parameters' zonation in the model corresponds to the description of the lithostratigraphic column of the recharge site (Table 2 and Figure 4). Values are based on reference literature (Carsel & Parrish, 1988;Freeze & Cherry, 1979;Johnson, 1992;Wösten et al., 1999). According to expert judgment and the aquifer system's modeling scale, plausible anisotropy ratios of the vertical to horizontal hydraulic conductivity (K sZ /K sX ) were chosen, assuming either isotropy for some layers (K sZ /K sX = 1) or strong anisotropy for others (K sZ /K sX = 0.1) ( Table 2).

Boundary and initial conditions
Recharge. Two different recharge conditions are imposed on the top of the model ( Figure 5). In one case, a constant specific recharge flux of 8.6 × 10 −5 m d −1 represents rainfall infiltration (Table 3). The other modeled case was a periodic variable recharge in which each year during the first 212 d, the specific recharge flux is 8.6 × 10 −5 m d −1 (wet season), and 6.8 × 10 −5 m d −1 in the remaining 153 d, representing irrigation return flow (dry season).

Specified pressure head.
A new approach to model multilayer systems with perched aquifers using Richards' equation is presented. At the lower limit of the model at 55.5-m depth, the regional deep aquifer level is represented with a null pressure head BC, h = 0 ( Figure 5, Table 3). For the disturbed system simulation, ponding conditions with a maximum ponding depth of 0.3 m are assumed as the MAR source in the streambed. A constant pressure head BC is used to simulate water inflow at the surface. During wet years, it is assumed that surface water surplus in the basin allows to input an amount of water close to the maximum ponding depth of the streambed (h = 0.27 m), whereas in average years, it would be less than this value (h = 0.20 m). Water deficit in dry years prevents MAR actions, and the input flow equals 0 (Table 1).
No-flow BC. No-flow lateral BCs are defined at both sides of the model domain (Table 3). In the MAR disturbed system, the right-side boundary is placed far enough (e.g., 4,000 m) in terms of the perched aquifer's hydraulic diffusivity to avoid interferences with the disturbance. Besides, in the MAR simulation, the no-flow BC in the left side represents the symmetric model's vertical flow line.
Specified total head. In Step 1 of the methodology (see Section 2.1), the piezometric level of the top perched water table, H p = −2 m, is imposed as a preset level. This is imposed at all cells located immediately below −2 m ( Figure 5). This BC is the same for both the constant recharge model and the variable recharge model (Table 3).
Specified flux. In Step 4 (see Section 2.1), at the same cells in which H p was imposed as a BC ( Figure 5), a flow rate (Q p ) is specified. Its value is determined in Steps 3 and 4 (see Section 2.1). In the constant recharge model, the specific flow rate is q pc = 2.207 × 10 −3 m d −1 , and in the case of variable recharge, it is q pv = 2.215 × 10 −3 m d −1 (Table 3), where q pv is the weighted average of the q p values that are generated with the H p condition by recharge period (i.e., 212 and 153 d).
Initial hydraulic conditions. Initial hydraulic conditions were specified as initial moisture content (θ 0 ) contours in the domain. Initially, each soil layer was assigned a soil water content close to field capacity, except for the top perched aquifer and two clayey aquitards that are saturated. These values were then replaced with the average per layer soil moisture content (θ 0 values in Table 2) calculated by the undisturbed system model (Step 4 in Section 2.1). Some soil layers were further subdivided to generate smoother moisture gradients, particularly the alluvial soil in the shallowest vadose zone (Table 2).

Study of the BCs for the representation of the perched aquifer
The validity of Steps 1-4 of the methodology presented in Section 2.1 is analyzed. Validation is performed in terms of the time to stabilization of the state variable, total hydraulic head, H, and the time to stabilization of the global water balance (WB g ) of the system. Simulations to validate the application of the different BCs (H p and Q p ) were run for 100 yr to allow the undisturbed system to reach stability with the imposed initial conditions.
The time to stabilization of H varies with depth. Also, these times are different from the time to stabilization of the WB g .
The disparity between these times indicates that once the system has stabilized globally, the different soil layers require more time to reach equilibrium in the internal distribution of moisture content. In the simulation, the time to stabilization depends on the initial moisture conditions and is influenced by the selected convergence criteria. Figure 6 shows the evolution in 100 yr of the total head, H, at the depth of the top perched water table (2 m), starting from arbitrary, although realistic, initial moisture conditions. The stabilization time of the water balance is practically similar in both models, t = 12,948 d in the case of constant recharge and t = 12,960 d for the variable recharge. This corresponds to 35 yr. In contrast, the times to stabilization of H in the entire domain differ markedly, t = 20,862 d (57 yr) for constant recharge and t = 29,990 d (82 yr) for variable recharge. The stabilized vertical profiles of H and h are equal for both BCs, H p , and Q p , in the models of the two types of recharge (Figure 7). Although the profiles of H and h give the same information, in that of h the effects of the degree of saturation and the characteristics of the different formations are observed more clearly. The unsaturated zone appears in the first 2 m (h < 0) and the linear increase of h between depths of 2 and 11 m defines the top perched aquifer (Figure 7b and 7d). Then, the decrease of h to negative values indicates the presence of the most superficial aquitard (11 to 18 m deep). It is interesting to see, in the thin sandy clay layer below the aquitard, the decrease in moisture content generated by the transition to the most permeable gravel package (18 to 35 m deep), which is at its base. Likewise, a small perched aquifer (approximately 2.5 m thick) is generated on these gravels favored by the existence of an aquitard of small thickness (35 to 37 m deep). Below, there is a very slight increase in humidity in the gravels underlying this aquitard (37 to 42 m deep) driven by the lower permeability of the sandy clay located below them.
Previous works focused on estimating the hydraulic conductivity of the perched aquifer (Flint & Ellett, 2004;Flint et al., 2012;Nimmo et al., 2017;Oostrom et al., 2013;Robinson et al., 2005) or its water level (Greco et al., 2018;Oostrom et al., 2013). Broadly, these studies modeled the properties of perched layers after rainfall or artificial recharge infiltration approximating steady flow driven by gravity simplified to Darcian formulations.
The methodological approach presented here provides a step forward in the characterization of flow through multilayer systems, as it allows defining an undisturbed system scenario prior to the application of any disturbance. This scenario accounts for both variable recharge terms (i.e., rainfall infiltration, irrigation return flow) and lateral flow maintaining the observed water table in the top perched aquifer. Furthermore, the condition is achieved by solving 2D Richards' equation (Equation 1) in the entire domain, instead of using 1D Darcian approximations to perched flow (Greco et al., 2018;Oostrom et al., 2013), which limit the applicability to uniform infiltration over an infinite plane such as widespread rainfall or flooding over a large area (Nimmo et al., 2017).

MAR simulations
The methodology validated in Section 3.1 has been applied to simulate a MAR disturbance of the Medina del Campo multilayer aquifer system as described in Section 2.2. This simulation requires a 2D modeling approach as the BCs now change along the X direction. Total simulation time was set to 70 yr. The recharge on the upper edge of the model in the first 60 yr is the periodic variable recharge described above This 60-yr run represents a spin-up period to allow the model to reach the undisturbed system with the specified flux condition ( Figure 6). For the following 10 yr, the input BCs were modified to simulate two different MAR scenarios: AD and MD. In the AD scenario, the conditions of MAR and periodic variable recharge corresponding to the average year are represented for the last decade (Table 1). In the MD scenario, MAR and periodic variable recharge follow the evolution of the MD (Table 1). Figure 8 shows the degree of saturation in the model domain at the end of the two 10-yr simulations (AD and MD). Artificial recharge does not manage to saturate the soil layers below the top confining aquitard in either of the two 10-yr scenarios. In the AD with constant MAR inflow, the topsoil layer ends up fully saturated up to a horizontal distance of 50 m (Figure 8a). On the contrary, in the MD scenario, the top layer does not end fully saturated because the last year of the simulation is of type "dry" (Figure 8b, Table 1). Both plots show a similar final distribution of soil moisture with depth. It is noteworthy that downward percolation from the top aquitard through the permeable gravel materials accumulates above the deeper aquitard up to a distance of 800 m in the horizontal direction. This last barrier hampers any effect of the MAR of the deep regional aquifer.
The lateral expansion of the artificial recharge input through the streambed in the top left corner of the model domain is approximately 1,800 m ( Figure 9). The lateral expansion of the MAR action is estimated as the horizontal distance at 2 m deep up to where the value of h remains equal in both MAR scenarios during the last day of simulation. Both scenarios show the same lateral expansion because this parameter is a function of the hydraulic diffusivity of the aquifer. The dry-type last year in the MD scenario produces F I G U R E 9 Distribution of the pressure head (h) at 2 m deep in the horizontal direction calculated by the model the last day of the simulation of a 10-yr managed aquifer recharge (MAR) action, under two climatic scenarios: average decade (AD) and modal decade (MD). The lateral expansion (LE) of the disturbance input is estimated as the horizontal distance up to where the value of h remains equal in both MAR scenarios due to hydraulic diffusivity. Note that the y axis is in logarithmic scale significantly lower values of h in the first 400 m controlled by differences in MAR influx and induces slightly lower values of h above 1,800 m in the X direction, reflecting the deviation in the perched water table from the average-type year (Figure 9).
Vertical profiles of h at the end of the MAR simulations highlight the distribution of saturated and unsaturated zones with depth ( Figure 10). The Q pv BC simulating the top perched aquifer allows for a dynamic response following the MAR disturbance. In the AD scenario, constant artificial recharge influx along the 10-yr period causes a larger increase of the water table in the top perched aquifer near the MAR influx point than in the MD case (2.5-m horizontal distance). At the other end of the model domain (3,850-m horizontal F I G U R E 1 0 Vertical profiles of pressure head (h) at three different horizontal distances for the last day of the simulation of a 10-yr managed aquifer recharge (MAR) action, under two climatic scenarios: (a) average decade and (b) modal decade. Light blue: perched aquifer; gray: aquitard distance) where there is a null effect of the artificial recharge, pressure head profiles in both scenarios are almost identical (Figure 10a, b). These profiles represent initial conditions for the MAR in the whole model domain after the 60-yr spin-up period. Maximum increases in h were measured at the bottom of the top perched aquifer (11 m deep) as the difference between the red line and the black dashed line: 2 m in the AD scenario and 1 m in the MD scenario. These pressure head variations represent the part of the upper unsaturated zone that has been saturated by the MAR (Figure 10a, b).
The small amount of recharge water that bypasses the top aquitard accumulates on top of the deeper aquitard as evidenced by the slight increases in pressure head at 35-m depth: 0.4 m in the AD scenario and 0.26 m in the MD scenario (Figure 10a, b). Therefore, the impact of the recharge is limited to the top perched aquifer and the site cannot be considered suitable for MAR.
The results highlight the layers of low permeability in MAR actions as restrictive elements of their effectiveness. In multilayer systems with perched aquifers, aquitards control the MAR process, regardless of the appropriateness of the aquifers' hydrodynamic characteristics.
The study by Flint et al. (2012) demonstrates the importance of iteratively monitoring and modeling the unsaturated zone in layered alluvial systems in the context of MAR actions. The authors show that adequate geologic and hydraulic-property data on perching layers are critical to success. Continuous monitoring in the unsaturated and saturated zones beneath a site provides data to develop and calibrate numerical models, better understand local unsaturated zone process, manage artificial recharge operations, and to determine the timing and volume of recoverable water for consumptive use.
The case study also confirms the suitability of unsaturated soil models for assessing MAR actions, as reported by Sallwey et al. (2018) . These authors point out that the full poten-tial of vadose zone models has not yet been exploited. The approach developed in Section 2.1 adds a valuable contribution to understand the hydraulic processes at a MAR site to determine the feasibility of the action. Moreover, VSD2TI can now be included in the list of vadose zone models used for MAR studies provided in the referred work.
To achieve broader applications, Nimmo et al. (2017) recommended the extension to the transient case as a further development of their model for quantification of flow through sedimentary layers. This would allow their model to predict "early-stage travel times and lateral travel times in perched layers, important for events of relatively short duration, and also estimate the time and conditions required to establish steadiness." With a totally different approach, the present study provides a basis for these developments. First, the response of the undisturbed multilayer system can be defined with the variably saturated flow equation. This response includes the lateral flow that supports observed water levels in perched aquifers and variations in recharge. Second, the initial time step for the simulations was 10 −6 d, and the model converged, which means that there is no limitation for prediction and modeling at shorter time scales. Therefore, based on the Neumann-type BC on the perched water table, the effect of any disturbance can be simulated, regardless of duration and intensity. Third, it allows representing the spatiotemporal variability of the water content in the system.
On the other hand, the methodology shows two main limitations: (a) the undisturbed situation corresponds to a dynamic system where all inputs, outputs, and values of the state variables have to be known; and (b) a 2D cross-section is simulated, assuming an infinite extension of the streambed to keep the approach quick and straightforward. Scaling the approach to a three-dimensional model could be difficult and tedious concerning modeling hypothesis, implementation of BCs, and computational cost. Although only heterogeneity in the Z direction (e.g., multilayer system) was considered in the simulations, there is no limitation to also model heterogeneities in the X direction (e.g., lenses) with the finite difference method used by VS2DTI, as long as suitable mesh refinement around the heterogeneity and model's domain dimensioning based on hydraulic diffusivity comparison are performed. Similarly, contaminant transport simulation will require an adequate spatiotemporal discretization of the model to avoid numerical instabilities. These issues must be further investigated in future research.

CONCLUSIONS
Perched aquifers developed on top of aquitards often show levels that respond to local recharge and lateral flow components. However, none of the existing approaches in the reviewed literature directly consider imposing a BC within the modeling domain to define an undisturbed system before simulating the effect of a given disturbance. A simple methodological framework to model the impact of a disturbance on multilayer aquifers with perched water tables using the variably saturated flow equation has been developed. The method is based on considering the undisturbed system and, once the system is disturbed (i.e., pond infiltration), the BC in the piezometric surface of the perched aquifer has to allow for a dynamic response. This is achieved by shifting from a Dirichlet specified head BC to a Neumann specified flux BC to represent the perched water table.
The method allows calculating transit times in multilayer systems in a simple and computationally economic and efficient way. A 2D simulation for MAR assessment is presented, observing the restrictive nature of the aquitards on the efficiency of the MAR action. Other possible applications of the framework include natural recharge evaluation, contaminant transport, groundwater exploitation planning, and longterm climate change impacts on groundwater resources. This methodology allows simulating and observing, with rigor, the effect of the characteristics of these particular hydrogeological systems on the possible disturbances that they suffer.