Inversion of interference hydraulic pumping tests in both homogeneous and fractal dual media

Advances in Water Resources 30 (2007) 314–334 www.elsevier.com/locate/advwatres

Inversion of interference hydraulic pumping tests in both homogeneous and fractal dual media Frederick Delay a

a,*

, Anne Kaczmaryk a, Philippe Ackerer

b

UMR 6532, CNRS, University of Poitiers, Earth Sciences Building, 40 Avenue du Recteur Pineau, F-86022 Poitiers, France b Fluid and Solid Mechanics Institute of Strasbourg, 2 Rue Boussingault, F-67000 Strasbourg, France Received 16 February 2006; received in revised form 22 May 2006; accepted 14 June 2006 Available online 22 August 2006

Abstract This work proposes a complete method for automatic inversion of data from hydraulic interference pumping tests based on both homogeneous and fractal dual-medium approaches. The aim is to seek a new alternative concept able to interpret field data, identify macroscopic hydraulic parameters and therefore enhance the understanding of flow in porous fractured reservoirs. Because of its much contrasted sensitivities to parameters, the dual-medium approach yields an ill-posed inverse problem that requires a specific optimization procedure including the calculation of analytical sensitivities and their possible re-scaling. Once these constraints are fulfilled, the inversion proves accurate, provides unambiguous and reliable results. In the fractal context inverting several drawdown curves from different locations at the same time reveals more accurate. Finally, hydraulic parameters drawn from inversion should be taken into account to improve in various situations the conditioning of up-scaled flow in fractured rocks.  2006 Elsevier Ltd. All rights reserved. Keywords: Interference pumping test; Inverse problem; Fractured porous reservoirs; Dual medium; Fractal media

1. Introduction Over the last decade, research in hydrodynamics and transport in underground reservoirs has been greatly developed, probably owing to a few strategic stakes such as the safety assessment of underground repositories, the prospecting of oil and water resources or the managing of existing reservoirs. Among topics related to this research, the simulation of synthetic (numerical) reservoirs has come up as a crux since it allows the mimicking at low cost of prospective exploitation scenarios. Whatever the method used for simulating numerical reservoirs, e.g., [27,28], the conditioning on available data is needed because it dedicates a generic tool to the specificity of the handled problem. This is why, despite the increasing place taken by numerical simulations in the toolbox of reservoir engineer-

*

Corresponding author. E-mail address: [email protected] (F. Delay).

0309-1708/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2006.06.008

ing, field studies such as interference pumping tests remain a necessity. Basically, an interference pumping test is carried out by pumping a well at a prescribed flow rate and by observing the hydraulic head drawdown at one or several wells located in the vicinity of the pumped well as a function of time. The interpretation of drawdown curves yields averaged values of the hydraulic conductivity and storage capacity of the reservoir. Although the spatial support over which the averaging is performed is not fully defined [30], interference tests condition the flow at a scale that is compatible with that of the reservoir exploitation. Very often, simple models prevail over the interpretation of drawdown curves because they yield values that are easily comparable to each other (which would not be the case if the complex heterogeneity of a reservoir was to be described) and they provide quality criteria, e.g. for deciding whether or not a portion of the reservoir is exploitable. In this context, the Cooper–Jacob analytical method [15] is the simplest tool for interpreting interference pumping tests. It assumes a

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

Darcian 2D convergent flow with cylindrical symmetry around the pumped well in a porous homogeneous medium. The drawdown in time (t) of the hydraulic head h [L] at a location r from the pumped well is written as: H0  h(r, t) = (Q/4peK) ln(2.25Kt/r2Ss) with H0 [L] the initial steady-state head, Q [L3 T1] the pumped flow rate, e [L] the reservoir thickness, K [L T1] the hydraulic conductivity, Ss [L1] the specific storage capacity. Due to its simplicity, this method has been widely used even in heterogeneous reservoirs and has become an inescapable reference, at least in the first stages of any field study. It may fail, however, in fractured rocks for which drawdown curves often increase more rapidly than linearly with ln(t) [10]. Note that fractured rocks have been increasingly exploited for oil and water in the past twenty years and that new analytical or semi-analytical models have been developed to account for their specific behavior, e.g., [2,7,14,17,26]. All are based on the statement that a noninteger dimension prevails for flow. Alternately, the dual-medium approach depicts homogenized flow (and transport) in a fracture field including relationships with the porous matrix. The concept was developed in the early sixties [6,42] and has been widely re-used, yielding an impressive abundance of works and ideas to grasp the macroscopic behavior of porous fractured rocks (e.g., spatial averaging [33,34], discrete fracture networks [20]). Surprisingly, no specific tool based on the dual-medium approach for interpreting interference pumping tests is reported in the literature. Nevertheless, a few manipulations of the equations rapidly show that the dual-medium approach yields drawdown curves that are similar to those observed in fractured rocks at least for short and intermediate pumping times. For very long times, dual media are stressed as a single medium yielding drawdowns linear with ln(t) for the homogeneous case and convex with ln(t) for the fractal case. This clear separation of behaviors between early and late pumping times justifies the alternative of the dual-medium approach compared with a single medium be it homogeneous or fractal. To our opinion, two reasons may explain the lack of tool based on dual media. First, the dual-medium approach has been mostly used in oil engineering and few attempts have been made with water resources. In the oil industry, interference tests (needing initial steady-state conditions) are scarce because they are too expensive. Most tests are performed at the well upon drilling or during exploitation. In view of the subsequent financial costs if exploitation is stopped, it is almost impossible to perform interference tests in oil fields. On the other hand, interference testing is a standard in hydrology and, since a clever managing of water resources has become a priority, the need for advanced interpretation models has increased. Second, the dual-medium concept is not a straightforward approach. For a model to be useful, one must obviously be able to fit the data by adjusting its parameters. As shown later in this work, the dual-medium approach is very sensitive to parameters, which does not

315

help to fit the data unambiguously. The problem amplifies if automatic inversion is contemplated, and this is probably why possible previous attempts have not been continued. This paper proposes to develop a complete approach to the interpretation of interference pumping tests in both homogeneous and fractal dual media. The fractal assumption with simple scaling power laws for hydrodynamic parameters (see Section 2) is just a choice which enlarges the capabilities of data inversion while keeping a handy tool. This is not an universally agreed framework (see e.g. [38,41,43]) but discussing on relevance of the fractal approach compared with others is beyond the scope of the present work. A direct calculation based on a specific finite-volume code (see Section 2) is coupled with an automatic inversion procedure, which avoids (or lowers) subjectivity in interpreting field data. As mentioned above, because the dual-medium approach is very sensitive to its parameters, a specific inversion procedure has been set up. It is based on the derivation of the analytical sensitivities of the direct model to its parameters and on a conditioning of the optimizer using re-scaled sensitivities (see Section 3). Section 4 deals with a few synthetic and field cases to address the accuracy and relevance of the inversion method. The dual-medium approach and its inversion are shown to provide reliable results, which could improve in the future the understanding of fractured reservoirs as regards their macroscopic hydrodynamic properties. Before going into details, a few arguments must be given on the assumption prevailing over the type of flow simulated by the direct model. Interference pumping tests are not designed to point out 3D flow effects. Tests are mostly carried out in full-penetrating open boreholes and then the reservoir is stressed over its whole thickness and over areas of large horizontal extension. In addition, the control parameters of the test (i.e. the pumped flow rate Q(t) and the observed head h(r, t) at a distance r from the pumped well) are in fact measured as homogenized values along the boreholes. As a result, the observation of 3D effects is made difficult except at very early pumping times [1,4,5]. Consequently, a confined radial flow with cylindrical symmetry can be assumed. The confined flow assumption is justified in most cases, even in unconfined reservoirs, since the drawdowns are weak compared with the wetted thickness of the reservoir. In the following section, isotropy is also assumed, as is the case for most practical studies in which drawdown curves are interpreted one at a time, even when several observation points are available for the same pumped well. Note, however, that the model could take a possible anisotropy into account, which would not change the techniques proposed hereinafter to solve accurately the inverse problem. The main difficulty would then be the discussion of the results and their comparison for instance with those from other interpretation methods, e.g., the Cooper–Jacob analytical solution.

316

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

2. The discrete direct model This section details the direct simulator of radial flow with cylindrical symmetry in a dual medium corresponding to a pumping at a constant flow rate in a single well. In the following expressions, indexes f and m refer to the fracture and matrix continua, respectively. Basically, the couple of equations at the Darcy scale describing flow in a dual medium is written as ohf ¼ r  ðK f  rhf Þ þ aðhm  hf Þ þ qf ot ohm ¼ aðhf  hm Þ Ssm ot

Ssf

ð1Þ

with h [L] the hydraulic head, K [L T1] the hydraulic conductivity (which is basically a tensor in multidimensional flow), Ss [L1] the specific storage capacity, a [L1 T1] the exchange rate coefficient between fractures and the matrix, and qf [T1] a sink–source term per unit volume. Expression (1) assumes that flow in the matrix is negligible (term in Km dropped) and that the pumping well (via term qf) is located in the fracture continuum. It must be pointed out that Expression (1) has a form that is similar to a transient diffusive motion for hf coupled with a first-order kinetics involving hf and hm. More generally, several works in the ongoing literature have reported the endemic difficulties encountered in solving solute transport by coupling mechanisms with very different characteristic times, e.g., [11,21]. As shown later when dealing with inversion, these difficulties are comparable, or maybe greater, for flow in a dual medium owing to the concurrent effects of parameters, much contrasted sensitivities, weak convergence rate as the consequence of the well-known robustness to parameters of the diffusion equation, etc. This is probably why, to our knowledge, the previous works do not mention inversion of hydraulic tests in a dual medium with an automatic optimization. Expression (1) can be reformulated in radial coordinates. As justified in the introduction, an isotropic confined radial flow with cylindrical symmetry is assumed. With parameters stationary in time and using (r, t) as radial and time coordinates, respectively, Expression (1) can be rewritten as   ohf ðr; tÞ 1 o ohf ðr; tÞ Ssf ðrÞ ¼ rK f ðrÞ ot r or or þ aðrÞðhm ðr; tÞ  hf ðr; tÞÞ þ qf ð0; tÞ Ssm ðrÞ

ð2Þ

ohm ðr; tÞ ¼ aðrÞðhf ðr; tÞ  hm ðr; tÞÞ ot

The discrete numerical form is based on a finite volume formalism with a concentric meshing at constant space step Dr and time step Dt. For a mesh of index i, the mesh walls are indexed (i  1/2) and (i + 1/2), and correspond to the surfaces 2pe(i  1)Dr and 2peiDr, respectively, with e the reservoir thickness. The numerical scheme is cell-centered

and the state variables (hf, hm) are assumed to be averaged values over the mesh i and centered in (i  1/2)Dr: hnfi ¼ hhf ðði  1=2ÞDr; nDtÞi. Given these notations, the finite volume scheme is built by writing a classical mass (volumetric) balance equation for each mesh ohfi ¼ Qfiþ1=2  Qfi1=2 þ ai DV i ðhmi  hfi Þ þ Qpd1;i ot ohmi ¼ ai DV i ðhfi  hmi Þ Ssmi DV i ð3Þ ot Ssfi DV i

