Flow in double-porosity aquifers: Parameter estimation using an adaptive multiscale method

Advances in Water Resources 73 (2014) 108–122

Contents lists available at ScienceDirect

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

Flow in double-porosity aquifers: Parameter estimation using an adaptive multiscale method Philippe Ackerer a,⇑, Nicolas Trottier a,b, Frederick Delay a a b

LHyGeS, Université de Strasbourg/EOST – CNRS, 1 rue Blessig, F-67000 Strasbourg, France Laboratoire de Modélisation des Transferts dans l’Environnement, CEA, Bât. 225, F-13108 Saint Paul lez Durance cedex, France

a r t i c l e

i n f o

Article history: Received 17 March 2014 Received in revised form 1 July 2014 Accepted 2 July 2014 Available online 16 July 2014 Keywords: Groundwater flow Double porosity Inverse method Multiscale parameterization

a b s t r a c t The double porosity model (DPM) concept is a widely used approach to simulate flow in fractured aquifers. Parameter estimation methods to calibrate standard groundwater flow models have been rarely applied to the DPM. In this paper, we make use of a multiscale parameterization which defines the spatial distribution of parameter by an adaptive discretization refined within the iterative inversion procedure. The mathematical developments, including a presentation of the adjoint state equations, and the adaptive multiscale algorithm are described in detail. The parameter estimation procedure is applied to a heterogeneous synthetic case, where each type of parameter is spatially distributed according to a different variogram. Without any prior knowledge of theses variograms conditioning the adaptive multiscale method, the latter is able to retrieve with good precision the reference model for both its local parameter values and their associated spatial distribution. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The quantification of flow in natural fractured reservoirs is an important issue for water resources and remains a very challenging topic. As for other groundwater reservoirs, monitoring water pressure heads and their evolution in time is a key to obtain data that can feed the inversion of a hydrodynamic model. In this study, it is dealt with the inversion of a flow model for fractured aquifers by focusing on two major questions. (1) How to develop an efficient inversion technique that can deal with an important number of hydrodynamic parameters because the aquifer is heterogeneous, and (2) how to propose a parameterization able to retrieve the heterogeneity of the aquifer while being parsimonious on the number of handled parameters? The existing conceptual models for fractured reservoir representations can usually be sorted into two categories [1]: Discrete Fracture Models (DFM), where the fractures are treated explicitly, or continuum models, where fractures and matrix are represented by a single or several separate but overlapping interacting equivalent porous continua (EPC). Reference works on DFM can be found in [2–6], among others. The EPC approach is very popular because it requires a less detailed knowledge of the fracture network (geometry, density, etc.) than DFN [7,8]. Since the pioneering ⇑ Corresponding author. E-mail address: [email protected] (P. Ackerer). http://dx.doi.org/10.1016/j.advwatres.2014.07.001 0309-1708/Ó 2014 Elsevier Ltd. All rights reserved.

works of Barenblatt et al. [9], Warren and Root [10], Moench [11], Long et al. [12] and Neuman [13], the EPC approach has been extended to multiple continua representations within both deterministic and stochastic frameworks, e.g., [14,15]. Reviews and discussions can be found, for instance, in Berkowitz [1] and Neuman [16]. When applicable, the main drawback of EPC approaches is the definition of the equivalent hydrodynamic parameters (porosity and permeability) of each continuum and of the parameters ruling the fluxes exchanged between different continua at the scale of the grid cell. The deterministic estimation of these parameters relies upon the use of an adapted up-scaling technique [17]. However, depending on the fracture network characteristics (density, percolation threshold, characteristic lengths of the fractures compared with the size of the investigated domain, etc.), up-scaling local properties at the scale of the grid cell may be flawed. The stochastic framework then becomes an interesting alternative by considering ensemble averages of the domain properties [1,4,18]. Another alternative is the computation of the local equivalent properties using local discrete fracture networks, e.g., [5,19–21]. Among these various model conceptions, we focus our work on the saturated double porosity model (DPM) in which fluid flow following Darcy’s law only occurs in fractures and in which the matrix is considered for its storage properties. Double porosity models are widely used for large-scale simulations and in engineering. The objective of this paper is to propose a methodology

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

for parameter estimation based on an inverse approach. The main parameters used to describe flow in a double porosity system are the hydraulic conductivity (or transmissivity) of the fracture continuum, the specific storage capacity of both fractures and matrix and the exchange rate coefficient describing water transfer between both continua. It is expected that the level of heterogeneity is different between these parameters (e.g., different variances and correlation lengths for a stochastic parameterization) as well as the impact of the parameters on water head fluctuations due to various stresses (in the case of interference pumping tests, for example). It must be acknowledged that very few metrics based on systematic studies inform on how a heterogeneous DPM would behave regarding its main parameters. Only some partial sensitivity analysis exist in the ongoing literature, but they discuss simplified configurations such as homogeneous and fractal DPM, e.g., [22–24]. Though DPM models have been widely used and enhanced, especially in the domain of oil engineering, e.g., [25,26], surprisingly, the inversion of a DPM has not been abundantly addressed in the last decade, except in a few recent publications, such as Kaczmaryk and Delay [23,24], Le Goc et al. [27], Ray et al. [28], Fahs et al. [29] or Trottier et al. [30]. Kaczmaryk and Delay [23,24] inverted interference pumping tests performed at the experimental site in Poitiers, France. Their DPM model assumed symmetric radial flow in homogeneous and fractal media. In the latter, the hydrodynamic parameters followed power-scaling laws regarding the lag distance between pumped and observed wells. The authors also added a hyperbolic wave propagation to account for the very rapid propagation of pressure depletions in karstic conduits. Le Goc et al. [27] proposed a hierarchical identification method relying on steady-state flow data and geometric information to estimate both the hydraulic properties and the hydraulic structures of channels embedded in a significantly less permeable medium. Ray et al. [28] reconstructed binary media (some kind of DPM) by inverting a Karhunen–Loève expansion of the covariance of the proportions of binary facies (highly permeable versus weakly permeable). The inversion of the Karhunen–Loève expansion resorted on an adaptive Markov Chain Monte Carlo. The parameter estimation approach in Fahs et al. [29] used a parameterization based on a zonation technique, and the inversion algorithm allowed for the identification of both the parameters assumed constant per zone and the geometry of the zones. Finally, Trottier et al. [30] applied an inverse procedure based on a multiscale parameterization to interpret interference pumping tests from the experimental site in Poitiers, France. This multiscale parameterization had already been applied to inverting a single continuum at the local (site) scale [31] as well as at the regional scale [32]. It offers the advantage of discriminating the subareas of a modeled domain that need to be highly parameterized, while others necessitate a few parameters. Differences in the local size of the parameterization can be the consequence of user-defined constraints or the consequence of the inverse problem itself. In the latter case, the number of spatially distributed parameters is automatically increased in subareas where the model does not fit data very well or in subareas where the convergence of the inverse problem is not optimal. The theoretical bases of the parameterization were never described in detail. Therefore, all mathematical derivations applied to transient DPM are properly presented in this paper. In the framework of descent methods for inverse problems, there is no rule grounded in rigorous metrics stating whether a sensitivity-based algorithm is better than a method approximating the gradient components of the objective function. On the one hand, a sensitivity-based method (e.g., a Levenberg–Marquardt algorithm) generates sharp descent directions in the parameter space when the model sensitivities to parameters are calculated by a direct differentiation of the forward problem. Nonetheless,

109

each iteration of convergence requires solving as much forward problems (sensitivity calculations) as the number of sought parameters. On the other hand, gradients approximations based for instance on the so-called adjoint state calculation are less demanding. Irrespective the number of sought parameters, calculating the adjoint state requires the resolution of a single problem (two for flow in a DPM, see hereafter) very similar in its form and computation needs to the forward problem. The downside is that only the gradient components of the objective function are calculated; the sensitivities are unavailable with the consequence of a poor evaluation of the uncertainty on parameters and their cross-correlation. In addition, Quasi-Newton algorithms seeking descent directions in the parameter space only on the basis of gradient components may be inaccurate and require many iterations of convergence. As a rule of thumb, Quasi-Newton algorithms never need more than fifty times the number of iterations used in the sensitivitybased methods. This makes the adjoint-state technique well-suited to highly-parameterized problems. Because the multiscale parameterization discussed hereafter can manipulate a large number of parameters, the automatic inversion of a DPM will be handled by means of a gradient-based method using the adjoint state equations. Given the approximations concealed in these methods, it is important that the adjoint state equations are developed rigorously to obtain precise trajectories. The work is presented as follows: Section 2 addresses the mathematical model used to simulate saturated flow in a DPM; the details of the parameter estimation procedure are described in Section 3; and finally, a numerical test inverting a complex synthetic flow scenario is proposed in Section 4.

