Research papersA hybrid inverse method for hydraulic tomography in fractured and karstic media
Introduction
The hydrodynamic prediction of fluid flow and solute transport in subsurface aquifers requires the spatial distributions of hydraulic properties, e.g. transmissivity (T) and storativity (S), to be accurately defined. Recently, hydraulic tomography (HT) has been proven to be an efficient technique for characterizing site-specific, spatial distributions of transmissivity (e.g. Gottlieb and Dietrich, 1995, Butler et al., 1999, Yin and Illman, 2009, Castagna et al., 2011, Klepikova et al., 2013, Illman, 2014). This method jointly analyses several sets of cross-hole pumping test data, and therefore integrates more hydraulic information contained in the recorded data (Illman, 2014). The efficiency of HT has been demonstrated by numerical (Hao et al., 2008), laboratory (Illman et al., 2007, Liu et al., 2007, Sharmeen et al., 2012) and field experiments (Straface et al., 2007, Illman et al., 2009, Wang et al., 2016).
Numerous inverse methods have been developed for conditioning hydrogeology models to observed hydraulic data and assessing the associated uncertainties (see e.g. Carrera et al., 2005 and Illman, 2014). Among others, the most popular ones are the geostatistical inverse methods and the stochastic sampling methods. The geostatistical approaches (e.g. Kitanidis and Vomvoris, 1983, Kitanidis, 1995, Yeh et al., 1995, Yeh et al., 1996, Kitanidis and Lee, 2014), which formulate the inverse problem as a regularized optimization, are efficient for high-dimensional problems. They have been widely used in many hydrogeological studies (e.g. Zhang and Yeh, 1997, Hanna and Yeh, 1998, Zhu and Yeh, 2005, Soueid Ahmed et al., 2014), and has recently been successfully applied to fractured media (Hao et al., 2008, Illman et al., 2009, Castagna et al., 2011, Sharmeen et al., 2012, Zha et al., 2015, Zha et al., 2016). Even though these inverse practices turn out to be successful, the solutions given by the geostatistical approaches, in many cases, correspond to local minima, and often suffer from inadequate exploration of the parameter space (Bates and Campbell, 2001, Keating et al., 2010).
In recent years, stochastic approaches based on sample-based Bayesian inference have become popular in aquifer characterization (Fu and Gomez-Hernandez, 2009, Vrugt et al., 2008, Mondal et al., 2010, Cui et al., 2011). These methods, using Markov chain Monte Carlo (MCMC) simulations to generate samples from the posterior target distribution of the model parameters, allowing for inference of both the unknown parameters and the uncertainties associated with the analysis. Comparing to geostatistical approach, MCMC simulations allow to avoid local minima and generate global statistical solutions that do not depend on a priori models.
In a traditional MCMC simulation that relies on a single-site updating scheme, a large number of forward model runs (equal to the number of parameters) are required to perform a single parameter sweep (Metropolis et al., 1953). The required number of simulations becomes prohibitive for high-dimensional inverse models that include a large number of model parameter. A popular approach to make the sampling adaptive for high-dimensional models is to construct surrogate models or emulators based on polynomial approximations (e.g. polynomial chaos or Karhunen–Loève expansions) that reduce the dimension of the problem (e.g. Marzouk and Xiu, 2009, Laloy et al., 2013). The approach is often improved with a two-stage scheme to achieve better sampling efficiency (Liu, 2001, Dostert et al., 2006, Cui et al., 2011). However, because of the requirement of smoothness in the model response fields, the emulator-based methods cannot be applied to highly heterogeneous media, such as fractured and karstic aquifers where sharp changes in the spatial distribution of hydraulic head/flux are commonly observed. Another approach for improving the efficiency of MCMC sampler is to adopt a multivariate proposal density (e.g. Haario et al., 2006, Fu and Gomez-Hernandez, 2009). For instance, the delayed rejection adaptive metropolis (DRAM) MCMC resorts to constructing an approximation to the posterior covariance matrix to guide the sampling process. Recently, many authors have proposed hybrid MCMC algorithms, in which local derivative information (gradient and Hessian for the forward model) is used in the construction of the proposal density (Girolami and Calderhead, 2011, Flath et al., 2011, Martin et al., 2012). Because the proposal distributions are tuned using derivative information contained in the past chain, the proposal distribution evolves to the target density much more efficiently. Although using Hessian information in conditioning MCMC samplers has been previously considered, most of the applications are restricted to models with a relatively small number of parameters, or to high-dimensional linear problems where the derivative information (gradient and Hessian) is easy to obtain (Flath et al., 2011). The applicability of the method to complex, nonlinear, hydrogeological models remains to be demonstrated. Also, there are only a few MCMC applications in the context of hydraulic tomography (Oliver et al., 1997, Fu and Gomez-Hernandez, 2009, Jardani et al., 2012). Therefore, in the present study, we have attempted to apply a Hessian-based MCMC method, originally proposed in Martin et al. (2012), to infer the posterior distribution of the unknown and uncertain parameters in a high-dimensional groundwater model.
No matter which type of inverse method is applied, inverse problems have an essential difficulty: non-uniqueness (Tarantola, 2005). Consequently, there are many models can fit the data equally well. To overcome this difficulty, the inversion models have to be constrained by a priori information. Conventionally, variogram models are adopted to pose assumptions about the spatial pattern of the model parameter field, in order to constrain the solution. However, the classical variogram models based on kriging or cokriging processes generate smooth spatial patterns. Therefore, they are not suitable for modelling fractured and karstified media. To capture the prominent heterogeneous nature and sharp property variation of fractured aquifers, indicator geostatistical methods based on categorical classification are developed (Goovaerts, 1996, Day-Lewis et al., 2000). However, the traditional indicator approaches require a large amount of transmissivity measurements. In the case where only sparse transmissivity data are available, the emergent categorical realizations are often geologically unrealistic (Illman, 2014). An alternative and promising method to model highly heterogeneous media is the transitional-probability (TP) based categorical approach (Ritzi et al., 1995, Carle and Fogg, 1996). The approach employs measured transition probabilities between different rock types as the weights of their spatial cross-correlation. The method has been extensively used in soil surveys and successfully applied to alluvial deposits (Carle and Fogg, 1997, Weissmann et al., 1999, Weissmann and Fogg, 1999, Lee et al., 2007), glacial deposits (He et al., 2014), and recently to fractured media (Park et al., 2004, Blessent et al., 2011). This method has also been used in the present work to generate a priori transmissivity structure fields, where the salient heterogeneity and anisotropy of fracture and karst networks are captured.
Our analysis differs from previous studies (e.g. Park et al., 2004 and Blessent et al., 2011) in how the spatial patterns are generated from TP geostatistical approach, and are subsequently calibrated to match hydraulic data. In Park et al. (2004) and Blessent et al. (2011) a classical two-stage approach (e.g. McKenna and Poeter, 1995, Day-Lewis et al., 2000) is used. In the first stage, the TP spatial patterns generated by fitting geostatistical models (transiograms) and hard data are treated as stochastic continuum models (i.e. rock type bodies being transformed into stochastic transmissivity zones for a flow model). In the second stage, for each realization, a set of mean hydraulic properties of transmissivity zones is then calibrated to match hydraulic data. In other words, variation of hydraulic properties within each rock type is not allowed. In such cases, history matching of hydraulic data can only be achieved to a limited level where the general trend displayed by the data is preserved, while disregarding a large amount of valuable information that may be crucial to infer fine-scale heterogeneity.
In our approach, we still follow the two-stage strategy. But, in the first stage, instead of fitting available sparse transmissivity estimates to transiogram models, we specify the explicitly the embedded transition probability matrix, based on integrated hydrogeological interpretation, to derive Markovian categorical rock-type realizations (Section 3.1). The stochastic realizations are conditioned to transmissivity estimates and inter-borehole connectivity both are determined from field cross-hole pumping tests. In the second stage, the generated TP realizations is used as a priori structural model of transmissivity, and we then apply a hybrid MCMC method (Section 3.2) to statistically infer the spatially varying transmissivity (mean and variance) of grid cells in selected parameter zones. In this sense, our inversion approach is based on an equivalent porous medium representation.
Section snippets
Site description and main hydraulic dataset
The location of this study is the Terrieu experimental site, which is located in the Montpellier region, Southern France (Fig. 1). The field site has a surface area of 1500 m2 (30 m × 50 m), and contains 22 boreholes that were drilled through the Cretaceous marly limestones of the Lez aquifer (total thickness of about 30–40 m; Fig. 2). The geological and hydrogeological backgrounds of the site have been described in detail in previous studies (Jazayeri Noushabadi, 2009, Wang et al., 2014, Wang et al.,
Categorical indicator classification
Because the transition-probability geostatistical approach is essentially an indicator-based method, it is necessary to define the categorical indicators. A good criterion for distinguishing one rock type from another is hydraulic transmissivity. The three indicators defined for further stochastic, categorical simulations are:
1) ‘bedding plane’, which has a transmissivity on the level of , and corresponds to the flow areas characterized by some intact aperture of the bedding plane;
2)
Synthetic case study
Prior to real-world applications, the proposed stochastic Newton inverse approach is validated using a synthetic case study. A synthetic rock type field was firstly generated on a 60 × 36 grid using the geostatistical method described in Section 3.2. Three rock types, i.e. ‘karst’, ‘fracture’ and ‘bedding plane’, were considered in the model (Fig. 7a). The rock type field was then translated into a log-transmissivity field by assigning a uniform initial hydraulic transmissivity (Section 3.1.1) to
Discussions
The proposal density of the proposed stochastic Newton (SN) method is composed of a deterministic part and a stochastic part. The stochastic Newton approach outperforms the deterministic inverse methods, because it allows not only the mean but also the uncertain range of the transmissivity field to be quantified. On the other hand, the proposal density of the stochastic Newton method is dynamically adjusted during the evaluation of the MCMC chains. Hence, despite that the SN method is
Conclusions
To conclude, in this study, a new stochastic Newton approach has been implemented for imaging the spatial transmissivity field of a real fractured and karstic aquifer. By integrating available geological and hydrogeological information collected from outcrops and boreholes, a 2D conceptual flow model was firstly conceptualized. On the basis of a regular parameterization, a transition-probability-based geostatistical method was employed to generate a series of equally plausible geological models
Acknowledgements
This work was performed within the framework of the MEDYCYSS observation site, part of the KARST observatory network (www.sokarst.org) initiated by the INSU/CNRS, which aims to strengthen knowledge-sharing and to promote cross-disciplinary research on karst systems. We thank S. F. Carle for providing the TPROGS code, and E. Pujades and two anonymous reviewers for their useful comments and suggestions.
References (81)
- et al.
Efficient solution of nonlinear, underdetermined inverse problems with a generalized PDE model
Comput. Geosci.
(2008) - et al.
Coarse-gradient Langevin algorithms for dynamic data integration and uncertainty quantification
J. Comput. Phys.
(2006) - et al.
Uncertainty assessment and data worth in groundwater flow and mass transport modeling using a blocking Markov chain Monte Carlo method
J. Hydrol.
(2009) - et al.
Steady-state hydraulic tomography in a laboratory aquifer with deterministic heterogeneity: multi-method and multiscale validation of hydraulic tomography tomograms
Hydrol. J.
(2007) - et al.
Geostatistical inverse modeling of the transmissivity field of a heterogeneous alluvial aquifer under tidal influence
J. Hydrol.
(2012) - et al.
Influence of the observation scale on permeability estimation at local and regional scales through well tests in a fractured and karstic aquifer (Lez aquifer, Southern France)
J. Hydrol.
(2011) - et al.
Numerical simulation of Forchheimer flow to a partially penetrating well with a mixed-type boundary condition
J. Hydrol.
(2015) - et al.
Bayesian uncertainty quantification for flows in heterogeneous porous media using reversible jump Markov chain Monte Carlo methods
Adv. Water Resour.
(2010) - et al.
Characterisation of the transmissivity field of a fractured and karstic aquifer, Southern France
Adv. Water Resour.
(2016) - et al.
Multi-scale alluvial fan heterogeneity modeled with transition probability geostatistics in a sequence stratigraphic framework
J. Hydrol.
(1999)
What does hydraulic tomography tell us about fractured geological media? A field study and synthetic experiments
J. Hydrol.
A Markov chain Monte Carlo scheme for parameter estimation and inference in conceptual rainfall runoff modeling
Water Resour. Res.
Three-dimensional transient hydraulic tomography in a highly heterogeneous glaciofluvial aquifer-aquitard system
Water Resour. Res.
Inverse modeling of hydraulic tests in fractured crystalline rock based on a transition probability geostatistical approach
Water Resour. Res.
A scaled stochastic Newton algorithm for Markov chain Monte Carlo simulations
Tech. Rep.
Pumping tests in networks of multilevel sampling wells: motivation and methodology
Water Resour. Res.
Transition probability-based indicator geostatistics
Math. Geol.
Modeling spatial variability with one- and multi-dimensional continuous Markov chains
Math. Geol.
Inverse problem in hydrogeology
Hydrogeol. J.
Joint estimation of transmissivity and storativity in a bedrock fracture
Water Resour. Res.
Scaling limits for the transient phase of local Metropolis-Hastings algorithms
J. R. Stat. Soc. Series B (Stat. Methodol.)
Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis Hastings algorithm
Water Resour. Res.
Facteurs d'échelle dans la hiérarchisation des écoulements au sein d'un aquifère karstique: analyse multi-échelles des propriétés hydrodynamiques et de transport de l'aquifère du Lez
Identifying fracture-zone geometry using simulated annealing and hydraulic-connection data
Water Resour. Res.
Modelling: picture perfect or abstract art?
Ground Water
Application of Markov analysis to the Banff Formation (Mississippian), Alberta
Math. Geol.
Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial Hessian approximations
SIAM J. Sci. Comput.
Note on the sampling distribution for the Metropolis-Hastings algorithm
Comm. Statist. Theory Methods
Riemann manifold Langevin and Hamiltonian Monte Carlo methods
J. R. Statist. Soc. B
Complexity
Ground Water
Stochastic simulation of categorical variables using a classification algorithm and simulated annealing
Math. Geology
Identification of the permeability distribution in soil by hydraulic tomography
Inverse Probl.
DRAM: efficient adaptive MCMC
Stat. Comput.
Estimation of co-conditional moments of transmissivity, hydraulic head and velocity fields
Adv. Water Resour.
Hydraulic tomography for detecting fracture zone connectivity
Ground Water
Transition probability-based stochastic geological modeling using airborne geophysical data and borehole data
Water Resour. Res.
Hydraulic tomography offers improved imaging of heterogeneity in fractured rocks
Groundwater
Hydraulic tomography in fractured granite: Mizunami Underground Research site, Japan
Water Resour. Res.
Comparison of aquifer characterization approaches through steady state groundwater model validation: a controlled laboratory sandbox study
Water Resour. Res.
Cited by (38)
Mapping conduits in two-dimensional heterogeneous karst aquifers using hydraulic tomography
2023, Journal of HydrologyCitation Excerpt :By interpreting the head responses and assessing the associated uncertainties, different kinds of HT inversion schemes have been reported in literature to determine the spatial distribution of hydraulic parameters. These inverse methods include general geostatistical inverse methods (Yeh and Liu, 2000; Zhu and Yeh, 2005; Li and Cirpka, 2006; Klein et al., 2017; Schwede et al., 2014; Cardiff et al., 2009, 2013; Mao et al., 2013a, 2013b; Zhao et al., 2015; Zhao and Illman, 2018; Li et al., 2021), travel time-based inverse methods (Brauchler et al., 2003, 2010, 2013; Vasco et al., 2019), pilot-point-based approaches (Lavenue and De Marsily, 2001; Castagna et al. 2011), sparse nonlinear optimizer and stochastic Newton approach (Wang et al., 2016, 2017), cellular automata-based approaches (Fischer et al., 2017a,b), and discrete network deterministic inversion method (Fischer et al., 2018, 2020). In this study, we focus on the geostatistics-based and travel time-based inversion methods of HT.
Bayesian inversion of laboratory experiments of transport through limestone fractures
2022, Journal of Contaminant Hydrology