where DVi = 2pe(i  1/2)Dr2 is the volume of the mesh, Qfk [L3 T1] are the volumetric fluxes at the walls k of the mesh, Qp [L3 T1] is the pumped rate at the well (counted negative). In a heterogeneous medium, Ssfi ; Ssmi ; ai are supposed to be averaged values over the mesh i. These values for the heterogeneity model handled in this work will be discussed later on. Each flux Qfk is classically given by Darcy’s law   ohf Qfiþ1=2 ¼ 2peiDrK fiþ1=2 or iþ1=2   ð4Þ ohf Qfi1=2 ¼ 2peði  1ÞDrK fi1=2 or i1=2 Again, K fk are averaged values, also called intermesh conductivity. Several forms of intermesh values for radial flow have already been proposed, see, e.g., [18,19,37] but no model, be it empirically- or theoretically-based, is above the others as regards transient flow in the close vicinity of the well. Fortunately, interference pumping tests are carried out over relatively large lag distances between wells and simplifications can be done for the heterogeneity model used in this work (see hereafter). Using a first-order finite difference approximation of the spatial derivatives in (4), substituting (4) in (3) and assuming a full-implicit scheme in time, the following discrete equations are obtained: hfnþ1  hnfi ði  1=2ÞDr2 Ssfi i Dt     nþ1 nþ1 ¼ iK fiþ1=2 hfiþ1  hnþ1  h  ði  1ÞK fi1=2 hfnþ1 fi f i i1   Qpd 1;i nþ1 ð5Þ þ ði  1=2ÞDr2 ai hnþ1 þ m i  hf i 2pe   hnþ1  hnmi nþ1 ði  1=2ÞDr2 Ssmi mi ¼ ði  1=2ÞDr2 ai hfnþ1  h m i i Dt To complete (5), the initial conditions are set up to h0fi ¼ h0mi ¼ H 0 8i, i.e., a negligible steady-state flow and pressure equilibrium between the matrix and fractures. The infinite medium assumption is handled by setting hnfNmþ1 ¼ hnmNmþ1 ¼ H 0 . The number Nm of meshes discretizing the domain is selected large enough to prevent the depletion wave due to pumping and propagated during the simulation from bouncing against the limit and then influencing drawdown curves at observed wells. In mesh #1, where the pumped well is located, because of the symmetry of the problem, no specific boundary condition is required. The general finite-volume scheme depicted above

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

holds, the fluxes in (4) are expressed as: Qf1þ1=2 ¼ 2peK f1þ1=2 ðhf2  hf1 Þ; Qf1=2 ¼ 0, and the prescribed pumped flow rate Qp (negative) is added to the mass balance equation (see Expression (3)). Note that a varying space step could be easily implemented but its justification is questionable. On the one hand, it complicates the code and on the other hand, a local refining (e.g., near to the pumped well) instead of a fine meshing all over the domain does not save significant computation time since the problem is 1D. At this stage, a simple numerical tool is available to mimic radial flow in heterogeneous media, provided one is able to average the medium heterogeneity in terms of concentric zones. It must be reminded that the aim is to interpret and invert interference pumping tests in fractured media and then to be able (1) to compare the inverted hydrodynamic parameters with other observed data or with results from other approaches, and (2) to improve for instance the conditioning of more sophisticated reservoir models with dynamic data. Therefore, introducing any type of spatial heterogeneity may make inversion much more cumbersome and the resulting parameters hardly comparable with other approaches. Moreover, there is always a crucial need for references to partially or fully homogenized approaches [23,36,44]. This quite simple concept is still able to fix ideas or make them develop in view of the results. On the other hand, fractured rocks have often been referred to as statistically homogeneous media with fractal behavior [3,32,35]. These reasons have motivated the choice of limiting the inversion of interference pumping tests to two types of dual media: a homogeneous one and a fractal one (with hydrodynamic parameters ruled by scaling laws for both the matrix and fractures).

Af  hfnþ1 ¼ Df  hnf þ C  hmnþ1 þ B ¼ Dm  hnm þ C  hfnþ1 Am  hnþ1 m

317

ð7Þ

Af is a tri-diagonal matrix whereas matrices Am, Df, Dm and C are diagonal. B is a vector in which only the first and last terms are non-null. The first term stores the prescribed rate pumped in the well at mesh #1, and the last one the prescribed head H0 in mesh #Nm. Due to the presence of hm in equations for hf (and conversely), the dual linear system in (7) is coupled and can be solved using the so called fixed point method: Af  hfnþ1;kþ1 ¼    C  hmnþ1;k then Am  hmnþ1;kþ1 ¼    C nþ1;kþ1 with k the convergence iteration index and an initial hf condition hmnþ1;0 ¼ hnm . In the present case, the system (7) is relatively simple and another method can be proposed. From the  innþ1(7) it can be written: second setn of equations 1 hmnþ1 ¼ A1  D þ A  C  hf . Substitution in  h m m m m 1 the first equation in (7) yields: A  C  A f m   n 1 n  h C  hnþ1 ¼ D  h þ C  A  D þ B which does not f m f f m m  involve hmnþ1 . The matrix Af  C  A1 m  C is tri-diagonal and the system is solved by a direct (opposed to iterative) nþ1 method using the Thomas 1  hf nþ1, the 1 algorithm.  n Knowing nþ1 resolution of hm ¼ Am  Dm  hm þ Am  C  hf is straightforward. This method does not require any convergence iterations to ensure the coupling in (7) and the computation is very rapid, for instance, about 3 s on a common PC for 500 meshes and 1000 time steps. Note that this onestep calculation merely lumps terms on the diagonal of the matrix Af. With the direct Thomas algorithm, this does not increase the risk to get an ill-posed linear system even with huge contrasts between matrix and fracture hydraulic properties. 2.2. Discrete scheme for fractal dual media

2.1. Discrete scheme for homogeneous dual media Stating that the medium is homogeneous with four constant parameters Kf, Ssf, Ssm and a simplifies Expression (5). The latter can be rewritten for each couple of equation ‘‘attached’’ to mesh i as follows:    nþ1 2 Ssf  ½ði  1ÞK f hfi1 þ ði  1=2ÞDr þa Dt  nþ1 þð2i  1ÞK f hfi  ½iK f hnþ1 fiþ1    Ss f ¼ ði  1=2ÞDr2 hnfi þ ði  1=2ÞaDr2 hnþ1 mi Dt ð6Þ Qpd1;i þ 2pe    Ssm þ a hnþ1 ði  1=2ÞDr2 mi Dt    Ssm n ¼ ði  1=2ÞDr2 hmi þ ði  1=2ÞaDr2 hnþ1 fi Dt The notation of the dual linear system of equations in (6) can be simplified using the matrix form

This model is developed using power laws in space for the hydrodynamic parameters of both the matrix and fractures. The parameters decrease with the lag distance r between the pumped well and the observed ones. This type of scaling laws, which inherits from the percolation theory and theoretical works on fractals [31,39], has already been used for single media, e.g., [2,14,17,25]. In the present case, dealing with a dual medium, the fractal behavior can be interpreted as nesting two structural entities of different scales. At the larger scale, the fracture medium encloses features of wide extension draining the reservoir. At the smaller scale, the matrix medium is also fractured but with smaller features, more densely distributed over space. As a result, the matrix is unable to connect the reservoir over large areas (which is ensured by the fractures) but conceals most of the fluid content. This type of behavior has often been reported in petroleum engineering applied to calcareous fractured reservoirs and, to a lesser extent, in hydrology [9]. The fractal dual medium is depicted by 8 parameters according to the following definitions:

318

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

K f ðrÞ ¼ K f0 ra ;

Ssf ðrÞ ¼ Ssf0 rb ;

Ssm ðrÞ ¼ Ssm0 rc ;

¼

aðrÞ ¼ a0 rd

ð8Þ

All power-law exponents a, b, c, d in (8) are considered to be positive. A specific storage capacity has the general following expression: Ss / /cqg, with / [–] the porosity of the medium, c [M1 L T2] the compressibility of the medium, q [M L3] the fluid bulk density and g [L T2] the gravity acceleration. In fractal media, the porosity follows scaling laws in the form / / rDd with D the mass fractal dimension of the medium and d its Euclidean dimension (here set to d = 2). Common values drawn from the literature, e.g. [39], give D between 1.5 and 1.95 for both synthetic and natural fractal objects, and given the above relationships, one may expect exponents b and c to span the range [0.05–0.5] for Ssf(r) and Ssm(r), respectively. On the other hand, the hydraulic diffusion is supposed to evolve as DH / rh with h the so-called dimension of the random walker [31,35]. In two-dimensional fractal media, h may vary between 0.3 and 0.8. Thus, with DH = Kf/Ss, one gets Kf / rDdh and an exponent a for Kf(r) in the range [0.4– 1.3]. Basically, the exchange rate coefficient between fractures and the matrix should have the significance of a diffusion coefficient and thus be ruled by rh. In practice, it buffers fluxes between the water content of fractures (Ssf) and that of the matrix (Ssm). In the following, the exchange rate coefficient a is therefore assigned with a scaling law similar to that of a storage capacity, which makes d in expression of a(r) vary in the range [0.05–0.5]. The discrete equations in Expression (5) hold for the fractal dual medium, and therefore their general matrix notation in (7) as well as their numerical resolution principle remain valid. However, given the medium heterogeneity, the scalar components enclosed in matrices and vectors {Af, Am, . . . , B} have to be re-expressed. First, take the intermesh conductivity (e.g., K fiþ1=2 ) between each concentric mesh of the model. Its theoretical value for radial flow should be as close as possible to that deduced from Thiem’s expression [37]. This expression is based on the flow rate conservation between two adjacent points located at distances ri and rj from the well: Q = 2peKij(hj  hi)/ (ln rj  ln ri), with Q the prescribed flow rate pumped in r = 0, e the aquifer thickness, hk the hydraulic head in rk (k = i, j) for a steady-state flow and Kij the equivalent hydraulic conductivity between locations ri and rj. For a regular cylindrical and concentric meshing, a few algebraic manipulations (Appendix A) show that the value Kij between neighboring meshes tends rapidly toward the harmonic mean when the distance from the pumping well increases. For the fractal hydraulic conductivity of Expression (8), the harmonic mean hK fiþ1=2 iH at the interface i + 1/2 located in iDr between concentric meshes i and i + 1 is written as 

hK fiþ1=2 iH

1

¼

1 Dr

Z

ðiþ1=2ÞDr

ði1=2ÞDr



K f0 ua

1

du ! hK fiþ1=2 iH

ða þ 1ÞK f0 Dra ði þ 1=2Þ

aþ1

 ði  1=2Þ

aþ1

ð9Þ

