Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging

Advances in Water Resources 52 (2013) 62–77

Contents lists available at SciVerse ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging A. Jardani a,⇑, A. Revil b,c, J.P. Dupont a a

Université de Rouen, M2C, UMR 6143, CNRS, Morphodynamique Continentale et Côtière, Bât. IRESE A, Mont Saint Aignan, France Colorado School of Mines, Dept. of Geophysics, Golden, CO, USA c ISTerre, CNRS, UMR CNRS 5275, Université de Savoie, 73376 Cedex, Le Bourget du Lac, France b

a r t i c l e

i n f o

Article history: Received 15 May 2012 Received in revised form 16 August 2012 Accepted 18 August 2012 Available online 7 September 2012 Keywords: Inverse problem Electrical data Salt tracer test Hydraulic conductivity Self-potential Resistivity tomography

a b s t r a c t The assessment of hydraulic conductivity of heterogeneous aquifers is a difficult task using traditional hydrogeological methods (e.g., steady state or transient pumping tests) due to their low spatial resolution. Geophysical measurements performed at the ground surface and in boreholes provide additional information for increasing the resolution and accuracy of the inverted hydraulic conductivity field. We used a stochastic joint inversion of Direct Current (DC) resistivity and self-potential (SP) data plus in situ measurement of the salinity in a downstream well during a synthetic salt tracer experiment to reconstruct the hydraulic conductivity field between two wells. The pilot point parameterization was used to avoid over-parameterization of the inverse problem. Bounds on the model parameters were used to promote a consistent Markov chain Monte Carlo sampling of the model parameters. To evaluate the effectiveness of the joint inversion process, we compared eight cases in which the geophysical data are coupled or not to the in situ sampling of the salinity to map the hydraulic conductivity. We first tested the effectiveness of the inversion of each type of data alone (concentration sampling, self-potential, and DC resistivity), and then we combined the data two by two. We finally combined all the data together to show the value of each type of geophysical data in the joint inversion process because of their different sensitivity map. We also investigated a case in which the data were contaminated with noise and the variogram unknown and inverted stochastically. The results of the inversion revealed that incorporating the self-potential data improves the estimate of hydraulic conductivity field especially when the selfpotential data were combined to the salt concentration measurement in the second well or to the time-lapse cross-well electrical resistivity data. Various tests were also performed to quantify the uncertainty in the inverted hydraulic conductivity field. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The steady-state ground water flow and solute transport are mainly controlled by the spatial distribution of the permeability and dispersivity of an aquifer. Permeability can vary over 12 orders of magnitudes and can be very heterogeneous at various scales, which in turn implies a complex pattern for groundwater flow and contaminant transports [1]. At the opposite, porosity and dispersivity does not exhibit such a broad range of variation. The hydraulic conductivity is most commonly estimated from invasive hydrogeological techniques, such as pumping tests. Despite recent advances in hydraulic tomography [2–4] the resolution of the inverted hydraulic conductivity depends strongly on the density of piezometers [5]. The limited number of available piezometers makes the reconstruction of the hydraulic conductivity field of heterogeneous ⇑ Corresponding author. E-mail addresses: [email protected] (A. Jardani), arevil@mines. edu (A. Revil), [email protected] (J.P. Dupont). 0309-1708/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.advwatres.2012.08.005

aquifers a difficult problem in most practical cases [6]. The use of geophysical methods can provide additional complementary information as broadly acknowledged in the last decade [7,8]. Recently, geophysical tools have benefited from (i) the evolution of the efficiency of numerical methods (for instance the mixed finite element approach) for solving partial differential equations and parallel computing [9,10], (ii) the development of improved petrophysical models connecting the geophysical signature to the hydraulic properties [11–14], and (iii) significant improvements in the technology of various sensors and filtering techniques with the possibility to developed multitask sensors. These developments have therefore given birth to a new era of three-dimensional time lapse geophysical imaging for tracking the changes of variables of interest like the moisture content, the salinity, and the pore fluid pressure [7,9–12,15]. Along these lines, the electrical resistivity imaging (ERI) is sensitive to changes in pore water electrical conductivity and temperature and therefore it has been used to track the subsurface

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

migration of conductive tracers (saline or heat tracer tests) with the goal to image the hydraulic conductivity of heterogeneous aquifers [16]. Pollock and Cirpka [17] presented recently a new analysis of the ERI data to image the distribution of the hydraulic conductivity. Their work is based on the analysis of temporal moments of potential electrical disturbances recorded during saline tracer tests. Various geophysical methods with different sensitivity maps can be also used in concert. For instance, Direct Current (DC) resistivity data can be also jointly inverted with Ground Penetrating Radar (GPR) data during salt tracer tests in the vadose zone, to determine the hydraulic conductivity and petrophysical properties such as electrical formation factor, the water content, and the effective grain radius of the sediments [18,19]. In this paper, we are interested by looking at the value of adding self-potential (SP) measurements in the joint inversion of geophysical data and in situ salt tracer sampling. The goal is still the same: the inversion of the hydraulic conductivity of an heterogeneous aquifer between two wells. The self-potential signals are passively recorded electrical potential signals associated with the occurrence of natural (source) currents in the ground. In the case of a salt tracer test, the source current density is generated by two contributions (i) the gradient in the activity of the salt (diffusion current) and (ii) an electrokinetic coupling (streaming current) directly associated with the flow of the ground water. The occurrence of this second contribution has been recently used to non-intrusively assess ground water flow For instance, Jardani et al. [20] used SP data to reconstruct the hydraulic head variations associated with pumping tests conducted in an alluvial aquifer. Suski et al. [21] showed that SP signals can be used to monitor the variations of the piezometric levels of an unconfined aquifer associated with an infiltration experiment from a ditch. SP data have been combined with the hydraulic head variations recorded during pumping tests to estimate the transmissivity of an heterogeneous aquifer [22,23] and to image in 3D the hydraulic conductivity [24]. The SP responds associated with the diffusion source current is recognized as an efficient method to delineate contaminated areas [25] and salt tracer tests [26]. We feel that the time-lapse analysis of the self-potential signals associated with a salt tracer test could emerged as a powerful tool in determining quantitatively hydraulic properties. As a side note, other geophysical methods have been conducted in order to evaluate the hydraulic conductivity of heterogeneous aquifers. For instance, Linde et al. [27] used cross-well ground penetrating radar (GPR) for salt tracer tests. Hyndman et al. [28] presented recently a relationship between the seismic slowness and hydraulic conductivity, which has been successfully used to predict the hydraulic conductivity of an alluvial aquifer from seismic and tracer test data. Hördt et al. [29] proposed the use of spectral induced polarization measurements to image the hydraulic conductivity of a sandy/gravel aquifer. The big picture is that combining hydrogeophysical and hydrogeological information need to be somehow coupled to reduce the uncertainty associated with the estimation of the hydraulic conductivity. This implies to take into account the uncertainty associated with the in situ measurements as well as with the geophysical data. Two approaches can be used to combine hydrogeological and geophysical data. In the first approach, the hydrological information is used to weakly constrain the inversion of the geophysical measurements [27]. The second approach is to fully couple the inversion of the two types of data [7,17,30]. Hinnell et al. [31] devoted an approach using the coupling of hydrological and geophysical data to reconstruct the hydrological parameters using ERI for tracking the infiltration front in vadose zone in a synthetic case study. In this paper, we are interested by the fully coupled joint inversion of geophysical and hydrogeological data to monitor a salt

63

