Modeling Groundwater Flow in Heterogeneous Porous Media with YAGMod

1 Università degli Studi di Milano, Dipartimento di Scienze della Terra “A. Desio”, via Cicognara 7, I-20129 Milano, Italy; laura.cattaneo1@unimi.it (L.C.); alessandro.comunian@unimi.it (A.C.); giovanna.defilippis@unimi.it (G.F.); chiara.vassena@guest.unimi.it (C.V.) 2 CNR-IDPA (Consiglio Nazionale delle Ricerche, Istituto per la Dinamica dei Processi Ambientali), via Mario Bianco 9, I-20131 Milano, Italy 3 CINFAI, Consorzio Interuniversitario Nazionale per la Fisica delle Atmosfere e delle Idrosfere, Piazza Niccolò Mauruzi 17, I-62029 Tolentino (MC), Italy * Correspondence: mauro.giudici@unimi.it; Tel.: +39-02-5031-8478; Fax +39-02-5031-8489


Introduction
The flow of water in porous sediments is described in mathematical terms by joining the mass conservation principle with the phenomenological Darcy's law [1,2].This results in non-linear equations for variably-saturated media and linear equations for fully-saturated ones.Non-linearity in the forward problem, which is solved to predict the state of the system in response to some forcing (e.g., pumping rates), is introduced by source terms or boundary conditions.The forward problem can be solved analytically only for simple cases, i.e., usually for homogeneous media, domains with simple geometries and for simple source terms.Therefore, analytical solutions can be useful to interpret the results of field tests, e.g., pumping and tracer tests, and of laboratory experiments, where the boundary conditions can be effectively controlled.On the other hand, their use to predict the behavior of complex natural systems is limited by the assumptions that are introduced to simplify the equations and the boundary conditions.Therefore, in practical applications, a variety of computer codes, based on the classical numerical methods of the solution of partial differential equations (finite differences, finite elements, finite volumes, etc.), is used to solve the forward problem.
One of the most challenging problems in hydrogeology is modeling variably-saturated groundwater flow processes.A fully-rigorous solution of this problem requires knowledge of the non-linear relationship of conductivity and matric potential with soil water content, for all the lithologies recognized in the subsoil.Therefore, approximated approaches, which introduce relatively simple modifications of the classical equations for saturated groundwater flow, are often applied.They give rise to numerical difficulties in the presence of dry cells or elements, e.g., under the influence of an extraction source term.When a cell becomes dry, i.e., its calculated water level falls below the bottom of the cell, two main problems arise [3][4][5][6].First, the dry cell cannot receive external water, if it is declared as "inactive"; neither can it contain any extraction source term, unless it is rewetted, i.e., the water level rises above a prescribed threshold.Second, drying and rewetting functionality often yields difficulties for the convergence of iterative algorithms used for the solution of the algebraic equations of the discrete model.Doherty [3] proposes an asymptotically small transmissivity to avoid drained cells being deactivated, even if they actually become dry: this approach uses a function that prevents cell transmissivity from becoming negative.The innovative idea of Keating and Zyvoloski [4] is a weak scaling for vertical connectivity, from partially-saturated to dry conditions.On the other hand, Niswonger et al. [6] use a quadratic approximation of the function that relates horizontal conductance to hydraulic head, over small intervals close to the fully-dry and fully-saturated limits.
Within this background, the first goal of this paper is to propose a code, YAGMod (yet another groundwater flow model), developed in Fortran90, for the simulation of constant-density, groundwater flow under stationary conditions, which is the extension of the codes developed by our research team over the years [7][8][9][10][11][12][13][14][15][16][17].YAGMod is based on a conservative finite difference scheme for stationary conditions and is oriented to the simulation of flow in saturated media.It takes into account the possible drying of shallow blocks of the domain with an original approach, which limits the number of additional parameters that have to be assigned by the user and that have a weak physical significance.Notice that YAGMod considers both prescribed distributed sources and variable point sources.While the former can be used to simulate aquifer recharge, the latter can be used to simulate draining systems or the effects of the water head drawdown on a water-well discharge.
The second goal of this paper is to compare the algorithm used by YAGMod to handle desaturated cells with the above-mentioned algorithms proposed by Doherty [3], Keating and Zyvoloski [4] and Niswonger et al. [6].This is performed by means of a simple, but significant test case.
Moreover, YAGMod includes a module for the model calibration, namely the identification of transmissivity from the knowledge of the reference head and source fields all over the domain.The calibration is performed for 2D flow conditions implementing the comparison model method (CMM).This method was originally proposed by Scarascia and Ponzini [18], successively developed by Ponzini and Lozej [19] and cast in a more formal mathematical framework by Ponzini and Crosta [20].The CMM has been applied and implemented with success to study 2D hydraulic flow in regional aquifers [14,16,17,[21][22][23], and therefore, for the moment, its implementation within YAGMod covers these kinds of systems only.
The third goal of this paper is a thorough and rigorous extension of the CMM to the case of a phreatic aquifer.This is possible because YAGMod takes into account the variation of the saturated thickness of each cell for the solution of the forward problem.
Generally speaking, the "best" model should permit one to properly describe the physical processes and should run with limited computational resources.Unfortunately, these objectives are usually conflicting: physical processes in natural systems might be very complex, non-linear and dependent on a huge number of physical parameters; the computer codes that can handle such situations require supercomputers with high memory capacity and a lot of parallel processors.Simplified approaches are commonly adopted by scientists and engineers.Nevertheless, relatively simple models and tests, like those proposed in this paper, can be very useful to cope with complex natural systems with a computationally-frugal approach, which can provide first insights into the relevant natural processes, on the most sensitive parameters, etc., as recently shown, e.g., by Hill et al. [24].

Forward Model
This section is dedicated to a description of YAGMod, and particular emphasis will be given to its innovative features.

Mathematical Model and Discretization
The 3D flow of groundwater, considering constant fluid density and steady-state conditions, is described by the balance equation: which, using the phenomenological Darcy's law q = −Kgradh, becomes: where q is Darcy's velocity (LT −1 ); h is the water or hydraulic head (L); f is the source term, i.e., the volume of water injected per unit time and unit volume of the porous medium (T −1 ); K is the hydraulic conductivity (LT −1 ).The numerical solution of Equation ( 2) is found with the finite difference method.The continuous physical system is replaced by a finite set of cells or blocks, which are identified by three integer indices (i, j, k), 1 ≤ i ≤ N x , 1 ≤ j ≤ N z and 1 ≤ k ≤ N z , and could be rectangular in the horizontal plane.The side lengths of the cells along the x and y directions, ∆x and ∆y, are assumed to be constant for all of the grid; the cell thickness, denoted by ∆z (i,j,k) , can vary for each cell.The center of a cell is called a node and is denoted with the same indices as the corresponding cell.Values of the hydraulic head are referred to each node, and the spatially-varying hydraulic conductivities are considered to be effective parameters of a cell.
The saturated thickness of a cell is given by: where z (top) (i,j,k) and z (bot) (i,j,k) represent the height of, respectively, the top and bottom surfaces of a cell.For each cell, an integral balance equation can be written.The water discharge into or from a cell is calculated considering only the six adjacent (first-neighborhood) cells (Figure 1).
In order to handle complex aquifers' geometries, each cell is identified by a domain code (Figure 2): I identifies the internal cells, for which the hydraulic head can vary freely; D identifies the cells where Dirichlet conditions are assigned, i.e., the hydraulic head is prescribed; E identifies the cells external to the domain or where no flow takes place.The integral version of Equation (2) could be discretized for each cell as: where h (i,j,k) is the hydraulic head at a node (L); K (i−1/2,j,k) and K (i+1/2,j,k) (LT −1 ) are called internode (or interblock) hydraulic conductivities along the x direction (analogous definitions are used for the similar terms along the y and z directions); ϑ (i+1/2,j,k) = ϑ (i,j,k) + ϑ (i+1,j,k) /2 is the arithmetic mean of the saturated thicknesses of the cells (i, j, k) and (i + 1, j, k) (an analogous definition is used for cells along the y direction); ∆z (i,j,k+1/2) = ∆z (i,j,k) + ∆z (i+1,j,k) /2 is the distance between two adjacent nodes along the vertical; F (i,j,k) is the cell source term, i.e., the volume of water injected in the cell (negative if extracted from) per unit time (L 3 T −1 ).Each of the nine terms appearing in the left-hand side of Equation ( 4) represents the water flux through an interface separating two cells.
A single value of hydraulic conductivity is assigned to every cell, and the internode hydraulic conductivity is calculated as the harmonic mean of the hydraulic conductivities of adjacent cells.
Equation ( 4) can be synthetically written for the most general case, by introducing internode transmittances as follows: Then, Equation ( 4) becomes:

Boundary Conditions and Source Terms
Dirichlet, Neumann or Robin boundary conditions can be assigned.The cells where Dirichlet boundary conditions (prescribed head) are assigned are simply identified by using a D label for the domain code: in that case, the hydraulic head does not change during the computation of the solution.Neumann and Robin boundary conditions are implemented as specific types of source terms, as specified hereinafter.

Distributed Source/Sink Terms
This type of source or sink term simulates areally-distributed fixed source terms, such as rainfall recharge.For each contribution, an array of , represents the flow rate into each cell of the domain; it is independent from h (i,j,k) , and its dimensions are (L 3 T −1 ).The user must consider that this type of source remains constant, even if the hydraulic head of a single cell becomes lower than the bottom of the cell during the iterative search of a solution.

Local Source/Sink Terms
In the YAGMod code, local sources or sinks, i.e., those that are concentrated in a single cell, are modeled with the paradigmatic equation: where F (loc) is the local contribution of the individual source or sink.F (loc) depends on the hydraulic head in that cell, h, according to the difference (h − H (cal) ).F (+) and F (−) are fixed fluxes (L 3 T −1 ); K (+) and K (−) are conductances (L 2 T −1 ); H (cal) and H (act) are two reference head values (L).H (act) is a threshold, which establishes if a source or sink is active or which couple of fluxes and conductances, (F (+) , K (+) ) or (F (−) , K (−) ), should be used to compute F (loc) ; H (cal) , which could be equal to H (act)  in many cases, is used to calculate the contribution to the source term, which linearly depends on h.All of these parameters can be separately defined for each source or sink.Different combinations of F (+) , F (−) , K (+) , K (−) , H (act) and H (cal) allow the user to generate a great variety of source terms, some of which are listed and briefly described below.
• Drain: represents the drain elevation.

• Robin boundary conditions:
These conditions can be used if the aquifer interacts with another water body and water exchange is controlled by the difference of the water head in the aquifer and in the external water body.They are introduced through Equation ( 6), by assigning the following parameters: cal) are the reference hydraulic heads; K (+) and K (−) represent the conductances for flux out or into the cell.K (+) and K (−) depend on the conductivity of the materials that separate the aquifer from the water body at the reference water head and on the distance from this water body.Notice that for the simulation of limited domains of aquifers with a large extension, it is usually impossible to prescribe physically-based boundary conditions.In those cases, Robin boundary conditions are very useful to introduce fictitious boundary conditions, which are more flexible than prescribed head (Dirichlet) or flux (Neumann) boundary conditions.In these situations, H (act) = H (cal) should be close to the estimated water head far from the aquifer system, and the conductances K (+) and K (−) can assume different values to take into account the geometrical, geological and hydrological characteristics of the aquifer.• River/aquifer interaction: H (act) is the height of the bottom of the river: therefore, if h ≥ H (act) , the river and groundwater are in contact, whereas, if h < H (act) , then they are separated by a vadose zone, i.e., partly-saturated sediments or rocks.In the first situation, H (cal) is the river stage; then, the river drains the aquifer if h > H (cal) and recharges the aquifer if h < H (cal) .It is quite common to assume F (+) = 0 and to consider K (+) as a function of the conductivity of the river bed sediments, their thickness and the area of the contact surface between the river bed and the aquifer in the considered cell.In the second situation, namely if h < H (act) , the river bed is assumed to be composed of fine-grained materials, which could be almost saturated, but poorly permeable, whereas the vadose zone between the river bed and the water table could be more permeable than the river bed sediments and approximated as dry.Therefore, the water flows through the river bed under a gravity-controlled, unit hydraulic gradient and freely flows through the relatively permeable vadose zone: then, K (−) = 0, whereas F (−) depends on the conductivity, thickness and extension of the river bed sediments in the considered cell and on the river stage.