If i is large compared to 1/2 (say, i above 25), the first-order development of the denominator in (9) simplifies and yields: a hK fiþ1=2 iH  K f0 ðiDrÞ which corresponds to the local value of Kf(r) at the interface iDr between the meshes i and i + 1. In other words, for interference testing in which observed wells are assumed to be far enough from the pumped well to avoid high fluid velocities and for a refined spatial discretization (say, more than 30 meshes between the pumped well and the nearest observed one), the local value of Kf(r) at the interface can be taken without loss of accuracy as the intermesh conductivity between adjacent meshes. Second, parameters Ssf, Ssm, and a write in the form: s(r) = s0rk and their value in each mesh should be calculated as their arithmetic mean within the mesh. Storage capacity is linearly linked to the medium porosity (see above) which is an additive property and justifies the use of an arithmetic mean (a volume V made of V1 and V2 with porosity /1 and /2, respectively, as a porosity / = (/1V1 + /2V2)/V) Z iDr 1 hsi iA ¼ s0 uk du ! hsi iA Dr ði1ÞDr   s0 Drk i1k  ði  1Þ1k ð10Þ ¼ 1k A truncated first-order development of i1k = (i  1/2 + 1/2)1k and (i  1)1k = (i  1/2  1/2)1k in e with e = 1/2(i  1/2) allows to simplify the arithmetic mean k hsi iA  s0 ðði  1=2ÞDrÞ , i.e. the local value of s(r) centered on mesh i where the hydraulic head hi is calculated. Again, a few computations did not show any loss of accuracy by using these simplifications. Even for indexes i lower than 10, the relative errors on heads are less than 0.5% which is often far below uncertainty on head measurements from field data. Moreover, these simplifications facilitate the derivation of the analytical sensitivities of the model to parameters and, as shown below, they also allow simplifying the descent direction when seeking the optimal solution to the inverse problem. Introducing these simplifications of mesh and intermesh parameters in Expression (5) yields the equations of mesh i for a fractal dual medium subjected to convergent radial flow h i 1a  ði  1Þ K f0 Dra hfnþ1 i1   Ss f0 ði  1=2Þb Drb þ ði  1=2ÞDr2 Dt     d 1a d 1a a þ a0 ði  1=2Þ Dr þ i þ ði  1Þ K f0 Dr   i1a K f0 Dra hnþ1  hfnþ1 fiþ1 i   Ssf ¼ ði  1=2Þ1b 0 Dr2b hnfi Dt h i Qpd1;i 1d þ ði  1=2Þ a0 Dr2d hmnþ1 þ i 2pe

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334



 Ssm0 c ði  1=2Þ Drc Dt i d þ a0 ði  1=2Þ Drd hnþ1 mi   1c Ssm0 2c Dr ¼ ði  1=2Þ hnmi Dt h i 1d þ ði  1=2Þ a0 Dr2d hnþ1 fi

ði  1=2ÞDr2

ð11Þ

uses an exact derivation of the gradient vector $Fk and an approximation for the Hessian matrix $2Fk, yielding  1    J Tnk  nk dk ¼  J Tnk  J nk ð15Þ with J nk the n · p Jacobian matrix of errors ¼ onki =osj ; i ¼ 1; . . . ; n; j ¼ 1; . . . ; p . The approximation of the Hessian   matrix is the result of dropping terms in

3. Inversion procedure For the sake of simplicity, only drawdown data from a single observed well are discussed hereinafter. However, the following development remains valid for multiple observed points (but a single pumped well) inverted with the same set of parameters. As is classical in dual media, it is assumed that the observed heads in a well are those from the fracture medium. The possible questioning of this assumption n is beyond the o scope of this paper. Thus, let us O denote hfi ; i ¼ 1; . . . ; n the set of observed heads in the fracture medium at a location r during a pumping test at the flow rate Q. The flow can be mimicked by the discrete model depicted in Section 2 and, for a given set of parameters s = {sj, j = 1, . . . , p}, nerrors between simulation and o observation write nðsÞ ¼ ni ¼ ðhCfi ðsÞ  hO Þ; i ¼ 1; . . . ; n fi with hCf the vector of simulated heads. A classical way to identify (invert) the best set s is to minimize an objective function F(s) n n  2 1 1X 1X n2i ðsÞ ¼ hCfi ðsÞ  hO F ðsÞ ¼ nT ðsÞ  nðsÞ ¼ fi 2 2 i¼1 2 i¼1 ð12Þ with T the transposition operator. This method, also called a simple least-square minimization, requires an optimization algorithm to seek s. For small sizes p of vector s, say less than 20 parameters, several works in the literature have shown the iterative Gauss–Newton optimizer to be efficient [12,16,24, among others]. This algorithm is based on a second-order approximation of the objective function. Using k as the iteration index when seeking s, it can be written as T 1 T F ðskþ1 Þ ¼ F ðsk Þ þ dk  rF ðsk Þ þ dk  r2 F ðsk Þ  dk 2

319

ð13Þ

with dk = sk+1  sk, $F(sk) the gradient vector (of size p) of the objective function ={oF/osj, j = 1, . . . , p}, $2F(sk) the Hessian matrix (of size p · p) ={o2F/osiosj, i, j = 1, . . . , p}. If vector sk+1 is the solution to the minimization of F, thus gradient $F(sk + 1) = 0. Calculating this gradient from (13) and canceling it yields the expression for calculating sk+1 as a function of sk

1 dk ¼ skþ1  sk ¼  r2 F ðsk Þ  rF ðsk Þ ð14Þ Expression (14) is a generic form also called the Newton method. More specifically, the Gauss–Newton algorithm

oJ Tnk =osj  nk from the complete derivation, assuming that n is small near to the optimum. Note that the so-called Levenberg–Marquardt optimizer, e.g., in [12], can be derived in mixing the Gauss– Newton algorithm and a Steepest-Descent one. Limiting the development in (13) to the first-order, the quickest way to decrease F(sk+1) compared with F(sk) is to choose dk =  $F(sk), namely a Steepest-Descent along the gradient of the objective function. The Levenberg–Marquardt method simply replaces the matrix in (15) by  1  J Tnk  J nk þ gI which allows to tend to a Steepest-Descent method for large values of g. This algorithm is often used to accelerate convergence when the solution sk is far from the optimum or when the Jacobian matrix is not well calculated. As shown below, solving the inverse problem in the present case needs specific adaptations and the Levenberg–Marquardt alternative becomes of weak interest. Since oni =osj ¼ ohCfi =osj , the Jacobian matrix stores the sensitivities of the model to parameters and the convergence of the optimization problem is all the more rapid that sensitivities are calculated accurately [29,40]. The simplest way to calculate the terms of the Jacobian matrix is the ‘‘perturbation’’ method, which consists in a term-by-term approximation     ni skj þ dsj  ni skj onki  ; dsj oskj i ¼ 1; . . . ; n; j ¼ 1; . . . ; p

ð16Þ

with dsj a small perturbation added to parameter sj. This is probably the most commonly used approach because no other development than the direct model is required. The latter is merely used p times around the solution sk. However, the results are sometimes disappointing because the perturbation ds is hard to set up: too large, the finite difference approximation in (16) is rough and the convergence of the problem is weak; too small, numerical discrepancies such as round-off errors, ‘‘division by zero’’, etc. are often experienced. A second approach consists in finding the explicit equations of the Jacobian terms. This method, also called the ‘‘analytical’’ method [29], can be used for models with a fairly simple parameterization. 3.1. Analytical sensitivities of the homogeneous dual medium The matrix formulation of the direct problem in a homogeneous dual medium (Expression (7) ) can be rewritten as

320

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

Af  hfnþ1 ¼ Wf ;

Am  hnþ1 ¼ Wm m

ð17Þ

where Wf and Wm are the right-hand-side vectors of the calculation of hydraulic heads in fractures and the matrix, respectively. After differentiation of (17) with respect to any parameter sj, one gets Af 

ohfnþ1 oWf oAf nþ1 ¼   hf ; osj osj osj

Am 

ohnþ1 oWm oAm nþ1 m ¼   hm osj osj osj

ð18Þ

Thus, provided that hnþ1 and hnþ1 have been calculated and f m that the analytic derivatives of Af, Wf, Am, Wm are availnþ1 able, sensitivities ohnþ1 f =osj and ohm =osj are calculable with a linear system similar to the direct problem (17) but with a modified right-hand-side term. These derivatives of Af, Wf, Am, Wm for a homogeneous dual medium are given in Appendix B. In terms of computation costs, the ‘‘analytical’’ method requires as much calculations as the ‘‘perturbation’’ method, but the result is as accurate as the direct computation of the state variables hf and hm. As a by-product, the plot of analytical sensitivities versus time gives information on the behavior of the drawdown at one location as regard the four parameters of the dual medium. In the following example, a pumping test at constant flow rate Q = 50 m3 h1 is simulated. The parameters of the dual medium are: Kf = 105 m s1, Ssf = 106 m1, Ssm=3 · 105 m1, a=5 · 1010 m1 s1 and correspond to average values often encountered in fractured oil reservoirs, e.g., [22]. The sensitivities versus time are calculated at location r = 100 m from the pumped well, and their behavior is consistent with the ‘‘empirical’’ but commonly held rules (Fig. 1). The sensitivity of heads in both the matrix and fractures to the fracture hydraulic conductivity Kf increases with time and evolves following almost the same shape as head drawdowns. The head in fractures is sensitive to the fracture specific storage capacity Ssf at early pumping times, to the matrix specific storage capacity Ssm at long times and to the matrix-fracture exchange rate a at intermediate times. On the other hand, the head in the matrix is weakly sensitive to Ssf, whereas its sensitivity to Ssm increases regularly with time and its sensitivity to a is maximal at intermediate times. Note that for a, the sensitivities of hf and hm are opposite in sign, i.e., an exchange from the matrix toward fractures makes the pressure increase in fractures and decrease in the matrix. In the end, the analytical sensitivities are consistent with the assertion stating that during a pumping test, the conductive fracture medium is requested first. Considering its small storage capacity, draining is rapid and the system is fed progressively by the matrix with an important exchange between the matrix and fractures. For long pumping times, the system depends on the matrix water content but is drained toward the pumped well by the fractures. The macroscopic hydraulic diffusion coefficient is in Kf/Ssm (which can be verified by computations, not provided here).

At least as important as their evolution in time are the scalar values of sensitivities: in the range 104–105 for Kf, Ssf, Ssm and about 109 for a. This very significant difference makes Expression (15) numerically ill-posed; the

 matrix Pn J Tn  J n storing terms m¼1 ðonm =osi Þ onm =osj ; i; j ¼ 1; . . . ; p is poorly conditioned. Take for instance a high sensitivity to parameter sk as compared to the others. The row and the column k of the matrix are filled with huge values, which results in a singular system in Expression (15). One could try to seek the parameter ln(s) instead of s since on/o ln(s) = son/os, which would re-scale the sensitivities. The problem with this method is that re-scaling is dependent on ‘‘local’’ values of s at each convergence iteration. These values may strongly vary during the optimization procedure and also become very different between parameters. They may also be very different from one problem to the other, which makes it difficult to express simple rules to improve conditioning. In the end, J Tn  J n is again poorly conditioned and out of control by the user. Another track (discussed below) has been followed in this work to propose a conditioning of sensitivities able to solve the illposed problem in Expression (15). To our knowledge, nothing in the literature has been reported on analytical sensitivities and their huge contrast for flow in dual media. This lack of data is probably one of the reasons why automatic inversion of interference pumping tests in these media has not been addressed yet. 3.2. Analytical sensitivities of the fractal dual medium As for the homogeneous assumption, the calculation principle of analytical sensitivities remains valid with a discrete direct model in the form A Æ h = W yielding A Æ oh/os = oW/os  oA/os Æ h. In the present case, remember that the parameters follow the expression s(r) = s0rk and therefore, eight terms have to be identified (see Expression (8)). The details of derivations allowing calculating the analytical sensitivities to parameters are given in Appendix C and an example of these sensitivities and their evolution in time is reported in Fig. 2. The settings of the computations are as follows: pumped flow rate Q = 50 m3 h1, observation distance r = 100 m, Kf0 = 103 m1+a s1, Ssf0 = 105 mb  1, Ssm0 = 104 mc  1, a0 = 5 · 109 md  1 s1, a = 0.8, b = c = d = 0.1. In the range r = 25–200 m, these parameters give s(r) values that are consistent with those previously tested for the homogeneous dual medium (Fig. 1) and the power-law exponents are compatible with theoretical values from 2D regular percolation networks [31,35]. The comparison between Figs. 1 and 2 show that sensitivities to parameters s0 (in s(r) = s0rk) of the fractal dual medium are comparable to their equivalent s of the homogeneous dual medium. The head hf shows a regularly increasing sensitivity to K f0 with time and a strong sensitivity to Ssf0 at short pumping times and both hf and hm show an increasing sensitivity to Ssm0 and a maximal absolute value of sensitivity to a0 at intermediate pumping times. Interestingly, sensitivities to power-law exponents k are

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

