Coupled hydrogeophysical inversion to identify non-Gaussian hydraulic conductivity field by jointly assimilating geochemical and time-lapse geophysical data

Journal of Hydrology 578 (2019) 124092

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Coupled hydrogeophysical inversion to identify non-Gaussian hydraulic conductivity field by jointly assimilating geochemical and time-lapse geophysical data

T

⁎

Xueyuan Kanga, Xiaoqing Shia, , André Revilb, Zhendan Caoc, Liangping Lic, Tian Lana, ⁎ Jichun Wua, a b c

Key Laboratory of Surficial Geochemistry of Ministry of Education, School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China Université Grenoble Alpes, Université Savoie-Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France Department of Geology and Geological Engineering, South Dakota School of Mines and Technology, Rapid City 57701, USA

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Patrick Lachassagne, Associate Editor

Reliable inversion of spatial heterogeneity of hydraulic conductivity is crucial to understand subsurface fluids migration. The Ensemble Smoother – Direct Sampling method (ES-DS) has proven to be an effective method to identify non-Gaussian hydraulic conductivity distributions by incorporating a variety of traditional hydrodynamic measurements, e.g., piezometric head. However, inversion problems for non-Gaussian parameters often suffer from a sparsity of the available data from direct sampling in boreholes. As a non-intrusive, cost-effective, and high sampling density method, time-lapse geophysical technique has not yet drawn much attention as a useful source of information for delineating the underlying non-Gaussian heterogeneity. In this study, we integrated coupled hydrogeophysical modeling and the ES-DS algorithm to estimate non-Gaussian hydraulic conductivity field by assimilating both geochemical and time-lapse geophysical datasets. Four synthetic Cases for a salt injection experiment, monitored by both sampling analysis and electrical resistivity tomography, are conducted to assess the ability of the proposed approach to characterize hydraulic properties by assimilating different types of data. Results show that using geochemical or geophysical data alone only allow a rough reconstruction of subsurface heterogeneity of aquifers but might lose the fine structure. By incorporating multisource datasets, the main patterns of the non-Gaussian reference fields can be reflected with an improved resolution. The time-lapse geophysical methods open up new opportunities to accurately characterize nonGaussian aquifers and monitor the dynamic processes of subsurface fluids.

Keywords: Non-Gaussian Parameters estimation Multi-source data assimilation Electrical resistivity tomography Ensemble smoother

1. Introduction Accurate estimation of hydrodynamic parameters such as hydraulic conductivity is a key issue for the understanding of groundwater flow and contaminant transport process (e.g., Gómez-Hernández et al., 2003). However, it is still a great challenge to characterize the heterogeneous aquifer directly, due to a scarcity of data from traditional hydrogeological investigations (Oliver and Chen, 2011). Characterization of hydrodynamic parameters can be better obtained from stochastic inversion approaches, using the direct data (e.g., hydraulic conductivities measurement) and indirect data (e.g., observed head and tracer data) (Zhou et al., 2014a). Among these inversion approaches, data assimilation methods have become popular and widely used in hydrogeology and petroleum

⁎

engineering communities since they can incorporate multi-source of information to jointly estimate model parameters and states (e.g., Ju et al., 2018; Oliver and Chen, 2011; Zheng et al., 2019; Zhou et al., 2014a). Various strategies have been proposed including the Kalman Filter (KF, Kalman, 1960) and its variants like the Ensemble Smoother (ES, VanLeeuwen and Evensen, 1996), the Ensemble Kalman Filter (EnKF, Evensen, 2009), and the Ensemble Smoother with Multiple Data Assimilation (ESMDA, Emerick and Reynolds, 2013). All of these methods assume however a multivariate Gaussian distribution for the parameters and state variables, which might be a questionable assumption in complex channelized aquifers. This is for instance the case when the distribution of the hydraulic conductivity is controlled by features such as paleo-channels with highly permeable preferential flow conduits (Zovi et al., 2017). The same can be found in presence of

Corresponding authors. E-mail addresses: [email protected] (X. Shi), [email protected] (A. Revil), [email protected] (L. Li), [email protected] (J. Wu).

https://doi.org/10.1016/j.jhydrol.2019.124092 Received 16 April 2019; Received in revised form 2 July 2019; Accepted 31 August 2019 Available online 05 September 2019 0022-1694/ © 2019 Elsevier B.V. All rights reserved.

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

complex sand-clay interbedding caused by marine transgressions and regressions (Baeteman et al., 1989). To address the challenge caused by the non-Gaussian priors, two strategies have been proposed in the literature. In the first set of approaches, reparameterization methods are used to transform parameters (e.g., log hydraulic conductivity) from non-Gaussian distributions into Gaussian distributions, such as, the normal-score EnKF (Zhou et al., 2011), iterative normal-score ensemble smoother (Li et al., 2018), the Gaussian anamorphosis (Schoeniger et al., 2012), and the deep learning method (Canchumuni et al., 2017). Another set of approaches are based on geostatistical simulation, which use a post-processing step to make the posterior parameters consistent with the prior probability distribution. For example, Sarma et al. (2008) developed a kernel EnKF to reconstruct the channelized hydraulic conductivity field based on the multiple-point geo-statistics (MPG). Jafarpour and Khodabakhshi (2011) introduced a probability conditioning method (PCM), which uses the updated parameters probabilities as secondary conditioning data to regenerate non-Gaussian conductivities field in the MPG. Laloy et al. (2018) proposed a deep neural network approach based on the MPG simulation to deal with the non-Gaussianity. Lan et al. (2019) combined the PCM and the ESMDA to estimate non-Gaussian fields of lithofacies indicators and hydraulic conductivities. Finally, Cao et al. (2018) coupled iterative ensemble smoother with Direct Sampling (an MPG method proposed by Mariethoz et al., 2010) via pilot points to inverse the non-Gaussian parameters. This approach can effectively incorporate dynamic hydraulic head data into groundwater model, and preserve the non-Gaussian structures as well. Although a wide range of methods have been proposed to characterize the non-Gaussian hydraulic parameters, most researchers focused on assimilating only flow and transport measurements, such as piezometric head and solute concentration data (Llopis-Albert and Capilla, 2010; Zhou et al., 2011; Li et al., 2012; Xu and GómezHernández, 2016, 2018). However, directly sampling observations are relatively sparse and the acquisition cost is quite high compared with geophysical methods. A relatively high uncertainty of characterization and predictions are expected with limited data sources. This uncertainty can be reduced by using non-invasive geophysical methods (especially electrical resistivity tomography, ERT and induced polarization, IP) with low costs and high sampling density (Binley et al., 2015). Geophysical methods have been broadly used in the characterization of Gaussian heterogeneity (Irving and Singha, 2010; Pollock and Cirpka, 2012; Camporese et al., 2015; Ahmed et al., 2016; Kang et al., 2018, 2019). However, only a few studies have been conducted in applying geophysical methods to identify non-Gaussian fields (Zovi et al, 2017; Gottschalk et al., 2017). Zovi et al. (2017) and Gottschalk et al. (2017) used ERT constraints on channel shape and geological properties, and then assimilated the state variables (e.g., water levels) or lithological data to characterize the non-Gaussian subsurface heterogeneity. Comparing to static ERT, time‐lapse ERT has a significant advantage in identifying subsurface structure and quantifying the dynamic of subsurface fluids (Jardani et al., 2013). However, to our knowledge, time-lapse geophysical methods have not been used to characterize the non-Gaussian fields to date. In this study, we leveraged the latest non-Gaussian data assimilation approach proposed by Cao et al. (2018) and built on our previous work (Kang et al., 2018, 2019) where multiphase fluids and ERT/IP forward model were coupled. The non-Gaussian hydraulic conductivity was calibrated by jointly incorporating concentration data and time-lapse ERT data. In this work, we first developed a data assimilation framework based on the ESMDA-DS algorithm, which was adapted from the ES-DS method proposed by Cao et al. (2018), and coupled hydrogeophysical inversion of concentration and ERT data. A synthetic salt injection experiment was used to evaluate the performance of the proposed framework. Four Cases with different types of datasets were quantitatively compared, including (1) using a small amount of