Screened Wells
YAGMod considers a new kind of source term that takes into account the dependence of the wells' extraction rate on the aquifer water head.In particular, this kind of source term allows one to turn off the pumping if a cell becomes dry.Sources in this category are denoted as "screened wells", as the user has to give as input data not only the (x, y) coordinates of the well, i.e., the node indices i W and j W , but also the top and bottom elevation of the screened interval (top W and bot W ) and the maximum well extraction rate, q w .The maximum extraction rate is subdivided among the cells intersected by screened intervals, as: with k = n (min) , . . ., n (max) , where: n (min) and n (max) are the indices (along the vertical direction) of the cells corresponding to the top and the bottom of the screened interval of the well; L is the screened thickness of the well corresponding to a fully-saturated porous medium within the (i W , j W , k) cell and is computed as: where: and: ) is corrected, so that the contribution of the well to the source term of the cell (i W , j W , k) is given by: The latter equation implies that, at a given location, the water extracted from a cell reduces as the square root of the thickness of the screened interval that intersects a fully-saturated portion of the aquifer in that cell.

Solution of the Balance Equations
Equation ( 5) can be written for each internal cell, resulting in a system of possibly non-linear equations that can be written in matrix formulation as: where x includes the values of the water head in the internal nodes, b (fix) includes the fixed source/sink terms (Section 2.2.1), b (var) includes the source/sink terms that depend on the water head of the aquifer (Sections 2.2.2 and 2.2.3) and the terms appearing in the left-hand side of Equation ( 5) that involve the water head at D nodes.A is a sparse, symmetric, diagonally-dominant matrix, which is strictly diagonally dominant if at least one D node is present in the domain; its elements are built with transmittances and, therefore, depend on x, as shown by Equations ( 3) and ( 4).The solution to Equation ( 12) could be obtained with any of the methods of solution for non-linear equations that can be found in textbooks of numerical analysis.In YAGMod, a simple approach, based on a generalization of the relaxation methods for the solution of systems of algebraic linear equations, is proposed.This choice is optimal from the point of view of the memory requirement.Other approaches, e.g., Newton's or conjugate-gradient methods, could be more efficient in terms of elapsed running time, if the code is properly modified to profit from parallel computers.However, it should be noted that the specific problem addressed in this paper includes non-differentiable terms in the system of equations, like those introduced by Equation ( 6) and by the sequence of equations from Equations ( 7) to (11).Several tests showed that the generalization of relaxation methods is in general quite robust, in particular for complex physical situations.