321

5

h (m) 5

3

x 10

Fracture 4

dhf /dKf

2.5

Matrix

dhm /dKf

2 3 1.5 2 1 1

0.5

0 2 10

4

6

10 time (s)

10

0 2 10

x 10

5

x 10

dhf /dSsf

dhf /dSsm 4

4

dhm /dSsm

3

3

2

2

1

1

0 2 10

6

10

5

4

5

4

10 time (s)

4

6

10 time (s)

10

dhm /dSsf

0 2 10

4

10 time (s)

6

10

9

1.5

x 10

1 0.5 0 -0.5

dhf /dalpha dhm /dalpha

-1 2 10

4

6

10

10

time (s) Fig. 1. Example of evolution with time of drawdowns and sensitivity to parameters at a given distance (100 m) from the pumped well for the homogeneous dual medium. hf, hm = hydraulic head in fractures and matrix, respectively, Kf = hydraulic conductivity of fractures, Ssf, Ssm = specific storage capacity of fractures and matrix, respectively, a = exchange rate between fractures and matrix.

comparable to those of their attached parameter s0 given a negative re-scaling factor (see Fig. 2). The explanation is as follows. Take again the formalism A Æ h = W ! A Æ oh/os = oW/os  oA/os Æ h. A is tri-diagonal for the fracture equations and diagonal for the matrix equations. W is a composite vector made of a linear combination between vectors and scalar products of (tri)diagonal matrices and vectors (see Section 2 and Appendix C). Thus, the calculation of sensitivities to any parameter s0 can be expanded as

j¼iþ1 X j¼i1

aij

j¼iþ1 ohj owi X oaij ¼  hj os0 os0 j¼i1 os0

ð19Þ

and a similar form applies to ohj/ok. A close look at coefficients aij and wj shows that s(r) appears in these coefficients as a linear combination of sk = s0k1kDrk with k = i  1, i  1/2, i. This is a consequence of the simplifications evoked in Section 2 about averages for mesh and intermesh parameters. The ith equation of sensitivity to s0 can be rewritten as

322

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

h (m) 3

1500

15000

2

dhf /dSsf0

dhf /dKf0

Fracture

1000

Matrix

dhm /dKf0

dhm /dSsf0

10000

500 1

5000

0 0 102

104

106

-500 102

104

time (s)

106

0 102

time (s)

104

106

time (s)

7

10000 8000

6 dhf /dSsm0

x 10

2

4

0

6000

2

-2

4000

0

-4

2000

-2

dhm /dSsm0

dhf /dalpha0

-6

dhm /dalpha0

0 102

104

106

time (s) 0 -0.2

-4 102

dhf /da dhm /da

104

106

time (s)

-8 102

0

0.5

-2

0

104

106

time (s)

-0.4 -4

-0.6 -0.8 102

104

-0.5

dhf /db

dhf /dc

dhm /db

dhm /dc

106

-6 102

dhf /dd dhm /dd

104

time (s)

106

-1 102

time (s)

104

106

time (s)

Fig. 2. Example of evolution with time of drawdowns and sensitivity to parameters at a given distance (100 m) from the pumped well for the fractal dual medium. hf, hm = hydraulic head in fractures and matrix, respectively, Kf(r) = Kf0ra = hydraulic conductivity of fractures, Ssf = Ssf 0 rb, Ssm = Ssm0rc = specific storage capacity of fractures and matrix, respectively, a = a0rd = exchange rate between fractures and matrix.

( ) ( ) j¼iþ1 X osk X X osk ohj aij ¼  hj ð20Þ os0 os0 os0 k k j¼i1 j¼i1 wi aij P with k a the notation corresponding to the linear combination of s(r) enclosed in the scalar a. Knowing that an expression similar to (20) holds for ohj/ok and j¼iþ1 X

osk ¼ k 1k Drk ; os0 osk osk osk ¼ s0 lnðkDrÞk 1k Drk ) ¼ s0 lnðkDrÞ ok ok os0

j¼iþ1 X j¼i1

ohj  s0 ln ðði  1=2ÞDrÞ ok 2( 3 ) ( ) j¼iþ1 X osk X X osk 4  hj 5 os0 os0 k k j¼i1

aij

wi

The above expression can be reduced to   oh oW oA A  S   h ok os0 os0 ð21Þ

thus, the ith equation of the sensitivity to k will write ( ) j¼iþ1 X ohj X osk ¼ aij s0 lnðkDrÞ ok os0 k j¼i1 wi ( ) j¼iþ1 X X osk þ s0 lnðkDrÞ hj ð22Þ os0 k j¼i1 aij

Assuming that radii rk = kDr (k = i  1, i  1/2, i) are large enough, so that ln(kDr) is almost constant whatever k, and thus set up to ln((i  1/2)Dr), one gets from (22)

ð23Þ

aij

ð24Þ

with S a diagonal matrix storing the terms s0 ln ((i  1/2)Dr), i = 1,. . . , Nm (Nm is the total number of meshes discretizing the domain). Thus, when sensitivity to s0 has been calculated using (20), sensitivity to k is obtained without further computations by stating that ohi ohi  s0 ln ðði  1=2ÞDrÞ ; ok os0

i ¼ 1; . . . ; Nm

ð25Þ

To be complete, the above reasoning should account for the fact that the handled problem is transient, and some sensitivity calculations may write for instance as (see Appendix C)

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

ohnþ1 oW oA nþ1 ohn  h þD ¼ ð26Þ ok ok ok ok Expression (26) shows that sensitivities ohn/ok at time nDt should be known before computing those at time (n + 1)Dt. From the above development of Expressions (19)–(24), it is stated that oW/ok =  S Æ oW/os0, oA/ ok =  S Æ oA/os0. The linear system (26) for a transient problem is therefore rewritten as

A

A

ohnþ1 oW oA nþ1 ohn ¼ S  þS h þD os0 os0 ok ok

ð27Þ

The following recursive reasoning can be suggested. Assuming that ohn/ok =  S Æ ohn/os0, Expression (27) is rewritten as   ohnþ1 oW oA nþ1 ohn ¼ S   h þD A ð28Þ os0 os0 ok os0 Knowing that sensitivity to s0 obeys the following expression: A

ohnþ1 oW oA nþ1 ohn ¼  h þD os0 os0 os0 os0

ð29Þ

then comparing (28) with (29) obviously yields ohn+1/ok = S Æ ohn+1/os0, which completes the recursive reasoning. Errors in sensitivities to power-law exponents between an exact calculation (see Appendix C) and a straightforward approximation discussed above (Expression (25)) remain to be checked. Remember that the basic key to this approximation is the simplification of averaged mesh and intermesh parameters in the discrete equations of the direct problem. As stated earlier, these simplifications are valid

323

for a fine discretization in space and/or at locations far enough from the pumped well (see Section 2). Since the approximation of sensitivities to power-law exponents is also based on the radii iDr where they are calculated indexed by large i values (see above), accuracy cannot be expected close to the pumped well and with a coarse meshing. The example in Fig. 3 comparing exact sensitivities with their approximation is based on the same settings as in Fig. 2. The mesh step Dr is 0.25 m and, for illustration purposes, only the sensitivities of the heads in fractures with respect to exponents a and d (exponents of the fracture hydraulic conductivity and the exchange rate between the matrix and fractures, respectively) are contemplated. As expected, for a given time (105 s), errors between exact sensitivities and their approximations decrease with the lag distance between pumped and observed points. They tend to infinity in r = 0 (i = 0) but rapidly decrease with the index i of the observed mesh. Beyond i = 200, relative errors are less than 10%, which is very small as compared with errors stemming from calculations using a perturbation method. In other words and from a practical standpoint, there is no efficiency loss in seeking an inverse solution when the exact sensitivities to power-law exponents are replaced with their approximations. Several reference inverse problems were computed with both exact and approximated sensitivities (not reported here), and no difference whether in terms of convergence rate or in terms of distance between the sought solution and the reference one was observed. In addition, it can be verified that errors close to the pumped well do not propagate and do not amplify with iterations over time, which could result in

dhf /da at t = 100 000 s

dhf /dd at t = 100 000 s

6

0.15

4

0.1

2 0

0.05

-2 -4

0 -0.05

-6 -8

-0.1 0

200

400

600

800

1000

0

dhf /da in mesh #300 (75 m)

400

600

800

1000

dhf /dd in mesh #300 (75 m) 0.2

1

0

0

-0.2

-1

-0.4

-2

-0.6

-3

-0.8

-4 -5 100

200

Index of meshes

Index of meshes

-1 1000

10000

Time (s)

100000

-1.2 100

1000

10000

100000

Time (s)

Fig. 3. Comparison between exact analytical sensitivities (solid lines) to power-law exponents and approximated ones (dots) in a fractal dual medium. Sensitivity of fracture head hf to exponent of fracture hydraulic conductivity (a) and exchange rate between the matrix and fractures (d). Up: sensitivities at a given time versus distance from the pumped well. Down: sensitivities at a given location versus time.

324

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

biased or numerically unstable calculations. For example (Fig. 3), at a distance of 300 meshes from the pumped well, i.e., 75 m, which is a common distance for interference tests, there is neither significant increase nor oscillations in errors over 200 time steps (i.e. times between 103 and 105 s, a common range for interference tests). To sum up, regarding these calculations of analytical sensitivities to parameters of the fractal dual medium, the computation time for a fractal approach with 8 parameters can be divided by a factor two without efficiency loss and then become equivalent to costs of the homogeneous approach with 4 parameters. This is a non-negligible advantage when solving inverse problems, since seeking an optimal solution may require many convergence iterations. 3.3. Conditioning of the Jacobian matrix

1

system dk ¼ skþ1  sk ¼  J Tnk  J nk  J Tnk  nk becomes singular. The purpose here is therefore to condition the Jacobian

matrix of errors J nk ¼ onki =osj ; i ¼ 1; . . . ; n; j ¼ 1; . . . ; p by multiplying sensitivities onki =osp by a scalar c to make them of the same order of magnitude as onki =osj ; j 6¼ p. In the following applications (see Section 4) sensitivities are re-scaled to values comparable with sensitivities to the fracture hydraulic conductivity. Because pinpoint accuracy is not required in this re-scaling, a single scalar coefficient c is chosen for all errors nki ; i ¼ 1; . . . ; n (but this coefficient can be chosen different at each convergence iteration of index k, and of course, different for sensitivities to sj and to sl, j 5 l, in the case several sensitivities are rescaled). Introducing these re-scaled sensitivities to sp in the Gauss–Newton algorithm gives

m

i; j ¼ 1; . . . ; p  1 0 P on 1 nm osmj B m C k T J nk  n ¼ @ P on A c nm osmp m

m

j ¼ 1; . . . ; p  1

p1 X

aij d j þ caip d p ¼ bi

j¼1

"

c

p1 X

ði ¼ 1; . . . ; p  1Þ

# apj d j

þ

capp d p

ð31Þ ¼ cbp

j¼1

On the other hand, the linear system without modification of sensitivities to sp would be p1 X

aij d j þ aip d p ¼ bi

ði ¼ 1; . . . ; p  1Þ

j¼1 p1 X

ð32Þ apj d j þ app d p ¼ bp

j¼1

