An inverse problem methodology to identify flow channels in fractured media using synthetic steady-state head and geometrical data

Advances in Water Resources 33 (2010) 782–800

Contents lists available at ScienceDirect

Advances in Water Resources j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / a d v wa t r e s

An inverse problem methodology to identify flow channels in fractured media using synthetic steady-state head and geometrical data R. Le Goc a,b,⁎, J.-R. de Dreuzy a, P. Davy a a b

Géosciences Rennes, UMR 6118 CNRS, Université de Rennes 1, CS 74205, F-35042 Rennes Cedex, France Itasca Consultants SAS, 64 chemin des mouilles, F-69134, Écully Cedex, France

a r t i c l e

i n f o

Article history: Received 17 December 2009 Received in revised form 12 April 2010 Accepted 14 April 2010 Available online 26 April 2010 Keywords: Inverse problem Simulated annealing Fractured media Flow channels

a b s t r a c t We present a methodology for identifying highly-localized flow channels embedded in a significantly less permeable medium using steady-state head and geometrical data. This situation is typical of fractured media where flows are often strongly channeled at the scales of interest (10 m–1 km). The objective is to identify both geometrical and hydraulic characteristics of the conducting structures. Channels are identified in decreasing order of importance by successive optimizations of an objective function. The identification strategy takes advantage of the hierarchical flow organization to restrict the dimension of the solution space of each individual optimization step. The characteristics of the secondary channels are strongly determined by the main flow channels. The latter are slightly modified by the secondary channels through the addition of a regularization term to the main channel characteristics in the objective function. As the objective function is strongly non-convex with numerous local minima, inversion is performed using a stochastic algorithm (simulated annealing). We assess the possibilities of the hierarchical identification strategy on simple synthetic steady-state flow configurations where hydraulic data are made up of 25 regularly spaced heads and of the boundary conditions. Those flow structures that are dominated by at most two simple channels can be identified with these head data only. Configurations comprising up to three complex and interconnected channels can still be identified with additional geometrical information including the distances of piezometers to their closest channel. The capabilities of the hierarchical identification strategy are limited to flow structures dominated by at most three equivalent flow channels. We finally discuss the perspectives of application of the method to transient-state data obtained on a more restricted number of piezometers. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Flow channeling is widespread in fractured media, mainly in lowpermeability rocks where flows concentrate in highly transmissive flow paths within a restricted number of fractures [39,51]. Highly-channeled media refer here to media where the main channels, i.e. those carrying the largest part of the flow, are sparse at the scale of interest. The hydrological properties of highly-channeled media are mainly controlled by the continuity and connectivity of the channels [21,28,46]. In such media, characterizing the hydraulic parameters comes down to identify these main channels and their hydraulic properties. Existing approaches based on extensive resistivity [22] or hydraulic data [30] yield zones of higher permeability rather than highly-localized channels. The objective of this article is to explore the possibilities of identifying highly-localized flow channels using head data. This is a rather new approach to inverse

⁎ Corresponding author. Itasca Consultants SAS, 64 chemin des mouilles, F-69134, Écully Cedex, France. Tel.: +33 4 72 18 04 20. E-mail address: [email protected] (R. Le Goc). 0309-1708/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.04.011

problems in hydrogeology, where most of the developments have focused on continuous porous media with permeability fields parameterized by zones or geostatistical functions [8,17,42]. We thus aim at solving the inverse problem with a parameterization relying on discrete flow structures. Parameterization is known to be a key issue for handling inverse problems, since optimization algorithms can effectively identify only a limited number of parameters [57]. Parameterization is also required for adapting the model complexity to the information available in the data to avoid under- and over-parameterization. Applied to fractured media, continuous-like parameterizations have led to the definition of wide fracture zones within a stochastic less permeable continuous medium [35,54], whereas 3D Discrete Fracture Network (DFN) parameterizations have led to more localized flow structures but failed to identify the major flow channels [6,7]. Because of the complexity of the problem, either simplifications or introduction of additional constraints have been considered. Gwo [27] simplified the inverse problem by discretizing the channels on a regular grid, where grid elements were either highly transmissive fractures or an almost impervious matrix. Fracture geometries and transmissivities were identified by a genetic algorithm but flow path geometries could be correctly estimated for single

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

straight-line flow paths only. Mauldon et al. [40] modeled a fracture network with a regular lattice where each element was either conductive or not. The geometry of conductivity was calibrated with simulated annealing using steady-state heads and the resulting models could identify the most transmissive areas rather than main flow channels. Reductions of the parameter space have also been proposed by replacing the network structure by self-similar structures like Iterated Function Systems (IFSs) conditioned to transient flow tests [19]. The main advantage of IFSs is to yield a large variety of shapes with very few parameters. Optimization performed by a stochastic method (simulated annealing) yielded a broad range of highly different structures. Another simplification consisted in identifying in 2D only one or two single structures [44,49] with non-linear least-square regression or co-kriging. Identification was successful only in those cases where the observation network surrounds the transmissive structure and covers an area close to the size of the fracture zone. Considering the lack of constraints, Bruel [3] used hydraulic and stress data to calibrate fracture transmissivities in a DFN with a previously identified geological structure. For the same reasons, Day-Lewis [14] used inter-well connectivity data and hydraulic data to estimate successively the geometry and the transmissivity of four fractured zones. They showed that accounting for discrete highly transmissive structures significantly improves the model accuracy. Following the same spirit, the standard parameterizations used in the petroleum industry and in the storage of high-level nuclear wastes now rely first on the definition of a reference geometrical structure of the highly transmissive zones, and second on the identification of hydraulic properties of fixed fracture families [4,5,18,26].

783

In the previous approaches, as in continuous models, the flow structure is deduced from geological observations or is limited to a very small number of structures (1 or 2) and only the hydrological properties (i.e. the transmissivities) are calibrated. The purpose of this paper is to set up an identification strategy that identifies both the structure and its hydrological properties by focusing on preferential flow pathways. This is a rather new approach since, whereas continuous approaches or equivalent porous media calibrations aim at identifying the mean property of media, we aim at identifying the extreme properties caused by major discontinuities. The strategy relies on an iterative procedure that identifies first the major flow and then, if the data contain enough information, the second-order channels. We assess the method on a set of synthetic configurations with an increasing level of complexity. In addition, we first use only synthetic steady-state head data and then introduce, for more complex configurations, basic synthetic geometrical information. The purpose is not to remain close to natural media and realistic data but rather to use simple synthetic configurations to assess the feasibility of the approach in terms of structure complexity and data requirement. We thus propose a parameterization designed to identify both the channel geometry and permeability. The parameterization is based on the main flow channels rather than on the fracture zones. This leads to much more localized flows than those obtained in heterogeneous porous parameterizations, while requiring fewer parameters and details than DFN-like parameterizations. The identification strategy is based on a hierarchical analysis of the flow paths. Structures of decreasing importance are identified using successive objective function optimizations performed by a simulated annealing algorithm. A stepby-step procedure is developed where results from first-order channels

Fig. 1. (a) Synthetic fracture network. (b) Flows computed with a constant sub-horizontal hydraulic gradient of direction given by the arrow. All fractures have the same transmissivity. The grey color is proportional to the logarithm of flow normalized by the total flow entering the domain. (c) Simplified fracture structure carrying 70% of the flow. (d) Sketch of the possible main flow channels, with a main channel (black), a secondary channel (dark grey) and third-order structures (light grey).

784

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

are slightly modified when adding second-order links. The first objective of this paper is to define the parameterization in terms of flow channels and to setup the identification methodology on 2D synthetic flow channel structures using synthetic hydraulic and geometrical data. Meanwhile, the second objective is to assess which type of flow configurations can be identified with head data only and to which extent the use of geometrical data can help solving the problem. In other words, we aim at defining the information concealed behind hydraulic and geometrical data. We present in Section 2 the parameterization adapted to highly-channeled media underpinning the strategy adapted to solve the inverse problem defined in Section 3. We test the methodology in Section 4 on three increasingly complex fracture flow structures using first hydraulic and then geometrical data. Section 5 presents a further evaluation of the methodology on a larger set of 20 examples with both hydraulic and geometrical data. Advantages and limitations of the data used as well as possible extensions of the inversion methodology to more classical transient-state data are discussed in Section 6. 2. Definition and principle of the hierarchical channel identification We only consider highly-channeled flow cases where flows are focused in very few structures at the scale of interest. We use this assumption to reduce the parameter space. When flows are more evenly distributed in the medium, geostatistically-based parameterizations are advised since they are dedicated to this kind of media [8,17,42]. Section 5 will assess quantitatively the domain of applicability of these