Fig. 1. Flowchart illustrating the coupled hydrogeophysical model (red dotted box) and data assimilation framework. Note that Ne is the number of the ensemble, Na is the number of iterations of the ensemble smoother with multiple data assimilation, K is the hydraulic conductivity, cl is the clay content, ϕ is the porosity, C is the solute concentration, ρ and ρa are resistivity and apparent resistivity, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

concentration data alone; (2) using ERT data alone; (3) combining a small amount of concentration and ERT data; and (4) using a large amount of concentration data. Further, the influence of the number of pilot points on the inversion is discussed. 2. Methods 2.1. Coupled hydrogeophysical model To simulate the resistivity anomaly caused by salt injection in nonGaussian heterogeneity field with the presence of channelized clay layers, a coupled forward model is developed among a flow and solute transport model, petrophysical models and a geophysical forward model. The coupled model is illustrated in Fig. 1 (red dotted box) and summarized in following steps: Step1. Specify the hydraulic conductivity (K) field To simulate a heterogeneous and non-Gaussian distributed hydraulic conductivity, we construct the channelized K field from a multiple-point geostatistical method, Direct Sampling (Mariethoz et al., 2010). This method directly samples the training image for a given data event. It can apply multiple‐point geo-statistics to continuous variables and to multivariate problems. Besides, it is fast, computationally efficient and easy to parallelize (Mariethoz et al., 2010; Meerschman et al., 2

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

unit volume of the aquifer [T−1]; ϕ is the total porosity; and Cs is the concentration of the source of the non-reactive component [ML−3]. Note that the salt concentration in this work is relatively low (< 1 g/L). Therefore, we can safely neglect its effect on density and buoyancy. The flow and solute transport forward models are solved by MODFLOW (Harbaugh et al., 2000) and MT3DMS (Zheng and Wang, 1999), respectively. Step 5. Convert C and φwcl , to ρ using the petrophysical models proposed by Revil et al. (2018) and Sen and Goode (1992). When the subsurface media is clay-free, the well-known Archie’s law (Archie, 1942) can be used as an appropriate physical model connecting the hydrological state to electrical properties. However, in presence of clay-rich materials, such a model may lead to errors (Patnode and Wyllie, 1950). To overcome this weakness, the resistivity model developed by Revil et al. (2018) is used to quantify the influence of the clay grains on resistivity,

2013). Notably, the hydraulic conductivity at each cell depends on the proportions of clay and sand. For the convenience of subsequent calculations, we transform K into permeability using:

k=

Kμ γ

(1)

where k denotes the permeability, µ is the dynamic viscosity (µ = 1.01 × 10−3 Pa·s at 25 °C), γ denotes the unit weight for water (γ = 9.8 × 103N/kg). Step 2. Calculate the volumetric clay content (cl) First, the lower threshold permeability of clay-free sand, ksd, is calculated by Revil and Cathles (1999) according to:

ksd =

2 dsd (ϕsd )3msd

(2)

24

where dsd is the grain diameter of sand (m), ϕsd represents the porosity of sand, and msd denotes the cementation exponent of sand. In this work, dsd = 1 × 10−4 m, ϕsd = 0.32, msd = 2.0 (Power et al., 2013). The resulting ksd = 4.5 × 10−13 m2, which corresponds to a fine sand. Permeability values larger than ksd are assumed as clay-free sand, whereas others are indicative of sand containing a fraction of clay. Then, Eq. (3) is used to calculate the volumetric clay content, cl (ratio of clay, including water, to the whole soil volume) in terms of the permeability of clay-sand mixtures: 1 3mcs

cl =

ksd

1 3mcs

ksd

1 = (Sw ϕ)mσw + (Sw ϕ)m − 1ρs (B − λ )(φwsd CECsd + φwcl CECcl) ρ

and where the pore water conductivity is related to the salinity and temperature according to the model of Sen and Goode (1992):

2.36 + 0.099T ⎞ 2 σw = (5.6 + 0.27T − 1.5 × 10−4T 2)·C − ⎛ ·C 3 ⎝ 1.0 + 0.214C ⎠

1

1 − ϕcl

(3)

ϕsd

where mcs represents the cementation exponent corresponding to the clayey sand mixtures and ϕcl is the porosity of pure clay. In this work, mcs = msd = 2.0, ϕcl = 0.56 (Power et al., 2013; Revil et al., 2017a). Therefore, we assumed cl = 0 for the grid cells where k > ksd. When k < ksd, the corresponding clay content cl can be calculated through Eq. (3). For more details about this equation, one can refer to Revil and Cathles (1999). Step 3. Calculate the total porosity (ϕ) and clay mass fractions (φwcl ) The total porosity is calculated from the clay content (cl) based on the equation developed by Berg (2007):

ϕ = ϕsd (1 − cl) + ϕcl cl

(4)

The mass fraction of the clay and sand grains cell are given by Power et al. (2013):

φwcl =

(φwcl ,

φwsd )

for each grid

cl (1 − ϕcl ) ρcl cl (1 − ϕcl ) ρcl + (1 − cl)(1 − ϕsd ) ρsd

(5)

φwsd = 1 − φwcl

(6)

where ρcl and ρsd denote the density of the clay and sand grains (in this work, ρcl = ρsd = 2650 kg/m3, the same as in Revil et al., 2017a). Step 4. Forward flow and solute transport model with K and ϕ In this work, the steady-state flow field can be obtained from the following governing equation (Bear, 1972):

1 ∇∙ ⎜⎛ ∇V ⎟⎞ = −Iδ (r) ⎝ρ ⎠

−1

where K is the hydraulic conductivity [LT ], and H denotes the hydraulic head [L]. The governing equation for non-reactive component transport is defined as (Zheng and Wang, 1999):

q ∂C = ∇ ·(D∇C ) − ∇ ·(vC ) + s C s ϕ ∂t

(11)