tracer test and to assess the hydraulic conductivity field. Using salt tracer tests with ERI is a problem that has been received a lot attention recently [32–37]. The non-uniqueness of the ERI inverse problem (and its sensitivity map) makes this method insufficient by itself. In other words, ERI needs to be combined with other sources of information like in situ measurements of the salt concentration in wells [33,38,39]. Irving and Singha [38] introduced a stochastic joint inversion approach of time-lapse cross-well ERI and salt tracer concentration data. Our work is following this idea but adding an additional method, the self-potential method, to the inverse problem. We use the pilot points approach for the joint inversion of ERI, SP, and salt tracer data because this approach reduces the over parametrization of the problem. While previous authors used deterministic methods to choose the values of the pilot points, we place the pilot points on a regular grid and we use a stochastic method based on the Markov chain Monte Carlo, McMC, sampling approach to estimate the values of the model parameters at the pilot points. 2. Theory We present in this section the equations governing the physical processes of groundwater flow and saline tracer transport in a heterogeneous unconfined aquifer. We also introduce the semicoupled equations connecting the salt concentration to the electrical resistivity (to interpret ERI) and to the source current density used to interpret SP data. 2.1. 2D flow and transport equations In steady-state conditions, the governing groundwater flow equation in a saturated and heterogeneous porous material is given by the following elliptic partial differential equation:

r  ðK rhÞ ¼ 0;

ð1Þ

subject to the boundary conditions,

h ¼ h0

at CD

ð2Þ

^  K rh ¼ q0 at CN ; n

ð3Þ

where h denotes the hydraulic head (in m), K is the hydraulic conductivity (m s1, we assume that the aquifer is isotropic). Eqs. (2) and (3) correspond to Dirichlet and Neumann boundary conditions, respectively. The hydraulic head h0 denotes the head fixed at the boundary CD , q0 is the hydraulic flux (m2 s1) assumed to be known ^ is the unit vector normal to the at the Neumann boundary CN , and n boundary CN . The constitution of the transport equation of a salt tracer consists in the coupling of the Darcy’s law for the Darcy velocity u (in m s1), and Fick’s law for the flux of the salt jd (in kg m2 s1) [40]:

u ¼ K rh;

ð4Þ

jd ¼ qf /D  rc þ qf cu; 2

ð5Þ

1

where D (in m s ) denotes the hydrodynamic dispersion tensor, / (unitless) denotes the connected porosity, c (unitless) denotes the solute mass fraction, and qf (in kg m3) represents the solute bulk density. The transport of the salt due to the injection of a salt tracer in the aquifer follows the advection–dispersion equation derived from a combination of the continuity equation (mass conservation equation for the salt) and the generalized Fick’s law given by Eq. (5):

@ðqf /cÞ þ r  jd ¼ 0: @t

ð6Þ

64

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

Eq. (6) can be written as:

@c / þ r  ð/D  rcÞ þ urc ¼ 0; @t

V¼ ð7Þ

which is subject to the following initial and boundary conditions:

8x at t ¼ 0;

c ¼ 0;

r¼ ð9Þ

^  ðDrcÞ ¼ 0 at C2 ; n

ð10Þ

where C1 and C2 denote the upstream and downstream boundaries associated with the position of the upstream (injector) well and a second vertical plane downstream far enough from the downstream (observation) well. In these equations, c0 (unitless) denotes the solute mass fraction of the salt in the source term. In the Fickian model, the hydrodynamic dispersion tensor is given by

 D¼

Dm

a

 aL  aT þ a T v I3 þ v  v;

v

ð11Þ

where Dm is the molecular (mutual) diffusion coefficient of the salt (in m2 s1) (for a NaCl solution, Dm is 1.60  109 m2 s1 at infinite dilution and is 1.44  109 m2 s1 at high salinities at 25 °C), v ¼ jvj, a  b represents the tensorial product between vectors a and b, aL and aT are the longitudinal (along v) and transverse (normal to v) dispersivities (in m), and a denotes the tortuosity of the pore space and can be obtained as the product of the electrical formation factor F by the porosity (see below) [41]. The numerical solution of the hydrodynamic-transport problem is performed by using the package Comsol Multiphysics 4.3. In the first step, we determine the distribution of the Darcy velocity by solving the flow equation in steady-state conditions, and then we use the Darcy velocity to solve the advection–dispersion equation to determine the distribution of the salinity in space and time. 2.2. Electrical resistivity imaging ERI ERI consists in inverting apparent resistivity data obtained by injecting an electrical current with two (current) electrodes and measuring the difference in the electrical potential at a set of (voltage) electrodes (at least two electrodes). For an isotropic and heterogeneous electrical conductivity distribution in 2D can be e ðx; ky ; zÞ [42]: defined with the potential spectrum V

r  rrV~ 

2 ky

  I  rV~ ¼  s dðx  xþs Þ  d x  xs ; 2

ð12Þ

with the following boundary conditions:

~ ¼ 0 at Cd ; krk ! 1; V

ð13Þ

~ ¼ 0 at CN ; ^  ðrrVÞ n

ð14Þ

~ denotes the electrical potential (in V) due to dipole source where V modeled in 3D, r denotes the electrical conductivity (in S m1), x denotes the position vector of the source of current, Is corresponds  to the injected current (in A), xþ s and xs are the locations of the positive and negative current sources, respectively, and ky denotes the real wave number determined by the algorithm discussed by Pidlisecky and Knight [42]. The boundary conditions include a Neumann boundary condition CN imposed at the insulating air–ground interface and a Dirichlet boundary condition Cd at the other boundaries [43]. The electrical resistivity problem (Eqs. (12)–(14)) is solved numerically with the finite element approach for a set of wave ~ is transformed from the wave numbers ky [42]. The solution V number domain to the spatial domain following the approach of Dey and Morrison [43] using an inverse cosine-Fourier transform,

p

Z

1

V~ cosðky yÞdky :

ð15Þ

0

The electrical conductivity of the soil is connected to the pore water electrical conductivity according to the model of Waxman and Smits [44]:

ð8Þ

cðx; tÞ ¼ c0 at C1 ;

2

rf F

þ rs ;

ð16Þ

where F ¼ /m (unitless) denotes the electrical formation factor, rf denotes the pore water electrical conductivity (S m1), and the power-law exponent m (>1, unitless) denotes the cementation exponent, which ranges from 1.3 to 1.5 for clay-free unconsolidated sands and from 1.8 to 2.5 for consolidated sandstones [13]. The conductivity rs (in S m1) corresponds to the conductivity of the (insulating) mineral grains coated by the electrical double layer and is often called the surface conductivity. It can be related to the cation exchange capacity of the surface of the mineral grains (clays, silicates, and carbonate minerals have all an electrical double layer but generally the cation exchange capacity is dominated by clays, see [13]). We will consider that the surface conductivity is only weakly dependent on the solute concentration and we will neglect this effect below (considering such an effect would not change however the conclusions obtained in this paper). The surface conductivity will be considered to be a constant for the whole aquifer (only change in the conductivity with the solute concentration matters). The pore water conductivity is related to the solute concentration c by a simple linear equation:





rf ¼ ac bðþÞ þ bðÞ e;

ð17Þ

where c is salinity concentration of the pore water expressed in kg m3, e is the elementary charge of the electron (e = 1.6  1019 C), bðþÞ (25 °C) = 5  108 m2 s1 V1 and bðÞ (25 °C) = 7  108 m2 s1 V1 denotes the mobility of the cations and anions in the pore water, respectively, and a denotes a coefficient to convert the unit of concentration from kg m3 to molecules m3. Because ionic mobilities depend slightly on the solute concentration, the true relationship between the pore water conductivity and solute concentration is not, strictly speaking, linear. That said, a linear model is usually a good approximation, which is broadly used in the literature. From Eqs. (16)–(18), the macroscopic electrical conductivity measurements are therefore related to the salinity c by:





r ¼ ac/m e bðþÞ þ bðÞ þ rs :

ð18Þ

