Assessment of hydraulic conductivity distributions through assimilation of travel time data from ERT-monitored tracer tests

Advances in Water Resources 84 (2015) 23–36

Contents lists available at ScienceDirect

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

Assessment of hydraulic conductivity distributions through assimilation of travel time data from ERT-monitored tracer tests E. Crestani∗, M. Camporese, P. Salandin Department of Civil, Environmental, and Architectural Engineering, University of Padova, via Marzolo 9, 35131 Padova, Italy

a r t i c l e

i n f o

Article history: Received 17 February 2015 Revised 23 July 2015 Accepted 30 July 2015 Available online 6 August 2015 Keywords: Hydraulic conductivity Data assimilation Travel time Electrical resistivity tomography Ensemble Kalman filter

a b s t r a c t Assessing the spatial distribution of hydraulic conductivity (K) in natural aquifers is fundamental to predict the spatio-temporal evolution of solutes, a process that is mainly controlled by the heterogeneity of K. In sedimentary aquifers, the vertical variations of K are typically more relevant than the horizontal ones in controlling the plume evolution at the local scale; such K vertical distributions can be inferred by combining the Lagrangian formulation of transport with the assimilation of tracer test data via the ensemble Kalman filter (EnKF). In this work, the data for the assimilation procedure are provided by monitoring tracer tests with electrical resistivity tomography (ERT). Our main objective is to show the possibility of directly using ERT data by assimilating the solute travel times, instead of the concentration values, thus avoiding the need for a petrophysical law. The methodology is applied to both a synthetic and a real test case and gives a satisfactory retrieval of the K field distribution, as well as of the solute evolution. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction A detailed knowledge of the hydraulic properties of natural aquifers is fundamental for the prediction of groundwater flow and contaminant transport. However, in natural porous media such properties are usually characterized by a large spatial variability, making their assessment difficult to achieve. Many inverse models have been developed to estimate aquifer hydraulic properties (e.g., [7,20,37,55,58,59]) and, in recent years, advances in computational methods and resources have opened new possibilities in their use. Traditionally, hydraulic conductivity and head measurements have been used to provide the necessary data in the application of inverse models (e.g., [8,22,29,34,56]), but in the last years data derived by geophysical methods have been increasingly used for constraining groundwater inverse problems (e.g., [5,25,27,44,45,51,54]). Ground penetrating radar (GPR) (e.g., [1]) and electrical resistivity tomography (ERT) (e.g., [3,30]) are now popular tools for hydrological and environmental research and applications, due to their capability to monitor changes in soil electrical properties, which in turn are linked to variations of moisture content or solute concentration. As a result, several joint inversion approaches have been recently developed in order to reconstruct the hydraulic conductivity distribution using time-lapse data deduced from geophysical

∗

Corresponding author. Tel.: +39 049 827 5438. E-mail address: [email protected] (E. Crestani).

http://dx.doi.org/10.1016/j.advwatres.2015.07.022 0309-1708/© 2015 Elsevier Ltd. All rights reserved.

measurements. The main advantage of such approaches is that the numerical models for the geophysical and hydrological processes are linked together so that the geophysical data are inverted directly to obtain the hydrological properties of interest. Recent applications include one-dimensional or two-dimensional infiltration experiments in the vadose zone monitored by GPR (e.g., [17,26,36]) or ERT (e.g., [23,24,33,47]) and ERT monitoring of tracer test experiments in shallow aquifers (e.g., [25,44]). In order to retrieve the hydraulic conductivity distribution, ERT may also be used to assess solute temporal moments [43] and travel times [40]. Among the many mathematical tools available for solving inverse problems, Bayesian data assimilation methods such as the ensemble Kalman filter (EnKF) [14] have been increasingly used by groundwater modelers. One of the first uses of the EnKF to retrieve K distributions can be found in [8], who estimate the hydraulic conductivity in both two- and three-dimensional synthetic domains by assimilating hydraulic head and K data. In [35], the MADE site K distribution is obtained by assimilating, in two subsequent steps, hydraulic head and concentration data resulting from a tracer test and using a steady-state flow model. Hendricks Franssen et al. [21] apply the EnKF to calibrate the conductivity and leakage coefficient together in real-time in an unconfined aquifer, while in a more recent work [32] piezometric heads and groundwater temperatures are jointly assimilated to estimate hydraulic subsurface parameters. In [34], the EnKF is used to map the hydraulic conductivity and porosity fields in a two-dimensional domain by assimilating dynamic piezometric data and multiple concentration measurements. Tong et al. [53] use the EnKF in a synthetic two-dimensional domain to

24

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

retrieve the hydraulic conductivity by assimilating solute concentrations measured in observation wells. In two recent studies, by assimilating ERT images into a model of Lagrangian transport through EnKF, Camporese et al. [5,6] investigate issues affecting the hydraulic conductivity retrieval performance related to the uncertainty of the electrical data and the choice of the prior information, as well as the filter divergence (or inbreeding) [16,22]. In addition, several variations of the original EnKF have been developed to improve its parameter updating capabilities, such as the dual EnKF (e.g., [19,39]), to deal with non-Gaussian distribution of the variables (e.g., [57]), and to limit or prevent the filter inbreeding (e.g., [52,56]). In this paper we present a method to retrieve the vertical hydraulic conductivity distribution by using ERT data obtained by a geophysical inversion without estimating concentration values from resistivities, which would require the on-site calibration of a petrophysical relationship (e.g., Archie’s law). The proposed method is applied at the local scale, i.e., in a volume whose horizontal and vertical extents are of the same order of magnitude. Using the same Lagrangian transport modeling framework described by Camporese et al. [5] and Crestani et al. [9], solute travel time data deduced by ERT are used as measurements in the EnKF assimilation procedure to update the aquifer K field. The model is first tested in a synthetic test case and then applied in a real-world study in Northern Italy.

2. Material and methods 2.1. Transport model and travel time approach The Lagrangian approach [12] is adopted to describe the movement of a tracer plume through a saturated, spatially heterogeneous porous medium. As described in [5] and [9], the solute cloud is driven by the effective velocity field obtained at steady state by solving the groundwater flow equation with appropriate boundary conditions and assuming that in natural sedimentary aquifers the spatial variability of the porosity can be considered negligible with respect to that of the hydraulic conductivity K (e.g., [18]). The steady-state assumption is justified by the fact that the method is intended to be applied for controlled tracer tests, typically characterized by a limited duration and carried out in the absence of significant transient forcings such as rainfall or variations of the boundary conditions. Furthermore, we assume the tracer as non-reactive (e.g., sodium chloride), as commonly adopted in tracer experiments. The solute injection of a non-reactive tracer is described as a particle system. Each particle has mass m = C0 (a)a and moves along the trajectory x = Xt (t; a, t0 ), where a is the starting position of the particle, C0 (a) its initial concentration, a is the infinitesimal volume associated to a and t0 is the injection time, hereinafter assumed equal to zero. The flow field statistics depend on the spatial distribution of K, which controls the trajectory via the classical Eulerian–Lagrangian relationship [11]. We use the travel time approach to describe the evolution of the solute as related to a control plane CP, fixed at a distance x¯ from the injection and perpendicular to the mean flow direction (Fig. 1). From the analysis of the trajectories, it is possible to evaluate the position η(η, ζ ), among all the possible positions ψ = (y, z), where the particles cross the CP and the travel time τ needed by each solute particle to reach the CP from the injection point. The cumulative solute mass can then be determined, as well as the temporal moments of the concentrations (e.g., [10,13]). Defined with δ the Dirac distribution, the solute flux mass through the CP at time t for a particle which was at x = a at t = 0 is

