<figure><img src="image_1.png" /></figure>

INSTITUT DE FRANCE
Académie des sciences

# Comptes Rendus
## Géoscience
### Sciences de la Planète

Jeremy White and Marsh Lavenue

**Advances in the pilot point inverse method: Où En Sommes-Nous maintenant?**

Volume 355, Special Issue S1 (2023), p. 9-17

Online since: 24 January 2023

**Part of Special Issue:** Geo-hydrological Data & Models

**Guest editors:** Vazken Andréassian (INRAE, France), Valérie Plagnes (Sorbonne Université, France), Craig Simmons (Flinders University, Australia) and Pierre Ribstein (Sorbonne Université, France)

https://doi.org/10.5802/crgeos.161

<figure><img src="image_2.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>

The Comptes Rendus. Géoscience — Sciences de la Planète are a member of the
*Mersenne Center for open scientific publishing*
www.centre-mersenne.org — e-ISSN : 1778-7025
# 1. Introduction

Since the 1970's, the advancements in computer models and associated tools and techniques used to simulate groundwater systems have enhanced hydrogeologists' understanding of the uncertainty in important simulated outcomes, and ultimately, the risks of decisions related to water supply planning, groundwater remediation and other commercial or academic pursuits. These models, tools and techniques collectively make up an important facet of water resource management decision support which aspire to simulate groundwater flow over space and time and the associated uncertainty in the simulation results. Appropriately representing the salient aspects of complex natural groundwater system with a necessarily simple numerical representation and then combining that simplified representation with formal inverse and data assimilation techniques to extract meaningful information from available state observations is one of the drivers behind the advancements in predictive groundwater modeling. Transformation from data to decisions has therefore been the underlying motivation for research in the area of decision support systems. In this regard, Ghislain de Marsily's contribution to the advancement of the inverse methods, a key component of decision support systems, cannot be understated.

Prior to de Marsily's doctoral research published in French in 1978, modelers had resorted to very simplistic representations of the hydrogeologic data obtained from field investigations due to the lack of techniques to appropriately leverage information embedded in the data. These simplistic techniques typically included the use of zones of pre-determined shape and of large spatial extent, which were assigned constant values for hydrogeologic properties such as transmissivity. The zones could not represent the expected (and necessarily stochastic) patterns of spatial heterogeneity often observed in hydrogeologic system properties at many scales. Therefore, the solutions from these early models were plagued with large and unknown errors initiating from the blocky, unrealistic representation of hydrogeologic properties used in the models and propagated to the decisions these models were meant to inform. The earliest attempts to use formal optimization methods to invert the solution of the groundwater flow equation were thwarted when it was discovered that errors

in the hydraulic head measurements, and moreover in the hydraulic gradient, could lead to instabilities in the inverse solution [Emsellem and De Marsily, 1971, Kleinecke, 1971, De Marsily et al., 2000].

de Marsily's doctoral research presented an alternative inverse method which overcame these challenges. He posed the coupling of a new field of research referred to as geostatistics, or spatial statistics, to represent the groundwater system property fields as continuously varying surfaces, thereby overcoming the known limitations of zonation-based schemes. He parameterized the inversion using an indirect inversion technique which enabled him to overcome the numerical instabilities from errors in the measured head data. Lastly, he also adopted the recently published adjoint technique [Chavent, 1974] to efficiently compute sensitivities of the computed head field to changes in the parameter values enabling him to quickly obtain an inverse solution in high dimensions. de Marsily called his technique the “Pilot Point Method” (PPM) because of the way in which he used specific locations to impose changes in parameter field, locations he called “pilot points”. The pilot points would be added to the geostatistical kriging equations and their values assigned optimally to change or “warp” the modeled transmissivity field to reduce the error between the simulated groundwater levels and their measured counterparts. The PPM method fundamentally changed the way groundwater modelers approached the inverse method and is still ubiquitous in groundwater modeling today.

Herein, we provide a brief overview of the theory and history of the PPM, focusing on the important contributions of de Marsily. We finish with some brief remarks related to some of the more recent advances in the PPM and how it continues to serve the decision-support groundwater modeling community.

# 2. Concepts and theory of the PPM