where ρ is the spatial distribution of resistivity in the media, Ω·m. V is the electric potential. δ is the Dirac delta function. r represents a single current electrode, idealized as a point source at the origin with current strength I (A). In this work, the geophysical forward problem is solved by AGI EarthImager software (AGI, 2010).

(7)

∇ ·(K∇H ) = 0

(10)

where ρ denotes the bulk resistivity of sand-clay mixtures, Sw is the saturation of water phase (in this study, Sw = 1.0), and m is the cementation exponent. We assumed the cementation and saturation exponents are equal (see Revil et al., 2018). The quantity σw denotes the pore water conductivity (S/m), ρs represents the mass density of the solid phase. The quantities B and λ (m2V−1s−1) represent the apparent mobility of the counterions for surface conduction and the polarization associated with the quadrature conductivity, respectively. Here, B = 4.1 × 10−9 m−2s−1V−1 and λ = 3.460 × 10−10 m−2s−1V−1 (Revil et al., 2017a, 2018). CEC (C/kg) represents the total amount of cations which can be sorbed on the mineral surface. CEC is proportional to the clay content (Mao et al., 2016). The influence of solute concentration and temperature upon the electrical conductivity of the electrolyte is described by the model of Sen and Goode (1992). In Eq. (10), T is temperature (°C) and C represents the ionic concentration (mol/L) associated with salinity. In this study, we assume that the temperature remains constant (25 °C). Step. 6 Simulate ρa (apparent resistivity) through geophysical forward model With the electrical properties (ρ) converted from petrophysical models, the raw geophysical observation data (ρa) can be obtained through DC resistivity forward model. Given the space-time distribution of bulk electrical resistivity, the forward model can predict apparent resistivity value by solving the partial differential equation:

− k 3mcs

( )

(9)

2.2. Ensemble smoother with multiple data assimilation - direct sampling (ESMDA-DS) In this work, the ESMDA-DS is used as an inversion tool to characterize the non-Gaussian heterogeneous parameter field (Y, loge- hydraulic conductivity). This method borrows computational efficiency from the ESMDA (Emerick and Reynolds, 2013) and leverages the capability of the DS (Mariethoz et al., 2010) for handling the non-

(8)

where C is the aqueous concentration of the non-reactive component [ML−3]; t denotes time [T]; D represents the diffusion coefficient [L2T−1]; v = (−K ∇H )/ ϕ [L2T−1]; qs means the volumetric flow rate per 3

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

Gaussianity. Note that we adapt the iterative ES-DS proposed by Cao et al. (2018) to the ESMDA-DS, because the ESMDA has less computational expense and better data matches than the iterative ES (see Emerick and Reynolds, 2013). This method can be summarized in following steps: Step 1. Build the augmented state vector x. Parameters of interest p (e.g., Y) are augmented with the state variables h (e.g., C and ρa) into a joint state vector x = [p h]T. Step 2. Decide the number of data assimilation (Na) and the coefficient (αi with i = 1, ...,Na ) for different times data assimilation. For i = 1 to Na, run Step 3 to 5. Step 3. For each realization, run the coupled hydrogeophysical model from time zero with updated parameters:

x if = G (x ia),

i = 1, 2, …, Ne

Table 1 Model parameters for the synthetic experiment.

(12)

where i is the ensemble member index, superscripts f and a denote forecast and analysis, respectively. Ne represents the number of the realizations in the ensemble. Through Eq. (12), we can obtain temporal evolutions of the state variables based on the model parameters, which can be used for update in Step 4. Step 4. Calculate the error covariance for the predicted state:

CYD =

CDD =

Ne

Parameter

Value

Petrophysical parameters Sand porosity Clay porosity Cementation exponent (sand) Cementation exponent (clay) CEC (sand) (C/kg) CEC (clay-bentonite) (C/kg) Water resistivity (before the injection) (Ω·m)

0.32a 0.56b 2.0a 2.0b 7.2a 4.8 × 104 b 20c

Flow model parameters Stress period Grid spacing (m) Model length (m) Model width (m) Model depth (m) Starting head (m)

1 0.5 × 0.5 × 0.5 40 0.5 20 100

Solute transport parameters Total simulation time (days) Source constant concentration (g/L) Longitudinal dispersivity (m) Horizontal transverse dispersivity (m) Molecular diffusion coefficient

100 1 0.5 0.05 0

a

1 Ne − 1

∑ (x if

1 Ne − 1

∑ (di − d¯)(di − d¯)T

− x¯ f )(di − d¯)T

b

(13)

i=1

c

Ne

3. Synthetic case study (14)

i=1

3.1. Case description

where CYD is the cross-covariance matrix between the forecast state and the predicted data, CDD is the covariance matrix of the predicted data, and di is the predicted data for the i th realization. Step 5. Update the ensemble by:

(

(Power et al., 2013). (Revil et al., 2017a). (Steelman et al., 2017).

)

x ia = x if + CYD (CDD + αl CD)−1 d obsi − di ,

i = 1, 2, …, Ne

To test the performance of the proposed framework, we conduct a synthetic salt injection experiment in a non-Gaussian setting. The solute transport process is monitored by both hydrogeological sampling and time-lapse ERT. The synthetic salt injection experiment is developed for a 2D heterogeneous confined aquifer (40 m × 20 m). The study area is discretized into 80 × 40 = 3200 blocks (the length of each block is 0.5 m, Δx = Δz = 0.5 m). No-flow boundary conditions are imposed at the upper and lower side of the aquifer. Constant head boundary conditions are set at x = 0 m (hydraulic head, H = 102 m) and x = 40 m (H = 100 m) to ensure groundwater flow. For the solute transport model, a linear inert tracer (NaCl) source with a constant concentration (C = 1 g/L) is set along a line at x = 1.0 m. The initial concentration is set to zero in the whole domain. The total simulation time is 100 days. The flow and transport parameters are shown in Table 1. For the ERT forward model, the domain boundaries are placed 240 m from the zone of interest to minimize boundary effects on geoelectric simulations. The non-Gaussian Y fields are generated through the DS simulation. The DS can model both categorical and continuous variables, based on the training image (Fig. 2a). In reality, the training image can be roughly inferred from expert information or geophysical methods (e.g., seismic technique). The resulting reference Y field is shown in Fig. 2b. The corresponding distribution of clay content is demonstrated in Fig. 2c. Note that the histogram of the reference Y field (Fig. 2d) shows a strong non-Gaussian character. The initial Y fields are generated from the same training image. Based on the reference Y field and the coupled hydrogeophysical model, we can obtain the dynamic process of solute transport (Fig. 2e) and the corresponding resistivity distribution (Fig. 2f). The number of multiple data assimilation in the ESMDA algorithm is chosen to be 4 (Na = 4). Larger Na might improve the characterization, because more iterations are performed. However, the computation burden will also be significantly increased. Therefore, Na should be set to achieve a trade-off between estimation accuracy and computation efficiency. We have pre-tested the performance of different Na. Based

(15)