As suggested earlier, the huge differences between values of sensitivities to different parameters may generate a numerically ill-posed problem when seeking an inverse solution with a Gauss–Newton algorithm. For the sake of simplicity in equation writing, it will be supposed that within the set of parameters s = {sj, j = 1,. . . , p}, only sensitivities oh/osp are very different from the others. The reasoning developed below remains valid, however, in the case of several sensitivities oh/osk, k among 1,. . . , p of very different orders of magnitude. If sensitivities oh/osp = on/osp are very strong (or very weak), remember    thatthe linear

 T 1  T  k d ¼  J nk  J nk  J nk  nk with 0P 1 P onm onm onm onm c os os os os i j i p  T  B m C m C J nk  J nk ¼ B @ P onm onm 2 P onm 2 A c c osp osj osp

Using a general notation A Æ d* = b, the linear system to be solved can be rewritten as

The comparison between (31) and (32) shows that multiplying sensitivities to parameter sp by a factor c is equivalent to the re-scaling dp ! cdp of the pth component of the descent direction d. More exactly, this modification is not strictly limited to dp. Stretching out (or compressing) one component of the descent vector influences other components since all are calculated by solving a linear system of p equations with p unknowns. It remains, however, that Expression (31) has a numerical solution and allows to invert the problem despite the huge contrast between sensitivities to parameters. All is needed is to grasp what is identified using Expression (31). One can write: con/osp  on/or; on/or = (on/osp)(osp/or) and identification between both expressions gives r = sp/c. Therefore, once Expression (31) has been solved, the parameters are updated from convergence iterations k to k + 1 as follows: k

skþ1 ¼ skj þ fd j j skþ1 p c

¼

skp c

ðj ¼ 1; . . . ; p  1Þ ð33Þ

k

þ fd p

k

with f a scalar limiting the descent along direction d . Several algorithms have been proposed to choose this step of descent f, e.g., [13]. In the examples proposed hereinafter, f was merely adjusted for Kf or K f0 (in the homogeneous and fractal cases, respectively) to be moved of less than 15% between two successive convergence iterations. We found this choice to be a good compromise between rapidity of convergence and accuracy in seeking optimal parameters. 4. Examples of inversion

ð30Þ

4.1. Synthetic cases Two synthetic cases are provided below for illustration purposes. Two reference calculations are first carried out using the direct model of Section 2 with parameters set up by the user. Then, drawdown curves from these reference calculations serve as data to be fitted by inversion.

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

Since the true solution is known, the aim is to check the capabilities of the inverse procedure to seek the optimal solution. Three major points can be addressed: (1) the capabilities to find the optimal solution rapidly; (2) the accuracy of the optimizer in terms of distance between the optimal sought solution and the reference solution; (3) the stability around the optimal solution. The first point is straightforward: the more rapid the optimizer, the smaller the number of convergence iterations required. The second point can be discussed while keeping in mind that accuracy must also take into account the distance between the exact solution and the initial one from which the search was launched. If this distance is small, the optimized solution should be very close to the reference one. If the initial distance is very large, one may expect much difficulty in seeking the optimal solution and consequently less accuracy. The third point can be addressed pragmatically by letting the optimizer continue its search even after an optimal solution has been found. If there are only slight oscillations around the optimum, the solution is stable. A more rigorous approach is P to analyze the n T matrix m¼1 ðonm =osi Þ

J n  J n . Since it stores terms onm =osj ; i; j ¼ 1; . . . ; p (Section 3, Expressions (15) and (30)), it can be viewed as a measurement of the crosscovariances between sensitivities to parameters. If the matrix is flat, all terms are of the same order, there is no preferential component to follow along the descent direction and the solution is stable. On the other hand, if the optimizer is able to experience new preferential directions with the hope of improving its search, most often the problem becomes unstable, with large oscillations of both parameter values and objective function. In the present case, the re-scaling of sensitivities flattens out the matrix J Tn  J n and vector J Tn  n, even at the first iteration. Thus, the above criteria do not apply strictly, but one may check that the flattening is preserved throughout the convergence iterations. Moreover, the mean values of components in J Tn  J n and J Tn  n inform on stability. At the first convergence iteration, the ratio hnon/osji/hon2/osi osji should be high, yielding high values of descent components; around the optimal solution, this ratio should be small and remain as such to ensure small variations in the descent direction and therefore stability. The first test deals with a homogeneous dual medium and the reference problem is designed as follows: Q = 50 m3 h1, Kf = 105 m s1, Ssf = 106 m1; Ssm = 3 · 105 m1, a = 5 · 1010 m1 s1, 500 radial meshes of step Dr = 4 m, a time step Dt = 500 s, and a single observed well to be fitted located at r = 150 m. Fig. 4 shows the two drawdown curves (observed and simulated) for the best fit and the evolution with the iteration index of the objective function and parameter values. In this run, the initial solution from which the search was launched is (103, 104, 3 · 103, 5 · 108) for Kf, Ssf, Ssm and a, respectively. Only sensitivities to a were re-scaled, on average by a multiplying factor 104, to get a numerically well-posed problem for the descent direction (Section 3, Expressions (30)–(33)).

325

This factor makes sensitivities to a of about the same values as sensitivities to the hydraulic conductivity Kf. An indicator of the re-conditioning of J Tn  J n is given by its eigenvalues, they span over a narrower range when conditioning improves. Without re-scaling the sensitivities to a the eigenvalues are (2.34 · 1011, 6.22 · 103, 3.39 · 103, 2.56), with re-scaled sensitivities one gets (1.11 · 104, 4.43 · 103, 1.87 · 103, 2.55). Fig. 4 shows that re-scaling these sensitivities makes all the sought parameters evolve evenly with the iteration index. With an initial solution chosen deliberately above the reference for all parameters, the latter go down in almost the same manner. This reduces by a factor of two the number of convergence iterations by avoiding inefficient ‘‘wanderings’’ of the optimizer within the parameter space. Otherwise, one parameter is sometimes increased, which forces another one to decrease, and because nonmodified sensitivities can be very different, the optimizer may require several additional iterations to put again one parameter on the right descent direction. It is also observed that stability around the optimal solution is good. Oscillations of both the objective function and parameters are very weak once the optimal solution has been found at iteration #32. This is supported by the flatness of the matrix J Tn  J n and vector J Tn  n, which is kept throughout convergence iterations. The ratio hnon/osji/ hon2/osi osji of 4 · 102/2 · 103 at the first iteration discussed above, goes down to 8 · 105/3 · 1011 at iteration #32 and beyond. The optimal sought solution is stable. In terms of accuracy, the optimal solution reached at iteration #32 is (1.01 · 105, 1.01 · 106, 3.02 · 105, 5.05 · 1010) for Kf, Ssf, Ssm and a, respectively, which represents a relative error of 1%. When the starting point is moved away from the reference solution, accuracy remains unchanged, only the number of convergence iterations is increased. For instance, starting from (102, 103, 3 · 102, 5 · 107), i.e. three orders of magnitude above the reference for all parameters, yields an optimum after 50 iterations and relative errors of less than 1% on each parameter. The second test deals with the fractal dual medium. Before discussing the results, a few details must be pointed out. Let us come back to the calculation of analytical sensitivities to power-law exponents. As stated earlier, provided that the spatial discretization is fine enough and that drawdowns are not calculated too close from the pumped well, analytical sensitivities to k for a parameter defined as s(r) = s0rk can be derived without calculation. Doing this way makes each sensitivity on/ok dependent on on/os0 by the simple relationship on/ok =  s0ln(r)on/os0 with r the location where n is calculated (see Section 3). For simplification purposes, take now a Gauss–Newton algorithm with only two parameters s0 and k, and the fitting of a single observed drawdown curve (i.e. all errors ni come from a single location r). The matrix and vector of the linear system of the Gauss–Newton algorithm allowing seeking the descent direction write

326

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334 -3

50

100 Objective function (m²)

Hydraulic head (m)

99.5

40

x 10

1

Kf(m.s-1) 0.8

99 30

98.5

20

98

0.6 0.4

97.5 10 0

0

10

20

30

40

50

Iterations

96.5 2 10

4

6

10

10

0

0

10

-3

x 10

3

5

0.8

2

0.6

1.5

0.4

1

0.2

0.5

40

50

40

50

x 10

Ssm (m-1) 2.5

30

-8

x 10

Ssf (m-1) 1

20

Iterations

Time (s)

-4

1.2

0.2

97

alpha (m-1.s-1) 4 3 2

0

0

10

20

30

40

50

0

1

0

10

Iterations

20

30

Iterations

40

50

0

0

10

20

30

Iterations

Fig. 4. Inversion of a synthetic test case for the homogeneous dual medium. Best fit between simulated (dots) and observed (solid lines) heads in fractures at 150 m from the well. Evolution with the convergence iteration index of the objective function and parameter values. Kf = hydraulic conductivity of fractures, Ssf, Ssm = specific storage capacity of fractures and matrix, respectively, a = exchange rate between fractures and matrix.



 J Tn  J n  d ¼ J Tn  n 0 P  2 onm

B m os0 J Tn  J n ¼ B @ P on 2 m b os0 m 0 P on 1 m n os0 m B m C J Tn  n ¼ @ P on A m b n os0 m

P onm 2 1 b os0 C m  2 C A P onm b2 os0

b ¼ s0 lnðrÞ

m

m

ð34Þ It is obvious that this linear system is numerically singular but has a solution in the form d k ¼ d s0 =b. Extending the above reasoning to the 8 parameters of a fractal dual medium allows to state that only the descent components of the 4 parameters of type s0 need to be calculated when a single drawdown curve is inverted. If several drawdown curves are inverted at once, the principle does not apply even when sensitivities to k are derived from those to s0. In both cases, however, the conditioning principle of the Jacobian matrix remains useful, and in this test (designed to fit a single drawdown curve) sensitivities to Ssf0 ; Ssm0 ; a0 have been multiplied by 3 · 103, 0.3, 4 · 106, respectively, to scale as K f0 . Again, re-scaling sensitivities to parameters improves conditioning of J Tn  J n . Eigenvalues of the matrix are (1.60 · 1011, 3.13 · 103, 8.32 · 102, 0.11) without rescaling and (55.64, 14.93, 2.62, 0.06) with re-scaled sensitivities.

The reference problem is designed as follows: Q = 50 m3 h1, Kf0 = 103 m1+a s1, Ssf0 = 105 mb1; Ssm0 = 104 mc1, a0 = 5 · 109 md1 s1, a = 0.8, b = c = d = 0.1, 2000 radial meshes of step Dr = 1 m, a time step Dt = 103 s, and a single observed well to be fitted located at r = 150 m. Note that the above parameters were set up to give effective values at r  100–150 m of the same order as those from the homogeneous test case (see above). This does not mean however that observed drawdowns are equivalent since parameters of the fractal approach obey scaling power laws and decrease with the distance. Starting above the reference solution for parameters of type s0 and below the reference solution for parameters of type k (sensitivities are of opposite sign) is a deliberate choice. The initial starting point is 102, 104, 103, 5 · 108, 0.5, 0.05, 0.05, 0.05 for K f0 ; Ssf0 ; Ssm0 ; a0, a, b, c, d, respectively. Again, thanks to the Jacobian conditioning, the parameters evolve smoothly without oscillation during the search (Fig. 5). The power-law exponents increase with iteration indices: b, c, d in following an exponential shape whereas a increases linearly. On the other hand, K f0 decreases following a negative exponential shape and Ssf0 , Ssm0 , a0 decrease almost linearly after sometimes a slight increase during the first dozen iterations. Stability is excellent without divergence of parameter values and objective function when the optimizer is let to seek beyond the optimum. As for the homogeneous problem, the matrix J Tn  J n and vector J Tn  n are kept flat during

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

327

-4

2

0.01

100 Head (m)

Objective function

Kf0

10

Ssf0

7

0.006

6

1

0.6

5

0.004

99.4

0.4 0.5

0.002