2. Mathematical model Many inversion exercises performed at the scale of regional aquifers rely on data in the form of hydraulic heads corresponding to water pressures plus local elevation averaged over the screened portion of monitored wells. Generally, these data barely discriminate the vertical components of flow, which makes researchers assume that the flow is two-dimensional and mainly horizontal (the so-called Dupuit assumption). It is obvious that a two-dimensional approach is a priori not appropriate when the vertical component of flow is non-negligible, for example in the presence of channeled flow in sub-horizontal layers connected by opened vertical fractures. But in this case, it is likely that the continuous approach to fractured reservoirs would not be suited too. Recently, Trottier et al. [30] used a dual continuum in a two-dimensional model to invert interference data from a karstic limestone aquifer. There exist vertical components of flow in this aquifer as evidenced by Chatelier et al. [33]. But at the scale of the aquifer, these vertical components are evenly distributed and homogenize fairly well. The two-dimensional approach becomes applicable provided one is not interested in identifying the local properties of a specific set of fractures. On a technical standpoint, there is no special difficulty in extending the approach discussed below to three-dimensional flow. Both the multiscale parameterization and the adjoint state formulas would remain similar. A three-dimensional approach would only complicate the numerical implementation by manipulating discretized blocks (for the forward problem and the parameterization) instead of planar grids. However, a key question, beyond the scope of the present study, is to know which data provide valuable information to condition accurately three-dimensional flow at the scale of the aquifer. As told above, we doubt that classical head measurements are relevant, and additional information is needed. Incidentally, one should also know how to bridge this information with a flow model.

110

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

For the sake of simplicity, it is assumed in the following that a horizontal two-dimensional flow occurs either in a confined aquifer or in an unconfined system with slight hydraulic head variations compared with the wetted thickness of the aquifer. The flow equations are therefore linear, based on their parameters. Again, choosing a linear two-dimensional approach does not result in some loss of generality in the inversion presented hereafter. In the context mentioned above, the dual continuum approach is the most widely used representation of flow in a fracture-matrix system. Darcian flow occurs in the fracture continuum, and flow is negligible in the matrix, which has the role of a storing reservoir fed by (or feeding) flow in the fractures. This behavior is depicted by the following equations:

Sf

@hf  r:ðT f rhf Þ ¼ rðhm  hf Þ þ q @t

ð1Þ

@hm ¼ rðhf  hm Þ @t

ð2Þ

Sm

The indexes f and m refer to the fracture and matrix continua, respectively. The term hk [L] is the hydraulic head in the continuum k (k = f, m), Sf [–] is the storage capacity of the continuum k, Tf [L2 T1] is the transmissivity of fractures, r [T1] is the exchange rate between the matrix and the fracture continua, and q [L T1] is a sink-source term only concealed in the fracture continuum (volumetric flux per unit of aquifer surface area). The transfer function in (2) is the simplest approach to model the water fluxes exchanged between the matrix and the fractures, e.g., [25]. It corresponds to first-order kinetics which is well known to be a firstorder finite-volume approximation of a diffusion equation. In this framework, the exchange rate r may have different definitions, e.g., [34], but to our opinion, the one with the best physical sense states that r is equivalent to eK(C/V)(Dh/Dl) where, e is the wetted thickness of the aquifer (because two dimensional flow is averaged over the vertical direction), K is the hydraulic conductivity at the interface between fractures and matrix, C is the exchange surface area developed at the interface between fractures and matrix in an averaging volume V, Dh is the mean head difference between fractures and matrix, and Dl is a mean distance between fractures and matrix, often associated with the mean spacing of fractures in a matrix block. The following initial and boundary conditions are associated with the flow Eqs. (1) and (2):

hk ðx; y; tÞ ¼ hk;0 ðx; yÞ for t ¼ 0;

k ¼ f;m

hk ðx; y; tÞ ¼ hk;D ðx; yÞ for ðx; yÞ 2 Ck;D ;

k ¼ f;m

ð3Þ

T f rhf ðx; y; tÞ:gCf ¼ qf ;N ðx; y; tÞ for ðx; yÞ 2 Cf ;N

nþ1

  n n ¼ b hm ; hf

nþ1

nþ1

n

hm;i ¼ hf ;i ð1  es Þ þ hm;i es ; nþ1

s¼

rE jEj þ rE0 jE0 j Sm;E jEj þ Sm;E0 jE0 j

Dt

ð5Þ

nþ1

where hf ;i (resp. hm;i ) is the hydraulic head in the fracture continuum (resp. in the matrix) at time step n + 1 for node i located in the middle of the common edge shared between the adjacent elements E and E0 . It is reminded that with non-conforming finite elements, the element properties and parameters are assigned as uniform values over each element. This makes in (5) that |E| (resp. |E0 |) is the area of the element E (resp. E0 ); rE, rE’ are the exchange rate coefficients assigned to elements E and E0 ; Sm,E, Sm,E’ are the matrix storage capacities of E and E0 , respectively; and Dt is the time step. If we consider an implicit form of Eq. (1) with the right-hand side discretized at time step n + 1, the solution (5) can be reintroduced in the nþ1 discrete form of (1). The substitution of the terms hm by the terms nþ1 n hf and hm modifies the components of the matrix A in (4) but allows for the resolution of both Eqs. (1) and (2) in a single step at each time iteration. This procedure also reduces numerical errors due to a sequential (split) resolution of (1) and (2). 3. Parameter estimation with the adaptive multiscale triangulation 3.1. The objective function

where Ck;D are boundaries of Dirichlet type prescribing head values hk;D , Cf,N are boundaries of Neumann type imposing water fluxes qf,N [L2 T1], and gCf is a unit vector oriented outward and normal to the boundaries Cf,N. There is no need for Neumann boundary conditions for the matrix continuum. By construction, Eq. (2) imposes null fluxes at the boundaries where the matrix hydraulic head is not imposed by a Dirichlet condition. Note also that in many cases, the Dirichlet conditions for the matrix continuum can be overlooked. The local (in space) kinetics in (2) can be solved at any boundary of the system by inheriting from the initial values hm,0 and the boundary conditions of the fracture continuum. The system of differential equations has to be solved numerically using finite volume or finite element methods for heterogeneous aquifers. Whichever the numerical method, the state equation (1) is transformed into a linear system of discrete equations in the form

Ahf

where A is the system matrix, b is the right-hand side vector, and n is the index of time steps. The matrix and vector components depend on the type of boundary conditions, on the sink-source terms, on the type of computation grid and obviously on the local values of the hydrodynamic parameters Tf, Sf, Sm and r. In this work, we used the two-dimensional nonconforming finite element method [35] as the discretization technique of Eqs. (1)–(3). The nonconforming finite element method is very close to the finite element method; the main difference is in the location of the seed nodes of the basis functions that interpolate the state variables hf and hm. These nodes are located at the element vertices for conforming finite elements and in the middle of the element edges for two-dimensional nonconforming finite elements (at the orthocenter of facets for tri-dimensional elements). The nonconforming finite element method associates the advantage of finite volumes (flux continuity at the elements’ edges and exact mass balance at the element level) with the advantage of finite elements (rigorous treatment of anisotropy and flexibility in geometry). More details on nonconforming finite elements are available in Crouzeix and Raviart [36], for example. Regarding Eq. (2) in the form of local first-order kinetics, a straightforward algebraic solution can be found, provided that the variable hf is assigned a constant value. By using an Euler implicit scheme for (2) with a constant hf, the solution to (2) becomes

ð4Þ

The formal differentiation of the adjoint state (explained below) is easier, at least for mathematical notations, when all simulated state variables for all time-steps are introduced in the objective function. In the common case where no reliable prior information on the parameters is available, the objective function will only include a measure of the distance between simulated state variables (heads) and observed variables:

Jðh; hðhÞÞ ¼

Nt X

^ k  hk ÞT Wðh ^ k  hk Þ ðh

ð6Þ

k¼1

where k is the time step, Nt is the number of time steps, ()T is the transpose operator, and h is the parameter vector (transmissivity, ^ k is the vector of size N with the storage capacity, etc.) of size Np. h c ^ k ¼ hk , if there is a measured head value hk in well components h o;j o;j i ^k ¼ hk , where j located in cell (or node) i at time step k. Otherwise, h i i k k ^k  h ¼ 0. N is hi is the computed head value that locally makes h c i i the number of cells (or nodes) of the discretized domain. W is a

111

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

Nc  Nc matrix weighting each error on heads in the objective function. Assuming that the measurement errors are not correlated in space and are invariant with time, this matrix is diagonal, and wii ¼ e1 m;j , where em,j is the variance of the measurement error at well j located in cell (or at node) i. In the absence of a measured value in i, ^k  hk ¼ 0. the weight wii can take any scalar value because h i i As stated earlier when justifying the assumption of two-dimensional flow, the head values observed in wells are averages. On the one hand, for wells connected to the fracture network of the aquifer, the common assumption is that the measured values correspond to heads hf of the fracture continuum. On the other hand, poorly connected wells could depend on head values hm of the matrix continuum. Recently, Delay et al. [37,38] suggested that head data from interference testing in a karstic aquifer might show non-reciprocal responses. In short, a general diffusion equation, irrespective the modeled mechanism, is usually reciprocal (the Lorentz reciprocity), which means that a stress at location B generates a response at location A equivalent to the response in B for the same stress in A. In dual porosity media, the violation of this reciprocity principle could be due to observations in some wells that mix the fracture (hf) and matrix (hm) responses [37]. To avoid any loss of generality in the inversion of data from a dual continuum, we consider that the head k signal monitored in wells ho;j could be the head value in the fractures or in the matrix or a linear combination of both. Therefore, after the discrete equations (5) and (6) have been solved, the vector of the computed head values hk is defined as k