two different approaches by using the channeling indicators previously developed in Le Goc et al. [36].

2.1. Hierarchical channel organization In fractured media, flows are highly channeled in a small number of fractures below the homogenization scale (i.e. Representative Elementary Volume, REV), which is roughly the distance between two main channels [37]. Below this scale, flows are focused in preferential flow channels, or “paths of least resistance” [52], and are organized hierarchically in channels of decreasing importance [15,16]. Fig. 1 gives an example of channeled flows obtained on a synthetic 2D fracture network composed of a square system where fracture positions and orientations are drawn from uniform distributions, and lengths are power-law distributed. Among all generated fractures (Fig. 1a), only a few carry a non-negligible amount of flow (Fig. 1b) and even fewer carry 70% of the total flow (Fig. 1c). Channels are visually identified and organized in decreasing order of importance in Fig. 1d with a dominant channel in black carrying 25% of the flow, a secondary channel in dark grey carrying 15% of the flow and higher-order structures in light grey carrying from 5% to 15% of the flow. The highly simplified network in Fig. 1d leads to a flow structure very similar to that of the simplified network in Fig. 1c, providing that the segment transmissivities are increased to keep the same total flow. All the remaining flowing structures of Fig. 1a are integrated in the background matrix.

Fig. 2. Principle of the identification strategy applied to the flow structure of Fig. 1d. Channels are identified in decreasing order of importance. At each step, an objective function with a restricted number of parameters is defined and optimized. Steps 1 to 3 illustrate the identification of the first channel with a shape of increasing complexity. Steps 4 to 5 illustrate the addition and refinement of a secondary channel. Step 6 shows how channel connectivity can increase with a channel connecting two previously identified channels. Steps 7 to 9 complete the identification strategy. The background color is correlated to the background permeability. It is expected to decrease with the number of structures identified.

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

2.2. Quantitative definition of a channel We define a channel as a succession of highly transmissive flow paths maintaining high flow rates over long distances. Thus, for twodimensional domains, a channel is a polyline, i.e. a series of connected line segments, defined by the location of its nodes and by its transmissivity. The channel extremities belong either to the boundary of the domain or to sinks/sources or to other channels. Fig. 1d displays five channels of varying complexities along which flows are at least one tenth of the total flow entering the domain. The channel parameterization is directly derived from this definition. The full fracture network

785

is parameterized by a restricted set of channels superimposed on a homogeneous background matrix of smaller permeability. This parameterization can be organized in three classes of parameters comprising channel geometry and topology (connectivity between channels), channel transmissivities, and matrix permeability. 2.3. Iterative channel identification The identification strategy relies on the hierarchical organization of the flow structure and consists in an iterative identification of the channels in decreasing order of importance. Fig. 2 illustrates this

Fig. 3. Flowchart of the flow channel identification strategy consisting of four nested loops underlined by different grey colors. The four loops are the calibration loop, the single channel identification loop, the hierarchical channel identification loop and the post-processing analysis loop. P is the last parameterization considered as valid, P ′ is the current parameterization, ns is the simulation number, nmax is the number of simulations to be performed, and the stopping criteria are given by (5) and (6).

786

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

algorithm. In this section, we present these four loops together with two additional sections dedicated to the definition of the objective function and to the method used for solving the direct problem.

strategy applied to the example in Fig. 1. At each step, we solve an inverse problem on a single channel while modifying only marginally the already identified ones (e.g. Fig. 2, steps 1 to 3). We thus restrict the complexity of each step of the inverse problem by reducing the parameter space. Variations in the previously identified parameters are strongly damped by the addition of regularization terms in the objective function. Only the newly added parameters are allowed to vary freely within their possible range. Because of the hierarchical channel organization, our strategy consists in refining the shape of the most important channels before introducing less important additional channels. This is illustrated in the simplified flow structure of Fig. 1d by an idealized succession of identification steps (Fig. 2). The channel refinement is performed by increasing the number of nodes in the polyline (Fig. 2, steps 2, 3, 5, 8 and 9). Once a channel has been identified, we add a new channel in the identification process (Fig. 2, steps 1, 4, 6, and 7). Each added channel enhances the connectivity between the system limits (Fig. 2, steps 1, 4, and 6) or between the channels (Fig. 2, step 7). The strategy consists in identifying first the main channel (Fig. 2, steps 1 to 3) and then the other channels in decreasing order of importance (Fig. 2, steps 4 to 9). One of the critical point of the algorithm is the criterion used to shift from the refinement of an existing channel to the addition of a new channel. This will be detailed in the following section.

3.1. Hierarchical channel identification We first describe the “hierarchical channel identification” loop in Fig. 3. This loop adds new channels as long as they yield a better agreement between the model and data. Each step of the loop consists in initializing a new channel, optimizing its structure parameters from the single channel identification loop and, after a stopping criterion (defined below) has been met, analyzing the benefits of its addition. The new channel is initialized by drawing its characteristics from uniform distributions between their extreme values given in Table 1. The positions of the channel extremities are drawn uniformly from the already existing channels and boundaries. The log10-transmissivity of the channel is chosen at least two orders of magnitude greater than the background permeability in a uniform distribution in the interval [min(log10Tmat)+2, max(log10Tmat)+2], where min(log10Tmat) and max(log10Tmat) are the limit values allowed for the decimal logarithm of the background matrix transmissivity. The improvement of the model by the possible addition of the channel is assessed by two criteria. The first criterion is given by a measure of the agreement between the model and data. The measure is the classical quadratic mismatch sum of squares RE between observed data d and modeling results d ′(nc) obtained with the number of channels nc of the current model:

3. Inverse problem methodology and implementation The algorithm is composed of four nested loops (Fig. 3). As will be shown in Section 3.3, the objective function is highly non-convex and possesses numerous local minima. We thus use a stochastic optimization yielding a distribution of sub-optimal solutions rather than a unique optimal solution. The external loop drives a complete identification process that ends with the final solution being derived from all the suboptimal solutions obtained using the identification strategy. This strategy, directly derived from the discussion in Section 2.3, is made up of three loops. The outer loop concerns the hierarchical channel identification that adds and identifies a new channel at each iteration. It calls the single channel identification loop that refines the structure of the channel by adding new nodes to the channel discretization. The inner parameter calibration loop identifies the characteristics of the polyline nodes. It is in this loop that the optimization of an objective function takes place. Once the distribution of solutions has been built up, we look for similarities and differences among the solutions with a clustering

Nw

REðnc Þ = ∑

!2 0 di −di ðnc Þ

i=1

σid

;

ð1Þ

where Nw is the number of observations, di and d′i are the ith observed and simulated data, respectively, and σid is the error associated with the ith observed data. The second criterion is given by the comparison of the parameter uncertainty PU with the structure discrimination SD. PU is based on the covariance matrix of the estimated parameters defined by Yeh and Yoon [55]. For a parameterization with channels nc and parameters NP(nc), PU writes [50,55]:

PU ðnc Þ =

∥

∥

 −1 Nw −NP ðnc Þ  T Jd ⋅Jd ; RE

ð2Þ

Table 1 Parameter list and possible values for the whole identification strategy. Parameter

Definition

Range of possible values

Value chosen for the method

Log10(Tmat) Log10(Tch) B1ch, B2ch

Log-transmissivity of the background matrix Log-transmissivity of channels Borders or channels connected by the current channel

Initialized at 10− 3 Drawn in a uniform distribution Drawn in a uniform distribution

C1ch, C2ch (x,y)ch σih

Position of the channel's extremity Position of points defining the channel shape Weight associated to the ith head value (Eq. (7))