99.2

Ssm0

8

0.8 99.6

x 10

9

1

0.008

99.8

1.5

-4

x 10

1.2

4 3

0.2

2 0

0

25

50

99 103

104

iterations

105

0 0

25

50

0

0

25

iterations

time (s)

50

1

0

25

iterations

50

iterations

-8

5

x 10

0.8

0.14

alpha0

4.5

0.75

4 3.5

0.14

0.14

c

b

a

d

0.12

0.12

0.12

0.1

0.1

0.1

0.08

0.08

0.08

0.06

0.06

0.06

0.7

3 0.65

2.5 2

0.6

1.5

0.55

1 0.5

0

25

iterations

50

0.5

0

25

50

0.04

0

iterations

25

iterations

50

0.04 0

25

iterations

50

0.04

0

25

50

iterations

Fig. 5. Inversion of a synthetic test case for the fractal dual medium. Best fit between simulated (dots) and observed (solid lines) heads in fractures at 150 m from the well. Evolution with the convergence iteration index of the objective function and parameter values. Kf(r) = Kf0ra = hydraulic conductivity of fractures, Ssf = Ssf0rb, Ssm = Ssm0rc = specific storage capacity of fractures and matrix, respectively, a = a0rd = exchange rate between fractures and matrix.

the convergence and the ratio hnon/osji/hon2/osi osji = 7/13 at first iteration is lowered down to 50/3 · 105 at iteration #38 and beyond. Again the problem is stable around the optimum. The optimal sought solution is reached after 38 iterations (but a very good result was already reached since iteration #25) and gives (7.4 · 104, 1.13 · 105, 1.14 · 104, 5.69 · 109, 0.74, 0.13, 0.13, and 0.13) for K f0 , Ssf0 , Ssm0 ; a0, a, b, c, and d, respectively. Thus, errors are larger in the fractal case than in the homogeneous one. This comment must be moderated, however, for the following reasons. First, parameters in the form s(r) = s0rk generate concurrent effects. For instance, to get a prescribed scalar value s(r), one may increase s0 and therefore increase k too, or conversely, decrease both s0 and k. Thus, there are many solutions to the prescribed value s(r), which does not help the optimizer as regard accuracy. Several other tests have been undertaken by making the inversion start from different initial points. Errors do not increase even if the initial point is moved away from the reference solution. In fact, the optimizer scans the space of possible solutions and the latter is wide enough to generate nonnegligible differences between solutions even if the objective function is minimized down to the same small value.

Second, matching a single drawdown curve from a single location is not a very well-posed inversion problem for parameters such as s(r) = s0rk that are supposed to vary in space. The objective function does not record the trace of any spatial variations and therefore, pinpoint accuracy cannot be expected. A few tests (not reported here) have been carried out by inverting at once 3–6 drawdown curves located at different distances r from the pumped well. The results are improved with errors of less than 20% on parameters s0, less than 0.1 for exponent a, and less than 0.02 for exponents b, c, and d. This proves that introducing any spatial reference into the objective function helps the optimizer to seek accurate solutions. Third, the parameters, in particular the power-law exponents, must be relocated in the fractal context of the model. Fractal networks are ‘‘statistical’’ objects in the sense that their properties such as their fractal dimension and the scaling-law exponents may vary over a non-negligible range, even within the same object. For instance, computations over 2D regular percolation networks with 2 · 106 bonds have shown the fractal dimension within a single network to vary between 1.55 and 1.99 for clusters of about 105 bonds. Similarly, this fractal dimension spans the range 1.65–1.95 when calculated for different networks whereas the theoretical value

328

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

should be 1.89 [39]. Note that the same behavior applies to scaling exponents of macroscopic permeability and diffusion coefficient of the networks. In the end, variations up to half an order of magnitude on parameters K f0 ; Ssf0 ; Ssm0 , a0, ±0.2 on a, ±0.08 on b, c, d can be considered as an accurate result compared with the variability expected from real natural media and uncertainty from field data. In the end, inverting the fractal model in analyzing a single drawdown curve yields less accurate solutions than the homogeneous approach, but this does not mean that the problem is ill-posed. The relationships in time and/or in space of the fractal parameters are associated with large variations versus time of the model sensitivities to parameters. All together, they bring an indirect spatial reference to the optimizer even in fitting a single drawdown curve. In addition, most concrete applications rest on very few observation points and do not ensure that information is consistent between different locations. In that case, inverting all curves at the same time may fail or at least result in a poor fitting between model and observations. It may appear more relevant to invert each drawdown curve independently, and then analyze sought parameters in terms of statistical spreading (as done in the field tests reported below). In practice, there is no reason to urge any user to fit all data at once. The eventual ill-posed character of inverting a fractal approach is probably more the consequence of concurrent effects of parameters in the scaling laws s(r) = s0rk than fitting a single curve. When a fractal behavior is foreseen, it remains however that inversion may possibly benefit from handling at the same time information from several locations. 4.2. A field case The following examples of inversion are based on data from the Hydrogeologic Experimental Site in Poitiers. This site has been recently developed by the University for the purpose of providing facilities to develop research and engineering on subsurface flow and water resources management. This site is still in the building phase but data from interference pumping tests are already available as well as a synthetic description of the geological context [8] (data soon accessible through the internet). In brief, the site is being developed over Jurassic fractured limestone constituting a confined aquifer about 100 m thick beneath 20 m of Tertiary clays. The recognized fracture network is not very dense and is made of sub-vertical fractures oriented N010 and N110, but the rock matrix is porous and karstic features have been observed on outcrops in the vicinity of the site. The first interference pumping tests carried out on the site (in 2004) have shown a typical behavior of fractured media, with drawdowns increasing more rapidly than linearly with ln(t). They have already been interpreted with a single fractal medium approach and discussed in terms of homogenization scale [9,17]. As stated in the introduction,

the dual-medium approach may be viewed as an alternative and enhances the exploitation of data. The aim here is not to undertake a complete reassessment of all available data with the dual-medium approach. This work should be undertaken in a near future. Here, the point is to show how inversion can be launched to get a few statistics on identified parameters. These statistics are unbiased by the user since inversion is solved automatically, which lowers subjectivity in the interpretation of the macroscopic behavior of the reservoir. To this end, a single pumping test with 21 observed drawdown curves at different locations was inverted. The identification of parameters was performed for each drawdown curve and the sought parameters are plot as a function of the lag distance between observed and pumped wells. Fig. 6 reports on inversion with the homogeneous dualmedium approach. The latter proved to be able to fit accurately the 21 observed drawdown curves as well as their convexity with respect to ln(t) up to t = 138 h of pumping. Identified hydraulic conductivity of fractures Kf is constant whatever the distance (and the orientation) between observed and pumped wells. On the other hand, Ssm, and Ssf, the specific storage capacities of the matrix and fractures respectively, decrease, following a negative exponential or a negative power-law shape with lag distances from 50 to 250 m. The exchange rate coefficient a between the matrix and fractures decreases similarly to storage capacities. This behavior is logical since a can be viewed as a buffer of the pressure depletion propagated by the fractures before it affects the matrix. This buffering effect depends on storage capacities too, i.e., if the matrix storage capacity is huge and the fracture storage capacity is small, the matrix is rapidly stressed but the depletion velocity in the matrix is controlled by the ratio a/Ssm. Thus, with an almost constant propagation in fractures (Kf = constant), but with decreasing storage capacities, it is not unusual for a to follow the same macroscopic behavior as storage capacities. These decreases of parameters with distance are in agreement with descriptions given for fractal networks in which the macroscopic water conduction and water content are assumed to decrease with the size of the system (via the network connectivity and its fractal dimension, which is less than the Euclidean one). This decrease of parameter values advocates for the development of the fractal dual-medium approach proposed in this work, and another series of inversions on the same data set with the fractal model was carried out. Again, the parameters fitted from inversion of the 21 observed drawdown curves are plotted with respect to the lag distance between observed and pumped wells (Fig. 7). Each drawdown curve was inverted independently despite the fact that fractal solutions are more accurate when inverting several curves at once (see Section 4). A rapid look at the raw data shows that they do not follow the basic fundamentals prevailing to the model of interpretation. All experimental curves have almost similar shapes whatever the distance from the pumped well, with the same

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

329

2.0E-05

1.2E-04 1.0E-04

Ssm (m-1)

Kf (m.s -1)

1.8E-05 1.6E-05 1.4E-05 1.2E-05

8.0E-05 6.0E-05 4.0E-05 2.0E-05

1.0E-05

1.0E-07 0

50

100

150

200

250

0

50

100

150

200

250

200

250

Distance (m)

1.0E-05

1.2E-08

8.0E-06

1.0E-08 -1. -1 α (m s )

Ssf (m-1)

Distance (m)

6.0E-06 4.0E-06 2.0E-06

8.0E-09 6.0E-09 4.0E-09 2.0E-09 1.0E-11

1.0E-08 0

50

100 150 Distance (m)

200

250

0

50

100

150

Distance (m)

Fig. 6. Evolution in space of fitted parameters for the homogeneous dual medium from 21 independent inversions of 21 observed drawdown curves at different distances from a single pumped well. Kf = hydraulic conductivity of fractures, Ssf, Ssm = specific storage capacity of fractures and matrix, respectively, a = exchange rate between fractures and matrix.

amplitude at the same time. This behavior cannot be mimicked by a diffusion equation even in a fractal dual medium. One may suggest the presence of karstic flow propagating very rapidly the pumping stress superimposed onto a classical draining of a fractured reservoir. The model will partly fail in inverting several curves at once but it can be used for one curve at a time and statistical interpretation of sought parameters. For all power-law exponents k, no relationship with the lag distance is observed (see Fig. 7). These exponents show statistical variations within acceptable ranges (a = 0.70– 0.72; b = 0.12–0.17; c = 0.14–0.15; d = 0.05–0.14) and below those encountered in calculations over synthetic fractal networks (see above). In other words, it can be considered that the fractal dual-medium approach identifies a single power-law coefficient per hydrodynamic parameter and therefore proves to be a valuable approach to homogenization of radial flow for distances between 50 and 250 m (distances experienced by the tests). Note that these power-law exponents are in ranges that are compatible with 2D regular percolation networks, which was also observed by inverting the same data with a single fractal medium approach [9]. Unfortunately, this invariance with the lag distance is not observed for parameters of type s0 that decrease between 50 and 250 m of about half an order of magnitude for K f0 ; Ssf0 ; Ssm0 and up to one order of magnitude for a0 (Fig. 7). This cannot be explained by statistical variations of power-law exponents. Inversion was re-launched by seeking only s0 with the following fixed values a = 0.71, b = 0.15, c = 0.145, and d = 0.1. A decrease of s0 with the lag distance was also observed with

a slightly narrowed variation range. As a result, the decrease in space provoked by the relation s(r)  rk is not sufficient to fully account for the behavior of the hydrodynamic parameters viewed by the interference test, and a decrease of s0 is observed too. However, these variations remain relatively small compared with possible ranges recorded in natural reservoirs, e.g., [22]. In view of these results, it could be questioned on the validity of the solutions: did the optimizer led to wrong solutions? This is unlikely since: (1) the direct model fits almost perfectly all experimental drawdown curves; (2) sought parameters are compatible with a priori known values from synthetic or natural fractal media (see Section 2, Expression (8)); (3) effective parameters values for distances experienced by the tests are of the same order as those found with the homogeneous approach. This first attempt is probably not enough to state that the fractal dual medium is unable to homogenize the behavior of the aquifer. However, this decrease of s0 suggests that the macroscopic behavior of the aquifer is more complicated than that of the fractal dual medium. Channeled flow due to karstic features may be responsible for this difference, but to date, the knowledge of the aquifer geometry is not sufficient to propose an argued answer and further works are required. Nevertheless, inversion of interference pumping tests in homogeneous and fractal dual media extends the set of tools for interpreting this kind of data and thus provide more information on the hydraulic macroscopic behavior of fractured reservoirs, even if in the present case some discrepancies may appear.

