Sensitivity of groundwater flow with respect to the drain–aquifer leakage coefficient

Throughout the world, groundwater supplies are hidden vital resources that are facing rising pressure owing to pollution and overconsumption from anthropogenic activities. Pollutions putting groundwater at risk notably include: discharge of waste and wastewater onto or into aquifers, use of chemicals such as fertilizers and pesticides, spreading of slurry, manure, and abattoir wastes, poor storage of solvents and petroleum products by using above or underground storage tanks. These last facilities represent a significant contamination risk for groundwater if appropriate mitigation measures, including surface and subsurface drainage, are not designed in their conception and construction (Environment Agency 2017). The drainage system should be designed to convey all potentially contaminated water and spills of fuel to suitable collection points for disposal or treatment. The Groundwater Daughter Directive (2006/118/EC) clarifies the requirements for measures to prevent or limit inputs of pollutants into groundwater. Factors to be considered when carrying out risk assessment for tank leaking into the ground are identified in Environment Agency (2017).

Groundwater modeling tools are of interest to investigate the risk of underground contamination by leakage from storage tanks and to evaluate the performance of drainage systems often designed to prevent groundwater pollution. This paper presents, from a modeling point of view, the impact of a drain set up in an unconfined aquifer to control water fluxes in order to minimize the risk of pollutant migration into the aquifer. The objective is to evaluate the influence of the exchange coefficient in the modeling of a drainage mitigation measure. The modeling of the drain–aquifer interaction is carried out by applying the model of classical river–aquifer interaction. This model is based on a linear relationship between the exchange rate and the difference between the river head and the groundwater head (Rushton & Tomlinson 1979). More precisely, the exchange flux per unit area of riverbed is modeled using the leakage concept presented by Kinzelbach (1986).

According to Bear (1979), mainly horizontal flows are observed at a distance from the river of about 1.5 times the aquifer thickness. In a steady-state unconfined aquifer, groundwater flows between the drain and the aquifer may be thus modeled under the Dupuit‒Forchheimer assumption by considering a horizontal two-dimensional spatial domain and by neglecting vertical flow components. The groundwater flow and leakage interactions between the drain and the aquifer are implemented by two numerical methods: a conforming finite element method (FEM) and a mixed hybrid finite element method (MHFEM). The FEM modeling of such exchanges is classical (Kinzelbach 1986). To the best of our knowledge, the implementation of waterbody–aquifer exchanges in a mixed formulation has not yet been published. A comparison between these two discrete approaches is then performed in terms of computed piezometric head, velocity field, and streamlines considering an engineering case study. The considered case deals with a regional aquifer located below an industrial platform with chemical storage tanks under rain leaching and water infiltration. Computations are performed with both an ideal homogeneous hydraulic conductivity over the domain and a realistic heterogeneous one. A pending question is: ‘how much are the simulation results sensitive to the value of the leakage coefficient?’ since neither an empirical formula exists to calculate this coefficient, nor a device to measure it directly. The uncertainty on the leakage coefficient in this specific context has thus to be evaluated.

As a measure of the model response to perturbations, sensitivity analysis is often acknowledged as a key element in the modeling process (Helton et al. 2006). Several methods exist in the literature. Among the statistical ones, a Monte Carlo method has been used to study the variability of the river–aquifer seepage flow to the spatial variability in the aquifer-saturated hydraulic conductivity (Bruen & Osman 2004). Sensitivities were also computed from the differentiation of groundwater flow equations with respect to the hydraulic conductivity, the aquifer thickness, or the aquifer recharge rate (Mazzilli et al. 2010). In our paper, sensitivities of head, velocity field, and streamlines are computed with respect to the leakage coefficient to provide a new insight into the groundwater modeling, notably to grasp the differences between FEM and MHFEM models in an accurate manner. From a technical implementation point of view, the automatic differentiation (AD) (Griewank & Walther 2008; Hascoët & Pascual 2013) of the numerical codes is performed in the so-called tangent linear (TL) model (Elizondo et al. 2002; Charpentier & Espindola 2005).

The paper is organized as follows. The section Drain–aquifer modeling and discretization methods presents the groundwater flow equation, the leakage modeling, the two-dimensional case study, and the FEM and MHFEM numerical methods, with special attention being paid to the mixed one. The section AD for sensitivity analysis presents the sensitivity analysis and AD in a general manner. It then describes a minimal effort strategy for the AD of existing computer codes. Numerical results, FEM and MHFEM, homogeneous and heterogeneous cases, are reported and discussed in the section Numerical results and sensitivity computations.

Rushton & Tomlinson (1979) studied the leakage between aquifers and rivers by using an idealized one-dimensional problem. More generally, groundwater modeling tools allow for the computation of the exchanges of water between rivers and aquifers and their influences on the quality and quantity of water within both domains (Fleckenstein et al. 2006; Doppler & Hendricks Franssen 2007) and the references therein, for instance.

