Threedimensional hydrogeological parametrization using sparse piezometric data
 Institut Terre et Environnement de Strasbourg, Université de Strasbourg/EOST/ENGEES, CNRS UMR 7063, 5 rue Descartes, Strasbourg 67084, France
 Institut Terre et Environnement de Strasbourg, Université de Strasbourg/EOST/ENGEES, CNRS UMR 7063, 5 rue Descartes, Strasbourg 67084, France
Correspondence: Dimitri Rambourg (d.rambourg@unistra.fr)
Hide author detailsCorrespondence: Dimitri Rambourg (d.rambourg@unistra.fr)
When modelling contamination transport in the subsurface and aquifers, it is crucial to assess the heterogeneities of the porous medium, including the vertical distribution of the aquifer parameter. This issue is generally addressed thanks to geophysical investigations.
As an alternative, a method is proposed using estimated hydraulic parameters from a 2D calibrated flow model (solely reliant on piezometric series) as parametrization constraints for a 3D hydrogeological model. The methodology is tested via a synthetic model, ensuring full knowledge and control of its structure. The synthetic aquifer is composed of five lithofacies, distributed according to a sedimentary pattern, and functions in an unconfined regime. The level of heterogeneity for hydraulic conductivity spans 3 orders of magnitude. It provides the piezometric chronicles used to inverse 2D flow parameter fields and the lithological logs used to interpolate the 3D lithological model. Finally, the parameters of each facies (hydraulic conductivity and porosity) are obtained through an optimization loop, which minimizes the difference between the 2D calibrated transmissivity and the transmissivity computed with the estimated 3D facies parameters.
The method estimates values close to the known parameters, even with sparse piezometric and lithological data sampling. The maximal discrepancy is 45 % of the known value for the hydraulic conductivity and 6 % for the porosity (mean error 26 % and 3 %, respectively). Although the methodology does not prevent interpolation errors, it succeeds in reconstructing flow and transport dynamics close to the control data. Due to the inherent limitations of the 2D inversion approach, the method only applies to the saturated zone at this point.
To simulate contamination transport in the subsurface and aquifers, it is crucial to assess and reliably describe the heterogeneities of the porous medium. The development of inverse methods in recent decades is mainly based on twodimensional flow models and focused on the horizontal structure of heterogeneities with the collection of piezometric data as a cornerstone (Poeter and Hill, 1997; Carrera et al., 2005; Hendricks Franssen et al., 2009). But the latter is less sensitive to the vertical structure of the aquifer, leaving its estimation dependent on complex and expensive field methods, for example, pumping tests (De Caro et al., 2020), tracer tests (Linde et al., 2006), electrical resistivity (Coscia et al., 2011; Priyanka and Mohan Kumar, 2019), radar tomography (Boni et al., 2020), selfpotential methods (Eppelbaum, 2021), crosshole testing (Klotzsche et al., 2013; Doetsch et al., 2010), hydraulic tomography (SanchezLeón et al., 2015; Luo et al., 2020; Fischer et al., 2020), and/or laboratory analysis, for example grainsize analysis from core samples (Marini et al., 2018) and ex situ permeability tests (Zhang and Brusseau, 2005). The collection of this information, describing the vertical heterogeneity of the aquifer, allows for the development of 3D inversion techniques. For example, some successful methods combine direct parameter quantification and stochastic geological modelling (Guadagnini et al., 2004; Fu and GómezHernández, 2008; Cardiff and Kitanidis, 2009), and others incorporate water head data and more advanced geophysical measurements to the (joint) inversion procedure (Straface et al., 2011; Lee and Kitanidis, 2014).
Hydrogeological models are usually twodimensional, and transmissivities are estimated through model calibration. Twodimensional models are easier to handle considering fieldwork, parameter measurements, data manipulation, and calibration than 3D models. We question here the possibility of transferring data from a 2D calibrated hydrogeological model to a 3D configuration. Viaroli et al. (2019) recently employed a simplified 2D model to specify the boundary conditions and recharge of a 3D model already designed by other means. Incidentally, the design and parametrization of the 3D model itself are even more seldom independent of geophysical methods. In this line, we propose an original method using data from a 2D calibrated flow model (solely reliant on piezometric time series) as parametrization constraints for a 3D hydrogeological model resulting from interpolation of borehole data. The use of 2D calibrated transmissivities allows our technique to be completely unrelated to geophysical methods and less heavily computational than 3D joint inversion approaches. Moreover, the articulation of the method also allows a preexisting 2D calibrated model to be taken advantage of, if any.
The method is tested on a synthetic test case constituted by five hydrofacies, distributed according to a sedimentary pattern, with a level of heterogeneity for hydraulic conductivity spanning 3 orders of magnitude. This work can be considered an improvement of the method proposed by Harp et al. (2008), who also tested the combination of 2D inversion and an interpolation method but on a twodimensional transect model composed of only two facies.
In order to evaluate our methodology's robustness, it is first carried out with a very profuse data sampling (piezometric for the 2D inversion and lithological for the 3D model interpolation), assessing the consistency between the different numerical codes. Second, a sparser sampling is tested to approximate more realistic field conditions.
The detail of the methodology is described in Sect. 2, including the synthetic data framework, the mathematical background of the tools used, and the link between them. The results for both samplings concerning the inversions, the facies interpolation, and the final model outputs (in terms of parameter, piezometric series, and contamination plumes) are discussed in Sect. 3.
The methodology we propose and analyse in this paper is the following:

We propose estimates of transmissivity and porosity from a 2D calibrated flow model based on piezometric heads. These transmissivities exist for each element/cell of the 2D mesh. It is a huge dataset constrained by the piezometric heads.

We analyse the aquifer lithology at boreholes. Lithology is usually described at each borehole. It provides a qualitative description of aquifer heterogeneity. This qualitative description can be interpreted in terms of facies. This description is used here to define an optimal number of facies that have been identified within the aquifer.

The 3D discretization of the aquifer is a vertical extension of the 2D model, which avoids interpolation of the 2D parameters. Facies are interpolated over the 3D domain based on the borehole local data.

The hydraulic conductivity and porosity for each facies are estimated through optimization using the 2D data, which are considered to be vertical integrations of the 3D data. Optimization is required because the number of unknowns is quite small (twice the number of facies) compared to the number of constraints (twice the number of elements of the 2D flow model at the most). Of course, the constraints are correlated through the flow model and cannot be considered independent.
To evaluate this approach, we built a synthetic test case (Fig. 1) generated by the following:

a 3D aquifer design;