Therefore electrical resistivity, through its dependence on the solute concentration, provides some information that is indirectly related to the hydraulic conductivity field, which controls the migration of the solute together with the dispersivity tensor. 2.3. Self potential tomography SPT The SP method consists in measuring passively the electrical potentials at the ground surface or in boreholes using a couple of unpolarizable electrodes and a multi-channel DC-voltmeter characterized by a high input impedance (>10 MX) and a sensitivity of at least 0.1 mV. The self-potential signal recorded during a salt tracer experiment is the sum of two distinct contributions in the source current density. The first component is associated with the flow of ground water, which drags the excess electrical charge contained in the pore water. This excess of charge is more precisely located in the electrical diffuse layer coating the surface of the minerals [20]. This contribution is known as the streaming current density. The second contribution is related to the gradient of the solute itself. It depends on the gradients of the activity of the charge carriers (ions) that are present in the pore water

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

[14,39,45]. This contribution is known as the diffusion current density. In an isotropic heterogeneous porous material, the total current density j (in A m2) is the sum of a conductive current density (given by Ohm’s law) plus the two source components described above (e.g., [14]):

j ¼ r  ru þ Q V u 

 kb T  r 2tðþÞ  1 r ln rf ; e

jS ¼ Q V u 

 kb T  r 2tðþÞ  1 r ln rf ; e

above implies that the self-potential field is related to the pathways of the salt plume migration within the aquifer and therefore to the hydraulic conductivity field. It is also controlled by the distribution of the electrical conductivity and therefore a joint inversion of self-potential and resistivity data improves the resolution of the resistivity tomography because the sensitivity maps of the two methods are complementary to each other.

ð19Þ

where u denotes the electrical self-potential field (in V, E ¼ ru denotes the corresponding electrical field, in V m1), r corresponds to the electrical conductivity of the porous material (in S m1, see discussion in Section 2.2), kb is the Boltzmann constant (1.381  1023 J K1), T is the absolute temperature (°C), t(+) = 0.38 (for NaCl) denotes the microscopic Hittorf number of the cations in the pore water [34], and Q V (in C m3) denotes the advective excess charge density (due to the diffuse layer) of the pore water per unit pore volume. The advective charge density Q V can be accurately predicted from permeability (in m2) according to log10 Q V ¼ 9:2  0:82log10 k [20]. The continuity equation for the electrical charge is r  j ¼ 0. Therefore the combination of the constitutive equation and the continuity equation yields a Poisson equation for the electrical potential [46]:

r  ðrruÞ ¼ r  jS ;

65

ð20Þ ð21Þ

with the following boundary conditions:

u ¼ 0 at Cd ; krk ! 1;

ð22Þ

^  rru ¼ 0 at CN : n

ð23Þ

The Neumann boundary condition CN is imposed at the insulating air–ground interface and a Dirichlet boundary condition Cd is imposed at the other boundaries. The physical model described

3. Forward modeling in an heterogeneous aquifer A salt tracer synthetic test involves injecting a known amount of a non-reactive (conservative) salt (e.g., NaCl) into an upgradient well and monitoring the perturbation of the electrical conductivity and self-potential associated with the migration of the salt tracer due to the natural hydraulic gradient plus dispersion and diffusion. As mentioned above, the salt tracer is assumed to be conservative. If this not the case, a sorption isotherm needs to be include in the model and the forward solver becomes a reactive transport model. The monitoring is performed until the salt tracer reaches a downstream well where the salt concentration is sampled at distinct depths. The geophysical data (DC resistivity and SP) are measured with a network of electrodes located at the ground surface as well as in the two wells. We built a synthetic model study based on this scheme where sampling ports are located at eight different sampling intervals in the observation wells and the electrodes are located at 15 depths (see Figs. 1–3 for details). The distance between the injection well and the observation well is 8 m and the depth of the (vertical) wells is 15 m (Fig. 1). In our synthetic model, we use a heterogeneous hydraulic conductivity field between the two wells. The heterogeneous model for the hydraulic conductivity field is generated by the SGSIM software [47] with an anisotropic exponential variogram with a maximum range of 10 m, a sill of 0.8, an inclination of 75° from the vertical axis and the anisotropy ratio between the minor range and the major range is the 0.45. The values of the hydraulic properties were chosen to approximate a sandy alluvial environment [13]. The porosity is assumed constant all over the domain

Fig. 1. True hydraulic conductivity field generated by the SGSIM code showing two areas between the two wells characterized by high hydraulic conductivity values. The salt movement is due to a constant hydraulic head gradient applied between the upgradient (injection) well (on the left side of the first panel) and the downgradient (observation) well (on the right). The figure on the left side shows the salt concentration distribution at t = 20 days. This figure shows the formation of preferential fingering associated with preferential flow pathways (the arrows denote the Darcy velocity). The open circles corresponds to the position of the electrodes.

66

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

Fig. 2. Four snapshots of the self-potential profiles for the 15 electrodes located in the downstream (observation) well. The curves illustrate that the SP anomalies exhibit substantial anomalies (6 mV) when the plume comes in the vicinity (few meters) to the electrodes. At the opposite, the change in the self-potential signals is only of 2 mV in areas characterized by lower hydraulic conductivity values.

(/ = 0.35). Note that the porosity is expected to vary inside a factor 2–3 in real aquifers while the permeability can vary over several orders of magnitudes. The generated hydraulic conductivity field is used to solve numerically, in steady state conditions, the flow problem and to determine the distribution of the hydraulic heads as well as the Darcy velocity distribution. The domain boundaries are placed 50 m from the area of interest located between the wells to minimize boundary effects on the simulations. The hydraulic heads imposed at the boundaries on the left and right sides of the domain were selected to impose a constant hydraulic head gradient between the two wells (value of 0.02). We use a no flow boundary condition at the top and bottom boundaries of the aquifer. The solution of the ground water flow problem is used to evaluate the advection term (see Eqs. (5) and (7)) and the streaming current density (see Eqs. (20) and (21)). The advection-dispersion equation of transport for the salt was simulated by assuming the injection of 1 kg m3 of a saline (conservative) tracer in the upgradient well (see the left side of Fig. 1). This well is assumed to be entirely screened. We assume that the values of the longitudinal and transverse dispersivities in Eq. (11) are 0.5 m and 0.05 m, respectively. They are assumed to be homogeneous and taken in agreement with the scale of the investigated area (see [37] for instance for a discussion regarding the scaling of these parameters). The transient transport equation is solved over a period of 100 days. This implies in turn that the characteristic time scale of the hydrogeological problem is much longer than the characteristic time required to take a snapshot in resistivity tomography (selfpotential, being passive, is real-time). The salt transport in the aquifer has a direct impact on the distribution of the electrical conductivity of the pore water and consequently on the bulk electrical conductivity of the porous medium. A series of simulated electrical conductivity measurements were carried out with six electrodes positioned on the surface and along the walls of the two boreholes (15 electrodes per

Fig. 3. The filled circles indicate the spatial locations of the pilot points where we assigned the hydraulic conductivity values selected randomly during each iteration of the McMC sampler. These values are used to build the field perturbation by using conditional simulation with a variogram assumed to be known. The open circles denote the locations of electrodes in both boreholes to monitor the variations in time of the apparent electrical resistivity and self-potential signals. The stars indicate the position of the monitoring of the changes of salt concentration (eight in situ sampling ports).

well, see Figs. 1 and 3). There is therefore a total of 36 electrodes used for the resistivity tomography. This corresponds to a reasonable array size in performing time lapse resistivity tomography (most of the resistivity meters work with over 64 electrodes). The configuration technique used for the voltage measurements is the dipole–dipole array. Five acquisitions were performed (simulated) at the following times t = [5, 10, 15 20, 30 days] which corresponds to a regular sampling in time with respect to the residence time of the salt between the injection well and the observation well. We also simulated the acquisition of self-potential measurements to record numerically the anomalies associated with the gradients of salinity and pressure head. In real field conditions, the electrodes used previously for ERI in the downgradient well can also be used to record (passively) the electrical potential field. The self potential measurements are assumed to be collected once a day during a period of 30 days. The self-potential distribution has to be recorded before the injection of the electrical current used to measure the apparent electrical resistivities to avoid any artifacts related to electrode polarization. More often, the electrical resistivity equipment records the self-potential data as well but this record is considered as noise by the practitioners of ERI and only used to correct the apparent electrical resistivity data. The measured self-potential data are reported with respect to a unique reference location located remotely (typically 100 m) from the study area. Fig. 2 shows the self-potential signals at the positions of the electrodes located in the second well. The values of the electrical potential increase by few millivolts (a typical accuracy is 0.1 mV in time lapse SP surveys) when the plume starts to be located in the vicinity of the observation well (few meters). Areas

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