10− 5 to 100 10− 3 to 102 0 to nb + nc where nc is the number of channels and nb the number of system borders 0 to the border/channel length Inside the system 10− 3 to ∞

σid σik α1 α2 λj,k

Weight associated to the ith distance value (Eq. (7)) Weight associated to the ith parameter (Eq. (7)) Stopping criterion associated to RE (Eq. (5)) Stopping criterion associated to PU and SD (Eq. (6)) Weighting term associated to previously identified parameters

N0 N0 ]0;1] ≥1 N0

ξi

Weighting term associated to channels during the post-processing

N0

Itmax min Fobj Tinit ns dc

Allowed number of iterations of the simulated annealing Aimed value of the objective function Initial temperature of the simulated annealing Number of single solutions used in the post-treatment Threshold of the clustering algorithm

≥1 ≥0 N0 ≥1 N0

Drawn in a uniform distribution Drawn in a uniform distribution 1 or equal to the difference between the channeled and the homogeneous configuration 1 1 1 2 ncur(k) − nk where ncur is the current step and nk the step where the kth parameter were introduced. Fobj/Si with Fobj the objective function value corresponding to the channel and Si its sensitivity. 5.104 10− 4 100, 101 or 102 50 Between 0.5 and 3.5 as a function of the configuration

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

with α1 ≤ 1 and that:

where Jd is the Jacobian matrix of RE:

Jd =

p ∂di σi ⋅ ∂pk σkd

! :

ð3Þ

i;k

In Eq. (3), di refers to the ith observation, σid to its associated error, pk to the kth parameter and σkp to its associated error. SD is used as a criterion to discriminate between two parameterizations and is based on the distance between the two corresponding solutions. It is given by Tsai et al. [50]: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u u 1 Nw d ðn + 1Þ−d ðn Þ 2 i c i c SDðnc + 1Þ = t ∑ : Nw i = 1 σid

ð4Þ

The addition of channel “nc + 1” significantly contributes to the improvement of the model providing that: REðnc + 1Þ bα1 ; REðnc Þ

787

ð5Þ

PU ðnc + 1Þ bα2 : SDðnc + 1Þ

ð6Þ

The first criterion α1, based on RE, indicates an improved agreement between simulated and observed data. The second criterion α2, based on PU and SD, ensures that the two models with channels nc and nc + 1 can be discriminated according to the parameter uncertainty. Overparameterization, leading to a higher parameter uncertainty [31,50], is avoided. These criteria have been used for the identification of zones in heterogeneous continuous media with α1 = 1 and α2 = 2 [10,50]. We take these same simple values in this study. 3.2. Single channel identification All channels are identified iteratively. At first, each channel is a single segment initialized in the hierarchical channel identification loop. Its parameters, i.e. the position of its extremities and its transmissivity, are optimized by the parameter calibration loop. Once optimized, the segment is the starting element of a polyline made up of a series of segments and nodes. At each further step, a new node is added and both its position and the transmissivity of the polyline are optimized. The stopping criteria of the refinement are those of the hierarchical channel

Fig. 4. Objective function versus model parameters in a simple test case comprising two channels. (a) Configuration used to compute the objective function; red dots indicate observation locations. Objective function values (Eq. (7)) with σhi = 1 and σri = ∞ as a function of (b) the first channel parameter, (c) the second channel parameters and (d) the matrix permeability.

788

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

Fig. 5. Sketch of the post-processing clustering algorithm identifying the prevalent characteristics (right) within a set of individual solutions (left). Colors in the individual solutions refer to the channels identified as prevalent characteristics.

identification loop (5 and 6). The full identification procedure stops when no improvement results from the last channel or from the addition of a new channel. 3.3. Definition of the objective function Determination of the polyline parameters is performed by the parameter calibration loop based on the optimization of an objective function made up of two terms. The first term is the mismatch between modeled and observed head data RE (1). The second term is a regularization term for current and previously identified channels [8]. The objective function writes: Nw   Fobj nc; ; NP ðnc Þ = ∑

i=1

!2 di −di′ ðnc Þ σid

nc −1 NPðjÞ

+ ∑ ∑ λj;k ⋅ j=1 k=1

!2 pk ðjÞ−pk′ ðjÞ σkp

;

(nc) is the sum of NP(j) for the first channels nc and λj,k is the strength of the regularization of the kth parameter of the jth channel. The regularization term (second term in Eq. (7)) restricts the possible variations in the parameters of the previously identified channels nc − 1. This strategy is critical to the success of the identification strategy. Neglecting it precludes the possibility to modify first-order channels when adding higher-order channels. On the opposite, searching all parameters simultaneously would exceed the abilities of existing optimization algorithms. The objective is to allow some degree of freedom to the previously identified parameters. In fact, the previously identified channels may be impacted by subsequent channel additions and modifications while the regularization term forces the former parameters to remain close to their previously identified values. We denote ncur the current step index and nk the step index where parameter k was introduced. The weighting coefficient λj,k is a function that linearly increases with the number of steps between ncur and nk:

ð7Þ λj;k = ncur ðkÞ−nk : where pk is the kth parameter introduced in the model with channels nc, pk′ is its value derived from previous steps, σkp is the associated error, NP(j) is the number of parameters describing the jth channel, NP

ð8Þ

Increasing the value of λj,k restricts the possible modifications of the already identified corresponding parameter, which means that the

Fig. 6. Flow rates ϕ of three tested configurations, ordered in increasing complexity. The simplest configuration (a) is composed of two straight-line channels. The intermediary configuration (b) is composed of two channels that are more complex. The complex configuration (c) has several second-order well-connected structures.

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

parameters added earlier in the model remain closer to the value identified in the previous steps. 3.4. Parameter calibration At each iteration, the single channel identification loop performs the optimization of the objective function (7) with a restricted number of parameters. In the case of the addition of a new channel, the number of parameters to be identified is 6: the position of the extremities of the channel (the connected structure and the curvilinear coordinate of the connection point), the log-transmissivity of the channel and the transmissivity of the background matrix. In the case of the modification of a channel, the number of parameters to be identified is 4: the position of the nodes of the polyline (2 parameters), the log-transmissivity of the channel and the log-transmissivity of the background matrix. The background “matrix” transmissivity is systematically calibrated as it is expected to decrease through the identification process. We note here that log-transmissivity as a model parameter was found to yield much better results than transmissivity. This is mainly due to the range of transmissivities varying over several orders of magnitude while geometrical parameters have a narrower variation range. The optimization algorithm has been constructed based on the properties of the objective function. Fig. 4 shows the objective function as a function of parameter values. Each curve represents the variations of the objective function according to the variations of a single parameter, the others remaining equal to their optimal value. The 1st and 2nd connected structures refer to the borders or channels

789

connected by the channel (between 0 and 4) and the position refers to the curvilinear coordinate of the intersection along the connected structure. All curves show that the objective function is not convex and displays numerous local minima (Fig. 4b and d). We thus choose a simulated annealing algorithm, which is among the simplest and most common stochastic inversion algorithms [48] used in groundwater calibration problems [14,38,43,56]. We give in Appendix A the chosen implementation of the simulated annealing. The algorithm is stopped either when a maximum number of iterations itmax is reached or when min the objective function falls below a minimum value Fobj . 3.5. Numerical method for solving the direct problem The direct problem consists in solving the flow equation in wellconnected channels embedded in a uniform homogeneous medium. We set up a fast flow simulation method by discretizing the matrix so that the discretization nodes of the matrix include the wells where the observed data are made available. Practically, the discretization grid is regular with a characteristic scale equal to the inter-well distance. The discretized matrix superimposed on the channels issues a network of 1D segments in which we solve the classical flow equation: ∇·ðT∇hÞ = 0;

ð9Þ

where T is the segment transmissivity. The boundary conditions are described in Section 4.1. The discretized equation leads to a linear

Fig. 7. (a) Objective function against iteration number during one particular optimization run out of the 50 optimizations performed for the simple configuration of Fig. 6a. Variation in the parameter values of (b) the first identified channel, (c) the second identified channel and (d) the background matrix against the number of iterations.

790

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

system solved by the multifrontal method implemented in the software UMFPACK [13].