In the PPM, the Bayesian viewpoint of the groundwater model inverse problem is adopted. Originally, the PPM focused on the estimation of the modeled transmissivity field. The log of measured transmissivity is used to transform the statistical distribution of transmissivity (considered to be logarithmic; [Law, 1944, Sánchez-Vila et al., 1996]) to be Gaussian. The unknown log transmissivities in between measurements are viewed as variables of a
Multi-Gaussian random function (RF). Log transmissivity at a point in space is considered a Gaussian random variable (RV) which may or may not have a distribution which is statistically dependent on, or correlated to, its neighboring transmissivities. A transmissivity measurement at location (x, y, z) is conceptually considered as a sample from the RV distribution at point (x, y, z). As presented by Matheron and Blondel [1962], Matheron [1963], kriging provides the best unbiased estimate of the RV at any point in space given observed data and the correlation between the data.

In de Marsily's work, the RV was log transmissivity and the estimation points were the grid block centroids of his finite-difference groundwater flow model. He employed ordinary kriging (OK) to estimate the initial log transmissivity field. OK provides estimates of a RV at any point (or area) in space through the use of measurements of the RV already available and their covariance [see Chilès and Desassis, 2018]. The OK estimator is given by:

$$Z_m^*(u) = \sum_{\beta=1}^{n} v_{\beta,m} Z(u_{\beta}) \quad (1)$$

where $Z_m^*(u)$ is the log of the transmissivity estimate at the centroid of grid block m (located at (u) = (x₁, x₂, x₃)), $Z(u_β)$ is the log of the measured transmissivity data at point $u_β$, and $v_{β,m}$ is the kriging weight assigned to measurement location $β$ given the estimation point (u). The covariance of the RV affects the magnitude of the kriging weights through the solution of the following system of equations used to determine $v_{β,m}$.