k

k

h ¼ ahf þ ðI  aÞhm

ð7Þ

where a is a diagonal weighting matrix (0 6 aii 6 1:0), and I is the identity matrix. Basically, the weighting matrix can vary in time and therefore conceal unknown values that become parameters of the inverse problem. For the sake of simplicity, it is assumed here that the matrix a is known and invariant in time. 3.2. The differentiation of the adjoint state The minimization of the objective function is performed using a gradient-based technique requiring the differentiation of the objective function with respect to all the components of the parameter vector h. A technique of optimization under constraints is used, e.g., [39,40], and we define the following constraints between the parameters h and the state variables hf and hm

uk ¼

k Ahf

  k1 k1  b hf ; hm

k

k1

ð8Þ

k

wk ¼ hm  es hm  ð1  es Þhf

The Lagrangian function associating the objective function (6) with the above constraints is

Lðh;hf ;hm ; k; lÞ ¼

Nt h X

k T

^ k  h Þ Wðh ^ k  h Þ þ uk;T kk þ wk;T lk ðh k

i

i dL d h d h k;T k i d h k;T k i k ¼ Jðh; h Þ þ u k þ w l dhp dhp dhp dhp k

ð11Þ

The differentiations of each of the three terms in (11) are detailed in Appendices A–C, respectively. All calculations completed, one obtains

i d h k ^ k  hk ÞT W½I  es ðI  aÞ dhf Jðh; h Þ ¼ 2ðh dhp dhp k  T ^ kþ1  hkþ1 W½es ðI  aÞ dhm 2 h dhp k

k

þ

@Jðh; h Þ @hp

ð12Þ

!  k T d @u @kk dkk ½uk;T kk  ¼ kk þ uk;T þ dhp @hp @hp dhp 2 ! !T ! !T k T k1 T dhf dhf @ uk @ uk 4 þ þ k k1 dhp dhp @hf @hf 3 ! !T k1 T dhm @ uk 5kk þ k1 dhp @hm 20 ! !T  !T 1 3 n T n T k k X dhf @k dh @k m 4@ Auk 5 ð13Þ þ þ n n dhp dhp @hf @hm n and

!T  k  d k;T k @wk @ l dlk ½w l  ¼ lk þwk;T þ dhp @hp @hp dhp 2 ! ! !T !T T T k k dhf @wk dhm @wk 4 þ þ k k dhp dhp @hf @hm 3 ! ! T k1 T dhm @wk 5lk þ k1 dhp @hm 20 ! !T  !T 1 3 n T n T k k X dhf @ l dh @ l m 4@ Awk 5 ð14Þ þ þ n n dh dh @h @h p p f m n Finally, the total differentiation of Lk with respect to any parameter k hp becomes the following, after reassembling all terms in dhf =dhp k and dhm =dhp : 2 ! !T X @ uk T k dL @kk dkk @wk 4 þ ¼ k þ uk;T þ dhp @hp @hp dhp @hp k # !T  k  k X dhkf @l dlk @Jðh; h Þ þ þ þ @hp @hp dhp dhp k 2 !T !T !T k @ ukþ1 @wk 4 @u kk þ kkþ1 þ lk k k k @hf @hf @hf # !T X dhk k T ^k s m  2½I  e ðI  aÞW ðh  h Þ þ dhp k 2 !T !T !T k kþ1 kþ1 @w @u 4 @w lk þ lkþ1 þ kkþ1 k k k @hm @hm @hm # XX ^ k  hk Þ þ  2½es ðI  aÞWT ðh

lk þ wk;T ð9Þ

k¼1

where ()T is the transposition operator, and kk and lk are vectors of size Nc at time iteration k of unknown scalar Lagrange multipliers. The vectors kk and lk are also called the adjoint states of the vectors k k of variables hf and hm , respectively. Very often, the adjoint state is calculated by stating that the state variables, the adjoint states and the parameters are independent vectors. This assumption greatly simplifies the algebraic developments but may be erroneous. In this work, elaboration of the equations ruling the adjoint state is performed rigorously without any prior guess. The total differentiation of the Lagrangian function in (9) with respect to a given parameter hp is k

Nt X dL dL ¼ dhp k¼1 dhp

with

ð10Þ

20

n

dhf 4@ dhp

!T

k

@k n @hf

!T

k



n T dhm

n k

@k n @hm

!T 1 3 Au k 5

þ dhp 20 ! !T  n T !T 1 3 n T XX dhf @ lk dhm @ lk A k 5 4 @ þ w þ n n dhp dhp @hf @hm n k

ð15Þ

112

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

We remind that the vectors k and l correspond to Lagrange multipliers that contribute to the problem of optimization as many unknowns as the number of constraints enclosed in u and w. The vectors k and l are calculated so that they allow for the elimination of the terms multiplying the head derivatives with respect to parameters in (15), i.e., 8  T  T  T k kþ1 k > ^ k  hk Þ ¼ 0 > < @ uk kk þ @ u k kkþ1 þ @wk lk  2½I  es ðI  aÞWT ðh @hf

@hf

@hf

@hkm

@hkm

@hkm

 T  T  T > > : @wk lk þ @wkþ1 lkþ1 þ @ ukþ1 kkþ1  2½es ðI  aÞWT ðh ^ k  hk Þ ¼ 0

for k ¼ 1; .. .; N t

ð16Þ Using (16), the definition of u and w from (8) leads to

8  T kþ1 > ^ k  hk Þ > < AT kk ¼ @b k kkþ1 þ ð1  es Þlk þ 2½I  es ðI  aÞWT ðh @hf

 T > > : lk ¼ es lkþ1 þ @bkþ1 kkþ1 þ 2½es ðI  aÞWT ðh ^ k  hk Þ k @hm

ð17Þ Both equations are solved backwards in time assuming an ‘‘initial’’ condition kNt þ1 ¼ lNt þ1 ¼ 0. The constraints u and w in L are expressions corresponding to null values when the direct problem is solved. In this context, dL/ dhp is no different than dJ/dhp, i.e., the component of the gradient of the objective function along the direction hp in the parameter space. By using the property in (16) and by canceling out u and w in (15) and using (8), one obtains

!T k X@ uk T dL dJ @wk @Jðh; h Þ ¼ ¼ kk þ lk þ dhp dhp @hp @hp @hp k " #  T k  T X @A k @b @es k @Jðh; h Þ k1 k k ¼ h  k þ hf  hm l þ @hp f @hp @hp @hp k ð18Þ As shown in (18), the calculation of each component of the gradient of the objective function sums up to scalar products. The adjoint state can also be derived to handle uncertain boundary conditions needing for inversion. As shown in Section 2, the dual continuum flow model is constrained by boundary conditions very similar to those of a single continuum. Even though we did not focus this study on the identification of boundary conditions, the specific adjoint state equations associated with this problem are developed in Appendix D. In the same vein, we derived above the adjoint state equations for a dual porosity model where no flow occurs in the matrix. The dual porosity model is the most widely used, but we could also extend the adjoint state equations to a dual permeability model with a flowing matrix. The flow equations in both fracture and matrix continua become similar. The adjoint state equations are derivable in the same way as for the above development. We note that the adjoint state equations for the matrix would be very similar to that of the fractures. The simplification due to the substitution of the matrix head variable in the flow equation of the fractures (see Eq. (5)) would no longer hold, making the terms in es in the adjoint equations to disappear. These terms would be replaced by linear combinations of the heads in the fractures and the heads in the matrix because the two flow equations are coupled and also because one can still consider that the observed heads are a combination of signals from the fractures and the matrix. 3.3. Parameterization and optimization procedure The Adaptive Multi-scale Triangulation (AMT) relies on the discretization of the modeled domain by triangular elements. The dif-

ferential equations of the direct problem are solved on a fine grid T sim with parameters hsim defined at the Nsim elements of T sim . On the other hand, the parameters hropt manipulated by the optimization procedure at the convergence iteration r are defined at the N rh nodes of a coarser grid T rh superimposed onto T sim . The parameter grid T rh can be locally refined at will between successive iterations. Notably, the grid T rh can be different for each type of parameter (e.g., transmissivity versus storage capacity). The multi-scale triangulation has many advantages because it can be refined at will and at specific locations. One can have for instance data distributed as patches with dense information separated by large areas where information lacks. The identification of hydraulic parameters is obviously weakly conditioned in regions lacking data, and it becomes useless to refine the parameter grid. In the same vein, one can be interested in refining the parameter grid locally because the prior knowledge on the system states that specific flow structures must be inferred. As shown later, the parameter grid can also be refined according to local convergence criteria of the inverse problem. This feature allows for a parsimonious parameterization increasing the number of sought parameters only in critical areas. The implementation of the AMT method is ruled by the following steps (see also Fig. 1). 1. An initial triangulation T 0h is defined. The initial values of each parameter are defined at the N 0h nodes of T 0h . 2. At each iteration r, the mapping of the parameters from T rh to T sim is performed by