where l is the times index of the ESMDA, l = 1,2,Na, CD is the covariance matrix of the measurements error, and dobs is the perturbed observations with noise of covariance αlCD. Step 6. Interpolate the updated Y at pilot points to un-sampled locations, using DS method. The extrapolated Y fields will be implemented for the next iteration. Notably, the pilot points play an important role in this approach. If there are no pilot points, the ESMDA-DS becomes the DS alone and therefore dynamic data cannot be incorporated into the MPG simulations. If all the nodes are pilot points, the proposed algorithm ends with the ESMDA, and the non-Gaussian patterns will not be preserved after assimilation. To handle the non-Gaussianity, we use the pilot points to bridge the ESMDA and DS for better flow and transport modeling.

2.3. Data assimilation framework Fig. 1 shows the flowing chart of the proposed data assimilation framework. First, an ensemble of non-Gaussian Y fields is generated by the DS. With the Y fields, the distribution of cl (clay content) and ϕ are obtained from Eqs. (3) and (4). Next, the groundwater model is propagated forward to get the evolution of solute concentration C. Using the petrophysical models, the distributions of resistivity is converted from C and cl. Through geophysical forward model, geophysical properties (ρ) are transferred to geophysical observations (ρa). Then, the ESMDA is used to update the ln K values at pilot points by incorporating geochemical and geophysical measurements (C and ρa). The updated ln K at the pilot points is used as hard data for the DS to regenerate the ln K fields. 4

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

Fig. 2. The properties of the reference field. (a) Training image for multiple-point geostatistics simulation. (b) ln K (loge hydraulic conductivity) distribution. (c) Distribution of clay content. (d) Histogram of the reference ln K field. (e) Reference concentration distribution at T = 50 d. (f) Reference resistivity distribution at T = 50 d.

acquired. Note that the resistivity data has a larger error, due to its lower accuracy. This level of noises is realistic in terms of real conditions (Jardani et al., 2013). Both the geochemical (C) and geophysical (ρa) measurements are acquired every 10 days, from t1 = 10 d to t10 = 100 d.

on the simulation results (not shown) and previous research (Emerick and Reynolds, 2013), Na is taken as 4 in this work. We set the coefficients of data assimilation as follow, α1 = 9.333, α2 = 7.0, α3 = 4.0 and α4 = 2.0, according to Emerick and Reynolds (2013). The realization size for ESMDA-DS is taken as 500. The number of the pilot points is 300. The spatial distribution of the pilot points is randomly generated (Fig. 3a) to demonstrate the robustness of the algorithm. Sensitivity to the number of pilot points will be further discussed in Section 5.

4. Results 4.1. Hydraulic conductivity estimation

3.2. Multiple observation data Fig. 4 demonstrates the ensemble mean, variance and histogram of the calibrated Y fields using different types of datasets. The initial fields (the first row in Fig. 4) are directly generated by the DS, without considering any dynamic data. Although each realization can preserve non-Gaussian feature, the ensemble mean only provides a relatively homogeneous result. Besides, the ensemble variance for initial fields is the largest, because no dynamic data are incorporated here. In Case 1, only limited amount of geochemical data (300 measurements in total) is assimilated. The limited data are very common in practice due to the high cost of traditional intrusive drilling method. After integrating these dynamic geochemical data, the ensemble mean can roughly delineate main preferential flow conduits of the reference. However, suffering from a scarcity of data, the estimated field fails to characterize the true heterogeneity structure in detail, e.g., the connectivity of the highly permeable preferential flow conduits is relatively poor. In addition, compared with the results from initial fields, overall ensemble variance in Case 1 declines considerably. Besides, the Y at the boundaries of channels have a larger uncertainty than those within the channels, which is consistent with the results of Cao et al. (2018) and Li et al. (2018). The histogram from Case 1 is closer to the reference, due to the integration of dynamic geochemical data. In Case 2, the low-cost and non-invasive ERT method is introduced to reconstruct the complex non-Gaussian heterogeneity. The ensemble mean can roughly recover the reference Y distribution. Compared with the results of Case 1, in Case 2, the preferential flow conduits are more connective and closer to the reference. It is because ERT can provide sufficient data (5610 measurements in total), which are highly sensitive to the conductive solute and clay material in subsurface. However, suffering from the low resolution of ERT, the results in Case 2 only produce a relatively smooth distribution of conductivities, which loses the fine structure of the reference. Also, lower measurement accuracy

To demonstrate the advantage of multi-source data assimilation, four Cases (Table 2) are designed to analyze the trade-off between the different types of datasets. The non-Gaussian ln K field is estimated using only a few concentration data in Case 1, only ERT data in Case 2, both concentration (small amount) and ERT data in Case 3, and a large amount of concentration data in Case 4. To characterize the heterogeneous aquifer, three vertical boreholes (the red dotted frames in Fig. 3b) are excavated at x = 10, 20 and 30 m for drilling sample and cross-borehole ERT measurements. To obtain the concentration data, 30 concentration observation points (the black triangles in Fig. 3b) are evenly distributed in three boreholes for Case 1 and 3. In Case 4, a large amount of concentration data is collected in 100 observation points (the red circles in Fig. 3b). To implement ERT measurements, 30 electrodes are installed in the vertical boreholes, and 21 electrodes are set at the top aquifer/aquitard contact (Fig. 3c). These electrodes are placed 2 m apart both in vertical boreholes and horizontal line. The dipole-dipole array (a commonly used electrode array in field investigation, e.g., Johnson et al. 2015) is used for the horizontal line, while a bipole-bipole configuration (e.g., Seferou et al., 2013) is used for the cross-well measurements. The resulting geophysical measurement vector includes 561 apparent resistivity data. Besides, the ERT data acquisition for this synthetic case is assumed to be taken as instantaneous snapshots. In other words, the duration required to record the geophysical data is small with the time required for the salt to move in subsurface. It follows that the solute transport can be neglected during each geophysical data acquisition. With the known electrode configuration and the reference Y distribution, we can obtain the reference C and ρa data through the coupled hydrogeophysical model. Adding 3% and 15% relative errors to the reference C and ρa, respectively, the noisy measurements are 5

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

self-potential and ERT for Gaussian permeability estimation via the pilot points and Markov chain Monte Carlo sampling methods. Their results also demonstrate that when incorporating multi-source datasets, the hydraulic aquifer properties can be better constrained. To further illustrate the advantage of integrating multi-source datasets, a large amount of concentration data (1000 measurements in total) is considered in Case 4 for comparison. Note that the estimated fields in Case 3 is even comparable to those in Case 4. However, a huge amount of geochemical data like Case 4 is impractical in reality. Therefore, combining sufficient geophysical data and a small amount of high-precision geochemical data is a good alternative. It can be easily implemented and can accurately depict the non-Gaussian fields. The results can also be investigated by comparing the cumulative distribution functions (CDFs) of the non-Gaussian Y ensembles for the initial fields and four Cases, with the CDF of the reference Y field (Fig. 5). The CDFs of the initial fields roughly illustrate the non-Gaussian character. However, their uncertainty is quite high and the ensemble mean is very different from the reference. Considering the geochemical data alone (Case 1), the dispersion of ensembles considerably declines, but the ensemble mean CDF cannot match the reference well. When using resistivity data alone (Case 2), the results suffer from an unacceptable huge uncertainty caused by the low resolution of ERT. By contrast, the results in Case 3 and 4 can better reconstruct the non-Gaussian structure with a lower uncertainty. To quantitatively analyze the performance of these four Cases, the root mean square error (RMSE) which is widely used to evaluate the performance of the inversion method (e.g., Chen and Zhang, 2006; Jardani et al., 2013) is applied as a criterion to measure the agreement between the estimated and reference distribution of Y:

RMSEY =

Table 2 Parameter sets for the four Cases of the synthetic salt injection experiment.

Number of C observation Number of ρa observation

Case 2

Case 3

Case 4

300 0

0 5610

300 5610

1000 0

Ng

∑ (Yr − Ya)2 i=1

(16)

where Ng is the number of grids, Yr and Ya represent the reference and the estimated (ensemble mean) distribution of loge- hydraulic conductivity, respectively. Optimally, RMSE should decline as the characterization improves. Fig. 6 displays scatterplots, the coefficient of determination (R2) and the RMSE of true versus estimated Y for initial fields and four Cases, respectively. These results indicate when using both geochemical and geophysical data sets (Case 3 and Case 4), the updated fields can match the reference better than when just a single dataset is considered (Case 1 and Case2).

Fig. 3. (a) The distribution of the 300 pilot points (black dots). (b) The distribution of the observation points for solute concentration. The black triangles represent observation points in Case 1 and 3, while the red circles correspond to Case 4. The red dotted frames correspond to the boreholes. (c) Sketch of the electrode (black dots) distribution for ERT measurements. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Case 1

1 Ng

4.2. Concentration data reproduction Besides directly assessing the characterization of hydraulic conductivity, we also evaluate the inversion results with the reproduction of concentration data. Based on the estimated Y fields from these four Cases, the solute transport model is rerun to obtain the corresponding C distributions. Fig. 7 demonstrates the ensemble mean, ensemble variance and the CDF of plume distribution (T = 50 d) for the initial fields and four Cases. Note that the C distribution calculated from the initial fields does not exhibit any decent features shown by the reference field. Besides, the variance and CDF maps demonstrate its uncertainty is quite high. By contrast, the results in Case 1–4 can reflect the reference C distribution to some extent. When a single dataset (i.e., geochemical or geophysical data, corresponding to Case 1 and 2, respectively) is used, the ESMDA-DS fail to characterize the fine architecture of the plumes. Moreover, variance and CDF maps show both two Cases still suffer from a relatively large uncertainty despite that the realizations slightly converge. Nevertheless, considering both geochemical and geophysical datasets, the results in Case 3 can better reflect the plumes distribution with a lower uncertainty. In addition, the estimation accuracy achieved by combining multi-source data is similar to the one achieved by using a large amount of drilling concentration data.

Note that C and ρa are solute concentration and apparent resistivity.

results in larger ensemble variance. To improve the resolution of ERT, the survey design optimization (e.g., Loke et al., 2014) and the image guided method based on GPR reflections (e.g. Zhou et al., 2014b) could be used to alleviate this problem in the future. In Case 3, we incorporate both the concentration and resistivity data to update the conductivities. The results show that the ensemble mean map can reconstruct the reference field with a high resolution. Both the connectivity of the channel and the spatial distribution of conductivities can be preferably recovered. This can be explained that ERT data can act as a good complement to concentration data. Despite its low resolution, ERT can provide sufficient data covering the entire study area. On the other hand, constrained by the high-accuracy concentration data, the estimated field can successfully converge to the reference with a low uncertainty. Similarly, Jardani et al. (2013) combined tracer test, 6

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

Fig. 4. Ensemble mean, variance and histogram (ln K) for the initial fields and four Cases. Note that the non-Gaussian ln K field is estimated using only a few concentration data in Case 1, only ERT data in Case 2, both concentration (small amount) and ERT data in Case 3, and a large amount of concentration data in Case 4.

discussed its influences on the inversion results. Fig. 8 demonstrates the ensemble means of updated hydraulic conductivity using different number of pilot points with ESMDA-DS, i.e., Npilot = 100, 300, 900 and 3200. Note that the ensemble mean is improved from Npilot = 100 to Npilot = 300, because the increasing number of pilot points help integrate more information from dynamic datasets. Also, the uncertainty of the estimated Y declines. When the number of pilot points is further increased from 300 to 900, the ensemble variance further decreases. However, when Npilot = 900, the histogram shows less non-Gaussianity since the ESMDA will turn those points into Gaussian distribution. When all the nodes are pilot points (Npilot = 3200), the ESMDA-DS turns into ESMDA alone, i.e., no DS is used for interpolation. The ensemble mean map from ESMDA shows overfitting and the uncertainty of the conductivities is underestimated. Besides, the histogram shows that after ESMDA assimilation, the estimated Y fields have almost become a Gaussian distribution, losing the bimodal distribution of the reference. Meanwhile, based on the estimated Y field by using different

To quantitatively evaluate the capability of C reproduction for these four Cases, we calculate the RMSE for CT=50d distributions (Table 3). It indicates the estimated C distribution in Case 3 is more accurate than those in Case 1 and 2. Obviously, integrating multi-source datasets can preferably reproduce the concentration data. All the computing work are run on a workstation with eight-core 2.6 GHz CPU and 32 GB RAM. The computation time for the four Cases are shown in Table 3. The main cost of the total CPU time is spent by the forward groundwater model and geophysical model. Note that incorporating ERT data (Case 2 and 3) requires additional calls of the forward ERT model, and therefore, results in a longer computation time.

5. Discussion One important issue for the ESMDA-DS algorithm is the optimization of the number of the pilot points. Therefore, we performed a sensitivity analysis on the number of pilot points based on Case 3 and 7

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

Fig. 5. Cumulative distribution function (CDF) of the initial Y fields, the final ensemble of retrieved Y fields for four Cases. Green, blue and red lines indicate the realizations, the reference and the ensemble mean for CDF, respectively. Note that the non-Gaussian ln K field is estimated using only a few concentration data in Case 1, only ERT data in Case 2, both concentration (small amount) and ERT data in Case 3, and a large amount of concentration data in Case 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

increase from 100 to 300, both the estimations for Y and CT=50d are improved. On the other hand, further increase in the number of pilot points might provide a worse estimation due to the loss of non-Gaussianity of the reference. In summary, a small amount of pilot points will not be enough to assimilate indirect data, while a redundant number of pilot points will

number of pilot points, the solute transport model is rerun to obtain the corresponding C distributions (Fig. 9). Fig. 9 illustrates that the uncertainty of the resulting concentration fields gradually declines with the elevated number of pilot points. To make a quantitative comparison, we calculate the RMSE for Y and CT=50d distributions (Table 4). When the number of pilot points

Fig. 6. Scatterplots of true versus estimated ensemble mean ln K for initial fields and different Cases. Note that the non-Gaussian ln K field is estimated using only a few concentration data in Case 1, only ERT data in Case 2, both concentration (small amount) and ERT data in Case 3, and a large amount of concentration data in Case 4. RMSE and R2 represent the root mean square error and the coefficient of determination, respectively. 8

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