computation of the 3D flow using the software TRACES (Hoteit and Ackerer, 2004);

selection of representative head data, used as constraints for the 2D inversion after vertical averaging;

estimation of transmissivities and mean vertical porosity by a 2D flow model calibration based on the selected head data using PINOGRI (Rambourg et al., 2020);

selection of representative boreholes for lithological data and facies definition;

design of the 3D facies distribution using a geostatistical interpolator (GemPy; de la Varga et al., 2019) or a deterministic interpolator (splines; Lee et al., 1997);

estimation of each facies' hydrodynamic parameters (hydraulic conductivity, porosity) using an optimization procedure constrained by the 2D calibrated values;

comparison of local hydrodynamic parameters, simulated water heads, and concentrations between the “true” aquifer and the reconstructed (estimated) aquifer.
The computations are run on a PC with Intel(R) Core(TM) i76700 CPU @ 3.40 GHz processor and 16 GB RAM.
2.1 Synthetic threedimensional dataset
2.1.1 The aquifer model
The synthetic aquifer consists of five hydrogeological facies (also referred to as hydrofacies) distributed along a sedimentary pattern over a 10×10 km area and 20 m depth (Fig. 2).
Each hydrofacies is characterized by a hydraulic conductivity 5 times higher than the underlying facies (Table 1). Their porosity is less heterogeneous as it is defined in the range of permeable sedimentary materials (10 %–30 %). In practice, a hydrofacies is defined by clustering lithofacies with comparable hydrodynamic properties. The limitations and pitfalls inherent in this step are not addressed in this study, which is assumed to be flawless.
Hydraulic boundary conditions are of nullNeumann type (no flux) except the northwest corner where a 21 m head is imposed (Dirichlet boundary), acting as the outlet of the aquifer. A constant pumping well (18 m^{3} h^{−1}) is positioned in the southeast part of the model, intercepting the whole thickness. For solute transport, a zero solute flux is prescribed at the boundary, except at the aquifer outlet, where the outflow is considered purely advective.
^{1} Proportion of facies at the scale of the whole domain. ^{2} Proportion of facies within the surface elements.
The aquifer is exclusively fed by rainfall (720 mm yr^{−1} on average), which is assumed homogeneous over the whole area. However, five different recharge patterns are imposed according to the most superficial facies, whose hydrodynamic parameters greatly influence the amount and dynamics of water infiltration. Thus, recharge zone 5 (formed by the most permeable surface facies) is subject to major and fast infiltration, in contrast to zone 1 (least permeable), where the seepage signal is very attenuated and spread over time (see Sect. 3).
2.1.2 Piezometric and contamination reference data generation
The behaviour of groundwater and dissolved contamination is computed using TRACES (Transport of Radioactive Elements in Subsurface) software (Hoteit and Ackerer, 2004), a numerical code written in FORTRAN 90 for the simulation of flow and reactive transport in saturated/unsaturated porous media.
The threedimensional flow model is the combination of the conservation of mass and Darcy's laws, generalized to also apply to the unsaturated zone (Darcy–Buckingham law), resulting in the Jacob–Richards equation (Eq. 1):
where θ and ϕ are the water content [–] and porosity [–], respectively, necessary to deal with the unsaturated zone. s and K are the specific storage coefficient [m^{−1}] and hydraulic conductivity tensor [m s^{−1}], respectively. h is the water head [m], and f is the sink–source term [s^{−1}].
To limit inconsistencies with the 2D inversion (where unsaturated flow is not addressed via a physical model), the 3D model is reduced to a fully saturated approach. Therefore, the flow equation (Eq. 1) is simplified and shown with the adequate initial and boundary conditions as Eq. (2).
where S is the storativity [–], equivalent to the effective porosity in an unconfined context. T is the transmissivity [m^{2} s^{−1}], the integration of the hydraulic conductivity over the vertical of the model. F [m s^{−1}] is the sink–source term, x is a position in Ω, the model domain, and h_{0}(x) represents the initial conditions. Γ_{D} and Γ_{N} are partitions of the domain boundaries that correspond to Dirichlet and Neumann conditions, respectively, and n is the unit vector normal to the boundary, counted positive outward. h_{D}(x,t) is the prescribed head value at the Dirichlet boundaries, and q_{N}(x,t) is the prescribed flux at the Neumann boundaries, both defined at each time t of the simulated period T.
The assumption of a locally constant transmissivity is satisfied, with water head variations of a maximum of 5.2 % (and 3.5 % on average) of the local mean water head.
TRACES addresses the migration of contaminants via an dispersion–diffusion equation, supporting adsorption, precipitation, and degradation (reactive transport) phenomena. However, the study considers one inert species, giving form to Eq. (3).
where C is the solute concentration [kg m^{−3}], D is the dispersion–diffusion tensor [m^{2} s^{−1}], q is Darcy's velocity [m s^{−1}], and Q is the solute sink–source term [kg m^{−3} s^{−1}]. C_{0}(x) is the initial concentration; and A(t), B(t), and q(t) are the parameters to define the boundary conditions (see Hoteit and Ackerer, 2004).
The dispersion and diffusion parameters are set identically for all the facies (Table 2).
These parameters are transferred into the dispersion–diffusion tensor as shown in Eq. (4).
where D_{m} is the molecular diffusion coefficient [m^{2} s^{−1}], and τ is the tortuosity factor of the porous medium [–] (according to Millington and Quirk, 1961). D_{T} and D_{L} are the transversal and longitudinal dispersivities [m], respectively, while δ is the Kronecker function, with i and j the position indexes in the tensors.
Flow and transport equations implemented in the code TRACES are solved under transient or steadystate computation in 2D or 3D heterogeneous domains. Mixed hybrid finite elements are used to solve the flow equation and the diffusive–dispersive components of the transport. A mass lumping formulation is used to limit the occurrence of numerical oscillations. The advective part of the transport is solved using a discontinuous Galerkin finite element method, which also prevents numerical oscillations in the simulations and strongly limits numerical diffusion. These numerical schemes ensure an exact mass balance at the element level and are very flexible in space discretization (triangular or quadrangular element in 2D, tetrahedral, prismatic of hexahedral elements in 3D).
As shown in Fig. 2, the model consists of 10 440 triangular prisms, 6193 nodes, and 27 544 facies. The horizontal edges have a characteristic length of 500 m, while the vertical edges are 2 m long.
The initial state of the water table is derived from a preliminary steadystate calculation involving averaged recharge. The aquifer is initially uncontaminated (C_{0}=0 over the whole domain) and undergoes a pollution episode from three surface sources (Fig. 2), each discharging 0.1 g s^{−1} for the first 24 h of the simulated time.
2.1.3 Structural and piezometric reference data sampling strategies
Structural (hydrofacies) and piezometric data are sampled following two subsequent strategies in order to validate the method (Fig. 3).
First, the methodology is conducted with a very dense dataset (400 control points) to assess its potential under ideal conditions and verify the numerical approaches' compatibility. The piezometric chronicles used for the inversion cover the 9 years of simulation (see Sect. 3).
Second, a sparser dataset is extracted to evaluate the method in more realistic conditions. The control points are reduced to 40, and the piezometric chronicles are shortened randomly (down to 2 years). Meanwhile, the hydrofacies logs extracted at the location of the control points are kept in their integrity.
The sparse sampling is pseudorandom so that each recharge area is covered. The points at each corner of the model are included in both sampling strategies to avoid extrapolation issues. The sparse sampling accounts for 10 % of the lithological information of the dense dataset and only 6 % of the water head data.
2.2 Twodimensional flow model inversion
The sampled piezometric data are used as inversion constraints for the PINOGRI (Parameter Inversion Numerically Optimized for Groundwater Issues; see Rambourg et al., 2020) software, developed at ITES (Strasbourg). As water heads are less sensitive to the vertical heterogeneities of the porous media, the inversion approach is restricted to a twodimensional scale, where heads are vertically averaged. This step results in the estimation of transmissivity and average porosity fields at the scale of each mesh of the model. The inversion procedure consists of minimizing an objective function (the quadratic difference between measured and computed piezometric heads) with parameter optimization guided by a gradient descent method.
Although piezometric data are subject to uncertainty in a field context, we do not address this aspect in the present study, and the water heads' measurement errors are considered negligible.
2.2.1 The flow model
Twodimensional groundwater flow in the aquifer is described by a diffusiontype equation, akin to the TRACES approach, but with a constant head over depth assumption (Dupuit–Forchheimer's hypothesis), reducing the problem's dimension.
The mathematical model is solved by a twodimensional nonconforming finite element method (Crouzeix and Raviart, 1973), ensuring flux continuity, mass balance (like the finite volume method), flexibility in geometry, and rigorous computation of full tensor transmissivity (like conforming finite elements, as stated by Ackerer et al., 2014). The time discretization scheme is implicit, giving the direct problem the form of Eq. (5).
where h is the water head vector at the calculation time t, and F is the sink–source term vector produced at the previous time step. A is the flow coefficient matrix, depending on the mesh geometry and the parameter vector.
2.2.2 The inverse problem
The groundwater flow parameters are estimated through the minimization of an objective function (Eq. 6) based on weighted least squares (Carrera and Neuman, 1986; Tarantola, 2005).
where J is the objective function, P represents the vector of the parameters to be estimated, h^{∗} is the measured piezometric head (obtained from vertical averaging of 3D sampled data), and h is the corresponding simulated values. T is the transpose operator, and W is the weighting matrix, depending on measurement errors, able to prioritize optimization effort on specific locations. Therefore, in this study, all data are weighted equally. Moreover, because no a priori hydraulic parameters information is added in the study, the objective function does not include the plausibility criterion of the maximum likelihood approach.
Due to the great number of parameters and measurements, the minimization of the objective function is led by a quasiNewton method, which is less timeconsuming compared to Gauss–Newton and other Jacobianbased approaches (Kitanidis and Lane, 1985). The gradient and an approximate Hessian of the objective function are calculated using the discrete adjoint state method (Carter et al., 1974) and the limitedmemory BFGS (Broyden–Fletcher–Goldfarb–Shanno) algorithm (Byrd et al., 1995), respectively. In our case, the adjoint state method is used to compute the gradient of the objective function (required by the LBFGS algorithm) as an optimization problem of a Lagrangian, constrained by the head values obtained from the direct calculation. On the other hand, instead of calculating the sensitivity coefficients for each parameter at each iteration required by Newton methods, the LBFGS algorithm (quasiNewton) approaches a Hessian approximation by converging an initial matrix (e.g. the identity matrix) according to the results from a limited number of previous iterations. As the parametrization of the inversion can lead to a high number of degrees of freedom, this set of techniques has been found more efficient than standard sensitivity approaches (Townley and Wilson, 1985). Finally, three stopping criteria are set to end the algorithm: (i) the objective function J or its gradient is sufficiently low, (ii) the adjustment of the parameters P or the decrease of J between two iterations is too small, or (iii) the number of iterations has reached a userset maximum value.
Incidentally, the inverse problems generally suffer from being illposed; i.e. the number of data (locally known piezometry) is too low compared to the number of unknowns (hydraulic conductivity and porosity at the scale of each mesh element). This leads to issues of nonuniqueness and instability of solutions. One way to limit these inconveniences is to reduce the number of unknowns via a parametrization technique. In PINOGRI, the parameter spatial pattern is inferred using an adaptive multiscale triangulation (AMT; Majdalani and Ackerer, 2011; Hassane and Ackerer, 2017). The parameters, borne by the vertices of the AMT mesh, are interpolated into each element of the calculus mesh (see Fig. 4). If during the inversion process, the minimization criteria are not met at the scale of a parameter cell, the latter is divided into four, increasing the optimization's degrees of freedom. Refinement is halted either when the objective function at the element level drops below a userdefined threshold, when the number of iterations reaches a userdefined maximum, or when the last iteration fails to produce a better optimization than the previous one. This adaptive approach allows for more flexibility and needs fewer preconceptions about the model structure than fixed parametrizations, such as zonation or interpolations. A detailed description of the mathematical developments and the algorithm can be found in Ackerer et al. (2014) and Hassane and Ackerer (2017).
A maximum of three adaptive multiscale iterations is set to ensure a satisfactory calibration while preventing overparametrization.
The boundary conditions of the 2D approach are exactly the same as for the 3D synthetic model. In contrast, the initial conditions cannot be integrally transposed, the knowledge of the water head being limited by data sampling. Thus, the initial water head for the 2D approach is derived from preliminary steadystate inversions constrained by timeaveraged water head data.
The parameter bounds of the inversion are $\mathrm{5}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}$ for the hydraulic conductivity [m s^{−1}] and 6 %–30 % for the porosity.
2.3 Facies interpolation
For comparison purposes, two interpolation methods are used to reconstruct the 3D facies distribution, both based on sampled lithological logs. The interpolation estimates the distribution between control points, using statistical dependency of measures in the case of geostatistical methods or drawing geometric surfaces independently of the spatial statistical repartition in the case of deterministic methods. For both approaches, the resulting 3D models contain only qualitative (indicator) data. The interpolation methods selected are not exclusive, as the general methodology presented can accommodate any interpolation tool that produces a 3D facies model.
2.3.1 GemPy (geostatistical interpolation)
GemPy (de la Varga et al., 2019) is an opensource 3D geomodelling package written in Python. It specializes in the reconstruction of stratigraphic series, with the possibility of modelling complex environments by adding faults, folds, plutonic intrusions, and other anomalies. GemPy's mathematical background is a development of the work of Lajaunie et al. (1997; Calcagno et al., 2008), using universal cokriging methods to interpolate potential fields (scalar fields).
Kriging covers a set of exact (unbiased) linear estimation techniques that minimize the estimation variance using the variogram, a function representing the correlation level of a random variable as a function of distance. Initially limited to stationary variables (simple and ordinary kriging), universal kriging has extended the use of this type of geostatistical methods to nonstationary variables. Eventually, cokriging not only uses the spatial correlation of a variable with itself, but also incorporates the crosscorrelations between two or more random variables.
In the software, the interpolation concerns two types of data: isosurfaces (including the interface between stacked lithologies and the boundaries of fault planes or unconformities) on the one hand and surface orientation (the gradient of the said isosurfaces) on the other hand. This last source of data allows for the computation of a very smooth and continuous sedimentary structure, which is rarely the case in other freeware geostatistical tools (dell'Arciprete et al., 2012; Langousis et al., 2017). In our case, the orientations (geological poles) are obtained by calculating the normal of the planes, defined by triangulation between the hydrofacies interfaces at the sampled data points.
Being specialized in geological modelling, GemPy handles the secondorder (weak) stationarity of universal kriging by assuming a linear trend in the mean value of the scalar field. In addition, the random function defined for universal cokriging does not bear any physical meaning as it only aims at ensuring equality at every point of the isosurface (no matter the value).
As a result, the crossvariogram, inherent to cokriging, cannot be empirically determined. The shape of the surfaces mainly depends on the orientations provided and on an arbitrary spherical covariance function that only balances the relative weight of the surfaces and their orientation in the cokriging. Hence, the variogram parameters do not bear any physical meaning as well and are arbitrarily chosen to ensure stability to the computation according to the GemPy's developers' guidelines (De la Varga et al., 2019): the nugget effect should be small (set to 10 in our case) and the range equal to the domain's extension (10 000 m in our case). As the variogram is not differentiated according to the search direction, the vertical component of the model must be exaggerated (×500 in our case) so that its dimension is compatible with the previously quoted values.
Finally, GemPy produces a 3D facies model made of $\mathrm{50}\times \mathrm{50}\times \mathrm{10}$ hexahedron elements. After rescaling on the z direction, it has the same extension as the mesh used for the other models and a finer resolution. Henceforth, the facies in the flow/transport mesh for TRACES are determined according to the majority facies of the GemPy elements intersecting each 3D prismatic element.
The calculation of both the scalar fields and their derivatives is handled by the Theano Python library, which also allows for developments toward stochastic modelling. A more precise description of the software is available in De la Varga et al. (2019).
2.3.2 Bsplines (deterministic interpolation)
Spline methods are also suitable for the construction of sedimentary models characterized by smooth surfaces. By definition, their interpolation adjusts continuous polynomial equations to the data, ensuring no discontinuities and exact fitting (Prautzsch et al., 2002). Splines can be assimilated to flexible surfaces constrained to fit the observation values while minimizing their bending energy. Contrary to a simpler deterministic method (e.g. trend surface) that operates via a single polynomial equation, splines represent the surface in pieces and therefore require the computation of a large number of equations.
However, this method is chosen for its ability to reproduce smooth surfaces, compatible with a sedimentary morphology, and can be easily carried out through a GIS procedure (QGIS/SAGA multilevel Bspline interpolation; Lee et al., 1997). To avoid anomalies in the stacking of the facies, the interpolation is carried on their thickness instead of their boundaries' z coordinates. In addition, the first underlying facies is not interpolated but considered to be the background (filling) lithology. The thicknesses of the four remaining facies are delivered in raster format, with integer values between 0 and 10 (i.e. the number of layers in the final 3D model) and a resolution of 200 m. Eventually, the facies stacking is transcribed for each column of prismatic elements in the 3D flow/transport mesh for TRACES according to the same majority analysis as for the GemPy procedure.
2.4 Hydrofacies parametrization and 3D simulations
The lithological models resulting from the interpolations do not have any assigned hydrodynamic parameters. Thus, an optimization procedure is implemented to find the hydraulic parameters of the five facies by minimizing the quadratic difference between 2D and 3D estimated transmissivities and porosities (Fig. 1, Eq. 7). Both previous steps of the methodology draw continuous data over the modelled domain (2D averaged parameters on the one hand and lithological structure on the other). Conceptually, the optimization could be performed with as many constraints as the number of elements of the mesh. However, the inversion and interpolation errors are expected to be minimal at the sampled data location. Therefore, the algorithm is carried out only with the parameter values and the lithological successions in these locations, minimizing uncertainties related to lack of sensitivity for transmissivity values or related to interpolation for lithological data.
The optimization is handled thanks to a Levenberg–Marquardt algorithm, whose unknowns are the hydraulic parameters (porosity and hydraulic conductivity) for each facies, i.e. 2×5 unknowns over the all domain. The constraints are the 2D mean values (transmissivity and porosity) at the sampled locations. In order to integrate the least number of preconceptions in the method, the bounds of values within which the algorithm can pick during the optimization are not differentiated by facies (the bounds are 10^{−6} and 10^{−2} m s^{−1} for the hydraulic conductivity, 2 % and 50 % for the porosity). The objective function of the optimization problem takes the form of Eq. (7).
where O is the objective function, i is the index for the constraint (i.e. the sampled location retained for the optimization), j is the index for each facies, and l_{i,j} is the thickness [m] of facies j at location i. p represents the parameters to be optimized (the hydraulic conductivity or the effective porosity of each facies) and P the 2D mean values calibrated during the inversion stage, weighted by the matrix σ representing this calibration uncertainty. We consider only the diagonal of the matrix, containing the inverse of the variance given at location i by the 2D calibration.
The final uncertainty of the optimized parameters is given by Eq. (8).
where ϵ_{p} is the uncertainty of the parameter p, $\widehat{O}$ is the objective function at end of the optimization, and m is the number of data. The coefficient φ is determined through Fisher's distribution, assuming a normal distribution of the uncertainty (for an estimation at 95 % of confidence, φ=1.96). C_{p} is the variance of the parameter p, derived from the Jacobian (sensitivity matrix) of the model.
Once the optimization estimated each facies' hydraulic parameters, the 3D model is parametrized. Flow and contamination simulations are carried out with TRACES, as described previously, with the new facies distribution and the new parameter set. The boundary conditions and the recharge distribution are kept unchanged from the reference synthetic model. However, the initial state data are directly taken from the 2D simulations.
3.1 Calibrated 2D models
The inversion algorithm is run 80 times for each sampling case to gather a set of possible solutions to the inverse problem, a model inversion lasting approximately 1 h on average in these conditions. Each set of solutions is called a batch.
To avoid overparametrization, the adaptive refinement of the parameter grid, initially composed of 21 nodes (i.e. degrees of freedom for the minimization), is limited to three refinements. The number of new parameter located at the parameter grid vertices is also constrained by the number and location of the piezometric control points. Therefore, the final average number of parameter nodes is 497 for the dense sampling (400 control points) and 74 for the sparser sampling (40 control points).
Each solution batch produces very stable parameter fields (Table 3) and piezometer chronicles with a mean absolute discrepancy of less than 5 mm compared to the sampled reference data. For the sparse sampling, the mean absolute error increases to 37 cm when all control points from the dense sampling are included in the evaluation.
The mean relative standard deviation of the parameter at the element mesh scale stays at a very low level in both cases. The variability of transmissivity is even lower for the sparse sampling (0.2 % vs. 0.5 %), as it has fewer degrees of freedom.
The consistency between the estimated and the reference parameters is described in Sect. 3.3.
3.2 Threedimensional interpolations
As the reference model is shaped according to a simple sedimentary pattern (absence of faults), both interpolation methods (geostatistical and deterministic) produce results of similar quality.
Differences in the facies composition of the models are marginal, even with a sparse distribution of the conditional data (Fig. 5).
The two percentages accompanying each model in Fig. 5 represent the proportion of elements incorrectly parameterized for the whole dataset and the conditional data, respectively. With a dense conditional data sampling, the deterministic approach yields slightly better results than the geostatistical one (2.6 % vs. 4.8 % of elements parametrized with the wrong facies). The GemPy algorithm handles sparser constraints slightly better (9.5 % vs 11.0 % of errors). However, the differences between the two interpolation methods can be considered small and assumed to be dependent on the case study.
3.3 Parameter comparison
The results of the inversions (mean 2D values and their associated variance) and the known facies distributions at the sampling location are used to optimize 3D hydrodynamic parameters as explained in Sect. 2.4. These optimized values and their uncertainties are shown in Fig. 6.
In all cases, the hydraulic conductivity of each facies stays in the same order of magnitude as the reference data (the mean errors account for 26 % of the known hydraulic conductivity values). The largest discrepancy (0.26 log units, 45.5 % of the reference value in a linear scale) affects the facies 4 in the sparsesamplingbased estimation. The gaps between the data and the estimated porosity are also very low, with mean and maximal absolute errors of 1.3 % and 2.9 % (i.e. 2.9 % and 6.4 % of the data values). These consistent results are due to the fact that the parametrization optimization is carried out with input at the sampled data location, where the 3D facies vertical succession is known and where the errors on the 2D mean parameters are low.
The uncertainty attached to the optimization varies with the density of the data. The values related to the dense sampling are negligible (the highest uncertainties are at 0.02 log units for the hydraulic conductivity and 0.5 % for the porosity); therefore only the ones of the sparse sampling are shown. At the most, the uncertainty extends over 0.44 log units for the hydraulic conductivity (facies 2) and 3.3 % for the porosity (facies 5). Consequently, the confidence intervals almost always include the corresponding reference value.
The conjunction of the parameter discrepancies and facies misplacements results in transmissivity errors, as shown in Fig. 7. The transmissivity discrepancies are always below 1 order of magnitude and mostly below 0.2 log units. The largescale heterogeneities are very well reproduced in every case, notably in the models based on sparse data. Comparisons between the 2D parameters and the 3D parameters show that the estimation errors in the former do not propagate entirely in the latter. Indeed, the errors of the inverse models are not only located at the interfaces between the largescale heterogeneities, but also at smaller scales, within the heterogeneities, when the density of the constraint data is reduced (in the absence of local piezometric data, the hydrodynamic parameters are less constrained; see bottom left of Fig. 7). The 3D interpolation techniques are not subject to these smallscale errors as they generate very smooth and continuous facies distributions (Fig. 5). Therefore, the final discrepancies are mainly at the transitions between the largescale horizontal heterogeneities where the facies interpolation generates errors (in particular, the overestimation of transmissivity is visible where the interpolations have extended facies 1 in excess).
3.4 Piezometric head comparison
In view of the few parametrization errors produced with the dense sampling, the results with respect to piezometry are shown only for the sparse sampling. In addition to the simulations incorporating the optimized hydrodynamic parameters per facies, four additional runs are conducted for each interpolation method in order to study the propagation of parameter uncertainties. For each additional run, the parameter values are set as follows: (i) at the upper bound of the confidence interval, (ii) at the lower bound, (iii) alternately at the lower (for the facies 1, 3, 5) and upper (facies 2, 4) bounds, and (iv) alternately at the lower (for the facies 2, 4) and upper (facies 1, 3, 5) bounds of the parameter confidence interval. These simulations generate chronicles and maps whose extreme values are retained to construct “envelope curves”, showing the final uncertainty of the piezometry (Fig. 8).
The final piezometric state for each model (water head averaged on the vertical) is shown in the central map of Fig. 8. The piezometric contours are well reproduced, especially for the 40 and 30 m isolines. Differences between the interpolation strategies are marginal. The main discrepancy between the models and the reference data is visible in the middle of the right border of the domain, where the interpolations underestimated the presence of facies 1, in favour of facies 2. However, the narrow confidence interval on the parameters of facies 1 is reflected by a low uncertainty on the piezometry in the lower right corner of the domain, where this facies predominates. In contrast, the uncertainty in piezometry is significant where facies 2, 3, and 4 predominate.
The water head fluctuations are also consistent with the reference chronicles over the entire period of simulation and for every recharge zone. Amongst the chosen piezometers, those identified with an asterisk were not used as constraint data in the sparse sampling. The mean absolute errors of each model are 44.6 cm for the GemPy model and 47.6 cm for the Bsplines. The deviations mainly take the form of a shift (by excess or by default) in the base level when the fluctuations are well reproduced. Overall, the parametrization discrepancies (Fig. 6) are too small to significantly modify the flow dynamics. Combined with the small differences in the model's composition (Fig. 5), the water head equilibrium is slightly shifted, as shown in the charts. Therefore, a more significant deviation would occur in a permanent flow simulation (with averaged recharge). As with the contour map, the highest levels of uncertainty are for piezometers intercepting significant thicknesses of facies 2 (P2, P3, P6, and P7), this one having the widest confidence intervals.
3.5 Contamination comparison
The study of contaminant transport is another prime use of hydrogeological models. The outputs of the transport simulations are shown in Fig. 9, in the form of breakthrough curves and isoconcentration maps (10 mg L^{−1}) at different times of the simulation. As for the piezometric data, envelope curves of uncertainty are deduced from the four simulations involving the parameters at the bounds of the confidence intervals.
The models' output errors are attributable to parametrization discrepancy on individual facies, interpolation misplacements, and the accumulation of these same errors upstream of each surveyed location. The results confirm that contaminant transport is much more sensitive to parametrization errors than piezometry. Indeed, the latter is governed by the transmissivity, where individual facies misparametrization can be buffered by the vertical integration. For contaminant transport, each voxel parametrization may influence the outcome.
For instance, source 1 is located on facies 1 in the reference model and on facies 2 in the interpolated models. Therefore, the dynamics of the breakthrough curve and the pollutant plume are significantly different (i.e. a lower spike due to a higher porosity and a faster depletion due to a higher conductivity). Conversely, the facies distribution is preserved in the location of source S2 (in facies 4), where the discrepancies are mainly due to hydrodynamic parametrization errors on this very same facies.
Overall, the results at the outlet E and the plume maps show that, in our case, the parametrization errors have a more significant impact on the fate of sources S2 and S3 than for source 1.
The confidence intervals on the facies parameters (especially on facies 2, 3, and 4) generate scenarios where the plumes disappear early on the one hand or extend in excess on the other hand.
3.6 Perspectives regarding uncertainties
In this study, only the uncertainty related to the calibration of the 2D parameters is calculated and propagated to the 3D models. However, it must be emphasized that in the context of distributed hydrogeological models, many other sources of uncertainty occur (Pechlivanidis et al., 2011).
First, there are data uncertainties: piezometric measurements, rainfall and radiative data (leading to recharge estimates), and lithological descriptions (both in terms of their categorization and altimetry) are all subject to error. Our approach via a synthetic case encouraged us to postpone this aspect and to concentrate on the analysis of the methodology. These uncertainties will be taken into account in future work dealing with real cases.
Second, all uncertainties related to the parameters have not been addressed in this study. Focused on the determination of hydrodynamic parameters, it did not integrate the uncertainty related to the dispersivity and molecular diffusion parameters inherent to the contaminant transport phenomena.
Third, the structure of the models and the way they are defined also involve uncertainties. In particular, interpolation methods produce uncertainty and propagate those inherent in the lithological data (Phillips and Mark, 1996; Lloyd and Atkison, 2001; White, 2017). Our data sampling strategy allowed this pitfall to be circumvented, but it becomes unavoidable when there are significant gaps in the data.
Monte Carlo simulations (Wagener and Kollat, 2007; Beven, 2009) allow for the joint estimation of these uncertainties, at great computational cost (the simulations must cover a range of value for each input or parameter at stake). Several random (Olsson and Sandberg, 2002) or Bayesian (Vrugt et al., 2009) sampling strategies are at hand in order to reduce the number of iterations needed to obtain a representative view of the uncertainties.
Integrating this type of analysis into the framework of our method is one the important improvements planned.
A method has been developed to assess 3D aquifer parameters by combining hydrodynamic parameters estimated by a 2D model calibration and 3D facies interpolation. While direct 3D parameter estimation is generally based on a heavy geophysical survey, the proposed methodology is based solely on piezometric series and geological logs (commonly available at the same locations). Both 2D flow model calibration and 3D interpolation parts of the algorithm are independent. Therefore, the approach is not restricted to the tools described in the article (i.e. other interpolation methods than GemPy and Bsplines can be used) and can potentially incorporate a preexisting 2D model.
The synthetic test carried out with a relatively sparse dataset yields a consistent hydrodynamic parametrization (highest discrepancy: 45.5 % of the initial value for hydraulic conductivity, 6.4 % for effective porosity) and quite low errors in facies distribution (11 % of misplaced facies at the most). Subsequently, the reconstructed piezometric series show very consistent dynamics, with a maximal mean difference of 47.6 cm, mainly due to shifting the base level, while the fluctuations and the hydraulic gradients are generally unaltered. The discrepancies concerning the transport simulations are more significant, the phenomenon being more sensitive to parametrization errors at the individual voxel scale.
Comparatively to joint inversion methods, the need of data acquisition and the computation efforts are lower. However, in a field context, the method is very dependent on the characterization of the hydrofacies and the quality of the piezometric survey. Indeed, if uncertainties related to the 2D flow calibration were propagated to the 3D parameter optimization and, thereafter, into the piezometry and contaminant transport simulation, other sources of uncertainty hindering the hydrogeological modelling process were not accounted for at the time of publication. Among them, we can cite the uncertainties related to input data, transport parameters, or the structure of the model itself (especially the lithology interpolation step).
Another pitfall of the method, inherent to the 2D step, lies in the fact that lowhydraulicconductivity facies may be masked by more permeable facies in the transmissivity term, making their parametrization somehow difficult. Also, some other important modelling points are not addressed in the study, for example, the vadose zone dynamics and the transport parameters, which require separate estimates. Following this synthetic case, we plan to test the method on a real case, in order to confirm its operational potential, completed by comprehensive sensitivity and uncertainty analyses.
TRACES, PINOGRI, and the synthetic case dataset can be provided on request to the corresponding author. GemPy is an opensource 3D geomodelling package (https://www.gempy.org, last access: 5 December 2022). The Bspline interpolation method used is an addin tool of the free and opensource geographic information system QGIS/SAGA.
DR and PA designed the method workflow, DR designed the synthetic case and carried out the simulations. RDC adapted the GemPy Python package, and PA developed TRACES (with coauthors cited in the references) and codeveloped PINOGRI with DR (and previous contributors to the code development, as cited in the references). DR prepared the manuscript with contributions from both coauthors.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This paper was edited by Alberto Guadagnini and reviewed by two anonymous referees.
Ackerer, P., Trottier, N., and Delay, F.: Flow in doubleporosity aquifers: Parameter estimation using an adaptive multiscale method, Adv. Water Resour., 73, 108–122, https://doi.org/10.1016/j.advwatres.2014.07.001, 2014.
Beven K.: Environmental modelling – An uncertain future?, Routledge, London, 310 pp., ISBN 9780415457590, 2009.
Boni, R., Meisina, C., Teatini, P., Zucca, F., Zoccarato, C., Franceschini, A., Ezquerro, P., BéjarPizarro, M., FernàndezMerodo, J. A., GuardiolaAlbert, C., Pastor, J. L., Tomás, R., and Herrera, G.: 3D groundwater flow and deformation modelling of Madrid aquifer, J. Hydrol., 585, 124773, https://doi.org/10.1016/j.jhydrol.2020.124773, 2020.
Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C.: A limited memory algorithm for bound constrained optimization, J. Sci. Comput., 16, 1190–1208, https://doi.org/10.1137/0916069, 1995.
Calcagno, P., Chilès, J.P., Courrioux, G., and Guillen, A.: Geological modelling from field data and geological knowledge: Part I. Modelling method coupling 3D potentialfield interpolation and geological rules, Phys. Earth Planet. In., 171, 147–157, https://doi.org/10.1016/j.pepi.2008.06.013, 2008.
Cardiff, M. and Kitanidis, P. K.: Bayesian inversion for facies detection: An extensible level set framework, Water Resour. Res., 45, W10416, https://doi.org/10.1029/2008WR007675, 2009.
Carrera, J. and Neuman, S. P.: Estimation of aquifer parameters under transient and steady state conditions: 1. maximum likelihood method incorporating prior information, Water Resour. Res., 22, 199–210, https://doi.org/10.1029/WR022i002p00199, 1986.
Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., and Slooten, L. J.: Inverse problem in hydrogeology, Hydrogeol. J., 13, 206–222, https://doi.org/10.1007/s1004000404047, 2005.
Carter, R. D., Jemp, L. F., Pierce, A. C., and Williams, D. L.: Performance matching with constraints, Soc. Petrol. Eng. J., 14, 187–196, https://doi.org/10.2118/4260PA, 1974.
Coscia, I., Greenhalgh, S. A., Linde, N., Doetsch, J., Marescot, L., Günther, T., Vogt, T., and Green, A. G.: 3D crosshole ERT for aquifer characterization and monitoring of infiltrating river water, Geophysics, 76, G49–G59, https://doi.org/10.1190/1.3553003, 2011.
Crouzeix, M. and Raviart, P.A.: Conforming and nonconforming finite element methods for solving the stationary Stokes equations, Revue Française d'automatique, informatique, recherche opérationnelle, Mathématique, 7, 33–75, 1973.
De Caro, M., Perico, R., Crosta, G. B., Frattini, P., and Volpi, G.: A regionalscale conceptual and numerical groundwater flow model in fluvioglacial sediments for the Milan Metropolitan area (Northern Italy), J. Hydrol., 29, 100683, https://doi.org/10.1016/J.EJRH.2020.100683, 2020.
de la Varga, M., Schaaf, A., and Wellmann, F.: GemPy 1.0: opensource stochastic geological modeling and inversion, Geosci. Model Dev., 12, 1–32, https://https://doi.org/10.5194/gmd1212019, 2019 (data available at: https://www.gempy.org, last access: 5 December 2022).
dell'Arciprete, D., Bersezio, R., Felletti, F., Giudici, M., Comunian, A., and Renard, P.: Comparison of three geostatistical methods for hydrofacies simulation: a test on alluvial sediments, Hydrogeol. J., 20, 299–311, https://doi.org/10.1007/s1004001108080, 2012.
Doetsch, J., Linde, N., Coscia, I., Greenhalgh, S. A., and Green, A. G.: Zonation for 3D aquifer characterization based on joint inversions of multimethod crosshole geophysical data, Geophysics, 75, G53–G64, https://doi.org/10.1190/1.3496476, 2010.
Eppelbaum, L. V.: Review of processing and interpretation of selfpotential anomalies: Transfer of methodologies developed in magnetic prospecting, Geosciences, 11, 194, https://doi.org/10.3390/geosciences11050194, 2021.
Fischer, P., Jardani, A., and Jourde, H.: Hydraulic tomography in coupled discretecontinuum concept to image hydraulic properties of a fractured and karstified aquifer (Lez aquifer, France), Adv. Water Resour., 137, 103523, https://doi.org/10.1016/j.jhydrol.2020.125438, 2020.
Fu, J. and GómezHernández, J. J.: A blocking Markov chain Monte Carlo method for inverse stochastic hydrogeological modeling, Math. Geosci., 41, 105–128, https://doi.org/10.1007/s1100400892060, 2009.
Guadagnini, L., Guadagnini, A., and Tartakovsky, D. M.: Probabilistic reconstruction of geologic facies, J. Hydrol., 294, 57–67, https://doi.org/10.1016/j.jhydrol.2004.02.007, 2004.
Harp, D. R., Dai, Z., Wolfsberg, A. V., Vrugt, J. A., Robinson, B. A., and Vesselinov, V. V.: Aquifer structure identification using stochastic inversion, Geophys. Res. Lett., 35, L08404, https://doi.org/10.1029/2008GL033585, 2008.
Hassane, M. M. and Ackerer, P.: Groundwater flow parameter estimation using refinement and coarsening indicators for adaptive downscaling parameterization, Adv. Water Resour., 100, 139–152, https://doi.org/10.1016/j.advwatres.2016.12.013, 2017.
Hendricks Franssen, H. J., Alcolea, A., Riva, M., Bakr, M., van der Wiel, N., Stauffer, F., and Guadagnini, A.: A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterisation of well watchments, Adv. Water Resour., 32, 851–872, https://doi.org/10.1016/j.advwatres.2009.02.011, 2009.
Hoteit, H. and Ackerer, P.: TRACES user's guide, IMFS, Strasbourg, France, 2004.
Kitanidis, P. K. and Lane, R. W.: Maximum likelihood parameter estimation of hydrologic spatial processes by the GaussNewton method, J. Hydrol., 79, 53–71, https://doi.org/10.1016/00221694(85)901817, 1985.
Klotzsche, A., van der Kruk, J., Linde, N., Doetsch, J., and Vereecken, H.: 3D characterization of highpermeability zones in a gravel aquifer using 2D crosshole GPR fullwaveform inversion and waveguide detection, Geophys. J. Int., 195, 932–944, https://doi.org/10.1093/gji/ggt275, 2013.
Lajaunie, C., Courrioux, G., and Manuel, L.: Foliation fields and 3D cartography in geology: Principles of a method based on potential interpolation, Math. Geol., 29, 571–584, https://doi.org/10.1007/bf02775087, 1997.
Langousis, A., Kaleris, V., Kokosi, A., and Mamounakis, G.: Markov based transition probability geostatistics in groundwater applications: assumptions and limitations, Stoch. Env. Res. Risk A., 32, 2129–2146, https://doi.org/10.1007/s004770171504y, 2017.
Lee, J. and Kitanidis P. K.: Largescale hydraulic tomography and joint inversion of head and tracer data using the Principal Component Geostatistical Approach (PCGA), Water Resour. Res., 50, 5410–5427, https://doi.org/10.1002/2014WR015483, 2014.
Lee, S., Wolberg, G., and Shin, S. Y.: Scattered data interpolation with multilevel Bsplines, IEEE T. Vis. Comput. Gr., 3, 228–244, https://doi.org/10.1109/2945.620490, 1997.
Linde, N., Finsterle, S., and Hubbard, S.: Inversion of tracer test data using tomographic constraints, Water Resour. Res., 42, W04410, https://doi.org/10.1029/2004WR003806, 2006.
Lloyd, C. D. and Atkinson, P. M.: Assessing uncertainty in estimates with ordinary and indicator kriging, Comput. Geosci., 27, 929–937, https://doi.org/10.1016/S00983004(00)001321, 2001.
Luo, N., Illman, W. A., Zha, Y., Park, Y.J., and Berg, S. J.: Threedimensional hydraulic tomography analysis of longterm municipal wellfield operations: Validation with synthetic flow and solute transport data, J. Hydrol., 590, 125438, https://doi.org/10.1016/J.JHYDROL.2020.125438, 2020.
Majdalani, S. and Ackerer, P.: Identification of groundwater parameters using an adaptive multiscale method, Groundwater, 49, 548–559, https://doi.org/10.1111/j.17456584.2010.00750.x, 2011.
Marini, M., Felletti, F., Beretta, G. P., and Terrenghi, J.: Three geostatistical methods for hydrofacies simulation ranked using a large borehole lithology dataset from the Venice Hinterland (NE Italy), Water, 10, 844, https://doi.org/10.3390/w10070844, 2018.
Millington, R. J. and Quirk, J. P.: Permeability of porous solids, T. Faraday Soc., 57, 1200–1207, https://doi.org/10.1039/TF9615701200, 1961.
Olsson, A. M. J. and Sandberg, G. E.: Latin hypercube sampling for stochastic finite element analysis, J. Eng. Mech., 128, 121–125, https://doi.org/10.1061/(ASCE)07339399(2002)128:1(121), 2002.
Pechlivanidis, I. G., Jackson, B. M., McIntyre, N. R., and Wheater, H. S.: Catchment scale hydrological modelling: a review of model types, calibration approaches and uncertainty analysis methods in the context of recent developments in technology and applications, Global NEST J., 13, 193–214, 2011.
Phillips, D. L. and Marks, D. G.: Spatial uncertainty analysis: propagation of interpolation errors in spatially distributed models, Ecol. Model., 91, 213–229, https://doi.org/10.1016/03043800(95)001913, 1996.
Poeter, E. P. and Hill, M. C.: Inverse models: a necessary next step in groundwater modelling, Groundwater, 35, 250–260, 1997.
Prautzsch, H., Boehm, W., and Paluszny, M.: Bézier and Bspline techniques, Springer, New York, USA, ISBN 9783642078422, 2002.
Priyanka, B. N. and Mohan Kumar, M. S.: Threedimensional modelling of heterogeneous coastal aquifer: Upscaling from local scale, Water, 11, 421, https://doi.org/10.3390/w11030421, 2019.
Rambourg, D., Ackerer, P., and Bildstein, O.: Groundwater parameter inversion using topographic constraints and a zonal adaptive multiscale procedure: a case study of an alluvial aquifer, Water, 12, 1899, https://doi.org/10.3390/w12071899, 2020.
SanchezLeón, E., Leven, C., Haslauer, C. P., and Cirpka, O. A.: Combining 3D hydraulic tomography with tracer tests for improved transport characterization, Groundwater, 54, 498–507, https://doi.org/10.1111/gwat.12381, 2015.
Straface, S., Chidichimo, F., Rizzo, E., Riva, M., Barrash, W., Revil, A., Cardiff, M., and Guadagnini, A.: Joint inversion of steadystate hydrologic and selfpotential data for 3D hydraulic conductivity distribution at the Boise Hydrogeophysical Research Site, J. Hydrol., 407, 115–128, https://doi.org/10.1016/j.jhydrol.2011.07.013, 2011.
Tarantola, A.: Inverse problem theory and methods for model parameter estimation, Society for Industrial and Applied Mathematics, Philadelphia, USA, ISBN 0898715725, 2005.
Townley, L. R. and Wilson, J. L.: Computationally efficient algorithms for parameter estimation and uncertainty propagation in numerical models of groundwater flow, Water Resour. Res., 21, 1851–1860, https://doi.org/10.1029/WR021i012p01851, 1985.
Viaroli, S., Lotti, F., Mastrorillo, L., Paolucci, V., and Mazza, R.: Simplified twodimensional modelling to constraint the deep groundwater contribution in a complex mineral water mixing area, Riardo Plain, southern Italy, Hydrogeol. J., 27, 1459–1478, https://doi.org/10.1007/s1004001819103, 2019.
Vrugt, J. A., ter Braak, C. J. F., Diks, C. G. H., Robinson, B. A., Hyman, J. M., and Higdon, D.: Accelerating Markov Chain Monte Carlo simulation by differential evolution with selfadaptive randomized subspace sampling, Journal of Nonlinear Sciences and Numerical Simulation, 10, 273–290, https://doi.org/10.1515/IJNSNS.2009.10.3.273, 2009.
Wagener, T. and Kollat, J: Numerical and visual evaluation of hydrological and environmental models using the Monte Carlo analysis toolbox, Environ. Modell. Softw., 22, 1021–1033, https://doi.org/10.1016/j.envsoft.2006.06.017, 2007.
White, D. R.: Propagation of uncertainty and comparison of interpolation schemes, Int. J. Thermophys., 38, 39, https://doi.org/10.1007/s1076501621746, 2017.
Zhang, Z. and Brusseau, M. L.: Characterizing threedimensional hydraulic conductivity distributions using qualitative and quantitative geologic borehole data: Application to a field site, Groundwater, 36, 671–678, https://doi.org/10.1111/j.17456584.1998.tb02842.x, 2005.