Check of the Physical Consistency of the Solution
The proposed model does not solve equations for variably-saturated conditions, but aims at finding a solution for fully-saturated groundwater flow: the cells that become dry during the iterative algorithm of solution are not eliminated from the domain, but are used as auxiliary cells in the sense to be specified below.
If h (i+1,j,k) < z (bot) (i+1,j,k) , then ϑ (i+1,j,k) = 0.If also the adjacent cells along the horizontal directions are dry, then the terms corresponding to horizontal fluxes in Equation ( 5) vanish, and therefore, the cell under examination is involved only for a balance along the vertical direction.This choice permits one to transfer the fixed source terms to deeper cells: this is necessary, e.g., to permit the aquifer recharge, which is assigned at the top active (i.e., internal) cells, to reach the water table.Instead, if the adjacent cells have a non-vanishing, possibly small, thickness, then the physical situation implies that there is a horizontal transfer of water.
When the solution procedure is completed (Section 2.3), the computer code checks the physical reliability of the solution reached.First, a recursive function searches for every continuous path connecting partially-or totally-desaturated cells with the uppermost active cells.At the end of this checking step, every totally-or partially-desaturated cell needs to be connected with the surface, in order to allow air to infiltrate into the porous medium.Second, another function searches for the totally-desaturated cells for which the sum of all of the source terms is negative: such cells actually contain some outflowing source terms, but extracting water from a dry cell cannot be physically acceptable.
If one of these two physically unacceptable conditions occur (desaturated cells not connected with the surface; net outflowing source term from dry cells), YAGMod prints a warning message to the standard output and in the summary output file.