3.6. Post-processing analysis of the solution distribution Both because of the existence of local minima and because of the use of a stochastic inversion method, the result of the identification is not necessarily the absolute minimum of the objective function. When running the identification strategy several times with different initial conditions, each solution depends on the initial conditions of the algorithm, represented here by the position of the newly added channels and points within the channels, and on the changes in the configurations made at each step of the simulated annealing algorithm. For a given problem, we obtain a distribution of sub-optimal solutions. We use these solutions to construct the distribution of solutions in the most external loop of the inversion scheme (Fig. 3). We set up a methodology to analyze the similarities and differences of the simulations based on hierarchical clustering algorithms [12,24]. It consists in finding the most common channels and in quantifying their probability of occurrence (Fig. 5). The method described in Appendix B enables a limited number of most probable channels also called meta-channels to be extracted from the channel characteristics and sensitivity.

3.7. Synthesis of parameters of the identification strategy The identification strategy defined in the previous paragraph is controlled by several parameters (Table 1). First, the results of the simulated annealing algorithm depend on the initial conditions, i.e. the initial parameterization and the initialization of the pseudo-random number generator. Second, the balance between the optimization duration and accuracy is controlled by the stopping criteria of the single channel identification loop, i.e. the maximal number of iterations itmax min and the expected objective function value Fobj The refinement level of the final parameterization depends on the stopping criteria of the hierarchical identification loop α1 (Eq. (5)) and α2 (Eq. (6)). The optimization algorithm itself is controlled by the initial temperature, the temperature schedule and the size of the searching neighborhood (Appendix A, 16). The shape of the objective function is a critical factor that depends on the weighting coefficients λj,k of the regularization term in (Eq. (7)). The final clustering algorithm based on the analysis of single solutions is controlled by the number of solutions ns used in the analysis, by the weighting coefficient ξi (Appendix B, 18) and by the stopping threshold dc. The parameters have been assigned values that are either those most commonly used in other studies or those resulting from tests in simple configurations containing one or two straight channels.

4. Results After a synthetic presentation of the tested configurations, we assess the identification strategy on synthetic flow structures of increasing complexity first by using head data only, and second by adding geometrical data.

4.1. Flow configurations and data density We test the identification strategy on a set of fracture flow patterns dominated by large fractures. The configurations are obtained from random fracture networks slightly above the percolation threshold from which we have extracted the elastic backbone. Transmissivities are all equal to 1.0. In these configurations, flow is channeled within major structures and present second-order flow paths, consistently with the definition of hierarchical flow. Three typical configurations of increasing complexity are illustrated in Fig. 6. The simple configuration (Fig. 6a) is composed of two channels, while the intermediary configuration further comprises some secondary flow (Fig. 6b) and the complex configuration has several second-order well-connected structures (Fig. 6c). For each tested configuration, we compute the steady-state head data and extract Nw values from Nw wells regularly distributed over the square domain according to a grid pattern (like in Fig. 4a). These Nw values will be used as the “observed” data di in (Eq. (7)) for the inverse problem. We take the inter-well distance dw as the fundamental length unit. The domain size L is expressed as a function of Nw and dw as L = N1/2 w ·dw. We use Nw = 25 and a typical inter-well distance of the order of 100 m, which corresponds to a domain size of the order of 1 km. The boundary conditions are derived from a fixed regional gradient that is not aligned with the edge of the system, the head gradient being fixed at 1. Instead of (Eq. (7)), results are analyzed according to the squared difference between head data and results Δh2 over the whole domain and obtained with:

2

Δh =

2 1  ∫ h−h′ dS; SS

ð10Þ

where S is the system surface. Practically, this integral is computed by using a 100 × 100-discretization finer than the data density. The advantage of Δh2 compared to RE is that a more precise quantitative measure of the agreement between the observed and resulting heads is obtained.

Fig. 8. Distribution of the squared difference between observed and model results on the whole surface Δh2 (10) derived from the 50 individual solutions for (a) the simple configuration of Fig. 6a, (b) the intermediary configuration of Fig. 6b and (c) the complex configuration of Fig. 6c.

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

791

4.2. Hierarchical inversion with steady-state heads At first, we analyze the behavior of the inversion. Fig. 7a–d display the variation of both the objective function and the parameter values as functions of the iteration number during one particular optimization run out of the 50 optimizations performed for the simple configuration of Fig. 6a. The objective function (Fig. 7a) follows several staircase-like decreases that are typical of simulated annealing algorithms [33,45]. The addition of new parameters through the single and hierarchical channel identification loops triggers successive sharp increases and decreases. As a whole, the channel addition or modification can reduce either significantly or marginally the objective function as shown by the addition and the modification of the second channel (Fig. 7a). The identification strategy stops on a two-channel configuration, showing that the procedure correctly avoids over-parameterization. As expected, the channel and background log-transmissivity values fluctuate much more than the geometrical parameters (Fig. 7b-c) because head data are less sensitive to permeability than flow data. Moreover, transmissivities are only relative to each other, since heads are only sensitive to the difference between channel transmissivities, with the absolute flow being unknown. The channel importance hierarchy is respected since the first identified channel has a transmissivity around three orders of magnitude larger than the next. The background matrix plays efficiently its role of substitution of less important channels (Fig. 7d). It has indeed a high relative transmissivity value (Tmat = 100) when only one channel is involved, and a small value (Tmat b 10− 5) when the optimal number of channels is reached. For the three configurations of Figs. 6, 8 gives the distribution of the squared difference between head data and results Δh2 over the whole domain as derived from the 50 individual solutions. The distribution is given at the intermediary steps containing 1, 2 and 3 channels and at the final stage of the identification (Fig. 8). It shows that, for the simple and intermediary configurations, the addition of one and two channels significantly restricts the occurrence of large Δh2 values. With three channels, most of the models yield head maps that are close to the observed heads (Δh2 b 10− 1) and the distribution of Δh2 is only slightly shifted towards lower values by the introduction of additional channels. For the complex configuration, however, results are less significant and show that the optimal Δh2 value is reached with one or two channels and the resulting Δh2 remains higher (Δh2∼ 10− 0.2). Fig. 9 shows the mean and the standard deviation of the position of the identified channels resulting from the combination of all solutions. The percentage next to the channel indicates the proportion of channels containing this channel in all the configurations resulting from optimization. Channels are colored according to their relative transmissivity. For the simple case (Fig. 9a), the method provides a globally good agreement between the 50 solutions and the initial structure. Both channels are well located and have a transmissivity ratio equal to 2.8 of the same order as the real one equal to 1.0. For the intermediary configuration (Fig. 9b), the map shows the presence of the two main flow channels representing slightly less than 45% of all the determined channels while the remaining channels (∼55%) are not in the true model (Fig. 9b). Finally, for the complex configuration (Fig. 9c), the most transmissive structures are sub-vertical and located in the rightmost part of the domain (Fig. 9c). They identify only partly the right vertical crossing channel while the other channels are not resolved but replaced by an average diagonal channel. The hierarchical identification method yields good results for the simple case in terms of both location and transmissivity. However, the two main channels of the intermediary configuration are resolved but the solution contains two dummy additional channels. The solution obtained for the complex case is far from the true model. To explain these results, we study what we call the head deviation equal to the difference between the observed head and the head computed for a

Fig. 9. Illustrations of the results of the post-processing clustering algorithm for the configurations of Fig. 6 by using head data only. The objective function is defined with σih equal to 1.

homogeneous case having the same boundary conditions. The idea is to damp some of the influence of the boundary conditions and, conversely, to provide graphic illustrations of the information

792

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

Fig. 10. Interpolated maps of the head deviation from the 25 data (a) for the simple configuration of Fig. 6a, (b) for the intermediary configuration of Fig. 6b and (c) for the complex configuration of Fig. 6c.