hrsim ¼ Mr hropt

ð19Þ N rh .

r

where M is a mapping matrix of size Nsim  The components of this matrix depend on the interpolation technique between the nodes of T rh and the triangles of T sim . In this work, we use a linear interpolation to define the local parameters hsim,E of an element E of T sim . For an element E included in (beneath) an element G of T rh , we have

hrsim;E ¼

3 X

xGi hr;G opt;i

ð20Þ

i¼1

where hr;G opt;i is the parameter at iteration r defined at the vertex i of the element G of T rh . The term xGi is the linear function defined by

(

xGi ðx; yÞ ¼ ai x þ bi y þ ci xGi ðxc ; yc Þ ¼ dic for c ¼ i; j; k

ð21Þ

where i, j, and k are the vertices of element G and di,c, the Kronecker delta function. 3. After mapping the parameters from T rh to T sim , the flow model is conditioned based on its parameters. Then, the forward problem is solved in addition to the two adjoint states of hydraulic heads. 4. If the objective function is smaller than the user-defined threshold esim, the computation is stopped because a valuable solution to the inverse problem has been found. However, the gradient of the objective function is computed according to Eq. (18). 5. If the gradient of the objective function is greater than the userdefined threshold eopt, the optimization continues. A quasiNewton method is used in which the Hessian matrix needed to calculate the descent direction of the parameter vector rþ1 hopt  hropt is approximated by a BFGS algorithm (Broyden– Fletcher–Goldfarb–Shanno), e.g., [41]. The routines are provided rþ1 by HSL [42]. If necessary, the new values hopt of the parameters are modified to respect user-specified constraints. The latter simply rely here on a so called box-constraints projection P of ~ hopt defined as

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

113

Fig. 1. General flowchart of the optimization algorithm that seeks inverse solutions of spatially distributed problems using AMT parameterization and adjoint state calculations.

8 rþ1 > < hrþ1 opt;i ( hopt;i > : hrþ1 opt;i

h i min max if hrþ1 opt;i 2 hopt;i ; hopt;i   h i max min min max 1 if hrþ1 ( hmin opt;i þ u0 hopt;i  hopt;i opt;i R hopt;i ; hopt;i ð22Þ

where u10 is a random number uniformly distributed between 0 and max 1, and hmin opt;i (or hopt;i ) is the a priori lower boundary (or upper boundary) to parameter hopt,i variations. 6. If the gradient of the objective function is smaller than the user-defined threshold eopt, the objective function may be considered to have reached a minimum. However, this minimum or the local errors between the simulated variables and the observed variables may remain too high. In this case, it is assumed that the parameter grid T rh is too coarse to correctly describe the spatial heterogeneity of parameters impacting local observed head data. Therefore, T rh is refined according to a local criterion (i.e., on each triangle G of T rh ) defined as the sum of the local components of the gradient of the

objective function

@J=@hr;G opt;i

[43,44] or the local value of the

objective function. If a triangle G of T rh has to be refined, it is divided into four smaller triangles by joining the middle points of the three sides of the initial triangle and then assigning parameter values to the nodes of the smaller triangles (see below). Because the refinement criterion is established at the element level, the continuity of the parameter field hsim mapped from hopt is not ensured at common edges between adjacent triangles of T opt . Stated differently, the linear interpolation in (20) and (21) is not continuous along an edge separating a single triangle of T h on one side from the two triangles on the other side. rþ1 7. Once the new parameter grid T opt is designed, one must assign the parameters located at the new vertices. For the vertices inherited from the previous grid T ropt , the values remain unchanged. For vertices c at the midpoint of edges, e.g., i-j, the new value is defined by

  8 r;G > hmin ¼ min hr;G > opt;i ; hopt;j > <   r;G r;G > hmax ¼ max hopt;i ; hopt;j > > : Lnðhopt;c Þ ¼ Lnðhmin Þ þ u10 ½Lnðhmax Þ  Lnðhmin Þ

ð23Þ

Notably, we handle here natural logarithms of parameter values simply because these values may vary over several orders of magnitude and the difference hmax  hmin may have no meaning because it is almost equal to hmax, whereas the difference between natural logarithm values clearly separates high orders of magnitudes from small ones. The random perturbation term (u10 ) in (22) provides some stochastic behavior to the overall inversion procedure. It helps movement within the parameter space along several directions while seeking solutions. This feature can be useful when the inverse problem is poorly posed and has a non-unique solution because of weak sensitivity to parameters or over-parameterization. 4. Numerical experiments The inversion technique relying on the AMT parameterization and the calculation of adjoint states is tested regarding a synthetic problem of horizontal two-dimensional flow in a confined fractured reservoir represented by a dual fracture-matrix continuum. The modeled domain has a rectangular shape with length 2000 m and width 1000 m. The boundary conditions are of Neumann type and prescribe null-fluxes through the limits of the domain, except for two portions (the northwest and southeast boundaries) assigned with Dirichlet type conditions that prescribe constant head values of 100 and 90 m (see Fig. 2). For the sake of simplicity, we did not define uncertain boundary conditions to be fitted by the inversion procedure. Nevertheless, the adjoint state can handle this type of unknowns as shown in Appendix D. The initial head values are uniform over the whole domain at 100 m in both the fracture and matrix continua. Starting at time t = 0, two pumping wells, denoted #1 and #2 in Fig. 2, then extract water from the aquifer at the constant rate of 4.16  103 m3 s1

114

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

Fig. 2. The two-dimensional domain of 2000 m length and 1000 m width handled by the inversion exercise. Thin lines of the frame correspond to Neumann null-flux boundary conditions; thick lines correspond to Dirichlet boundary conditions. Empty circles, denoted as wells no. 1 and no. 2, correspond to the locations of pumping wells, and the small black dots are the locations of observed drawdown data. Specifically, observations points referenced by P1 to P4 are locations where the simulated and observed heads are compared in Fig. 5.

and 3.33  103 m3 s1, respectively. This type of pumping stress does not correspond to a classical interference testing where a single well is pumped. Pumping simultaneously in several wells will result in a more variable head field that helps to identify the aquifer heterogeneities. In fractured media, increasing the pumping stress also induces a quicker response of the matrix and diminishes the total pumping time needed to observe the complete behavior of the aquifer. Solving the flow model depicted by Eqs (1)–(3) allows the calculation of the drawdown responses to the pumping stress. These responses are monitored over 1.16 days (i.e., 105 s) at 62 wells randomly distributed over the domain (Fig. 2). These wells are assumed to be part of a piezometer network for water level survey, each borehole being of small diameter and unable to house hydraulic tests. This case is frequent in concrete applications and justifies the absence of relevant prior information on hydrodynamic parameters. The usual response of a homogeneous infinite single continuum results in a linear evolution of the drawdowns with the logarithm of time. This is also the case for drawdowns at early and late pumping times in an infinite homogeneous dual continuum. At early times, the hydraulic diffusion expresses as the ratio of the fracture transmissivity to the fracture storage capacity. Because this diffusion is large, the drawdown responses become rapidly linear with the logarithm of time. At late times, the dual continuum behaves as a single medium with a hydraulic diffusion as the ratio of the fracture transmissivity to the matrix storage capacity. At mid pumping times when the matrix of a dual continuum takes over for the fractures to feed the extracting well with water, the drawdowns are almost flat with time (and its logarithm). Knowing these features, the drawdown sampling in the present exercise is chosen to obey a logarithmic time step, which avoids the oversampling of long-term behaviors compared with short-term behaviors. To our knowledge, the only example where specific samplings at very early times of pumping and with high frequency were useful is one of the cases depicted in Bodin et al. [45]. The problem was to invert interference data in response to an extraction-injection doublet in a karstic aquifer. The very rapid propagation of the pressure depletion through open conduits justified to capture the very early time responses. As is often the case when dealing with real drawdown responses in fractured media, we assume that the hydraulic head

monitored in open boreholes is the response of the fracture continuum hf e.g., [30]. This means that the possibility evoked in Section 3 involving data mixing both fracture and matrix signals will not be explored here. One can mention that such mixtures of signals may induce reciprocity gaps in the drawdown responses (see Section 3). Up to now, it is still not clear how to handle both reciprocal and non-reciprocal information in the same inverse problem and what is the incidence on the inverse solutions. In the end, 43 measurements of hf at each monitored well are available (for a total of 2666 head measurements to condition the inverse problem). To mimic the common real case where head data are associated with measurement errors, each sampled drawdown value is corrupted with a non-correlated uniform random noise ranging between 5 and +5 cm. This value may appear slightly high, but our experience in manipulating data from fractured limestone aquifers [30,31] shows that such head measurement errors may appear. Usually, they are not stemming from flawed pressure gauges but from boreholes due to the partial collapsing of their walls. Our goal is to address, in a single exercise, how the AMT parameterization and the adjoint state technique reliably retrieve different types of heterogeneity in a reservoir. The flow problem discussed above is therefore solved by considering, for the same flow scenario, that (i) the ratio s = r/Sm of the exchange rate coefficient between fracture and matrix to the matrix storage capacity is uniform at the value log (s) = 3, with log() as the decimal logarithm; (ii) the fracture storage capacity is mildly heterogeneous (r2log Sf = 0.16) and obeys a covariance function of Gaussian type with a large effective correlation length of 900 m; and (iii) the fracture transmissivity is also mildly heterogeneous (r2log10 T f = 0.2) but with an exponential covariance of shorter effective correlation length of 300 m. The sampled variograms of the reference parameter fields log (Sf) and log (Tf) are reported in Fig. 3. As shown by the settings just above, the parameter fields are not very heterogeneous and the mean lag distance between observation wells is far smaller than the parameter correlation lengths. We could have duplicated the flow scenarios over different type of parameter fields, but the flow equations in dual continua are not very sensitive to parameters. In our opinion, it is easier to roughly identify contrasted heterogeneities than to depict accurately a smoothly varying random field. Notably, AMT parameterization and the adjoint state technique have already been used for concrete case applications in highly heterogeneous aquifers e.g., [30,31].

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

Fig. 3. Sampled variograms of the reference parameter fields used to build the inversion exercise. These variograms handle the decimal logarithm of the fracture transmissivity (log Tf) and the fracture storage capacity (log Sf).

The inverse problem is solved without any prior guess on the spatial structure of the parameters or of their local values. The objective function is that of Eq. (6), which merely handles errors on heads hf between model simulations and ‘‘noisy’’ data from the synthetic reference problem (see above). Because head measurement errors are known a priori as being uncorrelated and of constant variance, irrespective of the space–time location of the measurement, the weighting matrix W in (6) is simply taken as the identity matrix. Notably, introducing only data from the fracture continuum does not mean that one addresses a simplified inversion in which the adjoint state associated with the matrix head hm could be overlooked. The state variables hf and hm are coupled by the flow problem, which requires an inversion based on the complete calculation of the adjoint states, as depicted in Section 3. The inversion resorting in the AMT parameterization seeks three types of parameters: Tf, Sf, and s. The choice of inverting s (and not independently seeking the parameters r and Sm) is motivated by the fact that r and Sm appear as a ratio in the flow equations linking the variables hm and hf (Eqs. (2) and (5)). The inverse problem is therefore mainly sensitive to s = r/Sm. To be more precise, the flow equations are not very sensitive to Sm and sensitive to r, as they are to Tf and Sf [22]. The matrix storage capacity Sm it therefore roughly tuned independently in the first iterations of the inverse problem. Then, when its values no longer evolve, (usually, when there is enough water available at the scale of the whole

115

system), the ratio s = r/Sm is inverted by modifying r and letting Sm unchanged. We also investigate the possible impact on inverse solutions of the initial parameter triangulation T 0h . The latter could influence the successive refinement processes (see Section 3.3) and therefore modify the inverted parameter fields interpolated from seeds located at the vertices of the triangulation T h . Six initial types of parameter grid are tested (geometry in Fig. 4a and three refinement steps in Fig. 4b), and the same grid is used when seeking the three types of parameters. This choice of a single grid complicates the task assigned to the AMT parameterization because in the present example, the same grid geometry must retrieve very different parameter spatial distributions (uniform, short correlation length, and large correlation length). For each AMT initial geometry, the same flow scenario is inverted 50 times by varying the initial parameter values preset at the nodes of T 0h . This multi-start procedure is not intended to provide a Monte-Carlo approach to the inverse problem. It simply addresses the robustness of the inversion by identifying the variability of optimal solutions when starting from different locations in the parameter space. One may expect that accurate trajectories in the parameter space calculated with rigorously developed adjoint state equations are able to steer the inverse procedure toward the same solution or at least, a narrow set of plausible solutions. In the present exercise where the inverse problem is conditioned by an important amount of data, a robust inversion technique should identify solutions close to the reference problem. The results are depicted in Figs. 5–8. Most of the 300 (50  6) solutions show simulated heads very close to the data. An example from four selected observed wells (locations in Fig. 2) with various drawdown behaviors is plotted in Fig. 5. With reference to the general drawdown behavior with time in a dual continuum (see above), on can check that the logarithmic time step samples evenly the whole drawdown curves (dots in Fig. 5). In all cases, the errors between simulations and observations are always within the measurement error of ±5 cm. The threshold esim (Section 3.3) used to stop the optimization and prescribed here at 2666 (i.e., the total number of data points) was reached after 4 or 5 grid refinement steps, irrespective the initial geometry of T 0h . This means that the number of parameters (by counting all the parameters at all the vertices of the AMT grid) required to obtain a relevant solution ranges between 363 and 606 (Table 1). One cannot, under any circumstances, consider that a good solution is obtained by overparameterizing the inverse problem. Although the number of parameters is not huge, it is however large enough to justify an

Fig. 4a. The six types of initial AMT grids used in the inversion exercise. Their numbering from 1 to 6 corresponds the references used in Figs. 6 and 8. These initial grids are dimensioned to exactly overly the flow field in Fig. 2.

116

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

Fig. 4b. Three grid refinements for initial grid 4.

adjoint state calculation and a Quasi-Newton algorithm to seek the optimal pathways in the parameter space. Both methods are reputed to be useful when dealing with highly parameterized problems, though the downside sometimes involves a lack of or

inaccurate convergence. In the present case, a complete development of the adjoint state associated with a parsimonious parameterization are the keys to a solution very close to the reference problem (see below). With regards to the sought parameter fields, the spatial distribution of s is almost uniform over the whole domain and very close to the constant value of the reference solution (i.e., log (s) = 3), with decimal logarithm values that range between 3.095 and 2.859, irrespective of the initial grid T 0h . Fig. 6 illustrates the sought fields of Tf and Sf. Independent of the geometry of the initial grid T 0h , the inverse solutions look very similar to the reference fields (top plot of Fig. 6), albeit occasionally with slight smoothing effects as the consequence of the linear interpolation between the seeds of the AMT grid. At first glance, there is no special discrepancy between the reference and the inverse solution appearing, for example, at locations that lack head data. Though the initial grids T 0h can be of different geometry, their successive refinements conditioned on local values of the objective function or its gradient (see Section 3.3) yield similar results. This observation reinforces the idea that the AMT parameterization is well suited to reveal the parameter heterogeneity because it comes from the inversion of the data, while ignoring any prior guesses of the spatial structure (covariance) of parameters. The absence of conditioning onto a spatial structure of parameters does not forbid the AMT parameterization from identifying spatially structured solutions, as shown in Fig. 6. Incidentally, obtaining the same inverse solutions, irrespective the initial AMT grid, argues on the robustness of the methodology. If a parameterization technique is aimed at identifying parameter heterogeneity, the latter must be identified on the basis of available data and not on artifacts due to the initialization of the inverse procedure. We analyzed the 50 equiprobable solutions for each type of initial AMT grid. For each AMT grid, the envelope of the variograms, described by the mean value of the variogram at a given distance plus/minus its standard deviation, is plotted in Fig. 7 (shaded area) and is compared to the exact solution. The envelopes of each AMT are very similar, and no significant effect of the initial grid has been detected. The exact solution is included within the envelopes, and one could draw from the calibrated fields reasonable uncertainty ranges for both the variance and the correlation length of the parameters. The uncertainty concerning the transmissivity field is

8

8

Drawdown (m)

P2 10

Drawdown (m)

P1 10 6 4 2

6 4 2 0

0 1

100

1

10,000

P4 10

8

8

Drawdown (m)

P3 10 Drawdown (m)

100

10000

Time (s)

Time (s)

6 4 2

0

6 4 2 0

1

100

10000

1

Time (s)

100

10000

Time (s)

Measured

Simulated

Fig. 5. Four examples of fittings between simulated and observed heads. The observations points, denoted P1 to P4 (location in Fig. 2), are chosen for the variability that they show in the time behavior of drawdowns.

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

AMT grid Scale

log(Tf) -3.25 -3.16 -3.06 -2.97 -2.88 -2.78 -2.69 -2.59 -2.50

117

log(Sf) -6.80 -6.70 -6.60 -6.50 -6.40 -6.30 -6.20 -6.10

Solution

1

2

3

4

5

6

Fig. 6. Reference fields and examples of inverse solutions for the decimal logarithm of the fracture transmissivity (log Tf) and the fracture storage capacity (log Sf). The indexes 1 to 6 refer to the initial geometry of the AMT parameterization grid (see Fig. 4). Dots on the reference solution correspond to locations of observed head data. At first glance, the inverse solutions appear very similar to the reference, irrespective of the initial AMT grid geometry.

quite low, whereas the uncertainty concerning storage capacity is significant. This may be due to its sensitivity, which is significant only over a short period after drawdown starts. Though not tested in the present case study, one could try to improve the evaluation of the storage capacity by giving more weight to the early pumping time data in the objective function. Another method would consist of artificially stretching or shrinking the gradient components of the objective function. This would provide some anamorphosis of the parameter space to improve descent directions for some

parameters, as done in Delay et al. [22] with model sensitivities to parameters. If we calculate the local coefficient of variation (the ratio of the standard deviation to the mean) of the parameters through the 50 equiprobable solutions for each type of initial AMT grid, these coefficients are very small for both Tf and Sf. The values are on the order of 0.05–0.15 for Tf, 0.1–0.25 for Sf, with a few values reaching 0.4 at locations close to the northwest corner of the domain. As mentioned above, the NW corner is assigned a constant head boundary

118

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

Fig. 7. Exact variograms (line and symbols) of the log fracture transmissivity (right) and log fracture storage capacity (left) compared to estimated variograms for AMT 5 (grey area represents the mean value of c ± standard deviation for each inversion based on AMT5).

Fig. 8. Plots of the coefficient of variation for the log fracture transmissivity (left) and the log fracture storage capacity (right). The coefficients of variations are calculated over 50 equiprobable inverse solutions corresponding to 50 different starting points in the parameter space. Indexes 1 to 6 refer to the initial AMT grid geometries reported in Fig. 4. Back dots correspond to locations of observed head data. In general, the coefficients of variation are very small except for the storage capacity in the northwest corner of the domain deeply influenced by the prescribed head boundary condition.

119

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

Table 1 Number of refinements of the parameter grid and total number of parameters (summing all the variables at all the vertices of the parameter grid) required to fit the observed heads.

Nb. Ref. Nb. Param.

Grid no 1

Grid no 2

Grid no 3

Grid no 4

Grid no 5

Grid no 6

5 606

4 363

4 582

4 558

4 564

4 447

condition. It is well known that transient drawdowns from interference testing in dual continua are sensitive to the fracture storage capacity, especially for short pumping times [22]. With head values in the vicinity of a Dirichlet condition not really allowed to follow the time variations resulting from pumping tests, the inverse procedure generates more erratic Sf values, resulting in larger coefficients of variation. Reinforcing this argument, the coefficients of variation of Sf are the smallest in the vicinity of the pumped wells (see Figs. 8 and 2), i.e., where the drawdown amplitudes are the largest, and result in a better sensitivity of the flow model to the fracture storage capacity. Concerning Tf, high and low values of the coefficient of variation are evenly distributed over the domain and not influenced by the location of the pumping wells or the density of available head data. In the end, the overall inversion algorithm appears robust, independent of the initial parameter grid geometry and able to retrieve simple to complex spatial distributions of parameters in the same inversion exercise. 5. Conclusion Even though Darcian flow in a dual continuum appears to be a valuable assumption for simulating hydrodynamics at the scale of a fractured reservoir, automatic inversion of this approach has not received much attention to date. Previous works showed that inversion of a dual continuum model is not straightforward because model sensitivities to parameters should be calculated rigorously to obtain reliable results, especially with very transient flow, such as those associated with interference well testing, e.g., [22–24]. To the best of our knowledge, this present study represents one of the first times that a complete methodology inverting flow data from a heterogeneous dual continuum approach has been proposed. This methodology relies on an adjoint state calculation of the gradient components of the objective function and a parameterization that maps the parameters onto the computation grid of the forward problem.

of parameter by refining the AMT grid is performed only where needed, i.e., where the model does not fit the local data or where one can expect better precision because the local gradient component of the objective function is not minimal.  The AMT is free from any prior guess on the eventual spatial structure of the parameters. According to the method used to refine the AMT grid (see just above), the heterogeneity of the inverted parameter fields is simply revealed by the inversion of flow data. The AMT technique just identifies the local parameter heterogeneity required to fit the data; it does not modify a prior guess of this heterogeneity for the model to fit the data. This feature may become interesting in the frequent actual cases where no reliable information on the medium heterogeneity is available. The AMT can also easily handle non-stationary heterogeneous fields by setting side-by-side the mildly homogeneous and highly heterogeneous subareas. The AMT technique will simply generate both coarse and refined zones in the same parameter grid.  Associated with a good identification of the pathways to be followed in the parameter space (by means of a complete adjoint state), the AMT technique will retrieve the various heterogeneity levels and spatial structures associated with each type of parameter involved in the same inversion exercise. This search of the best spatial field for each parameter is independent of the initial parameter grid used to launch the inversion procedure. The successive refinements of the parameter grids will rub out the eventual bias generated by the initial grid. This result reinforces the argument stating that the AMT seeks only what is needed for the parameters to invert the problem. Though the parameterization is free from any prior guesses, it will retrieve spatial structures (e.g., covariance) in the sought parameters only if this structure is required for the model to fit the data.

Acknowledgments  The adjoint state is a technique often used in highly parameterized inversions because it allows for rapid calculations of the gradient components of the objective function, though at the cost of a loss of precision. The complete development proposed here guarantees that the adjoint state will not suffer from oversimplification and that the gradient components are calculated with the same accuracy as the ones characterizing the numerical resolution of the forward problem. We also develop an adjoint state dealing with data in which the monitored signal could mix both responses from fracture and matrix continua. This option is not tested here but could be interesting for further works inverting flow in natural fractured media, which frequently do not obey the Lorentz reciprocity principle.  The AMT parameterization of a spatially distributed problem is based on mapping parameters from a coarse parameter grid onto the refined computation grid of the forward problem. As for other techniques, the AMT strongly reduces the parameterization of the inverse problem. However, the AMT adds the possibility to locally increase the parameterization for a depiction with higher resolution of the medium heterogeneity. Some type of parsimony principle is used because increasing the number

This work is supported by the French Research National Agency (project RESAIN no ANR-12-BS06-0010-02) and N. Trottier benefits from a CEA grant. Appendix A The objective function is defined by

Jðh; hðhÞÞ ¼

X ^ k  hk ÞT Wðh ^ k  hk Þ ðh

ðA1Þ

k k

k

k

with h ¼ ahf þ ðI  aÞhm , k being the index of iteration in time. The total differentiation of the objective function with respect to any parameter hp is

dJðh; hðhÞÞ @Jðh; hðhÞÞ @Jðh; hðhÞÞ dhf @Jðh; hðhÞÞ dhm ¼ þ þ dhp @hp @hf @hm dhp dhp

ðA2Þ

We refer to the following discretized water balance equation for the head in the matrix continuum k

k1

k

hm ¼ es hm þ ð1  es Þhf

ðA3Þ

120

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

The definition of hk is therefore rewritten as: k

k

X

k1

h ¼ ½I  es ðI  aÞhf þ es ðI  aÞhm

ðA4Þ

n

n

dhf @ k;T k n ½u k  dhp @hf

!

In view of the definition (A1), differentiating J with respect to hp as done in (A2) would yield for any iteration in time n = 1, . . . , Nt n n k X @J dhf ^ k  hk ÞT W @h dhf ¼ 2 ðh n n @hf dhp @hf dhp k n

k

ðA5Þ

n

X @J dhm ^ k  hk ÞT W @h dhm ¼ 2 ðh n n @hm dhp @hm dhp k k

n

k

! !T ! !T 3 k T k1 T k k dhf dhf @ u @ u 5kk þ ¼4 k k1 dhp dhp @hf @hf 0 ! !T 1 n T k X dhf @k @ þ uk A ðB4Þ n dhp @hf n

A development similar to that in (B3) and (B4) can be applied to the fourth term of the right hand side of (B2). Considering n @ uk;T =@hm ¼ 0 except for n = k  1, one obtains

X n

n

@ k;T k dhm n ½u k  dhp @hm

!

n

Noting that @h =@hf ¼ 0; k – n, @h =@hm ¼ 0; k – n þ 1, and using (A4), Expression (A5) becomes n dhf

n dhf

@J ^ n  hn ÞT W½I  es ðI  aÞ ¼ 2ðh n dhp @hf dhp n

n

By reintroducing (A6) in (A2), the differentiation of the objective function is finally rewritten as X dJðh; hðhÞÞ @Jðh; hðhÞÞ ¼ 2 dhp @hp k " # k k T k ^ k  h Þ W½I  es ðI  aÞ dhf þ ðh ^ kþ1  hkþ1 ÞT W½es ðI  aÞ dhm ðh dhp dhp

ðA7Þ We note that the first term of the right hand side of Equation (A7) is oJ(h, h(h))/ohp = 0 except for hp = s, where it becomes

h  i X @Jðh; hðhÞÞ ^ k  hk ÞT W ðI  aÞ hk  hk1 ¼ 2es ðh f m @s k

ðA8Þ

Appendix B We differentiate here the second term of L in Eq. (11) and corresponding to the scalar product at time iteration k of the adjoint state kk and the vector of constraints uk k1 uk ¼ Ahkf  bðhk1 ; hm Þ f

ðB1Þ k;T

!  k T n k X @ dhf @u k k;T @k k;T k k þu þ n ½u k  @hp @hp dhp @hf n ! n X @ k;T k dhm þ n ½u k  dhp @hm n n X @ k;T k dk þ n ½u k  @k dhp n

n

n

k

Since @ u

n =@hf

!

! n X @ X dkn T @ uk T k dkk k;T k dk ¼ k þ uk;T n ½u k  n @k dhp dhp @k dhp n n

ðB6Þ

Because the constraints u do not depend on k (see (B1)), Expression (B6) is simplified into n X @ dkk k;T k dk ¼ uk;T ½ u k  n @k dhp dhp n

ðB7Þ

In the end, Expression (B2) sums up to

!  k T d @u @kk dkk ½uk;T kk  ¼ kk þ uk;T þ dhp @hp @hp dhp 2 ! ! ! !T T T k k1 T dhf dhf @ uk @ uk 4 þ þ k k1 dhp dhp @hf @hf 3 ! !T k1 T dhm @ uk 5kk þ k1 dhp @hm 20 ! !T  !T 1 3 n T n T k k X dhf @k dh @k m 4@ Auk 5 ðB8Þ þ þ n n dhp dhp @hf @hm n

Appendix C

0 ! !T 1 n T X dhf @ uk @ ¼ kk A n dhp @hf n 0 ! !T 1 n T k X dhf @k @ þ uk A n dhp @hf n

¼ 0 except for n = k and n = k  1, one obtains

We differentiate here the third term of L in Eq. (11) corresponding to the scalar product at time iteration k of the adjoint state lk and the vector of constraints wk k

k1

k

wk ¼ hm  es hm  ð1  es Þhf

ðC1Þ

The total differentiation of the scalar product wk,Tlk with respect to any parameter hp is

ðB2Þ

The third and fourth terms of the right hand side of (B2) have a very similar form. Taking for instance the third term, one can write

dhf @ k;T k n ½u k  dhp @hf

ðB5Þ

k

The total differentiation of the scalar product u k with respect to any parameter hp is

X

! !T k1 T dhm @ uk ¼ kk k1 dhp @hm 0 !T 1 X dhn T @kk m @ þ uk A n dhp @hm n

The last term of the right hand side in (B2) can be rewritten as

ðA6Þ

@J dhm ^ nþ1  hnþ1 ÞT W½es ðI  aÞ dhm ¼ 2ðh n dhp @hm dhp

d ½uk;T kk  ¼ dhp

2

ðB3Þ

d ½wk;T lk  ¼ dhp

!T

 n @ lk X @ k;T k dh þ n ½w l  @hp dhp @h n ! n X @ k;T k dhm þ n ½w l  dhp @hm n  X @ d ln k;T k ½w l  þ @ ln dhp n @wk @hp

lk þ wk;T

ðC2Þ

A development similar to that reported in Appendix B can be carried n

out. By noting in view of (C1) that @wk =@hf ¼ 0 except for n = k, and k

@w

n =@hm

¼ 0 except for n = k and n = k  1, one obtains

121

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

!T  k  d @wk @ l dlk ½wk;T lk  ¼ lk þ wk;T þ dhp @hp @hp dhp 2 ! !T !T !T ! !T 3 k T k k1 T k dhf @w dhm @wk dhm @wk 5lk þ4 þ þ k k k1 dhp dhp dhp @hf @hm @hm 20 ! !T  !T 1 3 n T n T X dhf @ lk dhm @ lk A k 5 4@ þ þ ðC3Þ w n n dhp dhp @hf @hm n

In the end, the expression of the gradient components of the objective function will sum up to (see also Eq. (18))

X dL dJ ¼ ¼ vk dhp dhp k " # k X  @A k @b T k @es k @Jðh; h Þ k1 T k ¼ h  k þ ðh  hm Þ l þ @hp f @hp @hp @hp f k ðD4Þ

Appendix D Let us take the boundary conditions of a flow problem in a dual continuum as the parameters of the inverse problem. The aim is to obtain the gradient components of the objective function with respect to these specific parameters using the adjoint state method. It must be reminded that boundary conditions are prescribed entities in the numerical resolution of the (forward) flow problem. On the one hand, it does not make sense to compare these non-calculated entities with any measure of them. Obviously, if any measured information on boundary conditions is available, the model should prescribe its boundary conditions according to the measured values. On the other hand, it makes sense to evaluate the incidence of prescribed boundary values on the state variable (heads) of the flow problem and then, modify the boundary conditions to better fit the model onto available head data. In this context, the objective function remains unchanged compared with that of the classical inverse problem, and writes

Jðh; hÞ ¼

X ^ k  hk ÞT Wðh ^ k  hk Þ ðh

 T dJ @b ¼ kk dhp @hp

ðD5Þ

Except for the scalar product in (D5), obtaining these gradient values does not require any other calculations than that of an adjointstate for the general flow problem. We note also that these gradient components only depend on the adjoint state k which is associated with the constraints u derived from the calculation of the heads in the fractures hf. It does not mean that the heads in the matrix hm are not sensitive to boundary conditions. As shown in Eq. (D2), the vector b depends on hm and the resolution of the adjoint state k is also coupled with that of l.

ðD1Þ References

k

^ k is the vector of observed where k is the index of iteration in time, h heads, hk is the (composite) vector of simulated heads, and W is a weighting matrix. Because the boundary conditions are enclosed in the matrix system solving the flow problem, and also because the flow problem remains unchanged in terms of size and manipulated variables, the constraints attached to the Lagrangian function obey exactly the same definition as in Eq. (8), i.e. k1 uk ¼ Ahkf  bðhk1 ; hm Þ f k

k1

ðD2Þ

k

wk ¼ hm  es hm  ð1  es Þhf

A and b are the matrix and the vector of the discrete resolution of the flow problem, hf, hm are the vectors of calculated heads in the fractures and in the matrix, respectively. The vectors of constraints uk, wk are associated with the scalar vectors of Lagrange multipliers kk and lk, respectively, and corresponding to the adjoint state. A close look at the algebraic developments reported in Appendices A–C, shows that all the differentiations with respect to a parameter hp are independent of the type of parameter, be it for example a transmissivity value or a prescribed boundary condition. All calculations done, the differentiation of the Lagrangian function with respect to a parameter hp will obey to Eq. (15) in the main text. By simplifying the notations, this equation can be rewritten as

!T !T X X dhkf X dhk dL k k m ¼ v þ n þ fk dhp dhp dhp k k k i X Xh i XXh T T þ ð1k;n Þ uk þ ð0k;n Þ wk k

Knowing that: (1) the objective function and the parameter s do not depend on the boundary conditions; (2) in most discretization schemes of partial differential equations, the boundary conditions are enclosed in the vector b, the gradient of the objective function with respect to a parameter associated with a boundary condition will be written as

n

k

ðD3Þ

n

When the forward problem is solved, the constraints uk and wk correspond to null values that cancel out the fourth and fifth terms of the right-hand side of (D3). Solving the adjoint state consists in canceling out the vectors nk and fk, which will result in solving exactly the same linear system as in Eq. (17) of the main text.

[1] Berkowitz B. Characterizing flow and transport in fractured geological media: a review. Adv Water Resour 2002;25(8–12):861–84. http://dx.doi.org/10.1016/ S0309-1708(02)00042-8. [2] Dverstorp B, Andersson J. Application of the discrete fracture network concept with field data: possibilities of model calibration and validation. Water Resour Res 1989;25(3):540–50. http://dx.doi.org/10.1029/wr025i003p00540. [3] Cacas MC. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation. I: The flow model. Water Resour Res 1990;26(3):479–89. http://dx.doi.org/10.1029/89wr00041. [4] Dorn C, Linde N, Le Borgne T, Bour O, de Dreuzy J-R. Conditioning of stochastic 3-D fracture networks to hydrological and geophysical data. Adv Water Resour 2013;62(PA):79–89. http://dx.doi.org/10.1016/j.advwatres.2013.10.005. [5] Follin S, Hartley L, Rhén I, Jackson P, Joyce S, Roberts D, Swift B. A methodology to constrain the parameters of a hydrogeological discrete fracture network model for sparsely fractured crystalline rock, exemplified by data from the proposed high-level nuclear waste repository site at Forsmark, Sweden. Hydrogeol J 2013:1–19. http://dx.doi.org/10.1007/s10040-013-1080-2. [6] Berrone S, Pieraccini S, Scialò S. An optimization approach for large scale simulations of discrete fracture network flows. J Comput Phys 2014;256:838–53. http://dx.doi.org/10.1016/j.jcp.2013.09.028. [7] Follin S, Thunvik R. On the use of continuum approximations for regional modeling of groundwater flow through crystalline rocks. Adv Water Resour 1994;17(3):133–45. http://dx.doi.org/10.1016/0309-1708(94)90037-x. [8] Selroos J-O, Walker DD, Ström A, Gylling B, Follin S. Comparison of alternative modelling approaches for groundwater flow in fractured rock. J Hydrol 2002;257(4):174–88. http://dx.doi.org/10.1016/s0022-1694(01)00551-0. [9] Barrenblatt G, Zheltov I, Kochina I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J Appl Math Mech 1960;24:852–64. http://dx.doi.org/10.1016/0021-8928(60)90107-6. [10] Warren JE, Root PJ. The behavior of naturally fractured reservoirs. SPE J 1963;3(3):245–55. http://dx.doi.org/10.2118/426-pa. [11] Moench A. Double-porosity models for a fissured groundwater reservoir with fracture skin. Water Resour Res 1984;20(7):831–46. http://dx.doi.org/ 10.1029/wr020i007p00831. [12] Long JCS, Remer JS, Wilson CR, Witherspoon PA. Porous media equivalents for networks of discontinuous fractures. Water Resour Res 1982;18(3):645–58. http://dx.doi.org/10.1029/wr018i003p00645. [13] Neuman SP. Stochastic continuum representation of fractured rock permeability as an alternative to the REV and fracture network concepts. In: Farmer IW, Daemen JJK, Desai CS, Glass CE, Neuman SP, editors. Rock Mechanics. Proceedings of the 28th US symposium, Tucson, Arizona. Rotterdam: Balkema; 1987. p. 533–61. [14] Wu Y-S, Liu HH, Bodvarsson GS. A triple-continuum approach for modeling flow and transport processes in fractured. J Contam Hydrol 2004;73(4):145–79. http://dx.doi.org/10.1016/j.jconhyd.2004.01.002.

122

P. Ackerer et al. / Advances in Water Resources 73 (2014) 108–122

[15] Rutqvist J, Leung C, Hoch A, Wang Y, Wang Z. Linked multicontinuum and crack tensor approach for modeling of coupled geomechanics, fluid flow and transport in fractured rock. J Rock Mech Geotechn Eng 2013;5(1):18–31. http://dx.doi.org/10.1016/j.jrmge.2012.08.001. [16] Neuman SP. Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeol J 2005;13(1):124–47. http:// dx.doi.org/10.1007/s10040-004-0397-2. [17] Renard P, de Marsily G. Calculating equivalent permeability: a review. Adv Water Resour 1997;20:227–50. http://dx.doi.org/10.1016/s0309-1708(96) 00050-4. [18] Neuman SP. Universal scaling of hydraulic conductivities and dispersivities in geologic Media. Water Resour Res 1990;26(8):1749–58. http://dx.doi.org/ 10.1029/wr026i008p01749. [19] Jackson CP, Hoch AR, Todman S. Self-consistency of a heterogeneous continuum porous medium representation of a fractured medium. Water Resour Res 2000;36(1):189–202. http://dx.doi.org/10.1029/1999WR900249. [20] Jourde H, Fenart P, Vinches M, Pistre S, Vayssade B. Relationship between the geometrical and structural properties of layered fractured rocks and their effective permeability tensor. J Hydrol 2007;337(1–2):117–32. http:// dx.doi.org/10.1016/j.jhydrol.2007.01.027. [21] Fourno A, Grenier C, Benabderrahmane A, Delay F. A continuum voxel approach to model flow in 3D fault networks: a new way to obtain upscaled hydraulic conductivity tensors of grid cells. J Hydrol 2013;493(17):68–80. http://dx.doi.org/10.1016/j.jhydrol.2013.04.010. [22] Delay F, Kaczmaryk A, Ackerer P. Inversion of interference hydraulic pumping tests in both homogeneous and fractal dual media. Adv Water Resour 2007;30(3):314–34. http://dx.doi.org/10.1016/j.advwatres.2006.06.008. [23] Kaczmaryk A, Delay F. Interference pumping tests in a fractured limestone (Poitiers – France): inversion of data by means of dual-medium approaches. J Hydrol 2007;337(1–2):133–46. http://dx.doi.org/10.1016/j.jhydrol.2007. 01.025. [24] Kaczmaryk A, Delay F. Improving dual-porosity-medium approaches to account for karstic flow in a fractured limestone: application to the automatic inversion of hydraulic interference tests (hydrogeological experimental site, HES – Poitiers – France). J Hydrol 2007;347(3– 4):391–403. http://dx.doi.org/10.1016/j.jhydrol.2007.09.027. [25] Di Donato G, Huang W, Blunt M. Streamline-based dual porosity simulation of fractured reservoirs. In: Proceedings of SPE annual technical conference and exhibition, Denver, Colorado, USA, SPE 84036; 2003. [26] Huang W, Di Donato G, Blunt M. Comparison of streamline-based and gridbased dual porosity simulations. Soc Petrol Eng J 2004;43(2):129–37. http:// dx.doi.org/10.1016/j.petrol.2004.01.002. [27] Le Goc R, de Dreuzy J-R, Davy P. An inverse problem methodology to identify flow channels in fractured media using synthetic steady-state head and geometrical data. Adv Water Resour 2010;33:782–800. http://dx.doi.org/ 10.1016/j.advwatres.2010.04.011. [28] Ray J, McKenna SA, van Bloemen Waanders B, Marzouk YM. Bayesian reconstruction of binary media with unresolved spatial structures. Adv Water Resour 2012;44:1–19. http://dx.doi.org/10.1016/j.advwatres.2012. 04.009. [29] Fahs H, Hayek M, Fahs M, Younes A. An efficient numerical model for hydrodynamic parameterization in 2D fractured dual-porosity media. Adv Water Resour 2014;63(1):179–93. http://dx.doi.org/10.1016/j.advwatres. 2013.11.008. [30] Trottier N, Delay F, Bildstein O, Ackerer P. Inversion of a dual-continuum approach to flow in a karstified limestone: insight into aquifer heterogeneity

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41]

[42] [43]

[44]

[45]

revealed by well-test interferences. J Hydrol 2014;508:157–69. http:// dx.doi.org/10.1016/j.jhydrol.2013.10.039. Ackerer P, Delay F. Inversion of a set of well-test interferences in a fractured limestone aquifer by using an automatic downscaling parameterization technique. J Hydrol 2010;389(1–2):42–56. http://dx.doi.org/10.1016/ j.jhydrol.2010.05.020. Majdalani S, Ackerer P. Identification of groundwater parameters using an adaptive multiscale method. Ground Water 2011;4:548–59. http://dx.doi.org/ 10.1111/j.1745-6584.2010.00750.x. Chatelier M, Ruelleu S, Bour O, Porel G, Delay F. Combined fluid temperature and flow logging for the characterization of hydraulic structure in a fractured karst aquifer. J Hydrol 2011;400:377–86. http://dx.doi.org/10.1016/ j.jhydrol.2011.01.051. Landereau P, Noetinger B, Quintard M. Quasi-steady two equation models for diffusive transport in fractured porous media: large-scale properties for densely fractured systems. Adv Water Resour 2001;24(6):863–76. http:// dx.doi.org/10.1016/s0309-1708(01)00015-x. Vohralík M, Wohlmuth B. From face to element unknowns by local static condensation with application to nonconforming finite elements. Comput Methods Appl Mech Eng 2013;253:517–29. http://dx.doi.org/10.1016/ j.cma.2012.08.013. Crouzeix M, Raviart PA. Conforming and non-conforming finite element methods for solving the stationary Stokes equations. RAIRO Anal Numer 1973;7:33–76. Delay F, Ackerer P, Guadagnini A. Theoretical analysis and field evidence of reciprocity gaps during interference pumping tests. Adv Water Resour 2011;34:592–606. http://dx.doi.org/10.1016/j.advwatres.2011.02.006. Delay F, Ackerer P, Belfort B, Guadagnini A. On the emergence of reciprocity gaps during interference pumping tests in unconfined aquifers. Adv Water Resour 2012;46:11–9. http://dx.doi.org/10.1016/j.advwatres.2012.06.002. Sun NZ, Yeh WWG. A stochastic inverse solution for transient groundwater flow: parameter identification and reliability analysis. Water Resour Res 1992;28(12):3209–20. http://dx.doi.org/10.1029/92wr00683. Michalak AM, Kitanidis PK. Estimation of historical groundwater contaminant distribution using the adjoint state method applied to geostatistical inverse modeling. Water Resour Res 2004;40:W08302. http://dx.doi.org/10.1029/ 2004wr003214. Byrd RH, Lu L, Nocedal J, Zhu C. A limited memory algorithm for bound constraint optimization. Northwest Univ, Dpt Electrical Eng Computer Science, Tech Rep, NAM-08; 1994. HSL. A collection of Fortran codes for large scale scientific computation; 2011 http://www.hsl.rl.ac.uk. Ben Ameur H, Chavent G, Jaffré J. Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities. Inverse Prob 2002;18:775–94. http://dx.doi.org/10.1088/ 0266-5611/18/3/317. Hayek M, Ackerer P. An adaptive subdivision algorithm for the identification of the diffusion coefficient in two-dimensional elliptic problems. J Math Model Algorithm 2007;6:529–45. http://dx.doi.org/10.1007/s10852-006-9046-1. Bodin J, Ackerer P, Boisson A, Bourbiaux B, Bruel D, de Dreuzy J-R, et al. Predictive modeling of hydraulic head responses to dipole flow experiments in a fractured/karstified limestone aquifer: insights from a comparison of five modeling approaches to real-field experiments. J Hydrol 2012;454– 455:82–100. http://dx.doi.org/10.1016/j.jhydrol.2012.05.069.