67

characterized by low values of the hydraulic conductivity exhibit no (or very small) change of the electrical potential over time because the salinity gradient is very small in these areas.

4. Stochastic joint inversion 4.1. Description of the algorithm We choose a stochastic approach to determine the hydraulic conductivity fields m from three data sources (d denotes the data vector) coming from saline concentration, electrical resistivity, and self-potential SP (see Fig. 4). The pilot point method is used as a parameterization technique to reduce the number of the model parameters by identifying a few key-locations (called pilot points) in the aquifer where the values of the hydraulic conductivity can be assessed to reproduce the hydrogeophysical measurements. Of course, this can only be done by assuming certain spatial statistical characteristic of the hydraulic conductivity field [28]. The method consists in generating from the SGSIM code [47] a hydraulic conductivity field conditioned by the pilot point values that will be modified at each iteration to minimize the misfit between measured and simulated data (Fig. 4). Then, the spatial distribution of the hydraulic conductivity field is written as [30]:

KðxÞ ¼ 10f ðxÞ K;

ð24Þ

where K is the mean value of the hydraulic conductivity (in m s1), f ðxÞ denotes a space random function with zero as mean and the variance is equivalent as of the semi-variogram of the Log10 KðxÞ field (in other words f ðxp Þ denotes a vector of the 49 hydraulic conductivity values observed at the pilot point locations). The field is obtained by a conditional simulation from the values assigned at the pilot point locations f ðxp Þ. This approach requires a variogram cðpÞ. The variogram parameters (such as the sill, the anisotropy angle, the major range, the ratio between the major and minor ranges, and the variance) can be obtained from the analysis of available measurements and observations (possibly geological datasets, hydrogeological information (transmissivity), and/or geophysical information) or determined from the inversion approach itself [30]. The vector of model parameters to be explored during the inverse process to reconstruct the hydraulic conductivity is m ¼ K; f ðxp Þ; cðpÞ where K denotes the mean of the hydraulic conductivity, the values of the pilot point f ðxP Þ and parameters of the variogram. These parameters are perturbed at each iteration during the inversion process to honor the hydrological plus geophysical data (see Fig. 3). Note that the previously described geostatistical approach is not necessarily the correct one for real applications. There are some heterogeneous aquifers that cannot be modeled with a single semi-variogram or the semi-variogram approach needs to be coupled with a facies approach (or a ‘‘rock type’’ approach in carbonate reservoirs/aquifers). Our approach could handle these difficulties but other geostatistical descriptors of the heterogeneity may be required (for instance in fractured aquifers in which dual porosity models should be used). In a probabilistic framework, the inverse problem consists in exploring the posterior probability distribution of the model parameters. These probability distributions are based on combining the information from the hydrogeophysical data with some prior knowledge on the model parameters. We note P 0 ðmÞ the prior probability density of m and PðdjmÞ represents the probability corresponding to the data for m fixed, d denotes a Nx1-vector corresponding to the observed data. PðdjmÞ provides therefore a measure of the misfit between the predicted data from the model m and the measured data d. The posteriori probability density

Fig. 4. Sketch illustrating the use of the McMC algorithm to reconstruct the hydraulic conductivity field from the inversion of the data (tracer concentrations, apparent resistivity and self-potential data). The pilot points technique is used as a parameterization method to perturb the mean hydraulic conductivity K, which is an unknown parameter together with the variogram parameters. A prior model is introduced as inequality constraints to avoid the sampling of undesirable hydraulic conductivity values. The routine provides a set of realizations sampled from the posterior probability distribution.

pðmjdÞ of the model parameters m given the data d is obtained using Bayes formula:

pðmjdÞ / PðdjmÞP0 ðmÞ:

ð25Þ

We assume that observational uncertainties and correlations across the different data sets are statistically independent, an assumption which is commonly made in Bayesian joint inversion studies. The presence of the multiple and independent data sets requires rewriting the likelihood function as a product of the likelihood of each data type [48]:

68

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77 n Y Pi ðdi jmÞ;

Pðdi jmÞ ¼

ð26Þ

i

pðmjdi Þ /

" n Y

# Pi ðdi jmÞ P 0 ðmÞ;

ð27Þ

i

where i denotes the number of data set. In our case, we have three types of data (in situ salt concentration measurements, electrical resistivity tomography ERI, and SP). For each data type i, the likelihood function Pi ðdi jmÞ is considered to be Gaussian:

  1 1 T 1 Pi ðdi jmÞ ¼ h i1=2 exp  ðg i ðmÞ  di Þ C di ðg i ðmÞ  di Þ ; 2 ð2pÞN det C di ð28Þ where gi(m) denotes the forward modeling operator to simulate the data type i (see Section 2). The forward problem connects the generation of the data vector di to a given hydraulic conductivity model m. The (N  N) diagonal covariance matrix C di accounts for the noise-to-signal ratio of the data in the inversion. In the salt tracer test, the forward problem is therefore composed of three partial forward problems, the first forward problem g 1 ðmÞ corresponds to the flow and transport model that generates the distribution of the solute concentration (see Section 2.1), the second g 2 ðmÞ allows the simulation of the electrical resistivity data (see Section 2.2), and the third one g 3 ðmÞ allows the simulation of the self-potential signals (see Section 2.3). The prior probability distribution of the model parameters, if available, is also assumed to be Gaussian:

P0 ðmÞ ¼ h

1 ð2pÞM detCm

i1=2

  1 exp  ðm  mprior ÞT C 1 ðm  m Þ ; prior m 2

ð29Þ where mprior denotes the prior value of the model parameters and C m denotes the model diagonal covariance matrix incorporating the uncertainties related of the prior model of hydraulic conductivity. To avoid unrealistic values in the sampled hydraulic conductivities, inequality constraints of the type bmin 6 SðmÞ 6 bmax need to be included in the posterior probability density pðmjdÞ. Therefore, we rewrite pðmjdÞ as a constrained posterior probability distribution [49]:

h

i

pðmjdÞ / PðdjmÞP0 ðmÞ HðSðmÞ  bmin ÞT Hðbmax  SðmÞÞ ;

ð30Þ

where,

8 SðmÞ > bmin > <1 HðSðmÞ  bmin Þ ¼ 1=2 SðmÞ ¼ bmin : > : 0 SðmÞ < bmin

ð31Þ

In these equations, H denotes the Heaviside function, KðxÞ ¼ SðmÞ is a function used to generate the hydraulic conductivity field with the parameters m (K, f ðxp Þand cðpÞ), and bmax and bmin denote the lower and upper limits of the hydraulic conductivity field, respectively. We point out that the posterior probability pðmjdÞ described by Eq. (30) is non-Gaussian. In this paper, the pilot points are placed on a regular grid in order to sample uniformly the totality of the domain comprised between the injection and observation wells. For instance, a regular grid with a density of 2–3 pilot points per correlation range in each direction is a good choice (see [50]). Their unknown hydraulic conductivity values can be assess with the mean of the hydraulic conductivity K and variogram parameters through a Markov chain Monte Carlo (McMC) sampler via the exploration of the posteriori probability density pðmjdÞ, which is expressed by Eq. (30). McMC algorithms are based on the random generation routines of a set of