An Example
To show the relevance of these checks and the behavior of the code when dealing with screened wells, two simulations have been run using two different ways to manage an extraction source term.The domain is built up with 100 × 100 × 28 cells of 2.5 m × 2.5 m × 2.5 m dimensions.A low conductivity lens (K (L) = 10 −8 m/s) is included in a homogeneous, permeable medium (K (H) = 10 −4 m/s).In the first simulation, a single deep well, whose extraction rate is about 0.06 m 3 /s, is located beneath the lens at x = 125 m and z = 20 m; in the second one, the single deep well is replaced by a screened well, from 20 m to 45 m, whose extraction rate is the same as the other (obviously, distributed along the whole screened interval).In the latter case, the extraction rate is controlled by the hydraulic head.Results are shown in Figure 3.The simulation run with the deep fixed extraction rate yields physically inconsistent results: in fact, a group of completely desaturated cells is not connected with the surface.On the other hand, the simulation run with the screened well leads to a physically acceptable solution.

A Simple Test Case to Compare Different Approaches to Handle Dry Cells
A number of different approaches to manage dry cells has been proposed [3,4,6].Each approach calculates internode conductivities in a different way; in some cases, also effective extraction rates are calculated taking into account the saturated thickness.In this section, the algorithm implemented in YAGMod is compared to those approaches by means of a test case, which is very simple, but permits one to emphasize some significant properties of the different methods.In particular, the basic characteristics of the analyzed algorithms are briefly recalled using a simplified notation based on this example.A simple 2D domain has been constructed with a grid of 3 × 1 × 2 cells whose size is 100 m × 100 m × 20 m.This 2D domain is illustrated in Figure 4, together with the cell numbering used in the following for the sake of simplicity.At cells (1), ( 3), ( 4) and (6), the hydraulic head is prescribed in such a way as to generate a hydraulic gradient along the x direction: h (1) = h (4) = 40 m, h (3) = h (6) = 39 m.At cell (2), an extraction source term is assigned.
The balance equations for the internal cells ( 2) and ( 5) can be written as: and: In the numerical experiments conducted in this work, h (2) and h (5) vary from minimum to maximum values, i.e., in the interval from 0 m to 40 m.The balance errors ǫ (2) and ǫ (5) are calculated as: The study of the existence and uniqueness of the solution of the problem is based on the analysis of the total quadratic balance error: The method proposed by Doherty [3] uses for the horizontal interblock transmissivity an asymptotically small transmissivity function, in order to keep every cell active, even if it actually becomes dry.This approach uses a decay function that prevents the transmissivity of a dry cell from becoming non-positive: where T is the transmissivity; K is the hydraulic conductivity of a cell; ϑ is the saturated thickness, as for YAGMod; g and f , which are numerical parameters, and ϑ r , the residual saturated thickness, are parameters supplied by the user.To ensure that the function defined in Equation ( 18) is continuous and continuously differentiable, the following relationship must be satisfied: so that the user must specify only two parameters, f and ϑ r .The transmissivity, calculated with Equation (18) for every cell (i, j, k) of the domain, is used to calculate interblock transmissivity with harmonic average.For the vertical water balance, Doherty [3] considers that if any cell in the domain becomes dry, then water inputs from the upper layer remain active.To improve vertical water exchange with the lower layers, i.e., to permit recharge introduced at the model top to reach deeper cells, a linear reduction of vertical interblock resistance (reciprocal of conductance) is introduced using the following equations: where R (2,5) is the interblock resistance; h u is the user-supplied water level below which the linear reduction of resistance is activated; R (u) (2,5) is the "standard" interblock resistance given by: R where m is a user-supplied multiplier.In Figure 5, the results obtained for different values of m from one to 100 and F (2) = Q = 0.1 m 3 s −1 are plotted: no significant difference was noticed among the simulations, so that in further tests, m = 1 was assigned.The second investigated approach was proposed by Keating and Zyvoloski [4].Horizontal interblock transmittance is calculated as follows: ∆y ∆x For vertical transmittance, the change from intrinsic vertical internode conductivity in a partially-saturated cell to zero in a dry cell would lead to a discontinuity.To improve numerical stability, Keating and Zyvoloski [4] allow a weak scaling controlled by a user-specified parameter, ξ: Keating and Zyvoloski [4] suggest using a value of 10 for ξ, which provides both accuracy and stability for all problems.
The last approach considered for this comparison exercise is implemented in the UPW Package of the MODFLOW NTW model [6], a standalone version of MODFLOW 2005 [25].The UPW Package smooths the horizontal conductance function during wetting and drying of a cell.Using this method, horizontal interblock conductance for our test is calculated as follows, by making reference to the example of cells (1) and (2): where which is suggested to be very small (10 −5 ), and In this approach, vertical conductance is calculated as in standard MODFLOW 2005 [25]: The pumping rate is reduced as the head in the cell drops below a user-specified percentage of the cell thickness, as: where: and Φ is a user-specified fraction of the cell thickness, typically assigned as Φ = 0.25.
The results obtained with the application of the four described approaches (i.e., the three approaches described in this section and the one implemented in YAGMod), for three cases corresponding to an extraction rate Q varying from 0.1 m 3 s −1 to 0.2 m 3 s −1 and 0.3 m 3 s −1 , are represented in Figure 6.Notice that with the approach by Keating and Zyvoloski [4], h (5) cannot drop below the cell bottom; therefore, the results obtained with this method for extraction rates of 0.2 m 3 s −1 and 0.3 m 3 s −1 could not be significantly compared to those from other algorithms.Besides this remark, all of the methods give realistic results, even if the values that yield the least total balance error for alternative algorithms differ from each other.17), as a function of the hydraulic heads h (2) (x axis) and h (5) (y axis).From left to right, respectively, the results obtained with YAGMod and the approaches by Doherty [3], Keating and Zyvoloski [4] and Niswonger et al. [6].From top down, the results obtained with extraction rates of 0.1 m 3 s −1 , 0.2 m 3 s −1 and 0.3 m 3 s −1 .The blue circles point out the zone where the least value of the total quadratic balance error is located.
The hydraulic heads computed with YAGMod for an extraction rate of 0.1 m 3 s −1 are smaller than those obtained with the other approaches; on the other hand, for higher extraction rates, the behavior is more complex, and no systematic difference is shown.
The color scales of the plots of Figure 6 are normalized with respect to the minimum and maximum total quadratic errors, separately for each method.Therefore, the images show that all of the methods yield a single minimum, and for YAGMod, ǫ 2 tot increases from the least value more rapidly than for other methods.
Finally, each of the algorithms, which are here compared to YAGMod, is based on some auxiliary parameters.Such parameters are not related to physical processes or quantities, but are necessary to apply artifices to face the numerical problems arising from the simulation of drying cells.These auxiliary parameters, which are listed in Table 1, have to be assigned by the user.On the other hand, YAGMod does not require any additional parameter.This is useful to limit the arbitrariness associated with the assignment of non-physical, auxiliary quantities.
Table 1.Parameters to be assigned by the user for the algorithms with which YAGMod is compared.