$$\left\{ \begin{array}{l} \sum_{\beta=1}^{n} v_{\beta,m} C(u_{\beta} - u_{\alpha}) + \mu(u) = C(u - u_{\alpha}), \quad \alpha = 1, \dots, n \\ \\ \sum_{\beta=1}^{n} v_{\beta,m} = 1 \end{array} \right. \quad (2)$$

where the $v_{β,m}$ are the OK weights and $C(u_β - u_α)$ is the covariance of the RV between measurements points $β$ and $α$. $μ(u)$ is the LaGrange parameter associated with the constraint in the second expression in Equation (2) designed to “filter” out the local mean value to preserve stationarity.

The solution of the inverse problem using the PPM relies upon an expanded form of (1):

$$Z_m^*(u) = \sum_{\beta=1}^{n} v_{\beta,m} Z(u_{\beta}) + \sum_{p=1}^{N} v_{p,m} Z(u_p) \quad (3)$$

where $Z(u_p)$ is the log transmissivity at a pilot point, $v_{p,m}$ is its associated kriging weight and $N$ is the total number of pilot points. The location and number of the pilot points were subjectively chosen by de Marsily. He suggested the number of pilot points to be less than the observed transmissivity values and that the pilot points should be placed in areas with high hydraulic-head gradient. As illustrated in Equation (2), the kriging weights $v_{β,m}$ (and thus $v_{p,m}$) do not depend on actual values of the RF, $Z(u)$. Thus, they may be solved for, prior to assigning transmissivities estimates/values to the pilot points. de Marsily's PPM thus consists of estimating $Z(u_p)$ so as to minimize the weighted least-squares objective function of the measured-to-simulated residual vector, $\phi_Z$:

$$\phi_Z = (\underline{\mathbf{h}}^* - \hat{\mathbf{h}})^T \underline{\mathbf{V}}_h^{-1} (\underline{\mathbf{h}}^* - \hat{\mathbf{h}}) \quad (4)$$

where $\underline{\mathbf{h}}^*$ is the vector of simulated groundwater level, $\hat{\mathbf{h}}$ is the corresponding measured values of groundwater level, and $\underline{\mathbf{V}}_h$ is the covariance matrix of measurement noise.

The sensitivity derivatives required to assign the optimal value of pilot-point transmissivities are computed with the adjoint technique [Chavent, 1974] and the kriging equations (Equation (2)). The overall equation solved by the PPM may be presented as follows.

Let $Z_m$ represent the log transmissivity value assigned to grid block m. Using the chain rule, the sensitivity of the weighted least-squares objective function, $\phi$, to the log transmissivity value assigned to each pilot point may be expressed by:

$$\frac{d\phi}{dZ_p} = \sum_{m=1}^{M} \frac{\partial \phi}{\partial Z_m} \frac{\partial Z_m}{\partial Z_p} \quad (5)$$

where M is the total number of grid blocks in the flow model. The second term on the right-hand side of Equation (5) may be reduced by taking the derivative of Equation (3) with respect to pilot-point transmissivity. Therefore,

$$\frac{\partial Z_m^*(u_m)}{\partial Z_p(u_p)} = v_{p,m} \quad (6)$$

Substituting the right-hand side of Equation (6) into Equation (5) and simplifying yields,

$$\frac{d\phi}{dZ_p} = 2.3 \sum_{m=1}^{M} v_{p,m} k_m \frac{d\phi}{dk_m} \quad (7)$$

where $k$ is the permeability of grid block m. The sensitivity coefficient, $d\phi/dk_m$, can be obtained by
adjoint sensitivity analysis or by perturbation methods and used to determine the optimal values of the pilot-point's log transmissivity by modifying the initial guess of the pilot point transmissivities (taken as the kriged estimate at $v_p$ obtained from the measured data). de Marsily states that constraints to the pilot-point transmissivities could be applied to ensure that the pilot-point transmissivities lie within ±2 standard deviations of the kriged estimate.

As mentioned earlier, de Marsily's PPM was the first to solve many of the problems inherent in the earlier inverse techniques. By incorporating geostatistics and logarithms of transmissivity, he ensured a smooth but spatially varying transmissivity field even in the presence of significant errors in the head field. In addition, he stabilized the inverse solution by:

* adding prior information, e.g. regularization, in the form of initial kriged estimates at the pilot point locations,

* properly posing the inverse problem by reducing the number of unknowns to a small number of pilot points,

* implementing the transmissivity changes in the neighborhood of a pilot point in accord with the variogram, and

* utilizing efficient adjoint sensitivity techniques to obtain the necessary derivatives for the optimization routine.

### 3. Early PPM advances

The first mention of the PPM was in de Marsily's PhD dissertation [De Marsily, 1978], which describes using “pilot points” at discrete locations. One of the first applications of the PPM published in English was presented in De Marsily et al. [1984]. Here the PPM was used to obtain the permeability field which reproduced interference tests in a well field.

In the late 1980s and early 1990s, research of geostatistical methods by Ahmed and De Marsily [1987] and research of the inverse problem by Certes and De Marsily [1991] resulted in improvements to the PPM. Lavenue and Pickens [1992] coupled ordinary kriging and the PPM to history match a regional groundwater model at the Waste Isolation Pilot Plan (WIPP) site in Carlsbad, New Mexico. They also adopted an adjoint sensitivity method to optimally locate pilot points, which were sequentially

added to the model during an iterative inversion process. One of the significant features added to the model transmissivities during inversion was a high-transmissivity fracture zone which provided a fast offsite pathway to the WIPP-site boundary potentially threatening the WIPP site safety assessment.

A technical review panel comprising the leading groundwater inversion researchers of the 1990s was convened by Sandia National Laboratories to review the WIPP model construction, inversion and travel time calculations. This panel, called the Geostatistics Expert Group (GXG), was headed by Ghislain de Marsily. The GXG recommended that the PPM approach be expanded to calibrate an ensemble of conditionally simulated transmissivity fields in order to investigate groundwater travel time uncertainty to the WIPP-site boundary. de Marsily, Lavenue, and RamaRao did so by linking a conditional simulation front-end routine (using turning bands), a flow model (SWIFT II), an adjoint sensitivity routine to optimally located pilot points and a conjugate gradient optimization routine to assign the pilot-point transmissivity values. In 1995, the resulting new PPM method referred to as GRASP-INV and its application to the WIPP site produced an ensemble of 100 conditionally simulated transmissivity fields [Rama-rao et al., 1995, Lavenue et al., 1995].

This extension of the PPM method to produce an ensemble of solutions, all equally calibrated or “history matched” to the observed data, was a significant step forward in the quantification of posterior uncertainty, which ultimately allowed for a more complete risk assessment to support decision making. It also laid a foundation for the use of inversion and sensitivity techniques to assess data worth which guided future data acquisition campaigns at the WIPP site to further reduce uncertainties in our understanding of the hydrogeologic system.

Shortly afterward, Gomez-Hernandez and his fellow researchers at the University of Valencia, Spain, developed a non-linear pilot-point technique (or a “pilot-point-like” method as one reviewer pointed out) that differs from the original approach of De Marsily [1978] and Lavenue et al. [1995]. These works include Sahuquillo et al. [1992], Gómez-Hernánez et al. [1997], and Capilla et al. [1997]. The solution of the inverse problem using Gomez-Hernandez’ approach is significantly faster than that of Lavenue and Pickens [1992] or Certes and
De Marsily [1991] for several possible reasons. First, Gomez-Hernandez did not update the sensitivity derivatives as frequently as the other pilot-point approaches. The second major difference in Gomez-Hernandez's method is that he assumed the conductance between two finite-difference grid blocks is expressed by the geometric mean (as opposed to the harmonic mean) of the associated grid-block transmissivities, which significantly reduced the time required to compute sensitivities. The combined effect of these differences significantly sped up the inversion process using Gomez-Hernandez's technique which he called the sequential self-calibration (SSC) pilot point method.

Sandia National Laboratory decided to benchmark GRASP-INV against other inversion techniques represented by the GXG members and invited Gomez-Hernandez to participate. This unique bench-marking exercise, described in Zimmerman et al. [1998], evaluated the performance of virtually every inversion technique used by industry and academics at that time. The ultimate objective of the study was to determine which of several geostatistical inverse techniques is better suited for making probabilistic forecasts of the potential transport of solutes in an aquifer where spatial variability and uncertainty in hydrogeologic properties are significant. Seven geostatistical methods were compared on four synthetic data sets. The inverse methods tested were categorized as being either linearized or nonlinear [Carrera and Glorioso, 1991]. The linearized approaches are generally based upon simplifying assumptions about the flow field (e.g., a uniform hydraulic head gradient, a small $log(T)$ variance, etc.), that lead to a linearized relation between $T$ and head using a perturbation expansion of the head and transmissivity fields. This equation can then be solved analytically or numerically. The nonlinear approaches have no such restrictions placed on them and can, in principle, handle more complex flow fields or larger $log(T)$ variances. Linearized methods included the Fast Fourier Transform inverse method by Gutjahr et al. [1994], the Linearized Cokriging Method by Kitanidis and Lane [1985], and the Linearized Semianalytical method by Rubin et al. [1992], while the non-linear methods included the Fractal Simulation inverse method by Grindrod and Impey [1991], the Maximum Likelihood inverse method by Carrera [1994], the Pilot Point Method [Rama-

rao et al., 1995] and the SSC pilot point method by Gomez-Hernandez [Gómez-Hernánez et al., 1997]. While the non-linear inverse methods were superior to the linearized inverse methods, the most robust results came from the inverse methods that were able to incorporate rigorous geostatistics into the inversion process, that is, inversion suffered when a technique could not adequately represent nor modify the expected spatial (e.g. Prior) patterns of variability of the aquifer properties.

Lessons learned in Zimmerman et al. [1998] regarding the importance of robust geostatistics led de Marsily, Lavenue and RamaRao to enhance their PPM by adopting the GSLIB [Deutsch and Journel, 1992] routines for conditional simulation. They collaborated with Gomez-Hernandez in this effort which resulted in a new conditional simulator for the PPM which could produce an ensemble of multi-facies geostatistical fields with unique covariance structures using categorical indicator simulation followed by sequential gaussian simulation to overlay spatial variability within each facies. de Marsily, Lavenue and RamaRao linked this new capability to their PPM method and later applied it to the WIPP site regional model [Lavenue, 1998]. Unique to this enhancement was the ability of the PPM to optimize the spatial variability within each facies (fractured and unfractured zones) separately to match transient water levels. It also enabled the delineation of a sharp fractured/unfractured zone interface within the regional aquifer which reduced travel time uncertainty. Lavenue and de Marsily then extended this version of the PPM to three dimensions and demonstrated its utility by matching a series of three-dimensional interference tests [Lavenue and De Marsily, 2001].

Other research motivated by Zimmerman et al. [1998] included methods to enforce prior plausibility in a Bayesian sense, the work of Alcolea et al. [2006], Riva et al. [2010] and later Jiménez et al. [2016] worked to include expert knowledge into the PPM solution process so that the resulting property fields were adjusted in harmony with the prior estimates.

## 4. The tipping point

Notwithstanding the developments and advancements of the PPM cited above, its use was not widespread since the PPM was largely coupled to
groundwater model software not widely used by industry. This changed in the late 1990s, when John Doherty, adopted the PPM for his PEST software suite and produced a full suite of advanced software tools that allowed the PPM to be applied at scale in an unlimited number of settings. Unlike the existing adjoint-based PPM approaches, the PEST software suite relies on non-intrusive/model-independent approaches, such as finite-difference derivatives to estimate the first-order relation between parameters and simulation results, so that it can be applied to practically any forward numerical groundwater model. Doherty realized that like PEST, the PPM can be implemented in a non-intrusive/model-independent approach, so that the PPM could also be applied to practically any forward numerical groundwater model. This shifted the use of the PPM from the bespoke research and development projects cited above to the broader groundwater modeling community, a tipping point of scaling the PPM to industry at large.

Within the PEST framework, the PPM serves as an indispensable parameterization scheme, one that balances the need to express uncertain spatial heterogeneity in model inputs with the requirement to maintain a computationally tractable inverse problem and represent broad-scale heterogeneity. We generalize from transmissivity to “model inputs” here because with the non-intrusive PEST framework, pilot points are commonly used to represent realistic spatial heterogeneity in a wide-range of uncertain properties and boundary condition elements that numerical models require. For example, Knowling and Werner [2016] and later Knowling and Werner [2017] used pilot points to represent uncertain spatial patterns of recharge, while McKenna et al. [2020] and later White et al. [2020a] used pilot points in a multi-scale parameterization scheme to explicitly represent different scales of heterogeneity and associated uncertainty.

The issue of inverse problem stability and posedness with the PPM were both overcome through the adoption of formal Tikhonov regularization [Tikhonov, 1963, Tikhonov and Arsenin, 1977] [see Tonkin and Doherty, 2005, Fienen et al., 2009, Doherty et al., 2010, Doherty, 2015], in combination with truncated subspace solution techniques [see Oliver et al., 2008, Aster et al., 2018, Doherty, 2015]. These two advancements, which were implemented

at scale to function in concert with each other for the first time in PEST software suite, allowed practitioners to forego the need to *a priori* reduce the number of pilot points to stabilize the inverse problem—Tikhonov regularization and truncated subspace techniques optimally and unconditionally stabilize the inverse problem. Now practitioners were free to express expected property heterogeneity (and associated uncertainty) at scales that are relevant to the predictive outcomes of the model.

## 5. Où En Sommes-Nous maintenant?

With the recent advancements in the PPM mentioned above, the combined use of the PPM and PEST has led to wide-spread investigation of predictive uncertainties in applied groundwater modeling, which is facilitated by the PPM allowance for a more complete and robust expression of expected patterns of spatial heterogeneity in subsurface properties [see for example Moore and Doherty, 2005, James et al., 2009, Dausman et al., 2010, Herckenrath et al., 2011, Tonkin and Doherty, 2009, Christensen and Doherty, 2008, Keller et al., 2021].

The PPM, as embodied within Doherty's PEST software suite, continues to serve an important role in decision-support environmental simulation around the world and new developments continue. Several recent advancements have focused on expanding the use of pilot points from 2-point geostatistics to multipoint geostatistics [see Li et al., 2013, Mariethoz and Caers, 2014, Ma and Jafarpour, 2018, Khambhammettu et al., 2020, Liu et al., 2021]. The capability of these pilot point methods to facilitate the use of multi-point geostatistical property representations within formal data assimilation is a major advancement in settings where high-order multi-point geostatistical property representations are important, such as settings where connected permeability are important for representing the fate and transport of dissolved-phase constituents.

So “Where are We Now?” Groundwater modelers have more computational power than ever before to build highly complex models and invert them with PEST, and more recently, the PEST++ software suite [White et al., 2020b]. With the related advancements in history matching using the iterative ensemble smoother approaches [e.g. Chen and Oliver, 2013, White, 2018], practitioners are able to more
efficiently employ high-dimensional inversion and
uncertainty quantification in non-intrusive ways and
at scales not previously possible except for those
with access to high-performance computing. This
"gilded age" of inversion however comes with the
responsibility of parsimony. It cannot be overstated
that sound groundwater model analyses are built on
the foundation of a solid hydrogeologic conceptual
model that is focused on the predictive purpose of
the modeling, and a firm grasp of the sources and
magnitudes of model input (i.e. parameter) and con-
ceptual model uncertainty. Decisions based on the
associated numerical model can only then take full
advantage of the decision-support simulation tools
discussed herein.

A final comment: When asked about his upcom-
ing conference presentation in the early 1980s, de
Marsily mentioned that his talk was on applying geo-
statistics in hydrogeology, which at the time was vir-
tually unknown in the US. He went on to say that
he felt it was important to explain to the audience
how powerful geostatistical tools are to hydrogeol-
ogists and that it "was never too late to learn new
things". He was right on both points. Countless hy-
drogeologists have been positively impacted by the
"new things" Ghislain de Marsily taught us and by
his example of leading with humility, style and grace.
Thank you, Ghislain.

Conflicts of interest

Authors have no conflict of interest to declare.

References

Ahmed, S. and De Marsily, G. (1987). Comparison of geostatistical methods for estimating transmissivity using data on transmissivity and specific capacity. *Water Resour. Res.*, 23(9), 1717–1737.

Alcolea, A., Carrera, J., and Medina, A. (2006). Piv- lot points method incorporating prior information for solving the groundwater flow inverse problem. *Adv. water Resour.*, 29(11), 1678–1689.

Aster, R. C., Borchers, B., and Thurber, C. H. (2018). *Parameter Estimation and Inverse Problems*. Elsevier, Amsterdam.

Capilla, J., Gómez-Hernández, J. J., and Sahuquillo, A. (1997). Stochastic simulation of transmissivity fields conditional to both transmissivity and

piezometric data 2. Demonstration on a synthetic
aquifer. *J. Hydrol.*, 203(1–4), 175–188.

Carrera, J. (1994). Invert-4. a fortran program for solving the groundwater flow inverse problem. User's guide. Technical report, CIMNE, Barcelona, Spain. 160 pp. þ appendices.

Carrera, J. and Glorioso, L. (1991). On geostatistical formulations of the groundwater flow inverse problem. *Adv. Water Resour.*, 14(5), 273–283.

Certes, C. and De Marsily, G. (1991). Application of the pilot point method to the identification of aquifer transmissivities. *Adv. Water Resour.*, 14(5), 284–300.

Chavent, G. (1974). Estimation de paramètres dis-
tribues dans des equations aux derivees partielles.
In *Computing Methods in Applied Sciences and En-
gineering Part 2*, pages 361–390. Springer, Berlin.

Chen, Y. and Oliver, D. S. (2013). Levenberg–marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. *Comput. Geosci.*, 17(4), 689–703.

Chilès, J.-P. and Desassis, N. (2018). Fifty years of kriging. In *Handbook of Mathematical Geosciences*, pages 589–612. Springer, Cham.

Christensen, S. and Doherty, J. (2008). Predictive error dependencies when using pilot points and singular value decomposition in groundwater model calibration. *Adv. Water Resour.*, 31(4), 674–700.

Dausman, A. M., Doherty, J., Langevin, C. D., and Sukop, M. C. (2010). Quantifying data worth toward reducing predictive uncertainty. *Groundwater*, 48(5), 729–740.

De Marsily, G. (1978). *De l'identification des systèmes hydrogéologiques*. PhD thesis, Pierre and Marie Curie Univ., Paris.

De Marsily, G., Delhomme, J., Coudrain-Ribstein, A., and Lavenue, A. (2000). Four decades of inverse problems in hydrogeology. In Zhang, D. and Winter, C. L., editors, *Theory, Modeling and Field Investigation in Hydrogeology: A Special Volume in Honor of Shlomo P. Neuman's 60th Birthday*, Special Paper 348, pages 1–17. Geological Society of America, Boulder, Colorado.

De Marsily, G., Lavedan, G., Boucher, M., and Fasanino, G. (1984). Interpretation of interference tests in a well field using geostatistical techniques to fit the permeability distribution in a reservoir model. In Verly, G., David, M., Journel, A. G.,
and Marecha, A., editors, *Geostatistics for Natural Resources Characterization 2nd NATO Advanced Study Institute, South Lake Tahoe, CA, September 6–17, 1987*, pages 831–849. D. Reidel, Hingham, MA. Pt. 2.

Deutsch, C. and Journel, A. (1992). *GSLIB: Geostatistical Software Library and User's Guide*. Oxford University Press, New York, NY.

Doherty, J. (2003). Ground water model calibration using pilot points and regularization. *Groundwater*, 41(2), 170–177.

Doherty, J. E. (2015). *PEST and its Utility Support Software, Theory*. Watermark Numerical Publishing, Brisbane.

Doherty, J. E., Fienen, M. N., and Tonkin, M. J. (2010). Approaches to highly parameterized inversion: Pilot-point theory, guidelines, and research directions. U.S. Geological Survey Scientific Investigations Report 2010-5168, https://doi.org/10.3133/sir20105168.

Emsellem, Y. and De Marsily, G. (1971). An automatic solution for the inverse problem. *Water Resour. Res.*, 7(5), 1264–1283.

Fienen, M. N., Muffels, C. T., and Hunt, R. J. (2009). On constraining pilot point calibration with regularization in pest. *Groundwater*, 47(6), 835–844.

Gómez-Hernández, J. J., Sahuquillo, A., and Capilla, J. (1997). Stochastic simulation of transmissivity fields conditional to both transmissivity and piezometric data—i. Theory. *J. Hydrol.*, 203(1–4), 162–174.

Grindrod, P. and Impey, M. (1991). Fractal field simulations of tracer migration within the wipp culebra dolomite. *Rep. IM2856-1, Intera Information Technologies*, 2, 62.

Gutjahr, A., Bullard, B., Hatch, S., and Hughson, L. (1994). Joint conditional simulations and the spectral approach for flow modeling. *Stoch. Hydrol. Hydroaul.*, 8(1), 79–108.

Herckenrath, D., Langevin, C. D., and Doherty, J. (2011). Predictive uncertainty analysis of a saltwater intrusion model using null-space Monte Carlo. *Water Resour. Res.*, 47(5), 1–16.

James, S. C., Doherty, J. E., and Eddebbarh, A.-A. (2009). Practical postcalibration uncertainty analysis: Yucca mountain, Nevada. *Groundwater*, 47(6), 851–869.

Jiménez, S., Mariethoz, G., Brauchler, R., and Bayer, P. (2016). Smart pilot points using reversible-jump

markov-chain monte carlo. *Water Resour. Res.*, 52(5), 3966–3983.

Keller, J., Franssen, H.-J. H., and Nowak, W. (2021). Investigating the pilot point ensemble kalman filter for geostatistical inversion and data assimilation. *Adv. Water Resour.*, 155, article no. 104010.

Khambhammettu, P., Renard, P., and Doherty, J. (2020). The traveling pilot point method. a novel approach to parameterize the inverse problem for categorical fields. *Adv. Water Resour.*, 138, article no. 103556.

Kitanidis, P. K. and Lane, R. W. (1985). Maximum likelihood parameter estimation of hydrologic spatial processes by the gauss-newton method. *J. Hydrol.*, 79(1–2), 53–71.

Kleinecke, D. (1971). Use of linear programming for estimating geohydrologic parameters of groundwater basins. *Water Resour. Res.*, 7(2), 367–374.

Knowling, M. J. and Werner, A. D. (2016). Estimability of recharge through groundwater model calibration: Insights from a field-scale steady-state example. *J. Hydrol.*, 540, 973–987.

Knowling, M. J. and Werner, A. D. (2017). Transient recharge estimability through field-scale groundwater model calibration. *Groundwater*, 55(6), 827–840.

Lavenue, A. M. (1998). *A new pilot point inverse method in hydrogeology: Generating an ensemble of conditionally-simulated transmissivity fields*. PhD thesis, School of Mines, Paris.

Lavenue, A. M. and Pickens, J. F. (1992). Application of a coupled adjoint sensitivity and kriging approach to calibrate a groundwater flow model. *Water Resour. Res.*, 28(6), 1543–1569.

Lavenue, A. M., Ramarao, B. S., De Marsily, G., and Marietta, M. G. (1995). Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields: 2. Application. *Water Resour. Res.*, 31(3), 495–516.

Lavenue, M. and De Marsily, G. (2001). Three-dimensional interference test interpretation in a fractured aquifer using the pilot point inverse method. *Water Resour. Res.*, 37(11), 2659–2675.

Law, J. (1944). A statistical approach to the interstitial heterogeneity of sand reservoirs. *Trans. AIME*, 155(01), 202–222.

Li, L., Srinivasan, S., Zhou, H., and Gómez-Hernández, J. J. (2013). A pilot point guided pattern matching approach to integrate dynamic data
into geological modeling. *Adv. Water Resour.*, 62, 125-138.

Liu, J.-b., Jiang, S.-m., Zhou, N.-q., Cai, Y., Cheng, L., and Wang, Z.-y. (2021). Groundwater contaminant source identification based on QS-ilues. *J. Groundw. Sci. Eng.*, 9(1), 73-82.

Ma, W. and Jafarpour, B. (2018). Pilot points method for conditioning multiple-point statistical facies simulation on flow data. *Adv. Water Resour.*, 115, 219-233.

Mariethoz, G. and Caers, J. (2014). *Multiple-Point Geostatistics: Stochastic Modeling with Training Images*. John Wiley & Sons, Hoboken, NJ.

Matheron, G. (1963). *Traité de géostatistique appliquée tome ii: Le krigeage, Memoires, Bureau du Recherches Geologiques et Minières*, volume 24, issue 14. Éditions Technip, Paris.

Matheron, G. and Blondel, F. (1962). *Traité de géostatistique appliquée, tome i*. Editions Technip, Paris. Á. Martínez-del-Pozo et al.

McKenna, S. A., Akhriev, A., Ciaurri, D. E., and Zhuk, S. (2020). Efficient uncertainty quantification of reservoir properties for parameter estimation and production forecasting. *Math. Geosci.*, 52, 233-251.

Moore, C. and Doherty, J. (2005). Role of the calibration process in reducing model predictive error. *Water Resour. Res.*, 41(5), 1-14.

Oliver, D. S., Reynolds, A. C., and Liu, N. (2008). *Inverse Theory for Petroleum Reservoir Characterization and History Matching*. Cambridge University Press, Cambridge.

Ramarao, B. S., Lavenue, A. M., De Marsily, G., and Marietta, M. G. (1995). Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields: 1. Theory and computational experiments. *Water Resour. Res.*, 31(3), 475-493.

Riva, M., Guadagnini, A., De Gaspari, F., and Alcolea, A. (2010). Exact sensitivity matrix and influence of the number of pilot points in the geostatistical inversion of moment equations of groundwater flow. *Water Resour. Res.*, 46, 1-15.

Rubin, Y., Mavko, G., and Harris, J. (1992). Mapping permeability in heterogeneous aquifers using hy-

drologic and seismic data. *Water Resour. Res.*, 28(7), 1809-1816.

Sahuquillo, A., Capilla, J., Gómez-Hernández, J., and Andreu, J. (1992). Conditional simulation of transmissivity fields honoring piezometric data. *Hydraul. Eng. Softw. IV, Fluid Flow Model.*, 2, 201-214.

Sánchez-Vila, X., Carrera, J., and Girardi, J. P. (1996). Scale effects in transmissivity. *J. Hydrol.*, 183(1-2), 1-22.

Tikhonov, A. N. (1963). Regularization of incorrectly posed problems. *Sov. Math. Dokl.*, 4, 1624-1627.

Tikhonov, A. N. and Arsenin, V. Y. (1977). *Solutions of Ill-Posed Problems*. V.H. Winston & Sons, Washington, DC: John Wiley & Sons, New York, NY.

Tonkin, M. and Doherty, J. (2009). Calibration-constrained monte carlo analysis of highly parameterized models using subspace techniques. *Water Resour. Res.*, 45(12), 1-17.

Tonkin, M. J. and Doherty, J. (2005). A hybrid regularized inversion methodology for highly parameterized environmental models. *Water Resour. Res.*, 41(10), 1-16.

White, J. T. (2018). A model-independent iterative ensemble smoother for efficient history-matching and uncertainty quantification in very high dimensions. *Environ. Model. Softw.*, 109, 191-201.

White, J. T., Foster, L. K., Fienen, M. N., Knowling, M. J., Hemmings, B., and Winterle, J. R. (2020a). Toward reproducible environmental modeling for decision support: A worked example. *Front. Earth Sci.*, 8, article no. 50.

White, J. T., Hunt, R. J., Fienen, M. N., and Doherty, J. E. (2020b). Approaches to Highly Parameterized Inversion: PEST++ Version 5, a Software Suite for Parameter Estimation, Uncertainty Analysis, Management Optimization and Sensitivity Analysis. U.S. Geological Survey Techniques and Methods 7C26, 52 p.

Zimmerman, D., De Marsily, G., Gotway, C. A., Marietta, M. G., Axness, C. L., Beauheim, R. L., Bras, R. L., Carrera, J., Dagan, G., Davies, P. B., et al. (1998). A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. *Water Resour. Res.*, 34(6), 1373-1413.