¯ a) = C0 (a) a δ(ψ − η) δ(t − τ ) q(ψ , t; x,

Fig. 1. Sketch of the Lagrangian approach for modeling of tracer tests: the tracer injection is simulated by a particle swarm released in the injection volume V0 . Travel time data are collected across the control plane, which is divided in sub-control planes (sCPs).

M that has crossed the CP until the time t is





¯ V0 ) = M(t; x,

∞ 0



dτ



×

∞

−∞

dη



∞

−∞

dζ





t

da V0

dt 

0

¯ a), (2) dψ C0 (a) δ(ψ − η) δ(t  − τ ) g1 (τ , η; x,

A

where V0 is the volume injection and A the CP area. The joint pdf g1 can be assumed as the product of two independent pdfs [13], i.e., ¯ a) = g2 (τ ; x, ¯ a) g3 (η; x, ¯ a), the former relative to the probg1 (τ , η; x, ability of a particle that starts from a to cross the CP at time τ and the latter relative to the probability of a particle that starts from a to cross the CP in η(η, ζ ). Therefore, Eq. (2) leads to:

 



¯ V0 ) = M(t; x,



t

dt 

da

V0



A

0

 =



t

da V0

0

¯ a) g3 (ψ ; x, ¯ a) (3a) dψ C0 (a) g2 (t  ; x,

¯ a), dt  C0 (a) g2 (t  ; x,

(3b)

which allows us to determine the temporal behavior of the solute mass independently from the particle crossing position on the CP (pdf g3 ), if all the possible positions ψ(y, z) are included in A. From the cumulative distribution function (CDF) given by Eq. (3b), the mean K value and the variance of the formations crossed by the solute can be obtained [12]. We describe now the CP as the sum of a finite number nsCP of sub control planes (sCPs) that do not overlap each other, such that nsCP A = A, where Ai is the area of each sCP. Eq. (3a) can then be i=1 i rearranged as follows:





¯ V0 ) = M(t; x, =



 da

V0

0

nsCP   V0

i=1

where

 Ai

t

dt 

nsCP   i=1

da C0 (a)

Ai



t 0

¯ a) g3 (ψ ; x, ¯ a) dψ C0 (a) g2 (t  ; x, ¯ a) dt  g2 (t  ; x,

 Ai

¯ a) dψ g3 (ψ ; x, (4)

¯ a) is the probability that the solute mass dψ g3 (ψ ; x,

injected in a crosses the ith sCP. In a statistically homogeneous Lagrangian flow field, assuming that the injection volume reduces to a line or a plane perpendicular to the mean flow direction x, we can write





¯ V0 ) = M(t; x,

nsCP  

V0

i=1



×

(1)

and is characterized by a joint probability distribution function (pdf) ¯ a). By definition, the expected value of the total solute mass g1 (τ , η; x,



=

Ai



t 0

dt  g2 (t  ; x¯ − ax , b)

dψ g3 (ψ − b; x¯ − ax )

nsCP   i=1

db C0 (ax , b)

V0

db C0 (ax , b)



t 0

dt  g2 (t  ; x¯ − ax , b) Wi

(5)

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

where b is the vector of coordinates (ay , az ) that, along with the constant value ax , define the position a(ax , ay , az ) of the solute particle at t = 0. The ith sCP is crossed by only a fraction of all the par ticles and the term Wi = A dψ g3 (ψ − b; x¯ − ax ) can be computed i on the basis of the knowledge of the pdf g3 . For example, if the transverse displacement of the trajectories is negligible with respect to the dimensions of each sCP, it can be assumed that g3 (ψ − b)  δ(ψ −b), where δ is again the Dirac function. This implies that Wi = nsCP AV0,i / i=1 AV0,i where AV0,i is the projection of the injection volume V0 on the ith sCP, according to the mean flow direction. A cumulative distribution functions is then defined for each sCP. This distribution retains the statistical characteristics of the velocity field controlling the transport of the solute particles from the starting position b to the ith sCP and, as a consequence, of the hydraulic conductivities crossed along the path. 2.2. The ensemble Kalman filter The ensemble Kalman filter (EnKF) [14] is an extension of the traditional recursive Kalman filter algorithm [28], which improves the estimation of transient model variables (system state) by assimilating measurements of a physical quantity related to the modeled system. The EnKF is based on a Monte Carlo approach that approximates the probabilistic information conveyed by the conditional pdf of the system state given the measurements [15]. The ensemble mean of the Monte Carlo realizations represents the estimation of the considered variables in terms of conditional expected value, while the second moment gives a measure of the estimation uncertainty. Compared to the traditional Kalman filter, the main advantage of the EnKF is that it can deal with nonlinear models and measurement operators, thanks to the approximation of the probability distributions of the variables with an ensemble of realizations [16]. In the basic EnKF formulation, a spatio-temporal discretization of the mathematical model describing a physical process is used. We define as Nr the number of Monte Carlo realizations and y j (t ) the system state vector that collects the model variables for the jth realization ( j = 1, . . . , Nr) at time t. Each y j (t ) is propagated in time by the model, according to a vector-valued discrete-time state equation:

y j (t ) = M[y j (t − ), β , t]; t0 ≤ t − < t; y j (t0 ) = y0j , j

(6)

where M is the function relating the system state at the current time t to the system state at the previous time t − ,β j is the set of model j parameters, and y0 is the initial condition at time t0 . The model describing how the measurements and system state are related is also expressed as a vector-valued discrete-time equation:

z(ti ) = Hytrue (ti )

(7)

where z(ti ) is the vector of observations at the time t = ti of the ith assimilation step, H is the operator that maps the model state to the measurements and ytrue (ti ) is the (unknown) true system state. The measurements are affected by errors and Eq. (7) can be re-written as z j (ti ) = Hytrue (ti ) w j (ti ), w j (ti ) being a log-normally distributed noise with mean equal to one and variance proportional to the measurement value, assigned in terms of coefficient of variation. When a set of measurements is available at time ti , it is possible to update the system state according to the Kalman filter equation [16]: j yupd

(ti ) = y (ti ) + Pe H (HPe H + Re ) (z (ti ) − Hy (ti )), j

T

T

−1

j

j

(8)

where yupd (ti ) is the updated system state at the assimilation time, j

Pe is the prior estimate of the system state error covariance matrix, which is computed by sampling the ensemble statistics, and Re is the measurement error covariance matrix.

25

The EnKF is used in this study as an inverse modeling tool to estimate the model parameters, by including them in system state vector and updating them on the basis of the cross-correlations with the measurements [16,41]. In particular, we use a “restart” EnKF, i.e., after each update the new parameters are used to simulate the tracer evolution from the beginning (t = 0) to the next assimilation time. Therefore, the system state does not need to be updated, being always consistent with the parameters, and techniques such as the dual EnKF are not necessary. Furthermore, the recursive nature of the restart scheme helps to overcome possible issues related to non-Gaussianity of the variables, as demonstrated previously by Crestani et al. [9]. The K spatial distribution is estimated by assimilating travel time data obtained as described in Sections 2.3 and 2.4. 2.3. ERT data The electrical resistivity tomography (ERT) is a geophysical method that provides the subsurface electrical conductivity distribution from a number of current and voltage measurements made from a series of electrodes located on the soil surface or in boreholes (e.g., [30,38]). ERT is widely used to monitor saline tracer tests, thanks to its capability of mapping the solute evolution, due to the contrast of electrical conductivity between fresh groundwater and the tracer. Cross-borehole ERT, in particular, is able to show the vertical heterogeneity of the tracer arrival times [42]. Thus, we define the control plane as the surface contained between two boreholes, each equipped with electrodes. Using a state-of-the-art inversion procedure (e.g., [6,42]), it is possible to compute the electrical conductivity (σ ) distribution across the CP for as many snapshots as required to fully capture the tracer evolution. The ERT inversion procedure seeks an electrical conductivity distribution that, using an electrical forward model, reproduces the measured raw resistance data to a specified level of uncertainty. This is normally accomplished by solving a regularized optimization problem where the objective function to be minimized can be expressed as a weighted sum of data misfit and a regularization term, usually constructed on the basis of a roughness matrix. The optimal weight between data misfit and regularization term is commonly made dependent on the error level in the data, and the corresponding solution is the smoothest compatible with the data within their error bounds (an Occam’s type solution), as defined by Binley and Kemna [4]. To obtain the solute concentration from σ , a calibrated petrophysical relationship, such as the Archie’s law [2], is usually needed. Independently of the calibrated parameters, the Archie’s law defines an empirical and monotonic relationship between the bulk electrical conductivity and the concentration. In this work, we avoid the need of a petrophysical relationship by using directly the σ evolution to calculate the solute travel time distribution for the assimilation with the EnKF. A discretization grid is defined to achieve a proper solution of the transport problem described in Section 2.1. At each ERT measurement time ti , the electrical conductivity value for each node k belonging to the CP, σ k (ti ), is known. Defined as σ k (t0 ) the background conductivity (i.e., the natural electrical conductivity of the soil), due to the properties of Archie’s law for saturated media, the electrical conductivity variations σk (ti ) = σk (ti ) − σk (t0 ) are proportional to the corresponding concentration variations, Ck (ti ) = Ck (ti ) − Ck (t0 ) (e.g., [30,54]). Therefore, the cumulative distribution function of the travel time at ti can be obtained as:

ti

t =t1

σ¯ (tl )

tl =t1

σ¯ (tl )

CDF (ti ) = tlm

(9)

where σ¯ (ti ) is the sum of the electrical conductivity variations over  all the nodes of the CP (nodCP), i.e., σ¯ (ti ) = nodCP k=1 σk (ti ), and tm is the time when the electrical conductivity distribution takes again the background values, i.e., when all the solute particles have crossed

26

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

Fig. 2. (a) Three-dimensional spatial distribution of hydraulic log-conductivity (Y = ln(K )) in the reference aquifer. The color bar indicates Y values in ln(m/s). (b) Flow field representation in terms of streamlines.

the CP. In other words, it is implicitly assumed that all the tracer is recovered. With the same procedure, the travel time CDF of each sCP is computed by considering the electrical conductivity variations of the nodes belonging to the sCP.

where CDF(tj )k is the CDF value at the tj measurement time for the kth sCP. Note that the measurement vector can only be known completely at the end of the tracer test and therefore data assimilation cannot be carried out in real-time.

2.4. Travel time analysis 3. Test case set-up Given a two-dimensional (2D) control plane, the solute arrival time and its cumulative distribution function gives only information about the mean K and its variance within the domain contained between the injection and the CP. No further information can be deduced about the spatial distribution of K in such domain from the analysis of the travel time CDF. However, for a sedimentary aquifer the main K stratifications (i.e., vertical variations) can be identified by adequately dividing the CP in several horizontal layers (the sub control planes sCPs), each of them characterized by its own travel time CDF. In particular, assuming the main flow as horizontal, Eq. (5) becomes:





¯ l) = M(t; x,

nsCP   i=1

l 0

dz C0 (ax , z)



t 0

dt  g2 (t  ; x¯ − ax , z)Wi

(10)

where l is the total thickness of the injection volume, corresponding to the height of the injection well screening. To properly identify the vertical stratifications, i.e., to find an optimal subdivision of the CP in sCPs from the analysis of travel time data, we suggest using the distribution of t50 , the time at which half of the solute has crossed the CP. Specifically, we compute t50 for each element of the structured grid that discretizes the CP and then we take the horizontal average for each row of elements, the horizontal integral scale of K being significantly larger than the vertical one in sedimentary aquifers; in other words, the two-dimensional distribution of travel times across the CP is lumped into a one-dimensional vertical distribution of t50 . By definition, t50 is representative of the solute average velocity and, consequently, it is considered representative also of the K values that the particles “sample” moving from the injection position to the CP. For instance, high values of t50 indicate that a long time is required for half of the solute to cross the CP, i.e., that the average velocity and the average horizontal K are low, and vice versa. In particular, the subdivision of the CP in sCPs is carried out at the t50 inflection points, as they correspond to points where the solute average velocity reveals significant vertical variations. The travel time CDFs computed for each sCP provide the measurements for the EnKF assimilation procedure. The EnKF measurement vector of Eq. (7) is assembled as:

z(ti ) = [CDF (t1 )1 , CDF (t1 )2 , . . . , CDF (t1 )nsCP , . . . , CDF (ti )1 , . . . , CDF (ti )nsCP ]

(11)

3.1. Synthetic test case The proposed inverse model is first tested in a synthetic threedimensional (3D) heterogeneous aquifer. The model domain has dimensions 16 m × 8 m × 8 m, and is discretized along each direction into cells of side 0.25 m, for a total of 65 × 33 × 33 = 70785 corresponding nodes. Assuming a log-normal hydraulic conductivity distribution (i.e, Y = ln K is normally distributed), the Y reference distribution is obtained as a single unconditional generation [48] conforming to an anisotropic exponential covariance model, with a nominal spatial mean Y  = −4.60 ln(m/s), variance σY2 = 0.50 ln2 (m/s), vertical integral scale λv = 1 m, and horizontal integral scale λh = 8 m. The resulting field is shown in Fig. 2a. The steady state flow field is computed using a standard finite volume solver with the following boundary conditions. Dirichlet boundary conditions are applied at x = 0 m (hydraulic head h = 100.00 m) and at x = 16 m (h = 99.84 m), ensuring a 1% mean gradient, while Neumann no-flow boundary conditions are imposed along the remaining sides of the domain. The resulting Eulerian velocity field is used for the computation of the trajectories of a suitable number of particles, simulating a solute release (e.g., [9]). The flow field is represented in Fig. 2b in terms of streamlines, clearly showing the presence of preferential flow paths in correspondence of regions with large K values. The tracer test is simulated by assuming an injection volume of vertical size 6 m (from z = 1 m to z = 7 m) and longitudinal and transversal size of 0.5 m (centered in x = 0.875 m and y = 4.125 m). The solute is simulated by the release of 28,080 particles, distributed across the injection volume proportional to the fluid velocity, according to the flux concentration concept [31]. The trajectory of each particle is computed by means of the Pollock’s particle tracking post-processing algorithm [46]. Different scenarios, summarized in Table 1, are analyzed in order to investigate the effect of the distance between the CP and the solute source on the Y retrieval and the travel time identification. The CP is located at x = 4.125 m, x = 8.125 m, and x = 12.125 m (corresponding approximately to distances of λh /2, λh , and 1.5λh from the injection well), for scenarios 1, 2, and 3, respectively. Our goal is to assess the hydraulic conductivity values of an ideally layered

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

27

Table 1 Characteristics of the reference field, control plane (CP) location, number of subcontrol planes (sCPs), and prior geostatistical parameters of the synthetic experiments. CP location

Number of sCPs

(x coordinate) (m) Reference Scenario S1 Scenario S2 Scenario S3

4.125 8.125 12.125

7 7 6

Geostatistical parameters

σY2

λv (m)

λh (m)

−4.60 −4.60 −4.60 −4.60

0.5 0.75 0.75 0.75

1 1 1 1

8 ∞ ∞ ∞

Fig. 3. Geographic location and schematic of the tracer test study area with position of the boreholes (in light blue the wells drilled after the tracer test). The tracer was injected in N4, while N5 and N6, where ERT electrodes were located, define the CP (dotted line). The modeled area is denoted by the solid brown line. The solid black lines indicate the corrected CP position and the assumed mean trajectory. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

formation that captures the main heterogeneities of the synthetic 3D domain. The choice of a perfectly layered scheme (i.e., λh = ∞) is due to the use of a single CP and to the integral nature of the information contained in the CDF; in fact, travel times only provide information about the spatial average of the Y values crossed by the solute particles while moving from the source to the CP. Furthermore, the assumption of stratified formations is justified by the fact that in sedimentary aquifers the ratio λh /λv is typically much larger than one, and thus represents a commonly adopted scheme (e.g., [49]). Based on previous studies, an ensemble of Nr = 2000 Monte Carlo realizations are used in this test case, ensuring a good tradeoff between accuracy and computational effort. Each initial K field is generated with statistical properties different from the reference, in order to verify the capabilities of the model to converge to the true solution, despite wrong prior information. The mean value and the vertical integral scale λv are assumed to be known, for example, from pumping or slug tests and the analysis of drilling cores, respectively, while more uncertainties affect the definition of the variance σY2 . Therefore, the initial ensemble of realizations are generated with a 50% error on the true variance value. 3.2. Real field experiment To further evaluate the capabilities of the proposed approach, we apply it to a real field experiment where a saline tracer test monitored by ERT was carried out in September 2010 [42]. The study area is a 270 ha alluvial unconfined aquifer located in Settolo (North-eastern Italy) and recharged mainly by the Piave River. The water table is continuously monitored through a network of 25 observation wells arranged according to a multiscale setup, where a local-scale study area is nested within an intermediate-scale monitored region corresponding to the entire aquifer [60]. In particular, the area of interest for the present study, where the tracer test was carried out, is the local-scale portion of the aquifer (about 40 m × 80 m), equipped with a high-density

measurement network, as schematized in Fig. 3 [60] (http:// settolo.dicea.unipd.it/). The saline tracer was injected in borehole N4, while boreholes N5, N6, and N7 were equipped with electrodes for cross-hole ERT acquisitions. For the entire duration of the test, the flow field was forced by the water withdrawal from the well “base”, which was pumping at a constant rate of approximately 160 l/s (Fig. 3). The boreholes N5 and N6 delimit the CP surface, while boreholes N8, N9, and N10 were drilled after the tracer test, providing some core samples and hydraulic conductivity data collected by slug tests. The modeled area has dimensions of 25 m × 18 m × 25 m and the grid elements that discretize the domain have dimensions of 0.5 m × 1 m × 0.25 m, for a total of 97,869 elements. From the slug tests performed in wells N8, N9, and N10 at three different depths (about 10, 20, and 30 m), we estimated an average Y equal to −3.79 ln(m/s) and σY2 equal to 1.76 ln2 (m/s). Based on the drill cores, we assigned a value of 1 m to the vertical integral scale λv . Noflow boundary conditions are assigned to the top, bottom, and lateral boundaries of the domain, while Dirichlet boundary conditions are applied to the upstream and downstream boundaries using the data acquired from the hydraulic head measurements. The resulting gradient is about 0.74%, with an average velocity of approximately 10−4 m/s, consistent with a recent model calibration of the aquifer [60]. The results of the tracer test in terms of time-lapse resistivity images are described in detail in [42]. The analysis of these images brings forth two main considerations. First, a slower arrival of the solute in the upper part of the CP compared to the lower part is evident. Second, the solute appears to be skewed towards well N5, inconsistently with the expected mean flow direction, the latter being controlled by the withdrawal from well “base” (i.e, positioned along the direction N4-base, represented by the dashed line in Fig. 3). This is due to a private well just outside the study area that was pumping during the test and caused a local alteration of the mean flow field. In particular, the mean solute trajectory crossed the CP close to well

28

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

Fig. 4. Vertical distribution of t50 for the three synthetic scenarios. Dots indicate the inflection points used for the definition of the subcontrol planes.

N5 instead of the midpoint between N5 and N6. Because of the difficulties in accurately describing the real flow path at this small scale, we approximate the real trajectory with the continuous line of Fig. 3 and, to consider a CP consistent with the domain and perpendicular to the mean flow path, we rotated the CP of about 10° with respect to the plane between N5 and N6. This remodeling still ensures a mean travel distance consistent with the real one, with negligible errors in terms of travel time. 4. Results 4.1. Synthetic test case Each scenario of the synthetic test case is characterized by a different distance between the injection and the CP where the measurements are taken (Table 1). As this distance increases, the travel time CDFs are representative of Y values averaged over an increasing distance, producing some differences in the t50 distributions for each scenario, as shown in Fig. 4. The main features of the t50 distribution are similar for all scenarios, although the total vertical extension is reduced for scenarios 2 and 3 with respect to scenario 1, probably due to the presence of preferential flow paths in which the solute canalizes, as shown by the streamlines representation in Fig. 2b. By analyzing the t50 distribution and its derivatives, we identified the inflection points corresponding to the divisions of the CP in sCPs. Hereinafter, we use the convention of enumerating the sCPs from the bottom to the top of the domain. In the first scenario (S1), the distance between the CP and the injection is approximately half of the horizontal integral scale and the CP is divided into seven sCPs. For each sCP the travel time CDF is computed until the time when the whole solute has crossed the CP. In particular, the data are collected at the CP every 3600 s, for a total of 13 available measurement vectors, i.e., updates. For a first verification of the model capabilities, the true travel time CDF for the second sCP is compared with the Nr simulated CDFs computed prior to, during, and after the EnKF update procedure (Fig. 5), i.e., calculated with the prior K fields (open loop), with the K fields obtained after half of the available set of measurements have been assimilated, and with the K fields obtained at the end of the simulation, respectively. The simulated CDFs initially incorporate the true CDF with a large uncertainty, which is progressively reduced as the number of updates increases and the ensemble mean of the CDFs converge towards the true one. Fig. 6 a shows the portion of aquifer retrieved between the injection and the CP, while Fig. 6b reports the corresponding portion of the reference aquifer, where Y is averaged for each row of grid elements across 1.5 m in the transversal direction around the injection volume, i.e., in a volume of dimensions 0.25 m (height of each grid element) along the coordinate z, 1.5 m along y (centered around the injection), and the length between the injection and the CP (different for each scenario) along x. The comparison shows a good agreement

Fig. 5. Scenario 1, sCP2 : ensemble of the Nr simulated travel time CDFs (grey lines) and its average (black line) compared to the true CDF (red line) at times (a) t = 0 s (prior ensemble), (b) t = 26,100 s (intermediate time), and (c) t = 47,700 s (after the last update). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Spatial distribution of Y in the longitudinal cross-section at y = 4 m between x = 0 m and the CP location (x = 4.125 m). Panel (a) shows the Y distribution estimated in scenario 1, while panel (b) shows the average Y values computed in the reference aquifer for each row of grid elements. Thick black and dashed grey lines represent the sCP divisions and the injection upper and lower boundaries, respectively. The color bar indicates Y in ln(m/s).

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

29

Fig. 7. Scenario 1: comparison between the true travel time CDFs (red lines) and the CDFs computed in the retrieved aquifer for each sCP (black lines). Open loop results are shown in grey. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

between the simulated and the true average hydraulic conductivity values, with a satisfactory identification of the mean Y values on the different sCPs. The resulting Y field is further validated by simulating the solute injection in the retrieved aquifer and comparing the travel time CDFs thus obtained for each sCP with the true ones. Fig. 7 shows that the main behavior of the CDFs is reproduced for each sCP, with larger er-

rors associated to the sCPs where the true CDF variance is large, corresponding to regions where the dispersion coefficient, i.e., the variance of Y, is particularly high (e.g., sCP5). This is probably due to difficulties in the proper identification of the sCP subdivisions when the K vertical variations are smoothed. Fig. 7 also shows the results of the open loop simulation, i.e., a run with the prior Y distribution, demonstrating the improvement obtained by the proposed inversion model.

30

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

Fig. 8. Spatial distribution of Y in the longitudinal cross-section at y = 4 m between x = 0 m and the CP location (x = 8.125 m). Panel (a) shows the Y distribution estimated in scenario 2, while panel (b) shows the average Y values computed in the reference aquifer for each row of grid elements. Thick black and dashed grey lines represent the sCP divisions and the injection upper and lower boundaries, respectively. The color bar indicates Y in ln(m/s).

Table 2 Root mean square error of the travel time CDFs (RMSECDF ) and the relative value (percentage) with respect to the total time required by the solute to completely cross the CP (ttot ) at the end of the assimilation procedure (RMSECDF /ttot ) and for OL /ttot ). the open loop simulation (RMSECDF Scenario

RMSECDF (s)

RMSECDF /ttot [%]

RMSEOL CDF /ttot [%]

S1 S2 S3 Real field experiment

1911.7 7912.3 8769.6 18360

3.27 7.25 6.09 9.44

20.57 25.52 30.77 26.25

The error on the retrieved travel time CDFs is quantified through the expression of the root mean square error (RMSE):

 RMSECDF =

nsCP  i=1

((t10,r − t10,t )2 + (t50,r − t50,t )2 + (t90,r − t90,t )2 )i 3 · nsCP (12)

where t10 , t50 and t90 are the times at which 10%, 50%, and 90%, respectively, of the solute has crossed the CP in the retrieved aquifer (subscript r) and in the true one (subscript t). As reported in Table 2, RMSECDF for scenario 1 is approximately 1911 s, compared to a total tracer test duration of about 58,500 s. This relatively small error is proof that the layered schematization provides a good approximation of the reference 3D field. In scenario 2, the distance between the CP and the injection is approximately equal to the horizontal correlation length λh . This means that even though the Y values are still horizontally correlated within this distance (Fig. 2), appreciable variations in the hydraulic conductivity may occur in the longitudinal direction. For example, Y values around z = 4 m vary from about −5.5 ln(m/s) (i.e., K = 4.1·10−3 m/s) at the tracer source to −3.5 ln(m/s) (i.e., K = 3.0·10−2 m/s) at the CP. Overall, the subdivision of the CP in sCPs (Fig. 4) does not show significant differences from scenario 1, although the presence of non-horizontal preferential flow paths in the upper part is suggested by the reduced solute vertical extension. The comparison between the retrieved Y field and the true one horizontally averaged between the tracer source and the CP (Fig. 8a and b, respectively) shows that the Y distribution is well reproduced. The retrieved travel time CDFs (Fig. 9) are similar to the true ones, except

for sCP1, probably due to the large variations of the mean Y values in this sCP, demonstrated by the t50 distribution (Fig. 4, scenario S2), and to a significant difference between the travel time CDF in sCP1 and sCP2, i.e., between the Y values of the two adjacent sCPs. Overall, an improvement is still obtained with respect to the open loop simulation. The RMSECDF , as reported in Table 2, is about 7912 s, i.e., approximately 7% of the total time required by the solute to cross the CP. This error is higher than in scenario 1 but still relatively small. In scenario 3, the CP is located at a distance of approximately 1.5 λh from the tracer source. The vertical distribution of t50 is similar to the ones of scenarios 1 and 2, despite a slightly different distribution of the sCPs. Again, the true vertical distribution of Y is reproduced satisfactorily by the inverse model (Fig. 10) and the retrieved CDFs are comparable to the true ones (Fig. 11). In addition to sCP1, also sCP6 shows significant errors; as in scenario S2, this is likely due to significant t50 variations in both sCPs (Fig. 4, scenario S3). The RMSECDF is about 8769 s, i.e., 6% of the total cumulative travel time (Table 2), consistent with the results of scenario 2. The performance of the method is also evaluated in terms of RMSE of the retrieved Y fields:



RMSEY =

nod Nr 1 1  (Yt,i − Yr,i,mc )2 , nod Nr mc

(13)

i

where nod is the total number of nodes discretizing the aquifer, Yt,i is the hydraulic log-conductivity value at the ith node in the true aquifer and Yr,i,mc is the retrieved value at the ith node for the mcth Monte Carlo realization. The results are reported in Fig. 12, which shows that the time evolution of RMSEY is similar for the three scenarios. After an initial increase, probably due to the fact that only a small portion of the solute has crossed the CP at early times, the final values are comparable and confirm an improvement if compared to the prior estimates. Although this improvement is relatively limited (about 6% less than the prior value), it must be considered that RMSEY is computed across the whole domain and not just in the portion crossed by the plume. 4.2. Real field experiment The resistivity tomographies obtained as described in Section 2.3 are analyzed to compute the t50 vertical distribution and the CP

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

31

Fig. 9. Scenario 2: comparison between the true travel time CDFs (red lines) and the CDFs computed in the retrieved aquifer for each sCP (black lines). Open loop result are shown in grey. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

32

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

Fig. 10. Spatial distribution of Y in the longitudinal cross-section at y = 4 m between x = 0 m and the CP location (x = 12.125 m). Panel (a) shows the Y distribution estimated in scenario 3, while panel (b) shows the average Y values computed in the reference aquifer for each row of grid elements. Thick black and dashed grey lines represent the sCP divisions and the injection upper and lower boundaries, respectively. The color bar indicates Y in ln(m/s).

Fig. 11. Scenario 3: comparison between the true travel time CDFs (red lines) and the CDFs computed in the retrieved aquifer for each sCP (black lines). Open loop result are shown in grey. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

33

Fig. 12. Time evolution of RMSEY computed over the whole model domain.

Fig. 13. Vertical distribution of t50 for the Settolo aquifer. Dots indicate the inflection points used for the definition of the subcontrol planes.

subdivision. The former is reported in Fig. 13 and shows the presence of two main layers: the lower part of the domain is characterized by smaller travel times (i.e., larger hydraulic conductivity) than the upper part of the domain. Five sCPs are identified and for each of them the travel time CDF is calculated. The ERT data were collected approximately every 30 min providing 24 measurement vectors for data assimilation. In Fig. 14, the ensemble of the simulated travel time CDFs is compared with the measured CDF for the second and the fourth sCP after the last assimilation step. The convergence toward the measured CDF is evident, despite a larger uncertainty compared to the synthetic experiments. To further evaluate the performance of the inverse model, the values of hydraulic conductivity obtained from the model application in the real field experiment are compared with values estimated by applying the Rosetta software [50] with the grain size distribution data of the core samples collected in N11 (Fig. 3). Fig. 15a shows the qualitative comparison in terms of ratio K/Kav , where Kav is the vertical average of the K values. The stratification obtained by the proposed inverse approach matches reasonably well the one derived from the core samples. The same results are shown in a more quantitative way in Fig. 15b, which shows a scatter plot representation of the estimated K values against those obtained from the core samples. Compared to the prior values (open loop), the correlation between inverse model results and experimental data is evident (coefficient of determination R2 equal to 0.47),

Fig. 14. Settolo real-world experiment: ensemble of the Nr simulated travel time CDFs (grey lines) and its average (black line) compared to the true CDF (red line) after the last update in (a) sCP2 and (b) sCP4. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

despite significant differences in the order of magnitudes. However, it must be considered that the K values obtained by Rosetta are affected by large uncertainties, due to the lack of additional information such as the soil bulk density; therefore, the capabilty of the inverse model to capture the vertical trend of K is judged satisfactory. After simulating the injection a posteriori in the retrieved K field, we observe an overestimation of the travel times, especially for sCP4 and sCP5, located where Y values are smaller (Fig. 16). However, the overall trend of the simulated travel time CDFs is consistent with the measured CDFs and the RMSECDF is about 5.1 h, approximately 9% of the total time needed by the solute to cross the CP (Table 2). Overall, the results of this application are satisfying, considering the uncertainties associated to the field experiment and the schematization adopted. First of all, the prior ensemble of realizations is generated with statistical parameters estimated from a limited number of slug tests that are typically affected by considerable uncertainties. Second, the definition of the boundary conditions is more uncertain compared to the synthetic test case and the mean flow direction cannot be perfectly known. Finally, unexpected operational difficulties (i.e., as the private pumping well) forced us to introduce some approximations in the domain definition.

34

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

Fig. 15. (a) Vertical distribution of the ratio between local (K) and average (Kav ) hydraulic conductivity obtained from the inverse model (blue bars) and the Rosetta analysis (magenta bars) of the core samples collected in N11 (Fig. 3). The green lines show the value of K/Kav computed from the slug tests. (b) Scatter plot of the K values estimated by the inverse model and those obtained from Rosetta (blue dots) with the corresponding regression line (blue dotted line). Orange and green dots indicate the open loop (prior) values and the values estimated by the slug tests, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16. Settolo real-world experiment: comparison between the true travel time CDFs (red lines) and the CDFs computed in the retrieved aquifer for each sCP (black lines). Open loop results are shown in grey. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

5. Conclusions The purpose of the present work was to evaluate the possibility of retrieving the spatial distribution of hydraulic conductivity (K) and the solute travel times by monitoring tracer tests with the electrical resistivity tomography (ERT). In particular, to avoid the use of solute concentration values, whose estimation would require the calibration of a petrophysical law, we suggest exploiting directly the travel time information contained in ERT imaging. The proposed inversion approach combines a Lagrangian transport model with an Ensemble Kalman filter that uses the travel time cumulative distribution functions (CDFs) collected across a control plane (CP) as measurements. Starting from a state-of-the-art transport theoretical framework, we extended the approach to take advantage of a proper subdivision of the CP into sub-control planes. The method was applied to both a synthetic and a real-world test case, where the simulated aquifers were assumed as perfectly stratified; this hypothesis is justified for sedimentary aquifers and small travel distances and is consistent with the developed approach. In fact, travel time data allow the estimation of the vertical distribution of the average hydraulic conductivity between the solute source and the CP, while no information is available about the K horizontal variations. As a consequence, running the model with a finite value of the horizontal integral scale would not improve the retrieved K field, as with the proposed approach we can only know the time needed by the solute to reach the CP but we cannot recognize the different K zones crossed by the solute. In the synthetic test case we considered three scenarios, characterized by different distances between the solute source and the CP, in order to evaluate the reliability of the perfect layering schematization adopted. Despite an increasing distance between the solute source and the CP, the results were satisfying and the retrieved hydraulic conductivity distribution reproduced the travel time CDFs with an overall RMSE always lower than 10%, compared to the total time required by the solute to cross the CP. Not surprisingly, the best performance was obtained in scenario 1, where the distance between the injection source and the CP is minimum. Compared to the synthetic test cases, the real-world test case was affected by several issues, due to uncertainties in the identification of prior geostatistical parameters, mean flow direction, and proper boundary conditions, which caused the introduction of significant approximations in the definition of the model domain. As a result, the estimated hydraulic conductivity distribution was affected by a large initial error that produced a slower convergence to the solution. Nonetheless, the results obtained were satisfying, as the travel time CDF was reproduced with a RMSE still smaller than 10%, while the estimated K vertical distribution was qualitatively consistent with experimental data obtained from the analysis of core samples and from some slug tests. An additional tracer test is being planned with a more careful characterization of the mean flow direction and a more suitable location of the boreholes delimiting the control plane, such that the proposed inversion approach will be further validated with more accurate field data. Acknowledgment This work was supported by the University of Padova project CPDA115351/11 (“Identification of hydraulic parameters in sedimentary aquifers at the local and catchment scales.”). We thank Mohamad El Gharamti and one anonymous reviewer for their constructive remarks, which helped to improve the quality of the paper. References [1] Annan AP. Hydrogeophysics. Water Science and Technology Library, 50. New York: Springer; 2005.

35

[2] Archie GE. The electrical resistivity log as an aid in determining some reservoir characteristics. J Pet Technol 1942;146(01):54–62 SPE-942054-G. http://dx.doi.org/10.2118/942054-G. [3] Binley A, Cassiani G, Middleton R, Winship P. Vadose zone flow model parameterisation using cross-borehole radar and resistivity imaging. J Hydrol 2002;267(3):147–59. http://dx.doi.org/10.1016/S0022-1694(02)00146-4. [4] Binley AM, Kemna A. Hydrogeophysics. Water Science and Technology Library, 50. New York: Springer; 2005. p. 129–56. [chapter Electrical Methods]. [5] Camporese M, Cassiani G, Deiana R, Salandin P. Assessment of local hydraulic properties from electrical resistivity tomography monitoring of a threedimensional syntetic tracer test experiment. Water Resour Res 2011;47:W12508. http://dx.doi.org/10.1029/2011WR010528. [6] Camporese M, Cassiani G, Deiana R, Salandin P, Binley A. Coupled and uncoupled hydrogeophysical inversions using ensemble kalman filter assimilation of ert-monitored tracer test data. Water Resour Res 2015;51:3277–91. http://dx.doi.org/10.1002/2014WR016017. [7] Carrera J, Medina A, Galarza G. Groundwater inverse problem: discussion on geostatistical formulations and validation. Hydrogeologie 1993;4:313–24. [8] Chen Y, Zhang D. Data assimilation for transient flow in geologic formations via ensemble kalman filter. Adv Water Resour 2006;29(8):1107–22. http://dx.doi.org/10.1016/j.advwatres.2005.09.007. [9] Crestani E, Camporese M, Baú D, Salandin P. Ensemble kalman filter versus ensemble smoother for assessing hydraulic conductivity via tracer test data assimilation. Hydrol Earth Syst Sci 2013;17(4):1517–31. http://dx.doi.org/10.5194/hess17-1517-2013. [10] Cvetkovic V, Shapiro AM, Dagan G. A solute flux approach to transport in heterogeneous formations: 2. uncertainty analysis. Water Resour Res 1992;28(5):1377– 88. http://dx.doi.org/10.1029/91WR03085. [11] Dagan G. Solute transport in heterogeneous porous formations. J Fluid Mech 1984;145:151–77. http://dx.doi.org/10.1017/S0022112084002858. [12] Dagan G. Flow and Transport in Porous Formations. New York: Springer; 1989. [13] Dagan G, Cvetkovic V, Shapiro A. A solute flux approach to transport in heterogeneous formations: 1. the general framework. Water Resour Res 1992;28(5):1369– 76. http://dx.doi.org/10.1029/91WR03086. [14] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. J Geophys Res: Oceans 1994;99(C5):10143–62. http://dx.doi.org/10.1029/94JC00572. [15] Evensen G. The ensemble kalman filter: theoretical formulation and practical implementation. Ocean Dyn 2003;53(4):343–67. http://dx.doi.org/10.1007/s10236003-0036-9. [16] Evensen G. The ensemble kalman filter for combined state and parameter estimation. Control Syst IEEE 2009;29(3):83–104. http://dx.doi.org/10.1109/ MCS.2009.932223. [17] Finsterle S, Kowalsky MB. Joint hydrological-geophysical inversion for soil structure identification. Water Resour Res 2008;7:287–93. http://dx.doi.org/10.2136/ vzj2006.0078. [18] Gelhar L. Stochastic Subsurface Hydrology. Englewood Cliff, NY: Prentice Hall; 1993. [19] Gharamti M, Kadoura A, Valstar J, Sun S, Hoteit I. Constraining a compositional flow model with flow-chemical data using an ensemble-based kalman filter. Water Resour Res 2014;50(3):2444–67. http://dx.doi.org/10.1002/2013WR014830. [20] Hendricks Franssen H, Alcolea A, Riva M, Bakr M, van der Wiel N, Stauffer F, et al. A comparison of seven methods for the inverse modelling of groundwater flow. Application to the characterisation of well catchments. Adv Water Resour 2009;32(6):851–72. http://dx.doi.org/10.1016/j.advwatres.2009.02.011. [21] Hendricks Franssen HJ, Kaiser HP, Kuhlmann U, Bauser G, Stauffer F, Muller R, Kinzelbach W. Operational real-time modeling with ensemble Kalman filter of variably saturated subsurface flow including stream-aquifer interaction and parameter updating. Water Resour Res 2011;47(2):W02532. http://dx.doi.org/10.1029/2010WR009480. [22] Hendricks Franssen HJ, Kinzelbach W. Real-time groundwater flow modeling with the Ensemble Kalman Filter: Joint estimation of states and parameters and the filter inbreeding problem. Water Resour Res 2008;44(9):W09408. http://dx.doi.org/10.1029/2007WR006505. [23] Hinnell AC, Ferré TPA, Vrugt JA, Huisman JA, Moysey S, Rings J, Kowalsky MB. Improved extraction of hydrologic information from geophysical data through coupled hydrogeophysical inversion. Water Resour Res 2010;46(4):W00D40. http://dx.doi.org/10.1029/2008WR007060. [24] Huisman J, Rings J, Vrugt J, Sorg J, Vereecken H. Hydraulic properties of a model dike from coupled bayesian and multi-criteria hydrogeophysical inversion. J Hydrol 2010;380(12):62–73. http://dx.doi.org/10.1016/j.jhydrol.2009.10.023. [25] Irving J, Singha K. Stochastic inversion of tracer test and electrical geophysical data to estimate hydraulic conductivities. Water Resour Res 2010;46(11):W11514. http://dx.doi.org/10.1029/2009WR008340. [26] Jadoon KZ, Slob E, Vanclooster M, Vereecken H. Uniqueness and stability analysis of hydrogeophysical inversion for time-lapse ground-penetrating radar estimates of shallow soil hydraulic properties. Water Resour Res 2008;44:W09421. http://dx.doi.org/10.1029/2007WR006639. [27] Jardani A, Revil A, Dupont J. Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging. Adv Water Res 2013;52(0):62–77. http://dx.doi.org/10.1016/j.advwatres.2012.08.005. [28] Kalman RE. A new approach to linear filtering and prediction problems. Trans ASME–J Basic Eng 1960;82(Series D):35–45. http://dx.doi.org/10.1115/1.3662552. [29] Karahan H, Ayvaz M. Simultaneous parameter identification of a heterogeneous aquifer system using artificial neural networks. Hydrogeol J 2008;16(5):817–27. http://dx.doi.org/10.1007/s10040-008-0279-0.

36

E. Crestani et al. / Advances in Water Resources 84 (2015) 23–36

[30] Kemna A, Vanderborght J, Kulessa B, Vereecken H. Imaging and characterisation of subsurface solute transport using electrical resistivity tomography (ERT) and equivalent transport models. J Hydrol 2002;267:125–46. http://dx.doi.org/10.1016/S0022-1694(02)00145-2. [31] Kreft A, Zuber A. On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions. Chem Eng Sci 1978;33(11):1471–80. http://dx.doi.org/10.1016/0009-2509(78)85196-3. [32] Kurtz W, Hendricks Franssen H-J, Kaiser H-P, Vereecken H. Joint assimilation of piezometric heads and groundwater temperatures for improved modeling of river-aquifer interactions. Water Resour Res 2014;50(2):1665–88. http://dx.doi.org/10.1002/2013WR014823. [33] Lehikoinen A, Huttunen JMJ, Finsterle S, Kowalsky MB, Kaipio JP. Dynamic inversion for hydrological process monitoring with electrical resistance tomography under model uncertainties. Water Resour Res 2010;46(4):W04513. http://dx.doi.org/10.1029/2009WR008470. [34] Li L, Zhou H, Gomez-Hernandez JJ, Franssen HJH. Jointly mapping hydraulic conductivity and porosity by assimilating concentration data via ensemble kalman filter. J Hydrol 2012;428:152–69. http://dx.doi.org/10.1016/j.jhydrol.2012.01.037>. [35] Liu G, Chen Y, Zhang D. Investigation of flow and transport processes at the MADE site using ensemble Kalman filter. Adv Water Resour 2008;31(7):975–86. http://dx.doi.org/10.1016/j.advwatres.2008.03.006. [36] Looms MC, Binley A, Jensen KH, Nielsen L, Hansen TM. Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data. Vadose Zone J 2008;7(1):227–37. http://dx.doi.org/10.2136/vzj2007.0087. [37] McLaughlin D, Townley LR. A reassessment of the groundwater inverse problem. Water Resour Res 1996;32(5):1131–61. http://dx.doi.org/10.1029/96WR00160. [38] Monego M, Cassiani G, Deiana R, Putti M, Passadore G, Altissimo L. A tracer test in a shallow heterogeneous aquifer monitored via time-lapse surface electrical resistivity tomography. Geophysics 2010;75(4):WA61–73. http://dx.doi.org/10.1190/1.3474601. [39] Moradkhani H, Sorooshian S, Gupta HV, Houser PR. Dual state–parameter estimation of hydrological models using ensemble kalman filter. Adv Water Resour 2005;28(2):135–47. http://dx.doi.org/10.1016/j.advwatres.2004.09.002. [40] Müller K, Vanderborght J, Englert A, Kemna A, Huisman JA, Rings J, Vereecken H. Imaging and characterization of solute transport during two tracer tests in a shallow aquifer using electrical resistivity tomography and multilevel groundwater samplers. Water Resour Res 2010;46(3):W03502. http://dx.doi.org/10.1029/2008WR007595. [41] Nowak W. Best unbiased ensemble linearization and the quasi-linear kalman ensemble generator. Water Resour Res 2009;45(4):W04431. http://dx.doi.org/10.1029/2008WR007328. [42] Perri M, Cassiani G, Gervasio I, Deiana R, Binley A. A saline tracer test monitored via both surface and cross-borehole electrical resistivity tomography: comparison of time-lapse results. J Appl Geophys 2012;79:6–16. http://dx.doi.org/10.1016/j.jappgeo.2011.12.011. [43] Pollock D, Cirpka OA. Temporal moments in geoelectrical monitoring of salt tracer experiments. Water Resour Res 2008;44(12):W12416. http://dx.doi.org/10.1029/2008WR007014. [44] Pollock D, Cirpka OA. Fully coupled hydrogeophysical inversion of synthetic salt tracer experiments. Water Resour Res 2010;46(7):W07501. http://dx.doi.org/10.1029/2009WR008575. [45] Pollock D, Cirpka OA. Fully coupled hydrogeophysical inversion of a laboratory salt tracer experiment monitored by electrical resistivity tomography. Water Resour Res 2012;48(1):W01505. http://dx.doi.org/10.1029/2011WR010779.

[46] Pollock DW. Semianalytical computation of path lines for finite-difference models. Ground Water 1988;26(6):743–50. http://dx.doi.org/10.1111/j.1745– 6584.1988.tb00425.x. [47] Rings J, Huisman JA, Vereecken H. Coupled hydrogeophysical parameter estimation using a sequential Bayesian approach. Hydrol Earth Syst Sci 2010;14(3):545– 56. http://dx.doi.org/10.5194/hess-14-545-2010,2010. [48] Robin MJL, Gutjahr AL, Sudicky EA, Wilson JL. Cross-correlated random field generation with the direct fourier transform method. Water Resour Res 1993;29(7):2385–97. http://dx.doi.org/10.1029/93WR00386. [49] Savini F, Salandin P. How well the perfect layering model fits natural formations?. In: Hassanizadeh S, et al Editors RS, editors. Proceedings of the conference on XIV CMWR, 2; 2002. p. 1227–34. [50] Schaap MG, Leij FJ, van Genuchten MT. Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J Hydrol 2001;251(34):163–76. http://dx.doi.org/10.1016/S0022-1694(01)00466-8. [51] Scholer M, Irving J, Binley A, Holliger K. Estimating vadose zone hydraulic properties using ground penetrating radar: the impact of prior information. Water Resour Res 2011;47(10). http://dx.doi.org/10.1029/2011WR010409. [52] Tong J, Hu BX, Yang J. Assimilating transient groundwater flow data via a localized ensemble kalman filter to calibrate a heterogeneous conductivity field. Stoch Environ Res Risk Assess 2012b;26(3):467–78. http://dx.doi.org/10.1007/s00477-0110534-0. [53] Tong J, Hu BX, Yang J. Data assimilation methods for estimating a heterogeneous conductivity field by assimilating transient solute transport data via ensemble kalman filter. Hydrol Processes 2013;27:3873–84. http://dx.doi.org/10.1002/hyp.9523. [54] 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(6):W06013. http://dx.doi.org/10.1029/2004WR003774. [55] Vrugt JA, Stauffer PH, Wöhling T, Robinson BA, Vesselinov VV. Inverse modeling of subsurface flow and transport properties: a review with new developments. Vadose Zone J 2008;7:843–64. http://dx.doi.org/10.2136/vzj2007.0078. [56] Xu T, Jaime Gómez-Hernández J, Zhou H, Li L. The power of transient piezometric head data in inverse modeling: an application of the localized normal-score enkf with covariance inflation in a heterogenous bimodal hydraulic conductivity field. Adv Water Res 2013;54:100–18. http://dx.doi.org/10.1016/j.advwatres.2013.01.006. [57] Zhou H, Gómez-Hernández JJ, Hendricks Franssen H-J, Li L. An approach to handling non-Gaussianity of parameters and state variables in ensemble Kalman filtering. Adv Water Resour 2011;34:844–64. http://dx.doi.org/10.1016/j.advwatres.2011.04.014. [58] Zhou H, Gmez-Hernndez JJ, Li L. Inverse methods in hydrogeology: evolution and recent trends. Adv Water Res 2014;63(0):22–37. http://dx.doi.org/10.1016/j.advwatres.2013.10.014. [59] Zimmerman D, Marsily Gd, Gotway C, Marietta M, Axness C, Beauheim R, et al. A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. Water Resour Res 1998;34(6):1373–413. http://dx.doi.org/10.1029/98WR00003. [60] Zovi F. Assessment of heterogeneous hydraulic properties in natural aquifers at the intermediate scale. Civil and Environmental Sciences, University of Padua; 2014. Ph.D. thesis.