Inverse Modeling with the Comparison Model Method
YAGMod permits the application of the comparison model method (CMM) to estimate the conductivity field for a 2D flow field by solving an inverse problem.The CMM was proposed to identify T at every node of a discretization grid [18] and successively developed to directly compute internode transmissivities [19].Further modifications were proposed by Ponzini and Crosta [20], Ponzini et al. [26], Pasquier and Marcotte [27], Ponzini et al. [28].The CMM was applied to alluvial aquifers in Italy [14,16,21], Switzerland [22] and Canada [23,29,30] and to a carbonatic aquifer in southern Italy [17,31].

Fundamentals of the CMM Algorithm
In this section, we summarize the working fundamentals of the CMM for the computation of the conductivity field for 2D stationary flow.Therefore, the k index is omitted in the following equations.
The CMM requires the knowledge of a reference head field, h (ref) (i,j) , which is usually obtained from the interpolation of field measurements and any other relevant geological or hydrological information.Furthermore, a reference source term has to be estimated at every cell of the domain, (i,j) , and it coincides with the same term introduced in Equation ( 5) (for k = 1, as 2D flow is assumed).
Difficulties in estimating a realistic conductivity field from the application of the CMM could arise from the intrinsic instability of the inverse problem or the ill-conditioning of the discrete inverse problem [32]: these difficulties are mainly related to the behavior of the balance equation, which, when used for the forward problem, behaves as a K-to-h map, which smooths the high wavenumber (or short wavelength) components of the conductivity field.From the point of view of the inverse problem, errors on the head at high wavenumbers could be amplified, even if they have a small amplitude [33].Therefore, it is often useful to smooth the reference head and source fields.
An initial guess of the conductivity field, K (CM) (i,j) , must be assigned, for instance, by interpolating the values estimated from field tests.Then, the forward problem is solved for a comparison model (CM): the CM shares the same geometry and source terms as the model, which has to be calibrated.
In particular, the forward problem aims at finding the head field, h (CM) (i,j) , which solves Equation ( 5), for k = 1.For the CM, h (i,j) at the nodes where Dirichlet boundary conditions are assigned.From the results of the forward problem for the CM, the aquifer unit discharge, i.e., the water flow rate through the whole aquifer thickness per unit horizontal length, can be computed at each cell as: for the integral approach and: for the differential approach.