contained in the head. Fig. 10 displays the “head deviation” for the three tested configurations. In fact, head deviation maps are consistent with the position of most of the correctly identified channels except for the complex configuration. Furthermore, the two dummy additional channels obtained for the intermediary configuration are mainly located in poorly sensitive areas. Finally, as expected, Fig. 10 denotes that head data are less sensitive to each channel when the number of channel increases. These results show a strong influence of the boundary conditions on head sensitivities. 4.3. Use of head deviation To damp some of the influence of the boundary conditions, we introduce the head deviation information in the inversion process. The principle is to use the inverse of head deviations as head uncertainties (σid in (Eq. (7))) to increase the weight of the heads that are more sensitive to the channel structure than to the boundary conditions. We present the same synthetic map as in Fig. 9 by using hydraulic heads and their deviations (σid = 1/“head deviation”) (Fig. 11). For the simple configuration, the position of the channels remains accurate and the ratio between transmissivities is reduced to 1.4, closer to the true value of 1.0 (Fig. 11a). For the intermediary configuration, introducing the head deviation in the objective function moves the dummy channels to the corner of the domain crossing a zone of low sensitivity (Fig. 11b). We also note that the standard deviation of the channel location is inversely correlated to the head deviation. The two channel transmissivities are almost equal with a ratio of 0.9, which is close to the actual value of 1.0. Four of the fifty solutions obtained in this configuration are displayed in Fig. 12. They all contain the most important flow channel and three of the four contain the secondary channel, while embranchments are not resolved. Regarding the complex configuration (Fig. 11c), the result is not improved by the use of head deviation, since the two diagonal structures crossing the domain from the left to the right do not exist in the true model. However, their importance is limited because of their relatively small transmissivity (Fig. 10c). The inversion process yields relevant results for the simple and intermediary configurations of Fig. 6a and b. Head data are sufficient and the use of head deviation precludes the occurrence of channels in the areas where head data are poorly sensitive to the channels. However, the exclusive use of head data finds its limitations in the complex channeling case.

improvements with the addition of basic geometrical information. The objective is to check if the inverse problem can be solved with additional data of a different nature. We deliberately choose to add geometrical data in the same spirit as what has been done when adding a regularization term in more classical inverse problems for continuous permeability structures. An alternative not explored in this paper would be to use transient-state head data. 4.4.1. Geometrical data As basic geometrical data, we use the distance between a well and the nearest non-negligible channel. Although this kind of information is not at first designed to be realistic, some of it could be derived from well test interpretation [11] or from geophysical measurements like GPR [25,53]. Introducing explicitly distances in the objective function (7) leads to: !2  2 # hi −h′i ðnc Þ ri −r i′ ðnc Þ + σir σih i=1 !2 nc −1 NP ðjÞ p ðjÞ−p′k ð jÞ + ∑ ∑ λj;k ⋅ k ; σkp j=1 k=1

Nw   Fobj nc; ; NP ðnc Þ = ∑

"

ð11Þ

where h stands for observed heads, h′ for simulated heads, r for observed distances, r′ for simulated distances, and σih and σir are the uncertainties linked to head and distance data, respectively. The choice of σih and σir balances the relative importance of head and distance data. 4.4.2. Results with distance data alone We first check the possibility to identify the channels with the sole distance data excluding temporarily the head data by choosing σir = 1 and σih = ∞. With distance data only (Fig. 13a), the optimal results of the inversion are made up of a correctly identified vertical structure on the left. The other channels are not correctly identified. To explain this result, we analyze the distance map interpolated from the Nw = 25 data points (Fig. 13b). The dummy channels connect indeed small-distance zones but do not locate existing channels. The exclusive use of distance data is limited by the impossibility to identify the right connectivity of small-distance areas, the lack of sensitivity of distances to the hierarchical channel organization and the absence of information provided on transmissivities.

4.4. Addition of geometrical data Because of the impossibility to identify the channels in the complex structure of Fig. 6c with only hydraulic data, we test the possible

4.4.3. Results with heads and distances When used simultaneously, hydraulic and geometrical information are balanced by choosing σir = 2 and σid = 1/“head deviation”.

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

793

solution is more accurate (Fig. 13c) and only the location of the bottom channel is inaccurate. It may be due to the short distance between the flow structures and the wells at the bottom of the model domain. The diagonal structure already observed with only distance or head data is persistent when using both head and distance data. Its importance is however limited by its transmissivity that is twice to four times smaller than that of the other channels. Finally, the use of additional information represented here by the distance data lowers the dispersion of the retrieved structures. It is first observed in the error bars of the identified channels. It is also apparent from the increase in the cumulative percentage of the identified structures. This increase ranges from 70% in Fig. 11c to 80% in Fig. 13c. 4.4.4. Conclusion For the complex structure of Fig. 6c, the identification yields acceptable results only when using head and distance data simultaneously. Transmissivity still appears to be a regulating factor of the importance of the different structures. Resulting structures that are not present in the reference configuration have indeed a significantly lower transmissivity. 5. Systematic evaluation of the hierarchical identification strategy In this paragraph, we investigate the performance of the identification strategy on a larger set of 20 configurations. We use for the inverse problem the complete set of information, i.e. heads, head deviations and distances. 5.1. Tested configurations We apply the inversion strategy to the 20 configurations of Fig. 15. They are numbered in increasing order of apparent degree of complexity and, from C1 to C20, networks contain gradually more intricately connected channels. All transmissivities are equal to 1.0. The configurations analyzed in the previous sections are C2, C10 and C12. Even if this visual rating is subjective, it is closely consistent with the product of the channeling indicators Dcc and Dic introduced in Le Goc et al. [36] (Table 2). Dcc is a characteristic channel continuity scale taken as the distance over which flow rates remain consistently high. Dic is a characteristic inter-channel distance, with Dic/L N 0.5 denoting a hierarchical flow. Therefore, the channeling degree increases with Dcc·Dic. 5.2. Analyses of results

Fig. 11. Illustrations of the results of the post-processing clustering algorithm for the configurations of Fig. 6 by using head data and head deviation.

Even if distances cannot be used alone, they provide additional information that is different from that provided by hydraulic data. Thus, combining head and distance data yields a much better identification (Fig. 13c). Most of the principal channels are resolved and single solutions contain at most three channels including, for more than half of them, the main channels (Fig. 14). The optimal

The results of the identification strategy are compared to the reference configurations in Fig. 15. The agreement between the results and the true models globally tends to decrease with the apparent complexity as measured by the number of the configuration. This is confirmed for most of the tested configurations by the increase in RE as a function of the configuration complexity as measured by Dcc·Dic/L2 (black squares in Fig. 16). Three more complex configurations (C14, C15, and C20) depart from this tendency with still low RE values (white squares in Fig. 16). These smaller RE values, as compared to those of the other configurations having similar Dic.Dcc/L2 values, are characterized by the absence of the largest head deviations amplified by the squared difference in RE (Eq. (1)). If Dic and Dcc give a global rating of the flow structure complexity and the final value of RE gives a global rate of agreement between the final solution and the reference structure, we are also interested in the quality of the identification of the structural characteristics. Therefore, we use two additional visual indicators for the number of independent channels and for the agreement between the solution and the reference structure. These indicators are clearly subjective. However,

794

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

Fig. 12. Four individual solutions of the inversion method applied to the intermediary complex configuration of Fig. 6b. Fobj refers to the value of the objective function (7). The individual solutions identify the diagonal main flow structure represented in green. Three of them identify the secondary flow structures in red. Embranchments are structures too small to be resolved and are identified only in the third result (blue). The channeled structure is recalled in the background.

Fig. 13. Results for the complex configuration of Fig. 6c. (a) Results of the post-processing clustering algorithm with the objective function defined with distances alone. (b) Interpolated maps of the distances from the 25 data points. (c) Results with the objective function defined using both heads and distances.

for the simple structures handled here, we assume that the variability in their estimation remains limited. The number of independent channels ranges from 1 for C1 to 6 for C19. The visual quality of the agreement between the reference and resulting structures increases from 1 to 5 (Table 2). As previously noticed, the increase in the channel number, and hence in the flow structure complexity, induces on average a decrease in the solution quality (Fig. 17). We identify two different groups. The first group comprises the wellidentified structures characterized by a visual rating greater or equal to 3.