models. In the present case the model parameters are again the values of the hydraulic conductivity at the (fixed) pilot point, the mean of the hydraulic conductivity field K, and parameters of the variogram (for more details about the McMC algorithm with bound constraints used in our approach see Appendix A). Then, we select all the realizations of these models after a certain transient (burn-in) period required to avoid the influence of the starting realizations. This burn-in period depends on the initial guess of the model parameters. After this period, the standard deviation on the realizations of the model parameters stabilizes and the additional realizations can be used to assess the posterior probability distribution of the model parameters. More precisely after the burn-in period, we can infer the statistical properties of the ensemble of realizations such as the central tendency (median, mode, and mean) and the percentiles of dispersion to estimate the uncertainty of the solution of the inverse problem. A more complete but tedious approach way would be to plot the posterior probably density for each model parameter. 4.2. Results with a known semi-variogram To evaluate the effectiveness of the inversion with coupling hydrogeological and geophysical data, we compared seven cases with no noise added to the data and a perfectly know semivariogram (Cases 1–7 in Table 1). We first examined the individual inversions of each type of data, and then we proceed to combine two data sets with the parameters of the variogram supposed to be perfectly and independently known (the case of an unknown variogram is discussed later). At the end, we combined all the hydrogeological and geophysical data together (the variogram parameters are assumed to be known). In Section 4.3, we are going to investigate the case where the data sets are contaminated by an uncorrelated Gaussian random noise and we will estimate the parameters of the variogram assumed to be unknown (Case 8 in Table 1). The inversions were all performed with using 107 and 101.5 m s1 as the lower and upper bounds for the hydraulic conductivity to avoid sampling non-realistic value of K that would slow down the convergence of the chain. We use a total of 49 pilot points distributed on a regular grid in agreement with the correlation length scale of the semi-variogram (see their position in Fig. 3). Indeed, the density of the pilot point (about one pilot point per m2 is much smaller than the maximum range of 10 m used in the semi-variogram). The results of the inversion are presented from the computation of the median and the percentiles 25% and 75% of the sampled simulations (see Fig. 5). The discussion of the inversion is based on a comparison between the ‘‘true’’ and inverted hydraulic conductivity fields. Their values are calculated from the median of the sampled simulations of the pilot point values f ðxp Þ with the mean of the hydraulic conductivity K using the posterior probability of the hydraulic conductivity field. We also compare the spatial distribution of the true concentration of the solute and the concentrations resulting from the most likely hydraulic conductivity field inverted at time t = 20 days. Finally, we will also show the position of two saline fronts determined respectively from the 75% and 25% quartiles of the sampled hydraulic conductivity field. In the first seven cases, the variogram parameters and porosity are assumed to be perfectly known. We first evaluate the reconstruction of the hydraulic conductivity field from the inversion of the salt concentration data alone. These data are collected at the downgradient (observation) well at eight depths (see Fig. 3) during for the time lapse comprised between t = 0 and t = 100 days. A total of 20,000 simulations of the hydraulic conductivity field were generated using the McMC sampling process. The burn-in

69

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

Table 1 Summary of results of the inversion of the hydrogeophysical data. Case 1, Case 2, and Case 3 correspond to the inversion based on the in situ concentration data alone, the electrical resistivity data alone, and the self-potential data alone. Case 4, Case 5, and Case 6 denote the inversion based on in situ concentration sampling plus resistivity data, the concentration together with the self-potential data, and the resistivity and self-potential data, respectively. Case 7 corresponds to the joint inversion of the concentration, selfpotential, and resistivity data. Finally Case 8 corresponds to the joint inversion of the noisy concentration, self-potential, and resistivity data plus the semi-variogram. RMS C, RMS SP, and RMS ERT denote the root mean square errors on the concentration, self-potential, and resistivity data using the best inverted models.

Case Case Case Case Case Case Case Case

1 2 3 4 5 6 7 8

(A1) (A2) (A3) (B1) (B2) (B3) (C)

RMS error C (%)

RMS error SP (%)

RMS error ERT (%)

Median absolute deviation (m s1)

1.2 60 1.5 0.8 1.8 2.7 0.6 4.0

7 31 6.3 13 3.5 4.5 3.2 8.0

53 8.5 41.0 5.5 31 8.0 7.6 11.2

1.8 2.4 1.7 1.3 1.45 0.85 0.51 2.02