A Test
A simple test is used to show the application of the CMM.In particular, the reference K field is represented in Figure 7 (left image) for a domain with 10 × 10 cells and corresponds to a conductive inclusion in a less permeable domain.The reference hydraulic head is obtained by solution of the forward problem with Dirichlet boundary conditions and is represented by the image map in Figure 8.The results obtained from the CMM with a homogeneous initialization K (CM) = 10 −4 m s −1 are shown in Figure 7, both for the integral (central map) and the differential (right map) approaches.Simple test of the CMM: image plot of the reference hydraulic head and contour lines of the reference h field (purple) and of the h fields obtained for the K fields estimated with the CMM: yellow lines refer to the K field obtained with the integral approach (central map of Figure 7) and green lines to the K field obtained with the differential approach (right map of Figure 7).
The high conductivity region is quite well reproduced by the CMM, although it is a little bit faint, an effect that depends on the spectral response of the balance equation: this induces also some small differences between the water head computed from the K fields estimated with the CMM and the reference h field (see Figure 8), above all around the border of the high conductivity region.

Conclusions
Yet another groundwater flow model, YAGMod, has been described in this paper.The basic characteristics of this model are briefly summarized below.
1.When the water head in a cell is lower than the cell's top coordinate, the cell is considered as partially saturated, and its saturated thickness is used for the water balance.If the water head is lower than the cell bottom, then the cell is considered as dry, but it is not excluded from the water balance calculations.In fact, the water head in this cell is used to compute a vertical water balance that permits one to transfer source terms (namely, aquifer recharge) from the shallow layers down to the deep ones.2. A large number of different sources or boundary conditions, which depend on the water head in the cell, are simulated with a single prototype equation.3.If the position of the screened intervals of a water abstraction well is known, the extraction rate can be limited when the water head falls within the screened interval or turned off when the water head falls below the bottom of the screens.
The comparison of the algorithm used by YAGMod and those proposed by Doherty [3], Keating and Zyvoloski [4] and Niswonger et al. [6] shows that the different approaches yield different solutions, among which, it is not possible to select the best one on the basis of physical arguments.Nevertheless, the advantage of YAGMod is that it does not introduce any additional parameter whose values should be assigned by the user and which mostly have a limited physical significance.Differently, the other approaches require that the user defines one or more additional parameters.These parameters mostly have a limited physical significance, making the parameterization of the model more complex.
Moreover, YAGMod can be used also to solve the inverse problem with the CMM for a phreatic aquifer, with either an integral or a differential approach and taking into account the variability of saturation.This further strengthens the features of YAGMod, because the CMM is embedded in the source code and can be applied to estimate the hydraulic conductivity field, directly at the model scale.The present version of YAGMod is developed to apply the CMM with 2D flow fields, but it can be easily extended to 3D flows, if data with high quality and high density are available.
YAGMod is sufficiently flexible to be adapted to other situations, as was done to model groundwater flow in multi-layered coastal aquifers by De Filippis et al. [31], who modified YAGMod to cope with the variable thickness of the aquifer saturated with fresh water and to identify the spatial variability of a fractured and karst carbonatic aquifer with the application of the CMM to a single layer, while the parameters of other layers are fixed with estimates based on prior information.
The strategy used by YAGMod to cope with draining cells can be easily applied to integrated finite differences [34], but can be adapted also to other numerical techniques, which are based on the discretization of the integral balance equation over a given subdomain.