Apart from C13 and C18, it also corresponds to the networks made up of less than 4 independent channels where flow is hierarchical (Dic/L N 0.5). The other configurations including C14, C15, C17, C19, and C20 have more than 4 independent channels and lead to a visual rating smaller or equal to 2. These more complex configurations do not correspond to a hierarchical flow (Dic/L b 0.5), which explains the lack of efficiency of the proposed method. We finally note that the three complex configurations C14, C15 and C20 having a small RE value are among the worse ones in terms of visual rating, consistently with their small Dic/L values.

Fig. 14. Four individual solutions in red and their associated objective function value obtained for the identification strategy with the full objective function applied to the configuration of Fig. 6c. The reference configuration is recalled with gray colors in the background and Fobj refers to the complete objective function (11).

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

795

Fig. 15. Configurations used to test the identification strategy. Left and right columns show the reference configuration and the modeling result given by the post-processing analysis. Green dots stand for well locations.

796

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

Fig. 15 (continued).

Consequently, the mismatch between the resulting and reference heads gives only a very crude estimate of the inversion quality. 6. Discussion We discuss first the results obtained in Sections 4 and 5, and second the extension of the method to more classical transient-state data. 6.1. Discussion of the results To be able to solve the inverse problem with increasingly complex channel structures, we have progressively added more information. The simplest configurations (C1 to C4, Fig. 15) can be identified with the sole knowledge of hydraulic data. The wells surround the channels and the characteristic inter-well distance is smaller than the channel scale, a case for which identification has already been proved possible with only head data and a restricted number of nine wells [49]. For more complex configurations made up of three and four channels (C5 to C13, Fig. 15), the knowledge of the boundary conditions and distance from the well to

the nearest channels are sufficient to determine most of the channel organization. The most complex among the handled structures (C14 to C20, Fig. 15) cannot be identified. A possible explanation is that the data obtained from the 25 piezometers do not contain enough information. Consequently, we have applied the hierarchical identification strategy to 100 wells instead of 25 and found no significant improvement on the identified channels. The limitation may not come from the quantity of data but from the impossibility for the algorithm to identify configurations. In fact, a more detailed analysis shows that configurations C14 to C20 are not extremely channeled since transmissivities are all equal. Because the hierarchical flow organization assumption is not valid, configurations C14 to C20 cannot be identified with the proposed algorithm and should be treated with more classical inverse problem strategies. This shows that the nature of the flow structure between highly channeled and more evenly distributed should be determined a priori. This information might be given by other data like the distribution of flowing fractures or of water inflows in wells. It might also be deduced from a finer analysis of the hydraulic data, head and flow distribution, a possibility that we will explore in the future with transient-state data.

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

797

6.2. Potential extension to more classical transient-state data Model calibration based on more advanced hydraulic and geometrical data has been the topic of several studies. So far, most inverse problem methods seeking discrete flow structures have been performed on more complex data but with a simpler parameterization. Doughty et al. [19] sought IFS parameters using transient-state drawdown data. Day-Lewis et al. [14] determined inter-well connectivity structures also using transient-state drawdown data. Bruyelle et al. [5] worked on a known fracture network structure and restricted the parameterization to the identification of the fracture transmissivities. Gwo [27] determined the occupation probability of links on a discrete structure using transport solute breakthrough data. In this paper, we proceed in two steps. First, we identify complex structures with basic steady-state head and geometrical data, and second we replace these simple data with more advanced geometrical and hydraulic data. Geophysical data will likely not give such precise information as the channel closest to a well but will yield some similar information that should be used to constrain the hydraulic structure close to the wells. Concerning hydraulic data, the objective is to use transient-state data from a more restricted number of wells. The interest in using transient-state data is that they give more information on the flow structure and are less dependent on the boundary conditions than their steady-state counterpart. At not too late times, the hydraulic response is more influenced by the pumping location than by the far away boundary conditions. Pumping systematically in different wells provides a possible varying coverage of the sensitive zones while potentially modifying the main flow structures (e.g. hydraulic tomography [1,30,41]). Since our identification strategy relies on the hierarchical organization of the flow channels, it may lead to different and only partially consistent flow structures for varying pumping locations. Combining them in a consistent discrete transmissivity structure is the first critical issue in shifting from steady-state to transient-state flow data [20]. 6.3. Optimization of the hierarchical identification strategy A key issue in shifting to transient-state simulation is the increased computational demand. Choosing 10 discrete times in 10 pumping tests performed in 10 of the 25 previous wells leads to the necessity to

Table 2 Residual errors RE defined by Eq. (1), channeling indicators (Dic and Dcc) associated to the tested configurations (Fig. 15), threshold of the clustering algorithm (dc) and visual indicators of the channel number and solution quality. Dic/L and Dcc/L are in the interval [0;1]. The number of channel and the quality of the solution range between 0 and 5. Configuration

RE

Dic/L

Dcc/L

dc

Visual number of channels

Visual quality of the solution

C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 C16 C17 C18 C19 C20

0.0014 0.0080 0.0033 0.0080 0.0059 0.0083 0.023 0.0059 0.0080 0.0099 0.020 0.020 0.050 0.012 0.0068 0.029 0.037 0.034 0.044 0.013

0.97 0.65 0.94 0.72 0.85 0.62 0.86 0.79 0.65 0.65 0.64 0.75 0.40 0.34 0.46 0.63 0.42 0.30 0.26 0.43

1.0 0.97 0.98 0.92 0.89 0.88 0.96 0.94 0.92 0.97 0.98 0.78 0.98 0.56 0.65 0.77 0.76 0.71 0.62 0.76

3.5 1.0 1.0 1.5 1.0 2.0 1.5 1.5 0.5 0.7 1.0 0.5 1.0 0.5 1.3 1.0 0.4 1.4 1.0 1.5

1 2 2 2 2 2.5 2 2 2.5 2 2.5 3.5 3.5 4 4 3 3.5 5 6 4

5 5 5 5 4.5 4 4 4 3.5 4 3.5 4 3 1 1 4 2 3 1 1.5

Fig. 16. Quadratic mismatch sum RE between reference data and modeling results as a function of the product of channeling indicators Dic and Dcc normalized by the square of the system size L2. Black squares identify configurations for which RE increases with more complex configurations corresponding to smaller Dic·Dcc/L2 values. White squares stand for three complex configurations resulting in small errors RE.

solve 250 flow problems at specified times. The most appropriate method is the Laplace method with the Stehfest algorithm for the inverse Laplace transform [47]. To achieve a reasonable precision, eight full flow problems need to be solved at each time. In this example, the computational load increases by a factor of 2000 compared to a single steady-state flow simulation. The optimization algorithms should thus be improved to lower the number of flow simulations by one to three orders of magnitude. We argue that this may be possible, first because we took a very simple simulated annealing algorithm and did not optimize its performances, and second because transient-state flow data contain more information than steady-state flow data. The simulated annealing algorithm may be replaced by more advanced optimization methods like the Covariance Matrix Adaptation or Monte-Carlo Markov Chain methods [23,29]. The combination of the different terms of the objective functions may also be improved to minimize the required number of simulations [2,9,42]. The other interest in modifying the objective function is to increase its regularity and use gradient-like optimization methods. The optimization method may combine an initial stochastic search and a subsequent gradient-like scheme like what is done for the identification of permeability zones [50].

Fig. 17. Visual rating of the final post-processing solution versus the visual channel number.

798

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