330

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334 0.716

8.0E-04

a

Kf 0 0.712

7.0E-04

0.708

6.0E-04 0.704

5.0E-04

0.700

4.0E-04 0

50

100

150

200

250

0.696 0

50

Distance (m) 1.1E-05

100 150 Distance (m)

200

250

0.18

Ssf 0

9.0E-06

b

0.17 0.16

7.0E-06

0.15

5.0E-06 0.14

3.0E-06

0.13

1.0E-06

0.12

0

50

100 150 Distance (m)

200

250

0

50

100 150 Distance (m)

200

250

50

100 150 Distance (m)

200

250

50

100 150 Distance (m)

200

250

0.154

7.0E-05 Ssm0

6.0E-05

c

0.152 0.150

5.0E-05

0.148 4.0E-05

0.146

3.0E-05

0.144

2.0E-05

0.142 0.140

1.0E-05 0

50

100 150 Distance (m)

200

0

250

7.0E-09

0.14

α0

6.0E-09

d

0.13

5.0E-09

0.12

4.0E-09

0.11

3.0E-09

0.10

2.0E-09

0.09

1.0E-09

0.08

1.0E-11

0.07

0

50

100 150 Distance (m)

200

250

0

Fig. 7. Evolution in space of fitted parameters for the fractal dual medium from 21 independent inversions of 21 observed drawdown curves at different distances from a single pumped well. Kf(r) = Kf0ra = hydraulic conductivity of fractures, Ssf = Ssf0rb, Ssm = Ssm0rc = specific storage capacity of fractures and matrix, respectively, a = a0rd = exchange rate between fractures and matrix.

5. Conclusion This work proposes a new numerical tool based on a dual-medium approach in a radial convergent flow for automatically inverting interference hydraulic pumping tests. The direct problem is solved for both a homogeneous and a fractal dual media using a full-implicit finite-volume scheme. It has been shown, in particular for the fractal

medium, that replacing averaged mesh and intermesh parameters by local values did not affect accuracy but was a key to simplifications of model sensitivity to parameters. Inversion was performed via a Gauss-Newton algorithm in which sensitivities to parameters are derived analytically. These derivations yield a very rapid and accurate optimizer when fitting simulated drawdown curves to

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

experimental ones. In the case of the fractal dual medium for which parameters write as s(r) = s0rk, it was shown that sensitivity to power-law exponents k can be derived from that to s0 without further calculations. In terms of computation costs, this puts the inversion of a fractal medium down to the same level as the homogeneous approach. The exact calculation of sensitivities proved indisputably that dual media have a sequential response in time to the pumping stress. The fracture medium is solicited first and the system is sensitive to the fracture hydraulic conductivity Kf. Because of the limited storage capacity of fractures, the matrix is rapidly stressed to feed the pumping well and, at intermediate times, the exchange rate coefficient a between the matrix and fractures governs the system by buffering the pressure depletion between both media. At long times, the matrix has completely taken over from the fractures the feeding of the system and the matrix storage capacity Ssm becomes the most sensitive parameter. The whole system is drained by fractures but fed by the matrix and the macroscopic hydraulic diffusion tends to the ratio Kf/Ssm. Note that this general behavior drawn from evolution in time of analytical sensitivities is common to both the homogeneous and the fractal dual media. This sequential behavior in time as well as the much contrasted values of sensitivities to parameters are probably the main reasons why the previous works did not report on automatic inversion of pumping tests in dual media. This paper pointed out that the problem is numerically ill-posed and proposed a new method based on sensitivity re-scaling before calculating descent-directions to seek optimal parameters. This re-scaling proved relevant, allowing inversion and making all parameters evolve similarly during optimization. This last point avoids uncontrolled wanderings of the optimizer in the parameter space and therefore deeply reduces the number of convergence iterations to seek a good solution. Computations of synthetic problems have shown the whole inversion procedure to be accurate. For the homogeneous case, the optimizer finds the best solution with less than 1% of errors on parameters, whatever the initial location from which inversion is launched. For the fractal case, errors are more important due to concurrent effects between parameters s0 and k in s(r) = s0rk but they remain much lower than statistical variations of fractal synthetic networks. Moreover, results are improved when several drawdown curves at several locations are inverted at once. In other words, introducing elements of spatial variability in the objective function facilitates the capture of spatial evolution of parameters due to scaling laws. This statement is valid only if data from different locations are consistent between each other, i.e. stemming from the same medium. In the case of a heterogeneity very different from that encapsulated in the direct model (homogeneous scaling power laws in the present case), interpretation of several drawdown curves at once may fail. It is then more valuable

331

to interpret data independently and address statistical variability of sought parameters. A first attempt to invert real data from the calcareous fractured aquifer of the Hydrogeological Experimental Site in Poitiers (France) was performed. A rapid analysis of parameters identified with the homogeneous dual-medium approach on several drawdown curves at different distances from a single pumped well has shown the fracture hydraulic conductivity to be constant, and the other parameters to decrease with the lag distance between observed and pumped wells. Reassessed using the fractal dual-medium approach, the same data set yielded a better homogenization of flow with constant exponents of parameter scaling laws but a decrease of the dimensioning parameters (s0). This divergence between field data results and theoretical expectations suggests that more data processing and a better appraisal of the aquifer geometry are needed to propose an argued explanation. Acknowledgements The authors are grateful to the French National Program for Research in Hydrology that funded this work. Appendix A. Derivation of intermesh hydraulic conductivity The Thiem formula for a steady-state 2D radial flow expresses the pumped flow rate Q at the well as Q ¼ 2peK ij

hj  hi

 ln rj =ri

ðA:1Þ

with e [L] the thickness of the confined reservoir, hk [L] the hydraulic head in rk (r is the radial coordinate from the pumped well, located in r = 0). Kij [L T1] is the average hydraulic conductivity between ri and rj. At any location ri  j between ri and rj, the mass balance principle assumes that Q is preserved. Thus, it can be written from (A.1) hj  hi hj  hij ¼ 2peK j lnðrj =ri Þ lnðrj =rij Þ hij  hi ¼ 2peK i lnðrij =ri Þ

Q ¼ 2peK ij

ðA:2Þ

with hi  j the local head in ri  j and Kk the local hydraulic conductivity values in rk. From (A.2), hi  j is calculated as bj hj  bi hi bi þ bj Kj bj ¼ lnðrj =rij Þ hij ¼

with bi ¼

Ki ; lnðrij =ri Þ ðA:3Þ

Using, e.g., the expression of Q in Ki from (A.2) and substituting hi  j from (A.3) yields Q ¼ 2pe

bi bj ðhj  hi Þ bi þ bj

ðA:4Þ

Identification between (A.4) and the expression of Q in Kij from (A.2) yields

332

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

K ij ¼

bi bj lnðrj =ri Þ K i K j lnðrj =ri Þ ¼ K i lnðrj =rij Þ þ K j lnðrij =ri Þ bi þ bj

ðA:5Þ

For a cylindrical radial meshing at the constant step Dr, with rj = jDr, ri = (j  1)Dr; and ri  j = (j  1/2)Dr (i.e., the inter-node location), Kij is rewritten as K ij ¼

K i K j lnðj=ðj  1ÞÞ K i lnðj=ðj  1=2ÞÞ þ K j lnððj  1=2Þ=ðj  1ÞÞ

ðA:6Þ

For large values of j it can be written j 1 ¼1þ ¼ 1 þ e with j1 j1 1 ) lnðj=ðj  1ÞÞ  e e¼ j1 j 1=2 ¼1þ ¼ 1 þ e0 with j  1=2 j  1=2 1=2 ) lnðj=ðj  1=2ÞÞ  e0 e0 ¼ j  1=2 j  1=2 1=2 ¼1þ ¼ 1 þ e00 with j1 j1 1=2 ) lnððj  1=2Þ=ðj  1ÞÞ  e00 e00 ¼ j1 0

ðB:4Þ

Sensitivity to parameter Ssf is calculated by solving the system Af 

ohfnþ1 ohnf oDf n ohnþ1 oAf nþ1 ¼ Df  þ  hf þ C  m  h oSsf oSsf oSsf oSsf oSsf f

Am 

ohf ohnþ1 ohn m ¼ Dm  m þ C  oSsf oSsf oSsf

ðB:5Þ ðA:7Þ

the other terms in the derivation of (B.1) being zero. Knowing that matrix Df is diagonal, one obtains oAf oDf ¼ oSsf oSsf

ud ! 0 ld ! 0 2 d ! Dr ði  1=2Þ Dt

ðB:6Þ

Sensitivity to parameter Ssm is calculated by solving the system

00

2K i K j Ki þ Kj

ðA:8Þ

which is the harmonic mean between local conductivity values Ki, Kj in ri, rj, respectively. Appendix B. Analytical sensitivity to parameters of the homogeneous dual-medium approach Remember that for a radial convergent flow with cylindrical symmetry, the finite-volume discretized equations write (see Section 2) Af  hfnþ1 ¼ Df  hnf þ C  hnþ1 m þB Am  hmnþ1 ¼ Dm  hnm þ C  hnþ1 f

ðB:1Þ

Sensitivity to parameter Kf is calculated by solving the system Af 

oB ¼ NmH 0 oK f

nþ1

Assuming e  e  1/2e and replacing logarithms in (A.6) by values in (A.7) results in K ij 

for vector B, only the last term enclosing the flux at the external limit depends on Kf and writes: KfNmH0 with Nm the total number of meshes in space. Thus one gets

ohfnþ1 ohnf ohnþ1 oB oAf nþ1 ¼ Df  þC m þ  h oK f oK f f oK f oK f oK f

nþ1

Am 

ohf ohnþ1 ohn oDm n oAm nþ1 m ¼ Dm  m þ  hm þ C   h oSsm oSsm oSsm oSsm oSsm m ðB:7Þ

the other terms in the derivation of (B.1) being zero. Am, Dm are both diagonal matrices, which yields oAm oDm ¼ oSsm oSsm

d!

Dr2 ði  1=2Þ Dt

ðB:8Þ

Finally, sensitivity to parameter a is calculated by solving the system ohfnþ1 ohnf ohnþ1 oC nþ1 oAf nþ1 h  ¼ Df  þC m þ  hf oa m oa oa oa oa ohfnþ1 oC nþ1 oAm nþ1 ohnþ1 ohn h  þ  hm Am  m ¼ Dm  m þ C  oa f oa oa oa oa ðB:9Þ Af 

Matrices Am and C are both diagonal, which yields

nþ1

ohf ohnþ1 ohn Am  m ¼ Dm  m þ C  oK f oK f oK f

ohfnþ1 ohnf ohnþ1 Af  ¼ Df  þC m oSsm oSsm oSsm

ðB:2Þ

oAf oAm oC ¼ ¼ oa oa oa

ud ! 0 ld ! 0 d ! Dr2 ði  1=2Þ

ðB:10Þ

the other terms in the derivation of (B.1) being zero. Using in the following, the notations ud for upper diagonal, ld for lower diagonal and d for diagonal, one gets

Appendix C. Analytical sensitivity to parameters of the fractal dual-medium approach

ud ! i oAf ¼ ld ! 1  i oK f d ! 2i  1

Remember that for a radial convergent flow with cylindrical symmetry, the finite-volume discretized equations write (see Section 2)