The mathematical models describing the flow of water in an undeformable porous media are based on the Darcy's law and the mass continuity equation.

The leakage coefficient depends on the thickness and permeability of drain–aquifer interface. These parameters are generally unknown. Consequently, a major issue in the modeling of drain–aquifer exchanges is its parameterization. Classically, the leakage coefficient is chosen as being constant in space and/or time (Kinzelbach 1986; Doppler & Hendricks Franssen 2007).

The hydraulic conductivity tensor, the effective porosity, and the leakage coefficient can also be uncertain. This paper is only concerned with the sensitivity of groundwater flow to the leakage coefficient. More precisely, the objective is to evaluate the impact of a change in the leakage coefficient on the head, velocity field, and streamlines.

The considered case corresponds to a real site in France (Figure 2(a)). The exact location is not given for confidentiality reasons. It is based on a real chemical industrial platform with an area of 1 km² located between a navigation channel to the west and a river to the east. The water flows from the channel (height: 221.5 m) to the river (height: 217 m). The underlying aquifer is subject to drainage for environmental protection. The domain (Figure 2(a) comprises two drains laid out to an elevation of 217.1 m to 217.4 m in a north–south alignment to catch and to evacuate the groundwater flowing below the possible contamination sources. Two cases are investigated by considering the domain as homogeneous or heterogeneous, respectively (Figure 2(b)). In the homogeneous case, the hydraulic conductivity is set to 180 m/d and the barrier is not accounted for in the computation.

The 2D computational domain denoted by is discretized using an unstructured mesh containing M = 10,546 triangular elements. The size of the elements is adapted to handle the confined zone, the drains, and the barrier. Constant Dirichlet boundary condition along the channel and river (east and west boundaries) and Neumann boundaries (north and south boundaries), which are set to zero flux, are considered (see Figure 2(a)).

Equation (1) is now discretized by means of a mixed approach that computes the state variable and its gradient in a simultaneous manner (Meissner 1973). In the groundwater modeling case, the discrete state variable comprises the pressure head and the velocity field. Hybridization is applied to obtain a positive definite matrix (Chavent & Roberts 1991; Mosé et al. 1994; Fahs et al. 2009; Younes et al. 2010).

A classical issue in environmental sciences is the evaluation of the sensitivity of the model outputs to the input parameters (Shao et al. 2017). Parameters may have a physical or a geometrical meaning, or may be some initial or boundary conditions, for instance. This section discusses (i) the AD basics, (ii) sensitivity analysis, and finally, (iii) AD application to FEM and MHFEM numerical codes.

AD (Griewank & Walther 2008; Naumann 2012) is a set of techniques designed to reinforce computer codes with derivative computations. Basic principles of AD may be explained as follows.

Within AD, a computer code may be viewed as a sequence of statements run in a prescribed order. Given the input data, the execution flow includes information about how the code starts, the actual order of execution of the statements, and how it terminates. This execution flow is no more than a large composition of arithmetic operations and intrinsic functions. Being well-suited to sensitivity analysis, the differentiation in TL mode consists in the differentiation of this compound function by applying the chain rule and classical rules such as ‘the derivative of a sum is the sum of the derivatives’. AD may be applied to very large codes (Charpentier 2000). Note that the differentiation in the so-called adjoint mode, frequently used for identification purposes, is beyond the scope of the paper.

AD relies on two kinds of software. On the one hand, source transformation tools such as Tapenade (Hascoët & Pascual 2013) or Adifor (Bischof et al. 1991) are able to generate source codes containing derivative statements. On the other hand, operator overloading libraries such as Adol-C (Griewank et al. 1996) or Rapsodia (Charpentier & Utke 2009) may be used at compile time to propagate derivatives at runtime. In this paper, we use the Tapenade software (Hascoët & Pascual 2013). For small codes (less than 3,000 lines), the web interface gives access to the ‘Tapenade On-line AD Engine’. The user provides the source of its code, the name of the top routine to be differentiated, the ‘independent’ input variables (modeling parameters of interest) and, optionally, the ‘dependent’ output variables (state vector, for instance). When the provided information is coherent, Tapenade generates a TL code differentiated with respect to the independent variables. At runtime, this linear code propagates one (or several) direction(s) of perturbation to compute the dependent variables and their sensitivity(ies).

Given a code and differentiation instructions, Tapenade first proceeds to a syntactic analysis of the code structure and statements (parser phase). The parser checks the correctness of the statements. It notably finds out type incoherences, potential source of errors (equality test on real numbers), and I/O of active variables (independent and dependent), for instance. If these exist, parser errors and warnings must be cured (or at least understood) before the actual differentiation phase.

To emphasize the generality of the proposed method, let us consider the general model ⁠, the inputs and output of which are the modeling parameters P and the state variable x, respectively.

Given a model ⁠, a sensitivity analysis measures the impact of any small perturbation of its input parameter set P on the model response ⁠. Some qualitative information may be deduced from repeated evaluations of the model with different values of ⁠. However, when possible, sensitivity analysis based on derivative computations should be preferred as it generally provides more valuable qualitative and quantitative information on the model (Elizondo et al. 2002; Charpentier & Espindola 2005).

A minimal effort AD strategy to obtain a TL code for sensitivity calculation is needed. The present study takes advantage of two existing groundwater flow codes, namely, GW_FEM and GW_MHFEM, written in Fortran. Their general implementation scheme ‒ modeling parameters P, boundary conditions, mesh, simulation duration, FEM or MHFEM discrete problem, and writing output data … ‒ is similar (see Figure S1 in Supplementary material, available with the online version of this paper). As described in Figure S1, the dependent variable x (representing piezometric head, velocity, and streamline in this study) in the original user code GW_FEM may be automatically differentiated with respect to the independent variable P here restricted to the drain–aquifer leakage coefficient ⁠.

This differentiation method has been successfully applied to both the GW_FEM and the GW_MHFEM codes to generate the FEM_D and MHFEM_D codes, respectively. Some details about the codes (number of routines, number of lines) together with runtimes are provided for the original user codes and the differentiated ones in Supplementary material, Table S1 (available with the online version of this paper). One can notice that the MHFEM solver comprises an iterative bi-conjugate gradient method for the solution of the general linear system ⁠, the differentiation of which is more efficient and accurate by calculating the TL variable following ⁠. It can be observed (Table S1) that the GW_FEM code is faster than the GW_MHFEM code because the number of unknowns is greater in the MHFEM case. The ratios between the TL codes and their respective original code are lower than the theoretical bound of 4 related to the differentiation of the division operator (Morgenstern 1985). This is an excellent result since the TL codes comprise the original code statements for trajectory computations.

In the present case, the GW_FEM and GW_MHFEM codes are differentiated with respect to the leakage coefficient of the northern drain. Taylor tests are carried out considering the position of a particle at the groundwater upstream. Table 1 displays both GW_FEM and GW_MHFEM ratios thanks to the Taylor test performed for the streamlines. One can observe the expected behavior of which converges to 1 when is less than 10⁻⁶. Hence, the TL codes are correct. Moreover, similar convergences are obtained for Taylor tests performed on the head and the velocity field.

Computations are performed with the leakage coefficient m⁻¹ d⁻¹ for the two drains. Sensitivities are computed with respect to a perturbation of 50% of the leakage coefficient of the northern drain only.

Piezometric head, velocity field, streamlines and their sensitivity to the leakage coefficient value are computed in the homogeneous and the heterogeneous transmissivity cases using the GW_FEM and GW_MHFEM codes.

Computed piezometric head and streamlines are displayed in Figure 4. Piezometric head isolines are distributed between the channel (upstream head of 221.5 m) and the river (downstream head of 217 m).

In the homogeneous case, the groundwater smoothly flows from the channel to the river both in the southern and northern zones ⁠. In the middle zone ⁠, some of the fluid particles are caught by the drains. This effect is corroborated by the drawdown of the piezometric head to the drains. Although the two discretization methods yield the same flow pattern, the drawdown is greater in the FEM case. This is an expected result since the FEM and MHFEM formulations manage exchange fluxes between the aquifer and the drain in a different manner (see subsection Numerical methods).

Table 2 reports fluxes computed at the channel, the river, and the drain boundaries as well as the mass balance error (MBE) estimated as the sum of the inlet flow rate at the channel and the outlet flow rate at the river and the drains. Null values for the MBE prove that the two codes are conservative. Table 2 shows that the drains catch much more water in the FEM cases. This is less noticeable in the heterogeneous case.

A similar behavior, (Figure 4(c) and 4(d)), may be observed in the northern and southern zones in the heterogeneous transmissivity case, whatever the discretization method. Streamlines are slightly different in the middle zone. Near to the drain, the streamlines are confined by the almost impervious barrier. Flux values at the channel, river, and drain boundaries are reported in Table 2. MBE values are again 0.

Simulations of the flow with natural heterogeneous hydraulic conductivity without barrier, which is a management structure, have been carried out and compared with the one with the barrier in the Supplementary material (see Table S2, available with the online version of this paper). This confirms that the barrier increases the efficiency of drainage (see Figures S2 and S3, available online).

The sensitivity analysis is carried out with respect to a perturbation of 50% of the northern drain leakage coefficient, that is m⁻¹ d⁻¹, for a deeper understanding of the underground flow.

Figure 5 displays head sensitivity isolines computed in the homogeneous transmissivity case. As expected, a positive perturbation of 50% in the exchange capacity between the aquifer and the north drain results in a negative piezometric head sensitivity indicating a stronger drawdown. Sensitivity isolines accurately describe how and where the drain catches water. Although the drain catches less water using the MHFEM code, it can be seen that the piezometric head computed by MHFEM is the more sensitive near to the drain since the closed domain delineated by the isoline −0.04 m is larger in the MHFEM case. This is confirmed in Table 2 by ratios FEM_D/FEM and MHFEM_D/MHFEM between sensitivities and computed fluxes of 19.6% and 25.2%, respectively.

The same calculations were performed in the heterogeneous transmissivity case. A glance at Figure 6 suffices to see that the almost impervious barrier is efficient since a perturbation of the leakage coefficient has no impact outside the protected area. A comparison between Figures 5 and 6 shows that the presence of heterogeneities prevents the diffusion of the sensitivities (there is no propagation of the sensitivities beyond the barrier). Likewise, the computed flux sensitivity at the river is almost zero in Table 2. Even trends are identical; it can be noticed in Table 2 that the ratio between sensitivities and computed fluxes at the drain is twice larger in the MHFEM case (8.98%) than in the FEM case (4.42%).

Velocity fields and their sensitivity computed in the homogeneous transmissivity case are plotted in Figure 7. As expected from the difference in the exchange modeling (subsection Numerical methods), velocity vectors are larger close to the northern drain in the FEM case because more water is caught. One also observes a difference in the orientation of the velocity vectors at the east of the drain indicating that the FEM can catch water in this zone, contrarily to the MHFEM.

The sensitivities of the velocity fields are almost similar along the drain, but differ at the drain ends. This is in agreement with the two-drain modeling.

Velocities are smaller in the heterogeneous case (see Figures 7(a) and 8(a)), and sensitivities are even smaller. Two reasons explain this behavior. First, the hydraulic conductivity is subject to high variations in the domain (including small values). This heterogeneity has an important impact on the flow as it reduces its velocity in a significant manner. Second, the barrier mostly prohibits the flow down to the river. Important differences in the vector lengths of the two velocity fields are noticeable at the northern end of the drain because the MHFEM_D/MHFEM ratio is here twice as large as FEM_D/FEM ratio in Table 2.

Streamlines are computed by means of either a classical interpolation method (Cordes & Kinzelbach 1992) for FEM or the property for MHFEM. We compare streamlines computed using m⁻¹ d⁻¹ with approximated streamlines, both west/east oriented, built from heads and velocity fields and their sensitivity computed with the TL codes with parameter and perturbation values ⁠.

In the homogeneous transmissivity case (Figure 9) approximated streamlines are close to streamlines computed with m⁻¹ d⁻¹. One can notice some erratic behavior on the lower streamline of Figure 9. Although originating from the south-west of the drains, this particle is impacted by the north drain while running to the river. The erratic behavior occurs at a piezometric head maximum. This behavior is probably also favored by the small size of elements in this part of the computational domain. The particle hesitates between the drain and the river, finally reaching the river.

The same particles were used to initiate the streamline computation in the heterogeneous transmissivity case (Figure 10). The streamlines computed with m⁻¹ d⁻¹ are again well approximated by the approximated streamlines computed with ⁠.

The groundwater flow and drain–aquifer interactions are implemented in a conforming FEM and a MHFEM considering the drain as a river. The drain–aquifer exchange flux is modeled by means of the leakage concept and the Darcy's law where the leakage coefficient depends on the permeability, the thickness, and the leakage area of the drain bed. In this paper, the leakage coefficient per meter of drain is assumed to be constant. Numerical experiments have been carried out using an ideal homogeneous hydraulic conductivity over the domain, and an actual heterogeneous case study to measure the sensitivity of the computed piezometric head, velocity field, and streamlines to a perturbation of the leakage coefficient value. The proposed implementation of sensitivity analysis through AD is described with details as it may serve for other applications.

Fluxes computed at the channel, the river, and the drain boundaries have similar behavior but distinct values due to the difference in the modeling of waterbody exchanges. In any case, the drain catches more water using FEM. Discrepancies between the FEM solutions and the MHFEM solutions are difficult to evaluate from model computations. The accurate sensitivity analysis we propose provides much more information on piezometric head and velocity field, revealing the influence of a perturbation of the leakage coefficient in a clear manner, whatever the case study. The role of the impervious barrier is clearly stated since sensitivities are 0 outside the protected area. To the best of our knowledge, streamline reconstruction results are new in groundwater flow modeling. Although some erratic behavior is observed, the first-order approximated streamlines are very close to the computed streamlines.

Mohammad Moezzibadi was supported by the YEKAN Center Grant No. 139392351345000041.