7. Conclusion In this paper, we present a hierarchical identification method designed to identify both hydraulic properties and hydraulic structures using simple steady-state hydraulic and geometrical data. The hierarchical identification method is derived from the hierarchical organization of flow. The method identifies first the main flow channels before calibrating second-order channels. We test the method on synthetic configurations typical of low-permeability rocks where flows concentrate in highly transmissive flow paths within a restricted number of fractures. Configurations range from a couple of straight channels to configurations comprising up to six interconnected channels. The principal novelty of the proposed approach is to focus on preferential flow paths rather than on equivalent media properties. This is why all other flow settings and optimization algorithms have been taken as simple as possible. The boundary conditions are derived from a uniform head gradient and synthetic data are made up of steady-state hydraulic heads. For the most complex configurations, additional data including distances from the wells to the nearest channel are used. We use these distance data as highly-idealized data possibly provided by near-well geophysical measurements. Channels are identified in decreasing order of importance by using successive optimizations of an objective function. The identification strategy takes advantage of the hierarchical flow organization to restrict the dimension of the solution space of each individual optimization. Because of the successive optimizations, main flow channels strongly determine the characteristics of secondary channels. Additionally, main flow channels can be slightly modified by secondary channels through the introduction of a regularization term in the main channel characteristics in the objective function. The regularization term is weighting the rate of variations of the formerly identified structures compared to the newly added structures. Refinement of main channels is performed first before introducing additional less essential channels. The identification is stopped when both the improvements of the channel structures and the channel number become marginal. A matrix has been introduced in order to replace the channels that should be identified later. The matrix permeability sharply decreases as the identification proceeds. Modifications of the objective function have been introduced to allow marginal modifications of the main flow structures by secondary flow structures. The classical simulated annealing method has been chosen as the optimization algorithm because of the strong non-convex nature of the objective function. For each configuration, the identification strategy has been run 50 times yielding 50 solutions. The prevalent channels are extracted from the 50 solutions by a post-processing algorithm. We conclude that the simplest configurations comprising a couple of straight channels are identified with only the steady-state head data. The identification of similar configurations having complex rather than straight channels requires additionally the knowledge of the head deviations. For intricate interconnected configurations comprising up to three channels with complex shapes, head data are insufficient and the identification requires additional data. In this paper, we have introduced the distances from the wells to the closest channel as highly-idealized geometrical information in the channeling structure. The tested configurations comprising more than four complex channels were not identified even with both hydraulic and geometric data because of the absence of a strong enough hierarchy in the channel structure. Extension of the method to transient-state head data obtained from a more restricted number of wells is the natural outcome of this work. It will require improving the optimization method to lower the number of direct flow simulations by one to two orders of magnitude. Acknowledgements The French National Research Agency ANR is acknowledged for its financial funding through the MOHINI project (ANR-07-VULN-008)

and for its contribution to the development of numerical methods through the MICAS project (ANR-07-CIS7-004). Additional funding was provided by the French Association for Research and Technology ANRT (CIFRE-747/2006). The authors thank Olaf Kolditz and three anonymous reviewers for their helpful and insightful comments. Appendix A. Simulated annealing parameterization Simulated annealing performs a random walk in the parameter space directed towards the minimization of the objective function [34]. The acceptance of parameter sets leading to an increase in the objective function is necessary to get out of local minima. Its probability, however, decreases slowly in order to force the algorithm to converge to the global minimum. The simulated annealing algorithm is parameterized by the acceptance mechanisms of worse parameter sets, by the definition of the random walk in the parameter space and by its condition of termination. The probability of acceptance is handled by an energy criterion characterized by a temperature. The scheduled temperature has been chosen according to Ingber [32]:   1 Tk = Tinit ⋅exp −ci ⋅kD ;

ð12Þ

where Tinit is the initial temperature, D is the dimension of the parameter space, k is the number of accepted states and ci is a userdefined value that can be adapted to improve the algorithm performances, the default value of which is set to: ci = − log 10

−5

  1 : ⋅exp − log100 D

ð13Þ

The acceptance test defines whether a new parameter state is accepted or rejected. It writes:   F −F u b exp − i i−1 ; Ti

ð14Þ

where u is a random value drawn from a uniform distribution in [0;1], Fi is the value of the objective function with the new parameters, Fi − 1 is the value of the objective function at the last accepted state and Ti is the current temperature. Consequently, if the new parameters induce a decrease in the objective function, they are accepted for sure. Otherwise, they are accepted with a probability that depends on how they penalize the solution. This probability decreases as temperature decreases. The random walk within the parameter space is characterized by the parameter modifications from steps i to i + 1. For parameter p j:   j j j j pi = pi−1 + yi ⋅ pmax −pmin ;

ð15Þ

j j and pmax are the minimal and maximal possible values for where pmax j p , and yi is the following scaling factor:

yi = sign ðu−0:5Þ⋅Ti ⋅



 1 2⋅u−1 1+ −1 ; Ti

ð16Þ

where u is a random value drawn from a uniform distribution in [0;1]. With this formulation, the random walk shifts progressively from global to local thanks to the decrease in temperature. The optimization stops when the objective function becomes smaller than a fixed minimal value Fi ≤ Fmin or when the maximal number of iterations allowed is reached.

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

Appendix B. Clustering algorithm The clustering algorithm takes as input a set of channels and groups them together according to the distance between channels, accounting for their relevancy and their sensitivity. The degree of similarity between two channels i and j is given by the distance between them: 2

 2 312 yi ðxÞ−y j ðxÞ dx7 6 dij = 4 ∫ 5 Δ2ij Δ

;

ð17Þ

ij

where i and j are the indices of the two channels, y i(x) and y j(x) characterize the ordinate of the channels corresponding to abscissa x and Δij is the overlapping x-range between the channels. The clustering algorithm gathers the channels in meta-channels by minimizing the distances between their elements. It is controlled by the maximum distance dc between two channels in a meta-channel. Large dc values reduce the number of meta-channels but may lower their internal consistency. On the contrary, small dc values yield larger numbers of meta-channels with higher internal consistency. After investigating the clustering results for several different configurations, we fix dc between 5% and 35% of the characteristic domain scale L. Each meta-channel is then characterized by its mean position and transmissivity, and the uncertainty of its parameters. The mean location of a meta-channel is derived from the position of its components weighted by the coefficients ξi given by: ξi =

Si ; i Fobj

ð18Þ

i is the value of the objective function for the solution where Fobj containing channel i and Si is the weight of the channel within this solution. The transmissivity of the meta-channel is chosen as the geometric mean weighted by the coefficients ξi of its components. Larger weights ξi are obtained for smaller objective function values and for more sensitive channels. Si is more difficult to obtain than the final value of the objective function. First, we rank the channels within a solution by their relative weight. The weight of channel i (Si) is given by the difference in RE values (Eq. (1)) between the whole solution and the whole solution from which channel I has been removed. Second, from a system void of channels, we add progressively the channels within the solution according to their decreasing weight as defined above. We define Si as the difference in RE values before and after the addition of channel i. With this method, the channel with the largest Si is likely to be the main channel. The uncertainty on the meta-channel position is obtained as the standard deviation of its components weighted by the coefficients ξi.

References [1] Brauchler R, Liedl R, Dietrich P. A travel time based hydraulic tomographic approach. Water Resour Res 2003;39:1370. [2] Bruckner G, Handrock-Meyer S, Langmach H. An inverse problem from the 2D groundwater modelling. Inverse Probl 1998;14:835–51. [3] Bruel D. Using the migration of the induced seismicity as a constraint for fractured Hot Dry Rock reservoir modelling. Int J Rock Mech Min Sci 2007;44:1106–17. [4] Bruyelle J, Lange A. Automated characterization of fracture conductivities from well tests inversion, 71st EAGE Conference & Exhibition, 8-11 June 2009. Amsterdam, The Netherlands: Society of Petroleum Engineers; 2009. [5] Bruyelle J, Lange A. An extended evolution strategy for the characterization of fracture conductivities from well tests. Proceedings of the 11th Annual conference on Genetic and evolutionary computation. Québec, Canada: ACM, Montreal; 2009. [6] Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbeau A, Durand E, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation, 1, the flow model. Water Resour Res 1990;26:479–89. [7] Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbeau A, Durand E, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation, 2, the transport model. Water Resour Res 1990;26:491–500.

799