Fig. 5. The box plots showing the results of the realizations (upon convergence of the chain, in general about 15,000 realizations) for the hydraulic conductivity values at the position of the pilot points (#1 to #49) obtained from the McMC sampling to invert the concentration saline data in the second (downstream) well. The box plots represent inter-quartile ranges and the solid black line at the centre of each box denotes the median. The arms of each box extend to cover the central 95% of the sampled posterior distribution.

period corresponds to about 5000 realizations. This burn period was found to be pretty consistent for the various model parameters and simulations. After this period, the McMC sampling stabilizes and becomes independent of the initial model. Therefore, we infer the median field from the last 15,000 realizations. As a side note, we have checked that the inversion results does not depend on the initial starting value of the model parameters (see [7]). That said, the number of realizations corresponding to the burn-in period may depend on the initial guess for the model parameters, which means that the convergence of the chain has to be carefully monitored over time. For this reason, we started our algorithm with random values of the model parameters and we monitored the convergence of the chain in all our simulations. The first inversion is called inversion A1 below (see Figs. 6a–c). This inversion yields a hydraulic conductivity field in the range 107 and 101.8 m s1 with a median absolute deviation of 1.8 m s1. The median of the realizations sampled honors perfectly the saline concentration data (see Fig. 7a). However, the spatial distribution of the salty plume at t = 20 days reconstructed from this median field of the hydraulic conductivity and two quartiles are far from the true position of the saline front (Fig. 6b). We also computed what should be the SP data predicted from this inverted field. Fig. 7b shows that the self-potential signals in the observation well remains close to the SP observed (see Table 1) but have a tendency to overestimate the true SP data. That said, this inversion of the concentration data alone did not delineate the main

permeable areas. In other words, there is a lack of similarity between the inverted hydraulic conductivity field and the true one (see Fig. 6c). Consequently, the downstream salt concentration measurements in a well do not provide enough information to recover the hydraulic conductivity field. We now focus on inverting the electrical resistivity data alone (see Fig. 3). The results of this inversion (inversion A2), are shown in Fig. 6d–f. The hydraulic conductivity estimated by using the inversion of the resistivity data are far from being perfect as shown by the difference between the true and inverted hydraulic conductivity fields (Fig. 6f). In the top surface of the aquifer, where the hydraulic conductivity takes high values, the reconstructed hydraulic conductivity field generates an erroneous distribution of the saline tracer (see Fig. 7c and d and Table 1). This discrepancy seems to be due to the lower resolution of dipole–dipole array in this area. The inversion of the self-potential measurements alone corresponds to inversion A3. This inversion yields a hydraulic conductivity field in the range [107, 102 m s1] (see Fig. 6g–i) with a median absolute deviation of 1.7 m s1. This reconstruction seems to reproduces the salinity data (Fig. 6e and Table 1) and is better than the reconstructed hydraulic conductivity field obtained from the DC resistivity data alone. This is therefore demonstrating the high information content in the self-potential field. That said, the hydraulic conductivity field inverted from the SP data alone yields an area of lower hydraulic conductivity in the vicinity of the

70

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

Fig. 6. Uncoupled inversion of each data set to estimate the hydraulic conductivity field. The inversion A1, A2 and A3 denote the inversion of concentration data, resistivity, and self potential data, respectively. (a), (d), and (g) Inverted hydraulic conductivity fields. (b), (e), and (h) True and inverted positions of the salt plume at t=20 days (the inverted position is from the median and the quartiles 75% and 25% of all realizations sampled during the inversion of the observed data). (c), (f), and (i) Difference between the real and inverted hydraulic conductivity fields.

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

71

Fig. 7. The inversions A1, A2 and A3 denote the inversion of concentration data alone, the resistivity data alone, and the self potential data alone, respectively. For each inversion A1–A3, we illustrate in the first how the predicted concentration data compare with the true data at position P[x = 7.5, z = 8 m], A similar task is done in the third column for the self-potential data collected along the downgradient well C stands for solute concentration, ERI for electrical resistivity tomography and SP for self-potential.

upgradient well between depths 2 and 6 m that does not exist in the true conductivity field (see Fig. 6i). The reconstructed hydraulic conductivity field is also still far from the true field (Fig. 6i). We present now the joint inversion results with the data sets taken two by two (Fig. 8). The result of the joint inversion of

apparent resistivity and concentration data corresponds to inversion B1. This inversion improves considerably the mapping of the hydraulic conductivity field compared to the results obtained from the decoupled inversion data of the concentration and ERI (see inversion B1 in Fig. 8a–c and compare with inversions A1 and A2

72

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

Fig. 8. Results of the joint inversion of the data sets taken two by two to estimate the hydraulic conductivity field. The inversions B1, B2, and B3 denote the inversion of the ERI + C, SP + C and ERI + SP data, respectively. (a), (d), and (g) Hydraulic conductivity fields determined from the coupled inversion. (b), (e), and (h) True positions of the salt plume at t = 20 days and inverted position from the median and the quartiles 75% and 25% of all realizations sampled during the inversion. The third column represents, for each inversion, the difference between the true and inverted hydraulic conductivity fields.

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

73

Fig. 9. The inversion B1, B2 and B3 indicate the coupled inversions (ERI + C, SP + C and ERI + SP respectively). For each inversion, we represent in the first column the predicted and true concentration data at P[x = 7.5, z = 8 m] and in the second column the predicted and true SP data along the downgradient borehole at t = 20 days.

in Fig. 6a–f). In this joint inversion, we eliminated the occurrence of the areas of high conductivity hydraulic near the surface and resulting from the inversion of the apparent resistivity data alone

(Fig. 6d–f) and the all hydrogeophysical data were fairly well reproduced (see Fig. 9a and b) (also see RMS errors are reported in Table 1 for each type of data). However, the spatial distribution

74

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

of the salt front at t = 20 days reveals that the estimated hydraulic conductivity does not reproduce perfectly the salt motion in the aquifer (see Fig. 8b). The joint inversion of concentration data and self-potential (inversion B2) delineates pretty well the area of preferential flows located in the center portion of the aquifer (see Fig. 8d–f). This improvement is also observed in the joint inversion of the resistivity and self-potential data (see inversion B3 in Fig. 8g–i). In this inversion, both the spatial distribution of the hydraulic conductivity and the position of the saline front have been improved by the

joint inversion of the resistivity and self-potential data. The predicted salt concentration and self-potential data are also better fitted through the joint inversion procedure by comparison with the uncoupled inversions (see Fig. 9e and f). The joint inversion of the self-potential, resistivity, and in situ salinity sampling corresponds to inversion C (Fig. 10a–c). The joint inversion of all the data sets together is shown in Fig. 11. This inversion improves the reconstruction of the hydraulic conductivity field. It is characterized by a low of the median absolute deviation (0.51 m s1) due to the complementary

Fig. 10. Result of the fully coupled inversion of all the data together. (a) The hydraulic conductivity field estimated from the joint inversion of all data sets (concentration C, ERI and SP). (b) Spatial distribution of the saline plume at t = 20 days from the true conductivity field compared to the plume predicted from the inverted hydraulic conductivity field (with its confidence interval) using all the data (salinity, resistivity, and self-potential). (c) Difference between the true and inverted hydraulic conductivity fields.

Fig. 11. The predicted data obtained from the joint inversion of all data sets. (a) Predicted and true concentration at position P[x = 7.5, z = 8]. The predicted data is obtained from the joint inversion of the saline tracer test plus the geophysical information (electrical resistivity imaging and self-potential data). (b) Predicted and true self-potential profiles along the downgradient borehole at t = 20 days.

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

75

Fig. 12. Result of the fully coupled inversion of all the data together which are contaminated with uncorrelated Gaussian random noise having a standard deviation of 15% of the mean value of the measurements and the unknown parameters of the variogram. (a) The hydraulic conductivity field estimated from the joint inversion of all data sets (concentration C, ERI and SP). (b) Spatial distribution of the saline plume at t = 20 days from the true conductivity field compared to the plume predicted from the inverted hydraulic conductivity field (with its confidence interval) using all the data (salinity, resistivity, and self-potential). (c) Difference between the true and inverted hydraulic conductivity fields.

Fig. 13. The predicted data obtained from the joint inversion of all data sets with an unknown variogram and noise added to all the data. (a) Predicted and true concentrations at position P[x = 7.5 m, z = 8 m]. (b) Predicted and true self-potential profiles along the downgradient borehole at t = 20 days.

sensitivity of the methods to salinity and hydraulic head gradients. This fully-coupled joint inversion is clearly the best strategy to increase the resolution between the two wells in order to provide the most meaningful hydraulic conductivity distribution.

4.3. Inversion with an unknown semi-variogram In this section, we discuss the case of the joint inversion of all the hydrogeophysical data when they are contaminated by an uncorrelated Gaussian random noise. This noise is characterized

76

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

by having a standard deviation of 15% of the mean value of the measurements. Such a level of noise is realistic with respect to real conditions. In this case, we also assume no prior information regarding the variogram. In other words, we need to estimate the parameters of the variogram in the inversion. The noise associated with the measurements of the self-potential data may come to daily variations of temperature in the shallow subsurface [20]. In practice, self-potential measurements may be accompanied by the measuring changes of temperature to remove this spurious effect. Resistivity data may be also contaminated by noise especially with the dipole–dipole array as this array may be characterized by a mediocre signal-to-noise ratio. In real field conditions, the variogram can be estimated from available geological and hydrogeophysical information (e.g., downhole measurements). In our case, it can be deduced from the first DC resistivity tomogram taken before the injection of the salt injection. The spatial distribution of the electrical resistivity reflects the main geological electrofacies providing an idea of the variogram to use in the inversion. It provides also an idea of the number of pilot points to use. In this inversion, we performed 20,000 realizations with random starting values for the hydraulic conductivity at the pilot points. We started the inversion with an isotropic variogram with a sill of 1 and a range of the 6 m (here again we tested several choices and observed that the final result is independent on these choices). We emphasize that during this inversion, the constraints for the hydraulic conductivity were also maintained like in previous cases. Among the models sampled via posterior probability density, we used the last 14,000 simulations to determine the median value of the sampled hydraulic conductivities (K varies between 107 and 101.5 m s1). The result of the inversion is shown in Figs. 12 and 13. The reconstruction provides the main hydrofacies controlling the salt motion especially the permeable pathways in the center and at the bottom of the domain (Fig. 12). The inverted (anisotropic) variogram is characterized by a maximum range of 12 m, a sill of 0.5, an inclination of 70° from the vertical axis and an anisotropy ratio between the minor range and the major range of 0.30. This inversion with an unknown variogram is pretty good despite the fact that it corresponds to a larger uncertainty in the inverted model with respect to the other investigated cases (see Table 1), as expected. 5. Concluding statements We have presented a synthetic case study investigating the usefulness of adding self potential measurements to electrical resistivity imaging and in situ concentration measurements to reconstruct the hydraulic conductivity field of an heterogeneous aquifer during a salt tracer test. As most of the impedance meters used for ERI are also able to record the self-potential field before the injection of electric current, this approach does not require new equipment with respect to those used traditionally in hydrogeophysics for this type of problem. Despite the fact that self-potential is considered to be a source of noise in ERI studies, we show that it contains valuable information because of its influence to the solute concentration and Darcy velocity. In our study, we showed that the self-potential method could be therefore a complementary tool to investigate salt tracer tests and to reconstruct the hydraulic conductivity field of heterogeneous aquifers. In our synthetic experiment, the salt tracer was tracked directly by the sampling of the salinity concentration in the down-gradient well. The tracking of the salt was also carried out indirectly thanks to a combination of DC resistivity and self-potential measurements. The movement of the salt plume linearly perturbs the background electrical conductivity. Consequently, the ERI imaging permits to draw the limit of the salt front at different times. The hydraulic

and salinity gradients in the ground are both responsible for self potential signals with a similar magnitude that is easily measurable (>0.1 mV) both at the ground surface and in boreholes. Finally, we have applied a stochastic approach to jointly invert different data sets to reconstruct the hydraulic conductivity field. To reduce the numbers of the unknown parameters to map the hydraulic conductivity and the time of computation, we used the pilot point method as a technique of parameterization. To improve the efficiency of this approach; we choose to introduce constraints on the hydraulic conductivity values to avoid the sampling of unrealistic values of the hydraulic conductivity. In this paper, we focused on the reconstruction of the hydraulic conductivity, but we feel that the approach can be extended to estimate spatial distribution of the porosity from the hydraulic conductivity using a simple linear relationship is established through the regression equation between the core hydraulic conductivity and the core porosity [51]: K ¼ a log10 / þ b. The two coefficients a and b can be estimated also from the inverse process or from laboratory experiments realized on core samples assuming no sampling bias in getting the samples. Revil and Jardani [46] showed that the mean value of the dispersivity can be estimated from the self-potential data. Similarly, Fowler and Moysey [52] showed that the time-lapse electrical resistivity can be used to estimate the mean value of the dispersivity. These approaches can also be applied to the estimations of the transports parameters of heterogeneous systems. The McMC algorithm used above is demanding in terms of computing time. To overcome this issue, the parallelization of our code for 3D applications will be required, An alternative road is to apply deterministic techniques based on linearization of an objective function through the state adjoint operators. A future work will concern the inclusion of the induced polarization (IP) method because of its sensitivity to permeability (see [53–55]) and the possibility to record ERI, time-domain IP, and SP with the same equipment. Acknowledgements We thanks Seine Aval program for funding the project ‘‘TidHydrex’’. AR thanks the Environment Remediation Science Program (ERSP), US Department of Energy (DOE, award DE-FG0208ER646559) for funding. We deeply thank Jinsong Chen and three anonymous referees for very constructive comments and the time they have invested in helping us to improve our manuscript.

Appendix A A description of the Metropolis–Hastings algorithm with inequality constraint is given below. This algorithm is based on the following steps: Step 1: Specify an arbitrary starting value m0 which satisfies the constraints. Set i = 0. Step 2: Given the current value mi, use a symmetric transition density qðmi jmc Þ to generate a candidate for the next value in the sequence, mc. Step 3: Use the candidate value mc to evaluate the inequality conh i straints HðSðmÞ  bmin ÞT Hðbmax  SðmÞÞ . If any constraints is not respected, set aðmi ; mc Þ ¼ 0 and go to Step 5. h i c Þqðmc jmi Þ . Step 4: Calculate aðmi ; mc Þ ¼ min 1; ppðm ðmi Þqðmi jmc Þ Step 5: Generate an independent uniform random variable U from the interval [0, 1].

mc ; if U < aðmi ; mc Þ Step 6: Set miþ1 ¼ : mi ; if U P aðmi ; mc Þ Step 7: Set i = i + 1 and go to Step 2.

A. Jardani et al. / Advances in Water Resources 52 (2013) 62–77

References [1] Gelhar LW. Stochastic subsurface hydrology. New York: Prentice Hall; 1993. [2] Berg SJ, Illman WA. Three-dimensional transient hydraulic tomography in a highly heterogeneous glaciofluvial aquifer–aquitard system. Water Resour Res 2011;47:W10507. http://dx.doi.org/10.1029/2011WR010616. [3] Brauchler R, Doetsch J, Dietrich P, Sauter M. Derivation of site-specific relationships between hydraulic parameters and P-wave velocities based on hydraulic and seismic tomography. Water Resour Res 2012;48:W03531. http://dx.doi.org/10.1029/2011WR0108. [4] Cardiff M, Barrash W, Kitanidis PK. A field proof-of-concept of aquifer imaging using 3D transient hydraulic tomography with modular, temporarilyemplaced equipment. Water Resour Res in press;48:W05531. http:// dx.doi.org/10.1029/2011WR011704. [5] Butler J. Hydrogeological methods for estimation of spatial variations in hydraulic conductivity. In: Rubin Y, Hubbard S editors. Chapter 2 in Hydrogeophysics. Springer; 2004. [6] Li W, Englert A, Cirpka OA, Vanderborght J, Vereecken H. 2D characterization of hydraulic heterogeneity by multiple pumping tests. Water Resour Res 2007;43:W04433. http://dx.doi.org/10.1029/ 2006WR005333. [7] Chen JS, Hubbard SS, Gaines D, Korneev V, Baker G, Watson D. Stochastic estimation of aquifer geometry using seismic refraction data with borehole depth constraints. Water Resour Res 2010;46:W11539. http://dx.doi.org/ 10.1029/2009WR008715. [8] Daily W, Ramirez A, Labrecque D, Nitao J. Electrical-resistivity tomography of vadose water movement. Water Resour Res 1992;28(5):1429–42. [9] Kim J, Teixeira F. Parallel and explicit finite-element time-domain method for Maxwell’s equations. IEEE Trans Antennas Propag 2011;59(6): 2350–6. [10] Ha T, Pyun S, Shin C. Efficient electric resistivity inversion using adjoint state of mixed finite-element method for Poisson’s equation. J Comput Phys 2006; 214(1):171–86. [11] Knight R, Pyrak-Nolte LJ, Slater L, Atekwana E, Endres A, Geller J, et al. Geophysics at the interface: response of geophysical properties to solid–fluid, fluid–fluid, and solid–solid interfaces. Rev Geophys 2010;48:RG4002. http:// dx.doi.org/10.1029/2007RG00024. [12] De Barros L, Dietrich M, Valette B. Full waveform inversion of seismic waves reflected in a stratified porous medium. Geophys J Internat 2010;182: 1543–56. [13] Revil A, Cathles LM, Losh S, Nunn JA. Electrical conductivity in shaly sands with geophysical applications. J Geophys Res 1998;103:23’925–quo;936. [14] Revil A, Linde N. Chemico-electromechanical coupling in microporous media. J Colloid Interface Sci 2006;302:682–94. [15] Rubin Y, Hubbard S. Hydrogeophysics: water and science technology library. The Netherlands: Springer; 2005. p. 50. [16] Bowling JC, Zheng CM, Rodriguez AB, Harry DL. Geophysical constraints on contaminant transport modeling in a heterogeneous fluvial aquifer. J Contam Hydrol 2006;85:72–88. [17] Pollock D, Cirpka OA. Fully coupled hydrogeophysical inversion of synthetic salt tracer experiments. Water Resour Res 2010;46:W07501. http://dx.doi.org/ 10.1029/2009WR008575. [18] Binley A, Winship P, West LJ, Pokar M, Middleton R. Seasonal variation of moisture content in unsaturated sandstone inferred from borehole radar and resistivity profiles. J Hydrol 2002;267:160–72. http://dx.doi.org/10.1016/ S0022-1694(02)00147-6. [19] Linde N, Binley A, Tryggvason A, Pedersen LB, Revil A. Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data. Water Resour Res 2006;42:W12404. http://dx.doi.org/10.1029/2006WR005131. [20] Jardani A, Revil A, Barrash W, Crespy A, Rizzo E, Straface S, et al. Reconstruction of the water table from self-potential data: a bayesian approach. Ground Water 2009;47(2):213–27. [21] Suski B, Revil A, Titov K, Konosavsky P, Dagès C, Voltz M, et al. Monitoring of an infiltration experiment using the self-potential method. Water Resour Res 2006;42:W08418. http://dx.doi.org/10.1029/2005WR004840. [22] Rizzo E, Suski B, Revil A, Straface S, Troisi S. Self-potential signals associated with pumping-tests experiments. J Geophys Res 2004;109:B10203. http:// dx.doi.org/10.1029/2004JB003049. [23] Straface S, Fallico C, Troisi S, Rizzo E, Revil A. An inverse procedure to estimate transmissivities from heads and SP signals. Ground Water 2007;45(4): 420–8. [24] Straface S, Chidichimo F, Rizzo E, Riva M, Barrash W, Revil A, et al. Joint inversion of steady-state hydrologic and self-potential data for 3D hydraulic conductivity distribution at the Boise Hydrogeophysical Research Site. J Hydrol 2011;407:115–28. [25] Naudet V, Revil A, Rizzo E, Bottero JY, Bégassat P. Groundwater redox conditions and conductivity in a contaminant plume from geoelectrical investigations. Hydrol Earth Syst Sci 2004;8(1):8–22. [26] Martínez-Pagán P, Jardani A, Revil A, Haas A. Self-potential monitoring of a salt plume: a sandbox experiment. Geophysics 2010;75(4):WA17–25. [27] Linde N, Finsterle S, Hubbard SS. Inversion of tracer test data using tomographic constraints. Water Resour Res 2006;42:W04410. http:// dx.doi.org/10.1029/2004WR003806.

77

[28] Hyndman DW, Harris JM, Gorelick SM. Inferring the relation between seismic slowness and hydraulic conductivity in heterogeneous aquifers. Water Resour Res 2000;36(8):2121–32. [29] Hördt A, Blaschek R, Kemna A, Zisser N. Hydraulic conductivity estimation from induced polarisation data at the field scale – the Krauthausen case history. J Appl Geophys 2006;62:33–46. [30] Kowalsky MB, Finsterle S, Peterson J, Hubbard S, Rubin Y, Majer E, et al. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resour Res 2005;41(11): W11425. [31] Hinnell AC, Ferré TPA, Vrugt JA, Huisman JA, Moysey S, Rings J, et al. Improved extraction of hydrologic information from geophysical data through coupled hydrogeophysical inversion. Water Resour Res 2010;46:W00D40. http:// dx.doi.org/10.1029/2008WR007060. [32] Binley A, Henry-Poulter S, Shaw B. Examination of solute transport in an undisturbed soil column using electrical resistance tomography. Water Resour Res 1996;32(4):763–9. [33] Binley A, Cassiani G, Middleton R, Winship P. Vadose zone flow model parameterization using cross-borehole radar and resistivity imaging. J Hydrol 2002;267(3–4):147–59. [34] Slater L, Binley A, Versteeg R, Cassiani G, Birken R, Sandberg S. A 3D ERT study of solute transport in a large experimental tank. J Appl Geophys 2002;49(4): 211–29. [35] Kemna A, Vanderborght J, Kulessa B, Vereecken H. Imaging and characterization of subsurface solute transport using electrical resistivity tomography (ERT) and equivalent transport models. J Hydrol 2002;267: 125–46. [36] Vanderborght J, Kemna A, Hardelauf H, Vereecken H. Potential of electrical resistivity tomography to infer aquifer transport characteristics from tracer studies: a synthetic case study. Water Resour Res 2005;41:W06013. [37] Singha K, Gorelick SM. Saline tracer visualized with electrical resistivity tomography: field scale moment analysis. Water Resour Res 2005;41:W05023. [38] Irving J, Singha K. Stochastic inversion of tracer test and electrical geophysical data to estimate hydraulic conductivities. Water Resour Res 2010;46:W11514. http://dx.doi.org/10.1029/2009WR008340. [39] Pollock D, Cirpka OA. Fully coupled hydrogeophysical inversion of a laboratory salt tracer experiment monitored by electrical resistivity tomography. Water Resour Res 2012;48:W01505. [40] Oltean C, Buès MA. Infiltration of salt solute in homogeneous and saturated porous media – an analytical solution evaluated by numerical simulations. Transp Porous Media 2002;48:61–78. [41] Revil A. Ionic diffusivity, electrical conductivity, membrane and thermoelectric potentials in colloids and granular porous media: a unified model. J Colloid Interface Sci 1999;212:503–22. 10.1006 jcis.1998.6077. [42] Pidlisecky A, Knight R. FW2_5D: a MATLAB 2.5-D electrical resistivity modeling code. Comput Geosci 2008;34:1645–54. [43] Dey A, Morrison HF. Resistivity modelling for arbitrarily shaped twodimensional structures. Geophys Prospect 1979;27:106–36. [44] Waxman MH, Smits LJM. Electrical conductivities in oilbearing shaly sands. SPEJ Soc Pet Eng J 1968;8:107–22. [45] Maineult A, Bernabé Y, Ackerer P. Detection of advected concentration and pH fronts from spontaneous potential measurements. J Geophys Res 2005;110:B11205. http://dx.doi.org/10.1029/2005JB003824. [46] Revil A, Jardani A. Stochastic inversion of permeability and dispersivities from time lapse self-potential measurements: a controlled sandbox study. Geophys Res Lett 2010;37:L11404. http://dx.doi.org/10.1029/2010GL043257. [47] Deutsch CV, Journel AG. GSLIB: geostatistical software library and user’s guide. New York: Oxford University Press; 1992. [48] Aines R, et al. The stochastic engine initiative: improving prediction of behavior in geologic environments we cannot directly observe. Rep. UCRL-ID148221. Livermore, CA: Lawrence Livermore Natl. Lab.; 2002. 58 pp. [49] Michalak AM. A Gibbs sampler for inequality-constrained geostatistical interpolation and inverse modeling. Water Resour Res 2008;44:W09437. http://dx.doi.org/10.1029/2007WR006645. [50] Capilla JE, Gomez-Hernandez JJ, Sahuquillo A. Stochastic simulation of transmissivity fields conditional to both transmissivity and piezometric data. 2. Demonstration on a synthetic aquifer. J Hydrol 1997;203:175–88. [51] Soltani A, Le Ravalec-Dupin M, Fourar M, Rosenberg E. Three-dimensional characterization of permeability at the core scale. Transp Porous Media 2010;84(2):285–305. [52] Fowler D, Moysey S. Estimation of aquifer transport parameters from resistivity monitoring data within a coupled inversion framework. J Hydrol 2011;409:545–54. [53] Revil A, Florsch N. Determination of permeability from spectral induced polarization data in granular media. Geophys J Int 2010;181:1480–98. http:// dx.doi.org/10.1111/j.1365-246X.2010.04573.x. [54] Revil A. Spectral induced polarization of shaly sands: influence of the electrical double layer. Water Resour Res 2012;48:W02517. http://dx.doi.org/10.1029/ 2011WR011260. [55] Revil A, Koch K, Holliger K. Is it the grain size or the characteristic pore size that controls the induced polarization relaxation time of clean sands and sandstones? Water Resour Res 2012;48:W05602. http://dx.doi.org/10.1029/ 2011WR011561.