ðB:3Þ

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

Af  hfnþ1 ¼ Df  hnf þ C  hnþ1 m þB

ðC:1Þ

Am  hmnþ1 ¼ Dm  hnm þ C  hnþ1 f

Sensitivities to parameters s ¼ K f0 and a from K f ðrÞ ¼ K f0 ra are calculated by solving the system ohnþ1 ohnf ohnþ1 oB oAf nþ1 f  ¼ Df  þC m þ  hf os os os os os ohnþ1 ohn ohnþ1 Am  m ¼ Dm  m þ C  f os os os Af 

ðC:2Þ

the other terms in the derivation of (C.1) being zero. Using in the following the notations ud for upper diagonal, ld for lower diagonal and d for diagonal, one gets ud ! i1a Dra oAf 1a ¼ ld ! ði  1Þ Dra oK f0 1a d ! i1a Dra þ ði  1Þ Dra ud ! K f0 lnðiDrÞi1a Dra ld ! K f0 lnðði  1ÞDrÞði  1Þ oAf h ¼ d ! K lnðiDrÞi1a Dra f0 oa 1a

þ lnðði  1ÞDrÞði  1Þ

1a

Dra

ohfnþ1 ohnf ohnþ1 ¼ Df  þC m os os os n ohfnþ1 oAm nþ1 ohnþ1 oh oD m  hnm þ C    hm Am  m ¼ Dm  m þ os os os os os ðC:7Þ

Af 

the other terms in the derivation of (C.1) being zero. Am and Dm are both diagonal matrices, which yields oAm oDm Dr2 1c ði  1=2Þ Drc ¼ d! oSsm0 oSsm0 Dt oAm oDm Dr2 1c ¼ d ! Ssm0 ln ðði  1=2ÞDrÞ ði  1=2Þ Drc oc oc Dt ðC:8Þ Finally, sensitivities to parameters s = a0 and d from a(r) = a0rd are calculated by solving ohfnþ1 ohnf ohnþ1 oC nþ1 oAf nþ1 :h  ¼ Df : þ C: m þ :h os m os os os os f ohfnþ1 oC nþ1 oAm nþ1 ohnþ1 ohn :h  þ :h Am : m ¼ Dm : m þ C: os f os os os os m Af :

Dra

ðC:3Þ

ðC:9Þ

i

Matrices Am and C are both diagonal, which yields

for vector B, only the last term enclosing the flux at the external limit depends on Kf(r) and writes: K f0 Nm1a Dra H 0 with Nm the total number of meshes in space. Thus one gets oB oB ¼ K f0 lnðNmDrÞNm1a Dra H 0 ¼ Nm1a Dra H 0 oK f0 oa

oAf oAm oC ¼ ¼ oa0 oa0 oa0

oAf oAm oC ¼ ¼ od od od

ðC:4Þ Sensitivities to parameters s ¼ Ssf0 and b from Ssf ðrÞ ¼ Ssf0 rb are calculated by solving the system ohnþ1 ohnf oDf n ohnþ1 oAf nþ1 f ¼ Df  þ  hf þ C  m   hf os os os os os nþ1 ohf ohnþ1 ohn Am  m ¼ Dm  m þ C  os os os ðC:5Þ Af 

the other terms in the derivation of (C.1) being zero. Knowing that matrix Df is diagonal, one obtains oAf oDf ¼ oSsf0 oSsf0 oAf oDf ¼ ob ob

ud ! 0 ld ! 0 2

d ! Dr ði  1=2Þ Dt

1b

333

Drb

ud ! 0 ld ! 0 2

d ! Ssf0 ln ðði  1=2ÞDrÞ Dr ði  1=2Þ Dt

1b

Drb

ðC:6Þ Sensitivities to parameters s ¼ Ssm0 and c from Ssm ðrÞ ¼ Ssm0 rc are calculated by solving the system

ud ! 0 ld ! 0 d ! Dr2 ði  1=2Þ1d Drd ud ! 0

ðC:10Þ

ld ! 0 d ! a0 ln ðði  1=2ÞDrÞ Dr2 ði  1=2Þ

1d

Drd

References [1] Acuna JA, Ershaghi I, Yortsos YC. Practical applications of fractal pressure transient analysis of naturally fractured reservoirs. Soc Petrol Eng Form Eval 1995:173–9. [2] Acuna JA, Yortsos YC. Application of fractal geometry to the study of networks of fractures and their pressure transient. Water Resour Res 1995;31(3):527–40. [3] Adler PM, Thovert J-F. Fractures and fracture networks. Dordrecht Netherlands: Kluwer Academic Publications; 1999. [4] Aprilian S, Abdassah D, Mucharam L, Sumantri R. Application of fractal reservoir model for interference test analysis in Komojang geothermal field (Indonesia). Soc Petrol Eng 1993;26465:511–20. [5] Bangoy LM, Bidaux P, Drogue C, Plegat R, Pistre S. A new method of characterizing fissured media by pumping tests with observation wells. J Hydrol 1992;138(1–2):77–88. [6] Barenblatt GI, Zheltov IP, Kochina IM. Basic concepts on the theory of seepage of homogeneous liquids in fissured rocks. J Appl Math Mech (Engl Transl) 1960;24:1286–303. [7] Barker JA. A generalized radial flow model for hydraulic tests in fractured rocks. Water Resour Res 1988;24(10):1796–804. [8] Bernard S. Characterization of limestone fractured reservoirs: application to the Hydrogeological Experimental Site (HES) of the University of Poitiers. PhD thesis, University of Poitiers, France, 2005.

334

F. Delay et al. / Advances in Water Resources 30 (2007) 314–334

[9] Bernard S, Delay F, Porel G. A new method of data inversion for the identification of fractal characteristics and homogenization scale from hydraulic pumping tests in fractured aquifers. J Hydrol; in press. doi:10.1016/j.jhydrol.2006.01.008. [10] Bogdanov II, Mourzenko VV, Thovert J-F, Adler PM. Pressure draw-down well tests in fractured porous media. Water Resour Res 2003;39(1):1021. doi:10.1029/2000WR000080. [11] Carrayrou J, Mose R, Behra P. Operator-splitting for reactive transport and comparison of mass balance errors. J Contam Hydrol 2004;68:239–68. [12] Carrera J, Neuman SP. Estimation of aquifer parameter under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information. Water Resour Res 1986;22(2): 199–210. [13] Carrera J, Neuman SP. Estimation of aquifer parameter under transient and steady state conditions: 2. Uniqueness, stability, and solution algorithms. Water Resour Res 1986;22(2):211–27. [14] Chang J, Yortsos YC. Pressure-transient analysis of fractal reservoirs. Soc Petrol Eng J 1990;18710:631–43. [15] Cooper HH, Jacob CE. A generalized graphical method for evaluating formation constants and summarizing well field history. Am Geophys Union Trans 1946;27:526–34. [16] Delay F, Porel G. Inverse modeling in the time domain for solving diffusion in a heterogeneous rock matrix. Geophys Res Lett 2003;30(3):1147. doi:10.1029/2002GL016428. [17] Delay F, Porel G, Bernard S. Analytical 2D model to invert hydraulic pumping tests in fractured rocks with fractal behavior. Geophys Res Lett 2004;31(16). doi:10.1029/2004GL020500. [18] Desbarats AJ. Spatial averaging of transmissivity in heterogeneous fields with flow toward a well. Water Resour Res 1992;28(3): 757–67. [19] Desbarats AJ. Spatial averaging of hydraulic conductivity under radial flow. Math Geol 1994;26(1):1–21. [20] Dershowitz W, Miller I. Dual porosity fracture flow and transport. Geophys Res Lett 1995;22(11):1441–4. [21] Gerke HH, Van Genuchten MTh. Macroscopic representation of structural geometry for simulating water and solute movement in dual porosity media. Adv Water Res 1996;19(6):343–57. [22] Huang W, Di Donato G, Blunt MJ. Comparison of streamline-based and grid-based dual porosity simulations. Soc Petrol Eng J 2004;43(2):129–37. [23] 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 Res 2001;24(6):863–76. [24] Larocque M, Delay F, Banton O. Stochastic inverse modeling: field scale application and uncertainty analysis. Ground Water 2003;41(1):15–23. [25] Le Borgne T, Bour O, de Dreuzy J-R, Davy P, Touchard F. Equivalent mean flow models for fractured aquifers: insights for a pumping test scaling interpretation. Water Resour Res 2004;40(3). doi:10.1029/2003WR002436.

[26] Leveinen J. Composite model with fractional flow dimensions for well test analysis in fractured rocks. J Hydrol 2000;234:116–41. [27] de Marsily G, Delay F, Goncalves J, Renard P, Teles V, Violette S. Dealing with spatial heterogeneity. Special Issue: The future of Hydrogeology. Hydrogeology J 2005;13:161–83. doi:10.1007/s10040004-0432-3. [28] Marsily G, de Delay F, Teles V, Schafmeister A. On some current methods to represent the heterogeneity of natural media in hydrogeology. Hydrogeology J 1998;6(1):115–30. [29] McLaughlin D, Townley LR. A reassessment of groundwater inverse problem. Water Resour Res 1996;32(5):1131–61. [30] Neuman SP, Di Federico V. Multifaceted nature of hydrogeologic scaling and its interpretation. Rev Geophys 2003;41(3). doi:10.1029/ 2003RG000130. [31] O’Shaughnessy B, Procaccia I. Diffusion in fractals. Phys Rev A 1985;32:3073. [32] Ouillon G, Sornette D. Unbiased multi-fractal analysis. Application to fault patterns. Geophys Res Lett 1996;23(23): 3409–12. [33] Quintard M, Whitaker S. Transport in chemically and mechanically heterogeneous porous media. I: Theoretical development of regionaveraged equation for slightly compressible single-phase flow. Adv Water Res 1996;19(1):29–47. [34] Quintard M, Whitaker S. Transport in chemically and mechanically heterogeneous porous media. III: Large-scale mechanical equilibrium and the regional form of Darcy’s law. Adv Water Res 1998;21:617–29. [35] Sahimi M. Flow, dispersion and displacements processes in porous media and fractured rocks: from continuous models to fractals, percolation, cellular automata and simulated annealing. Rev Mod Phys 1993;65(4):1393–534. [36] Samardzioska T, Popov V. Numerical comparison of the equivalent continuum, non homogeneous and dual porosity models for flow and transport in fractured porous rocks. Adv Water Res 2005;28(3):235–55. [37] Sanchez-Vila X, Axness CL, Carrera J. Upscaling transmissivity under radially convergent flow in heterogeneous media. Water Resour Res 1999;104(3):613–21. [38] Stach S, Cybo J, Chmiela J. Fracture surface – fractal or multi-fractal. Mater Charact 2001;26:163–7. [39] Stauffer D, Aharony A. Introduction to percolation theory. London: Taylor & Francis; 1994. [40] Tarentola A. Inverse problem theory: methods for data fitting and model parameter estimation. New York: Elsevier; 1987. [41] Walsh JJ, Watterson J. Fractal analysis of fracture patterns using the standard box-counting technique: valid and invalid methodologies. J Struct Geol 1993;15:1509–12. [42] Warren JE, Root PJ. The behavior of naturally fractured reservoirs. Soc Petrol Eng J 1963;3:245–55. [43] Wilson TH. Scale transitions in fracture and active fault networks. Math Geol 2001;33(5):591–613. [44] Zimmerman RW, Chen G, Hagdu T, Bodvarsson GS. Numerical dual-porosity model with semi-analytical treatment of fracture-matrix flow. Water Resour Res 1993;29(7):2127–37.