Figure 1 .
Figure 1.(a) Plan view of a domain's layer; (b) vertical section view.Red arrows are examples of groundwater fluxes through first-neighborhood cells and considered in the discrete model.

Figure 2 .
Figure 2. Domain code example in a 2D domain: I, E and D codes correspond, respectively, to internal, external and prescribed-head (Dirichlet boundary conditions) cells; the blue line denotes the border of the domain.

Figure 3 .
Figure 3. Contour lines of the hydraulic head (2-meter contour interval) and map of the saturation field (fully-saturated cells are drawn in blue; partially-saturated cells are drawn in green; dry cells are drawn in orange) along a vertical section for a 3D problem solved with YAGMod (yet another groundwater flow model).(top) A single deep well, whose fixed extraction rate is about 0.06 m 3 /s, is located at x = 125 m and z = 20 m, beneath a low conductivity lens (dashed area).(bottom) The single deep well is replaced by a screened well, from 20 m to 45 m (yellow line).

Figure 4 .
Figure 4. Scheme of the 2D domain used for the comparison test.

5 )
is the modified interblock resistance, calculated as the reciprocal of the "enhanced interblock conductance" C

Figure 5 .
Figure 5. Map of the balance error given by Equation (17), for the approach of Doherty [3], as a function of the hydraulic heads h (2) (x axis) and h (5) (y axis), for different values of m ((left) m = 1; (center) m = 10; (right) m = 100).At the bottom right corner of each graph, the zone containing the minimum error value is enlarged.

and h ( 2
bottom and top level of the cell corresponding to h up ; η is a small value, usually taken as η

Figure 6 .
Figure 6.Map of the total quadratic balance error given by Equation (17), as a function of the hydraulic heads h (2) (x axis) and h (5) (y axis).From left to right, respectively, the results obtained with YAGMod and the approaches by Doherty[3], Keating and Zyvoloski[4] and Niswonger et al.[6].From top down, the results obtained with extraction rates of 0.1 m 3 s −1 , 0.2 m 3 s −1 and 0.3 m 3 s −1 .The blue circles point out the zone where the least value of the total quadratic balance error is located.

Figure 7 .
Figure 7. Simple test of the comparison model method (CMM): hydraulic conductivity field; the color scale refers to log 10 K, with K expressed in m s −1 .(left) Reference field; (central) field estimated from the CMM with the integral approach; (right) field estimated from the CMM with the differential approach.

Figure 8 .
Figure 8. Simple test of the CMM: image plot of the reference hydraulic head and contour lines of the reference h field (purple) and of the h fields obtained for the K fields estimated with the CMM: yellow lines refer to the K field obtained with the integral approach (central map of Figure7) and green lines to the K field obtained with the differential approach (right map of Figure7).