Fig. 7. Ensemble mean, ensemble variance and cumulative distribution function (CDF) of the final estimated ensembles of C fields (T = 50 d) for four Cases and initial fields. For the third column, green, blue and red lines indicate the realizations, the reference and the ensemble mean for CDF, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

fail to preserve the non-Gaussian structure of the reference field. Therefore, the number of pilot points should be carefully set to achieve a trade-off between data assimilation and non-Gaussianity preservation (see Cao et al., 2018).

data assimilation framework integrating the coupled hydrogeologicalERT forward model and the Ensemble Smoother with Multiple Data Assimilation – Direct Sampling (ESMDA-DS) method to inverse the nonGaussian parameter field. To test the performance of the proposed framework, we conducted a synthetic salt tracer injection experiment in a non-Gaussian hydraulic conductivity field. This injection experiment is monitored by simulating both hydrogeological sampling and timelapse ERT. Quantitative comparisons for four Cases are performed, using the following constraints: (1) a small amount of concentration

6. Conclusions In order to overcome the difficulties due to a sparsity of observations from traditional hydrogeological investigations, we developed a 9

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

number of pilot points on the results are further discussed. The results show that using a single dataset (concentration or resistivity), the heterogeneous aquifer can be roughly reconstructed. However, neither geochemical nor geophysical data alone can characterize the fine heterogeneous structure. Combining geochemical and geophysical datasets can provide more detailed information of the nonGaussian field. Indeed, these two datasets preferably complement each other. ERT provides low-accuracy but sufficient data covering the entire study area. As a supplement, a small amount of highly accurate concentration data can successfully reduce the uncertainty of the nonGaussian inversion. Incorporating multi-source datasets can even achieve a similar estimation accuracy as using a large amount of concentration data. Also, results from the solute transport predictions demonstrate that integrating geochemical and geophysical data can better delineate the morphology of the plumes. For the first time, the time-lapse geophysical datasets are used to identify the non-Gaussian pattern of the hydraulic conductivity field. As

Table 3 Performance of the synthetic Cases in terms of root mean square error of the retrieved Y fields and concentration distributions by forwarding the model with the inverted Y fields.

RMSEY RMSEC Computation time

Initial field

Case 1

Case 2

Case 3

Case 4

2.34 0.76 /

2.07 0.45 20.2 h

1.75 0.37 31.6 h

1.63 0.26 31.9 h

1.55 0.24 20.4 h

Note that RMSEY, and RMSEC represent the root mean square error for logehydraulic conductivity field and concentration distribution at T = 50 d, respectively.

data alone; (2) resistivity data alone; (3) both concentration and ERT data; (4) a large amount of concentration data alone. Then the solute transport predictions are performed to evaluate the reproductive capability of the estimated hydraulic conductivities. The influence of the

Fig. 8. Ensemble mean, variance and histogram (ln K) for different number of pilot points based on Case 3 (both concentration data and ERT data are assimilated). Note that Npilot is the number of pilot points. 10

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

Fig. 9. Ensemble mean, ensemble variance and cumulative distribution function (CDF) of the final estimated ensembles of C fields (T = 50 d) for different number of pilot points based on Case 3 (both concentration data and ERT data are assimilated). For the third column, green, blue and red lines indicate the realizations, the reference and the ensemble mean for CDF, respectively. Note that Npilot is the number of pilot points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

non-intrusive, low-cost and high sampling density approaches, timelapse geophysical methods will open up new opportunities to identify the complex non-Gaussian subsurface heterogeneity and characterize the contaminant plumes. Note that we used a simple synthetic 2D example with two lithofacies to illustrate effectiveness of proposed framework in this work. In reality, the aquifer could be 3D and multivariate non-Gaussian. As a result, solving 3D non-Gaussian inversion problem will be more timeconsuming. But, in principle, this framework can be easily applied to realistic cases where the subsurface heterogeneity is more complicated. Moreover, the ESMDA-DS supports parallel computing, which can effectively alleviate computational burden. In addition, the petrophysical models are assumed to be perfect in this study. However, these models may face some uncertainties in reality and the petrophysical uncertainty might have an impact on the

Table 4 Performance of different number of pilot points for Case 3 (both ERT and concentration data are incorporated) in terms of root mean square error of the retrieved Y fields and concentration distributions by forwarding the model with the inverted Y fields. Number of pilot points

100

300

900

3200

RMSEY RMSEC

1.73 0.43

1.63 0.26

1.8 0.28

1.98 0.24

Note that RMSEY, and RMSEC represent the root mean square error for logehydraulic conductivity field and concentration distribution at T = 50 d, respectively.

11

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

inversion result (Linde et al., 2017; Brunetti and Linde, 2018). For instance, the petrophysical parameters such as cementation exponent might show some heterogeneity just like the permeability (Linde and Doetsch, 2016). Therefore, further studies in this field should account for petrophysical uncertainty to obtain a more realistic uncertainty assessment of the coupled inversion. Note that the characterization of the hydraulic properties from tracing test only provides an accurate estimation of the highly permeable channels. However, it suffers from some uncertainty in the low permeable area. Thus, to improve the characterization, a future work will concern to further incorporate the induced polarization (IP) method due to its high sensitivity to the clay content and the lowpermeability area (see Revil 2012; Revil et al., 2017b).