[8] Carrera J, Alcolea A, Medina A, Hidalgo J, Slooten LJ. Inverse problem in hydrogeology. Hydrogeol J 2005;13:206–22. [9] Carrera J, Neuman SP. Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information. Water Resour Res 1986;22. [10] Chang L-F, Sun NZ, Yeh WW-G. Optimal observation network design for parameter structure identification in groundwater modeling. Water Resour Res 2005;41. [11] Chang YC, Yeh HD, Huang YC. Determination of the parameter pattern and values for a one-dimensional multi-zone unconfined aquifer. Hydrogeol J 2008;16:205–14. [12] Chelcea S, Bertrand P, Trousse B. A New Agglomerative 2–3 Hierarchical Clustering Algorithm, Innovations in Classification, Data Science, and Information Systems; 2005. p. 3–10. [13] Davis TA. Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method. ACM Trans Math Softw 2004;30:196–9. [14] Day-Lewis FD, Hsieh PA, Gorelick SM. Identifying fracture-zone geometry using simulated annealing and hydraulic-connection data. Water Resour Res 2000;36: 1707–21. [15] de Dreuzy J-R, Davy P, Bour O. Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 1. Effective connectivity. Water Resour Res 2001;37:2065–78. [16] de Dreuzy J-R, Davy P, Bour O. Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 2. Permeability of networks based on lognormal distribution of apertures. Water Resour Res 2001;37. [17] de Marsily G, Delhomme J-P, Delay F, Buoro A. Four decades of inverse problems in hydrogeology. In: Zhang DW, C.L., editors. Theory, modeling and field investigation in hydrogeology : a special volume in honor of Shlomo P. Neuman's 60th birthday. Colorado, USA: Geological society of America, Boulder; 2000. p. 1–17. [18] Donado LD, Sanchez-Vila X, Ruiz E, Elorza FJ, Bajos C, Vela-Guzman A. Calibration of hydraulic and tracer tests in fractured media represented by a DFN model. In: IAHS W, UK, editors. International Conference on Calibration and Reliability in Groundwater Modelling : From Uncertainty to Decision Making, IAHS, Wallingford, UK, Hague , The Netherlands; 2005. [19] Doughty C, Long JCS, Hestir K, Benson SM. Hydrologic characterization of heterogeneous geologic media with an inverse method based on iterated function systems. Water Resour Res 1994;30:1721–45. [20] Fienen MN, Clemo T, Kitanidis PK. An interactive Bayesian geostatistical inverse protocol for hydraulic tomography. Water Resour Res 2008;44:W00B01. [21] Fogg GE. Groundwater flow and sand body interconnectedness in a thick, multiple aquifer system. Water Resour Res 1986;22:679–94. [22] Francese R, Mazzarini F, Bistacchi A, Morelli G, Pasquare G, Praticelli N, et al. A structural and geophysical approach to the study of fractured aquifers in the Scansano-Magliano in Toscana Ridge, southern Tuscany, Italy. Hydrogeol J 2009;17:1233–46. [23] Fu J, Gómez-Hernández J. A blocking Markov Chain Monte Carlo method for inverse stochastic hydrogeological modeling. Math Geosci 2009;41:105–28. [24] Gordon AD. Hierarchical classifications. Classification, 2nd edition (Monographs on Statistics and Applied Probability, 82). Chapman & Hall / CRC; 1999. p. 69–111. [25] Grasmueck M. 3-D ground-penetrating radar applied to fracture imaging in gneiss. Geophysics 1996;61:1050–64. [26] Grenier C, Bernard-Michel G, Benabderrahmane H. Evaluation of retention properties of a semi-synthetic fractured block from modelling at performance assessment time scales (Aspo Hard Rock Laboratory, Sweden). Hydrogeol J 2009;17:1051–66. [27] Gwo J-P. In search of preferential flow paths in structured porous media using a simple genetic algorithm. Water Resour Res 2001;37. [28] Hanor JS. Effective hydraulic conductivity of fractured clay beds at a hazardous waste landfill; Louisiana Gulf Coast. Water Resour Res 1993;29:3691–8. [29] Hansen N, Ostermeier, A. Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation. IEEE International Conference on Evolutionary Computation 1996, IEEE Press, 1996, p. 312–7. [30] Hao YH, Yeh TCJ, Xian JW, Illman WA, Ando K, Hsu KC, et al. Hydraulic tomography for detecting fracture zone connectivity. Ground Water 2008;46:183–92. [31] Hill MC. The practical use of simplicity in developing groundwater models. Ground Water 2006;44:775–81. [32] Ingber AL. Adaptive Simulated Annealing (ASA). McLean, VA: Lester Ingber Research; 1993. [33] Ingber AL. Simulated annealing — practice versus theory. Math Comput Model 1993;18:29–57. [34] Kirkpatrick S, Gelatt Jr. CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220:671–80. [35] Lavenue M, de Marsily G. Three-dimensional interference test interpretation in a fractured aquifer using the pilot point inverse method. Water Resour Res 2001;37. [36] Le Goc R, de Dreuzy JR, Davy P. Statistical characteristics of flow as indicators of channeling in heterogeneous porous and fractured media. Adv Water Resour 2010;33:257–69. [37] Long JCS, Remer JS, Wilson CR, Witherspoon PA. Porous media equivalents for networks of discontinuous fractures. Water Resour Res 1982;18. [38] Marryote RA, Dougherty DE, Stollar RL. Optimal groundwater management; 2. Application of simulated annealing to a field-scale contamination site. Water Resour Res 1993;29:847–60. [39] Martinez-Landa L, Carrera J. An analysis of hydraulic conductivity scale effects in granite (Full-scale Engineered Barrier Experiment (FEBEX), Grimsel, Switzerland). Water Resour Res 2005;41:13. [40] Mauldon AD, Karasaki K, Martel SJ, Long JCS, Landsfeld M, Mensch A, et al. An inverse technique for developing models for fluid flow in fracture systems using simulated annealing. Water Resour Res 1993;29:3775–89.

800

R. Le Goc et al. / Advances in Water Resources 33 (2010) 782–800

[41] McDermott CI, Sauter M, Liedl R. New experimental techniques for pneumatic tomographical determination of the flow and transport parameters of highly fractured porous rock samples. J Hydrol 2003;278:51–63. [42] McLaughlin D, Townley LR. A reassessment of the groundwater inverse problem. Water Resour Res 1996;32. [43] Nakao S, Najita J, Karasaki K. Hydraulic well testing inversion for modeling fluid flow in fractured rocks using simulated annealing: a case study at Raymond field site, California. J Appl Geophys 2000;45:203–23. [44] Renshaw CE. Estimation of fracture zone geometry from steady-state hydraulic head data using iterative sequential cokriging. Geophys Res Lett 1996;23. [45] Sambridge M, Mosegaard K. Monte Carlo methods in geophysical inverse problems. Rev Geophys 2002;40. [46] Silliman SE. An interpretation of the difference between aperture estimates derived from hydraulic and tracer tests in a single fracture. Water Resour Res 1989;25: 2275–83. [47] Stehfest H. Remark on algorithm 368. Numerical inversion of laplace transforms. Commun ACM 1970;13:47–9. [48] Tarantola, A. Inverse Problems Theory, Methods for Data Fitting and Model Parameter Estimation. Elsevier: Netherlands; 1987. [49] Tiedeman CR, Shieh PA, Christian SB. Characterization of a high-transmissivity zone by well-test analysis: Steady state case. Water Resour Res 1995;31:27–37.

[50] Tsai FT-C, Sun N-Z, Yeh WW-G. A combinatorial optimization scheme for parameter structure identification in ground-water modeling. Ground Water 2003;41:156–69. [51] Tsang C-F, Neretnieks I. Flow channeling in heterogeneous fractured rocks. Rev Geophys 1998;36. [52] Tsang YW, Tsang CF. Flow channeling in a single fracture as a two dimensional strongly heterogeneous permeable medium. Water Resour Res 1989;25:2076–80. [53] Tsoflias GP, Van Gestel JP, Stoffa PL, Blankenship DD, Sen M. Vertical fracture detection by exploiting the polarization properties of ground-penetrating radar signals. Geophysics 2004;69:803–10. [54] Vesselinov VV, Neuman SP, Illman WA. Three-dimensional numerical inversion of pneumatic cross-hole tests in unsaturated fractured tuff 2. Equivalent parameters, highresolution stochastic imaging and scale effects. Water Resour Res 2001;37:3019–41. [55] Yeh WWG, Yoon YS. Aquifer parameter identification with optimum dimension in parameterization. Water Resour Res 1981;17:664–72. [56] Zheng C, Wang P. Parameter structure identification using tabu search and simulated annealing. Adv Water Resour 1996;19:215–24. [57] Zimmerman DA, de Marsily G, Gotway CA, Marietta MG, Axness CL, Beauheim RL, et al. A comparison of seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. Water Resour Res 1998;34:1373–413.