J.E., 2017. Integrating non-colocated well and geophysical data to capture subsurface heterogeneity at an aquifer recharge and recovery site. J. Hydrol. 555, 407–419. https://doi.org/10.1016/j.jhydrol.2017.10.028. Harbaugh, A.W., Banta, E.R., Hill, M.C., McDonald, M.G., 2000. MODFLOW-2000, the U. S. Geological Survey Modular Ground-water Model. U.S. Geological Survey, Branch of Information Services, Reston, VA, Denver CO. Irving, J., Singha, K., 2010. Stochastic inversion of tracer test and electrical geophysical data to estimate hydraulic conductivities. Water Resour. Res. 46 (11), W11514. https://doi.org/10.1029/2009WR008340. Jafarpour, B., Khodabakhshi, M., 2011. A probability conditioning method (PCM) for nonlinear flow data integration into multipoint statistical facies simulation. Math. Geosci. 43 (2), 133–164. https://doi.org/10.1007/s11004-011-9316-y. Jardani, A., Revil, A., Dupont, J.P., 2013. Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging. Adv. Water Resour. 52, 62–77. https://doi.org/10.1016/j.advwatres.2012.08.005. Johnson, T., Versteeg, R., Thomle, J., Hammond, G., Chen, X., Zachara, J., 2015. Fourdimensional electrical conductivity monitoring of stage-driven river water intrusion: accounting for water table effects using a transient mesh boundary and conditional inversion constraints. Water Resour. Res. 51 (8), 6177–6196. https://doi.org/10. 1002/2014WR016129. Ju, L., Zhang, J., Meng, L., Wu, L., Zeng, L., 2018. An adaptive Gaussian process-based iterative ensemble smoother for data assimilation. Adv. Water Resour. 115, 125–135. https://doi.org/10.1016/j.advwatres.2018.03.010. Kalman, R.E., 1960. A new approach to linear filtering and prediction problems. Trans. ASME-J. Basic Eng. 82, 35–45. https://doi.org/10.1115/1.3662552. Kang, X., Shi, X., Deng, Y., Revil, A., Xu, H., Wu, J., 2018. Coupled hydrogeophysical inversion of DNAPL source zone architecture and permeability field in a 3D heterogeneous sandbox by assimilation time-lapse cross-borehole electrical resistivity data via ensemble Kalman filtering. J. Hydrol. 567, 149–164. https://doi.org/10.1016/j. jhydrol.2018.10.019. Kang, X., Deng, Y., Revil, A., Shi, X., Shi, L., Wu, J., 2019. Combined induced polarization and electrical resistivity tomography for characterizing DNAPL source zone architecture in low-permeability media. J. Contam. Hydrol. (Under Review). Laloy, E., Hérault, R., Jacques, D., Linde, N., 2018. Training-image based geostatistical inversion using a spatial generative adversarial neural network. Water Resour. Res. 54 (1), 381–406. https://doi.org/10.1002/2017WR022148. Lan, T., Shi, X., Chen, Y., Wu, J., 2019. A comparison of cumulative density function and probability conditioning method based on the ES-MDA to estimate non-Gaussian parameters for groundwater modeling. Adv. Water Resour. (Under Review). Li, L., Zhou, H., Franssen, H.J.H., 2012. Jointly mapping hydraulic conductivity and porosity by assimilating concentration data via ensemble Kalman filter. J. Hydrol. 428–429 (1), 152–169. https://doi.org/10.1016/j.jhydrol.2012.01.037. Li, L., Stetler, L., Cao, Z., Davis, A., 2018. An iterative normal-score ensemble smoother for dealing with non-Gaussianity in data assimilation. J. Hydrol. 567, 759–766. https://doi.org/10.1016/j.jhydrol.2018.01.038. Linde, N., Doetsch, J., 2016. Joint inversion in hydrogeophysics and near-surface geophysics. ch.7 In: Moorkamp, M., Lelievre, P., Linde, N., Khan, A. (Eds.), Integrated Imaging of the Earth. John Wiley & Sons, Inc., Hoboken, New Jersey, pp. 119–135. https://doi.org/10.1002/9781118929063.ch7. Linde, N., Ginsbourger, D., Irving, J., Nobile, F., Doucet, A., 2017. On uncertainty quantification in hydrogeology and hydrogeophysics. Adv. Water Resour. 110, 166–181. https://doi.org/10.1016/j.advwatres.2017.10.014. Llopis-Albert, C., Capilla, J.E., 2010. Stochastic inverse modelling of hydraulic conductivity fields taking into account independent stochastic structures: a 3D case study. J. Hydrol. 391 (3–4), 277–288. https://doi.org/10.1016/j.jhydrol.2010.07. 028. Loke, M.H., Wilkinson, P.B., Uhlemann, S.S., Chambers, J.E., Oxby, L.S., 2014. Computation of optimized arrays for 3-D electrical imaging surveys. Geophys. J. Int. 199 (3), 1751–1764. https://doi.org/10.1093/gji/ggu357. Mao, D., Revil, A., Hinton, J., 2016. Induced polarization response of porous media with metallic particles – part 4: detection of metallic and nonmetallic targets in timedomain induced polarization tomography. Geophysics 81 (4), D359–D375. https:// doi.org/10.1190/geo2015-0480.1. Mariethoz, G., Renard, P., Straubhaar, J., 2010. The Direct Sampling method to perform multiple-point geostatistical simulations. Water Resour. Res. 46, W11536. https:// doi.org/10.1029/2008WR007621. Meerschman, E., Pirot, G., Mariethoz, G., Straubhaar, J., Van Meirvenne, M., Renard, P., 2013. A practical guide to performing multiple-point statistical simulations with the Direct Sampling algorithm. Comput. Geosci. 52, 307–324. https://doi.org/10.1016/j. cageo.2012.09.019. Oliver, D.S., Chen, Y., 2011. Recent progress on reservoir history matching: a review. Comput. Geosci. 15 (1), 185–221. https://doi.org/10.1007/s10596-010-9194-2. Patnode, H.W., Wyllie, M.R.J., 1950. The presence of conductive solids in reservoir rocks as a factor in electric log interpretation. J. Pet. Technol. 2 (2), 47–52. https://doi.org/ 10.2118/950047-G. Pollock, D., Cirpka, O.A., 2012. Fully coupled hydrogeophysical inversion of a laboratory salt tracer experiment monitored by electrical resistivity tomography. Water Resour. Res. 48 (1), W01505. https://doi.org/10.1029/2011WR010779. Power, C., Gerhard, J.I., Tsourlos, P., Giannopoulos, A., 2013. A new coupled model for simulating the mapping of dense nonaqueous phase liquids using electrical resistivity tomography. Geophysics 78, N1–N15. https://doi.org/10.1190/geo2012-0395.1. Revil, A., 2012. Spectral induced polarization of shaly sands: influence of the electrical double layer. Water Resour. Res. 48, W02517. https://doi.org/10.1029/ 2011WR011260. Revil, A., Cathles, L.M., 1999. Permeability of shaly sands. Water Resour. Res. 35 (3), 651–662. https://doi.org/10.1029/98WR02700.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work was financially supported by the National Key Research and Development Program of China (No. 2018YFC0406402) and the National Natural Science Foundation of China (No. 41672229 and U1503282). We are grateful to the Associate Editor and two anonymous reviewers for their insightful comments and suggestions, which significantly improve the quality of this work. References AGI, 2010. EarthImagerTM 2D and 3D, Resistivity and IP Inversion Software. Advanced Geosciences Incorporated. http://www.agiusa.com. Ahmed, A.S., Jardani, A., Revil, A., Dupont, J.P., 2016. Specific storage and hydraulic conductivity tomography through the joint inversion of hydraulic heads and selfpotential data. Adv. Water Resour. 89, 80–90. https://doi.org/10.1016/j.advwatres. 2016.01.006. Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Ind. Min. Metall. Pet. Eng. 146 (01), 54–62. https://doi. org/10.2118/942054-g. Baeteman, C., Bogemans, F., Paepe, R., 1989. The Upper Quaternary Deposits of the Changjiang Coastal Plain, Shanghai Area. Belgian Geological Survey. http://www. vliz.be/imisdocs/publications/ocrd/260065.pdf. Bear, J., 1972. Dynamics of Fluids in Porous Materials. Dover, New York. Berg, C., 2007. An effective medium algorithm for calculating water saturation at any frequency. Geophysics 72 (2), E59–E67. https://doi.org/10.1190/1.2432261. Binley, A., Hubbard, S.S., Huisman, J.A., Revil, A., Robinson, D.A., Singha, K., Slater, L.D., 2015. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resour. Res. 51, 3837–3866. https://doi.org/ 10.1002/2015WR017016. Brunetti, C., Linde, N., 2018. Impact of petrophysical uncertainty on Bayesian hydrogeophysical inversion and model selection. Adv. Water Resour. 111, 346–359. https://doi.org/10.1016/j.advwatres.2017.11.028. Camporese, M., Cassiani, G., Deiana, R., Salandin, P., Binley, A., 2015. Coupled and uncoupled hydrogeophysical inversions using ensemble Kalman filter assimilation of ERT-monitored tracer test data. Water Resour. Res. 51, 3277–3291. https://doi.org/ 10.1002/2014WR016017. Canchumuni, S.A., Emerick, A.A., Pacheco, M.A., 2017. Integration of ensemble data assimilation and deep learning for history matching facies models. OTC Brasil. Offshore Technology Conference. Cao, Z., Li, L., Chen, K., 2018. Bridging iterative Ensemble Smoother and multiple-point geostatistics for better flow and transport modeling. J. Hydrol. 565, 411–421. https://doi.org/10.1016/j.jhydrol.2018.08.023. Chen, Y., Zhang, D., 2006. Data assimilation for transient flow in geologic formations via ensemble Kalman filter. Adv. Water Resour. 29 (8), 1107–1122. https://doi.org/10. 1016/j.advwatres.2005.09.007. Emerick, A.A., Reynolds, A.C., 2013. Ensemble smoother with multiple data assimilation. Comput. Geosci. 55 (3), 3–15. https://doi.org/10.1016/j.cageo.2012.03.011. Evensen, G., 2009. The ensemble Kalman filter for combined state and parameter estimation. IEEE Contr. Syst. Mag. 29, 83–104. https://doi.org/10.1109/MCS.2009. 932223. Gómez-Hernández, J.J., Franssen, H., Sahuquillo, A., 2003. Stochastic conditional inverse modeling of subsurface mass transport: a brief review and the self-calibrating method. Stochastic Environ. Res. Risk Assess. 17 (5), 319–328. https://doi.org/10. 1007/s00477-003-0153-5. Gottschalk, I.P., Hermans, T., Knight, R., Caers, J., Cameron, D.A., Regnery, J., Mccray,

12

Journal of Hydrology 578 (2019) 124092

X. Kang, et al.

10.1175/1520-0493(1996) 124<2898:DAAIMI>2.0.CO;2. Xu, T., Gómez-Hernández, J.J., 2016. Characterization of non-Gaussian conductivities and porosities with hydraulic heads, solute concentrations, and water temperatures. Water Resour. Res. 52 (8), 6111–6136. https://doi.org/10.1002/2016WR019011. Xu, T., Gómez-Hernández, J.J., 2018. Simultaneous identification of a contaminant source and hydraulic conductivity via the restart normal-score ensemble Kalman filter. Adv. Water Resour. 112, 106–123. https://doi.org/10.1016/j.advwatres.2017. 12.011. Zheng, C., Wang, P.P., 1999. MT3DMS: a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentation and user's guide. AJR Am. J. Roentgenol. 169 (4), 1196–1197. Zheng, Q., Zhang, J., Xu, W., Wu, L., Zeng, L., 2019. Adaptive multifidelity data assimilation for nonlinear subsurface flow problems. Water Resour. Res. 55 (1), 203–217. https://doi.org/10.1029/2018WR023615. Zhou, H., Gómez-Hernández, J.J., Hendricks Franssen, H., Li, L., 2011. An approach to handling non-Gaussianity of parameters and state variables in ensemble Kalman filtering. Adv. Water Resour. 34 (7), 844–864. https://doi.org/10.1016/j.advwatres. 2011.04.014. Zhou, H., Gómez-Hernández, J.J., Li, L., 2014a. Inverse methods in hydrogeology: evolution and recent trends. Adv. Water Resour. 63 (2), 22–37. https://doi.org/10.1016/ j.advwatres.2013.10.014. Zhou, J., Revil, A., Karaoulis, M., Hale, D., Doetsch, J., Cuttler, S., 2014b. Image-guided inversion of electrical resistivity data. Geophys. J. Int. 197, 292–309. https://doi.org/ 10.1093/gji/ggu001. Zovi, F., Camporese, M., Hendricks Franssen, H., Huisman, J.A., Salandin, P., 2017. Identification of high-permeability subsurface structures with multiple point geostatistics and normal score ensemble Kalman filter. J. Hydrol. 548, 208–224. https:// doi.org/10.1016/j.jhydrol.2017.02.056.

Revil, A., Sleevi, M.F., Mao, D., 2017a. Induced polarization response of porous media with metallic particles — part 5: influence of the background polarization. Geophysics 82 (2), E77–E96. https://doi.org/10.1190/GEO2016-0388.1. Revil, A., Coperey, A., Shao, Z., Florsch, N., Fabricius, I.L., Deng, Y., Delsman, J.R., Pauw, P.S., Karaoulis, M., de Louw, P.G.B., van Baaren, E.S., Dabekaussen, W., Menkovic, A., Gunnink, J.L., 2017b. Complex conductivity of soils. Water Resour. Res. 53 (8), 7121–7147. https://doi.org/10.1002/2017WR020655. Revil, A., Qi, Y., Ghorbani, A., Ahmed, A.S., Ricci, T., Labazuy, P., 2018. Electrical conductivity and induced polarization investigations at Krafla volcano, Iceland. J. Volcanol. Geotherm. Res. 368, 73–90. https://doi.org/10.1016/j.jvolgeores.2018.11. 008. Sarma, P., Durlofsky, L.J., Aziz, K., 2008. Kernel principal component analysis for efficient, differentiable parameterization of multipoint geostatistics. Math. Geosci. 40 (1), 3–32. https://doi.org/10.1007/s11004-007-9131-7. Schoeniger, A., Nowak, W., Franssen, H.H., 2012. Parameter estimation by ensemble Kalman filters with transformed data: approach and application to hydraulic tomography. Water Resour. Res. 48 (4), W04502. https://doi.org/10.1029/ 2011WR010462. Seferou, P., Soupios, P., Kourgialas, N.N., Dokou, Z., Karatzas, G.P., Candasayar, E., Papadopoulos, N., Dimitriou, V., Sarris, A., Sauter, M., 2013. Olive-oil mill wastewater transport under unsaturated and saturated laboratory conditions using the geoelectrical resistivity tomography method and the FEFLOW model. Hydrogeol. J. 21, 1219–1234. https://doi.org/10.1007/s10040-013-0996-x. Sen, P.N., Goode, P.A., 1992. Influence of temperature on electrical conductivity on shaly sands. Geophysics 57 (1), 89–96. https://doi.org/10.1190/1.1443191. Steelman, C.M., Meyer, J.R., Parker, B.L., 2017. Multidimensional investigation of bedrock heterogeneity/unconformities at a DNAPL-impacted site. Groundwater 55, 532–549. https://doi.org/10.1111/gwat.12514. VanLeeuwen, P.J., Evensen, G., 1996. Data assimilation and inverse methods in terms of a probabilistic formulation. Mon. Weather Rev. 124 (12), 2898–2913. https://doi.org/

13