Advances in Water Resources 63 (2014) 179–193
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
An efficient numerical model for hydrodynamic parameterization in 2D fractured dual-porosity media Hassane Fahs a,b, Mohamed Hayek c, Marwan Fahs a,⇑, Anis Younes a a
LHyGes, Université de Strasbourg/EOST, CNRS, 1 Rue Blessig, 67000 Strasbourg, France Lebanese International University, School of Arts and Sciences, 146404 Beirut, Lebanon c AF-Consult Switzerland Ltd, Täfernstrasse 26, CH-5405 Baden, Switzerland b
a r t i c l e
i n f o
Article history: Received 23 April 2013 Received in revised form 21 November 2013 Accepted 21 November 2013 Available online 3 December 2013 Keywords: Dual-porosity media Adaptive multiscale parameterization Refinement indicator Adjoint-state method
a b s t r a c t This paper presents a robust and efficient numerical model for the parameterization of the hydrodynamic in fractured porous media. The developed model is based upon the refinement indicators algorithm for adaptive multi-scale parameterization. For each level of refinement, the Levenberg–Marquardt method is used to minimize the difference between the measured and predicted data that are obtained by solving the direct problem with the mixed finite element method. Sensitivities of state variables with respect to the parameters are calculated by the sensitivity method. The adjoint-state method is used to calculate the local gradients of the objective function necessary for the computation of the refinement indicators. Validity and efficiency of the proposed model are demonstrated by means of several numerical experiments. The developed numerical model provides encouraging results, even for noisy data and/or with a reduced number of measured heads. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Fluids flow through fractured porous media is a process encountered in many areas of the geosciences ranging from groundwater hydrology to oil production. Several conceptual models have been developed to describe the hydrodynamic of fractured porous media [17,19,29,33,39,50,52,55]. An attractive one is the dual porosity model which is firstly introduced by Barenblatt et al. [2] in order to quantify flow in fractured rock. Dual model has been commonly used to describe the preferential movement of water and solutes at the macroscopic scale, the variably-saturated fractured rock formations and structured soils, the leaching of pollutants into soils and groundwater, and the solute transport through clay soils [4,16,21–26,30–35,37,38,43,45–47]. The main idea of the dual-porosity model is to consider at each point of the computational domain two porosities, one associated with the matrix, and the other associated with the fractures system. For each point of space, the model attributes hydraulic properties and pressure heads for both the matrix and the fractures. It assumes a macroscopic Darcy flow over the reservoir for the fractured system and no flow within the matrix system. Mathematically, the model consists of two equations, the mass conservation in the matrix and in the fractures. These equations are coupled by means of an exchange linear term.
⇑ Corresponding author. Tel.: +33 3 68 85 04 48. E-mail address:
[email protected] (M. Fahs). 0309-1708/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.advwatres.2013.11.008
The dual-porosity model introduces different parameters such as hydraulic conductivity, specific storage capacity, and the exchange rate coefficient between the fractures and the matrix. Most of these parameters cannot be measured directly due to prohibitive costs and because the relevant scale of measurement might be unknown or incompatible with the addressed problem. Therefore these parameters should be estimated by model calibration. There are few studies concerning the estimation of the hydraulic parameters of the fractured porous media using the concept of dual-porosity model [7,16,46]. All these studies consider problems involving the well-defined geometrical distribution of the heterogeneity in the computational domain. Hence, they are limited to homogeneous domains or to cases where the identifiable geologic indicators and the hydraulic parameters are clearly correlated. However, when the distribution of the heterogeneity is not prescribed, the problem becomes more complicated. This situation can be encountered in the case where hydraulic parameters do not correlate well with lithology. In this case the parameters should be calculated in each cell of the computational mesh. This calculation requires measurement data in each mesh element, which is practically impossible due to the high cost of the experimental measurements. An alternative technique is the parameterization which consists in reducing the number of unknown parameters using different techniques such as the zonation, the point estimation, the heuristic interpolation function, the pilot points or the conditional simulation. A vast literature exists on these subjects (e.g., review of [9]) and references therein. The
180
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
parameterization by zonation [6,8,20,38,40,49,51] consists in defining sub-domains (zones) in which the parameters are assumed to be constant. By this way the number of parameters can be significantly reduced but the shape and the number of these zones should be a priori known [28]. Thus, the parameterization requires an algorithm (the so-called zonation algorithm) to determine clearly the number of zones and their geometrical distribution. The aim of this work is to develop a robust and efficient numerical model for the parameterization of the hydrodynamic in 2D dualporosity media. Advanced and appropriate techniques are used for each computational task of the developed model. The developed model is based on the adaptive multi-scale method which consists in solving the problem through successive approximations by refining the parameters at the next refiner scale all over the domain and stopping the process when the refinement no longer induces a significant decrease in the objective function any more [3,12,27]. In this work we combine the adaptive multiscale method with the refinement indicators as the zonation algorithm. Refinement indicators do not require significant additional calculations to define an adaptive progressing parameterization but require specific algorithms to define the shape of the zones. The adaptive multi-scale method with refinement indicators has been used for the parameterization of the hydrodynamic in 2D saturated porous media [3,12,27] and 1D unsaturated porous media [28]. It is extended here to treat 2D fractured porous media via the dual porosity model. To do so, special attention should be paid to some supplementary points of difficulty such as the consideration of two state variables (the pressure head in the fractures and the matrix), the treatment of multidimensional parameters, the zonation algorithm for an unstructured triangular mesh and the optimization of the computational efforts. For the optimization procedure, we use the Levenberg–Marquardt algorithm [36,41]. This algorithm interpolates between the Gauss–Newton algorithm and the method of gradient descent. It converges to the solution even if the initial values are very far from the final minimum. The Jacobian and the Hessian matrices (needed for the Levenberg–Marquardt algorithm) are calculated using the sensitivity method [48]. For the calculation of the local gradients of the objective function with respect to the mesh element parameters (necessary for the calculation of the refinement indicators) we use the adjoint-state method [10,48]. This method is more convenient for that calculation than the sensitivity method since it requires the resolution of only one supplementary system at each time step. The adjoint state method is modified here to take into account the two state variables (pressure heads in the matrix and the fractures). For the direct problem, we use the mixed finite element (MFE) method since it is locally conservative and produces accurate and consistent velocity fields even for highly heterogeneous domains [54]. For the zonation we describe an appropriate algorithm for treating the unstructured triangular mesh. The algorithm proceeds by using two different meshes for the spatial discretization and for the refinement indicators. Efficiency and robustness of the developed numerical model are evaluated by studying its capability to provide proper parameterization even in critical cases, such as complex zone distribution, lack of data and/or measurement errors. 2. The direct problem
Sm
@hm ¼ aðhf hm Þ @t
where indexes f and m refer to the fracture and matrix continua, respectively, h ½L is the hydraulic head, S ½L1 the specific storage capacity, a ½L1 T1 the exchange rate coefficient between the fractures and the matrix, f ½M3 T1 a sink-source term per unit volume and q is the macroscopic fluid flux density given by Darcy’s Law:
qf ¼ K f rhf ;
ð3Þ
where K ½LT1 is the hydraulic conductivity (which is a tensor in multidimensional flow). Eq. (2) assumes that flow in the matrix is negligible (term in K m dropped out) and that pumping wells are located in the fractures continuum. 2.2. Discretization The direct problem is solved using the MFE method. The method proceeds by approximating simultaneously the hydraulic head h and the velocity q. It is locally conservative and leads to accurate and consistent evaluation of the velocity field even in highly heterogeneous domains [18,42]. The hybridization technique is often used with the MFE method to reduce the unknowns to the pressure at the edges or faces and to obtain a positive definite linear system [5,11]. To remove non-physical oscillations when the MFE method is used with small time steps, we use the mass lumping formulation introduced in [53]. A detailed review of this method is given by Younes et al. [54]. Note that both rectangular and triangular meshes can be used for the spatial discretization. However the MFE method is developed here for unstructured triangular meshes since they are more suitable for practical problems with complex geometry and local mesh refinement. In the following, we recall the main steps of the lumped formulation of the MFE method for solving the flow equations in dualporosity media. Inside each mesh element E, the velocity is approximated with linear basis functions as follows
qf ¼
3 X Q Ef;j wEj
ð4Þ
j¼1
where Q Ef;j is the flux (in the fracture media) across the edge j of the element E and wEj is the Raviart–Thomas basis function given by
wEi ¼
1 2jEj
x xE;i
z zE;i
Z h E
i 1 E E ðK Ef Þ qf :wEi ¼ hf Thf ;i E
At the Darcy scale, the partial differential equations describing the flow in a dual-porosity medium are written as:
j¼1
ð1Þ
ð6Þ
where hf and K Ef are respectively the fracture hydraulic head and E the fracture conductivity at the element E and Thf ;i is the average fracture head on the edge i of E. Substituting Eq. (4) into Eq. (6), leads to the following matrix form 3 X E E BEij Q Ef;j ¼ hf Thf ;i
@hf ¼ r qf þ aðhm hf Þ þ ff @t
ð5Þ
where jEj is the area of element E and xE;i and zE;i are the coordinates of the vertex i of E (opposed to the edge i). With the MFE method, Darcy’s law (3) is treated separately from the mass conservation equation. The variational formulation of Darcy’s law is given by (see [1,53])
2.1. Mathematical model
Sf
ð2Þ
ð7Þ
R 1 T where the local matrix BEij ¼ E wEi ðK Ef Þ wEj can be expressed analytically for triangles. Superscript T indicates the transposition operator.
181
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
The flux is expressed in terms of the head using the mass lump E across the edge i of the eleing. We firstly approximate the flux Q f ;i ment E under steady state conditions and without sink/source terms by inverting Eq. (7) and using Eq. (1), we get
E Q f ;i
X E E ¼ Nij Thf ;j
ð8Þ
j
where N Eij ¼
detðK Ef Þ T 1 ri K f rj jEj
and ri is the edge vector face to the ver-
tex i of the element E. In the second step of the mass lumping procedure, the accumulation, the sink/source and the exchange terms are distributed over the edges to obtain the total flux. For the edge i of E, the total flux is given by E
Q Ef;i
E jEj SE dThf ;i aE ðThE ThE Þ f E ¼Q m;i f ;i f ;i f f 3 dt
! ð9Þ
Using Eqs. (13) and (14), the matrix form of the direct system can be written as follows
(
nþ1
A Thf
nþ1
I Thm
¼ Fnf
ð15Þ
¼ Fnm
where A is the mass matrix of size Ne Ne (Ne is the number of edges of the computational mesh). I (of size Ne Ne) is the identity nþ1 nþ1 matrix, Thf and Thm (of size Ne) are respectively the vectors representing the hydraulic head in the fracture and the matrix at the time step (n + 1) and Fnf and Fnm (of size Ne) are the vectors of second member. 3. The inverse problem
where SEf and ffE are respectively the fracture storage coefficient and the fracture sink-source term of the element E, aE is the exchange rate coefficient between the fracture and the matrix continua at E the element E and Thm;i is the average matrix head on the edge i of E. Using Eq. (8), the general expression of the flux can be written as follows,
From a practical point of view, the resolution of system (15) requires four types of parameter K f , Sf , Sm and a. However, parameters Sm and a appear as fraction in the system (14). Hence they are strongly correlated and cannot be estimated simultaneously. To avoid this problem, we assume that the storage coefficient Sm is known and we search to estimate the other corresponding parameters. In other words, we propose to estimate the triplet of parameters fK f ; Sf ; ag for a given value of Sm .
Q Ef;i ¼
3.1. Formulation of the inverse problem
! X E E jEj E dThEf;i E E Nij Thj Sf aE ðThm;i Thf ;i Þ ffE : 3 dt j
ð10Þ
Once the fluxes are obtained, they are used to write the continuities of flux between elements E and E’ sharing a common edge i 0 ðQ Ef;i þ Q Ef ;i ¼ 0Þ. Also continuities of the pressure heads E E0 E E0 ðThf ;i ¼ Thf ;i and Thm;i ¼ Thm;i Þ are used to reduce the number of unknowns. This leads to the following equation:
jEjSEf jE0 jSEf þ 3 3
0
!
E 0 X E0 E0 dThf ; i X E E jEjaE jE0 jaE ¼ Nij Thf ; j þ Nij Thf ; j þ þ dt 3 3 j j E
E
ðThm;i Thf ; i Þ þ
jEj E jE0 j E0 f þ f 3 f 3 f
ð11Þ
Now, the time derivative is discretized using an implicit scheme, which gives: 0! E;nþ1 X E E;nþ1 X E0 E0 ;nþ1 jEjSEf jE0 jSEf Thf ;i þ Nij Thf ;j Nij Thf ;j 3 3 Dt j j 0 jEjaE jE0 jaE E;nþ1 E;nþ1 Thm;i Thf ;i þ þ 3 3 0! E;n jEjSEf jE0 jSEf Thf ;i jEj jE0 j E0 þ þ ffE þ f ¼ 3 3 f 3 3 Dt
E;nþ1
Thm;i
E;nþ1
Thf ;i
ð12Þ
ð13Þ
By substituting (13) in (12), we get
jEjSEf
¼
3
þ
jE
0
0 jSEf
!
3
j
E;n 0 E E0 Thf ;i jEjaE jE0 jaE SaE þþSaE0 Dt E;n m m þ Thm;i e Dt 3 3
jEj E jE j E0 f þ f 3 f 3 f
3.2. Optimization of the objective function The optimization procedure is performed using the Levenberg– Marquardt algorithm [36,41]. This algorithm is a modification of the Gauss–Newton method which navigates between the Gauss– Newton algorithm and the method of steepest descent [44]. Combining the robustness of the latter with the computational efficiency of the Gauss–Newton method, the Levenberg–Marquardt algorithm is a highly valued optimization tool for nonlinear problems of least-squares fitting in many engineering fields. This algorithm proceeds to the solution iteratively starting from an initial set of parameters. It consists in adding a positive quantity g to the Hessian diagonal. Hence, the iterative process is given by: T
ð14Þ
nþ1 Thm;i
where N t is the total number of time steps, and are respectively the predicted hydraulic head in the fracture and the matrix on the edge (i) of the computational mesh at the time step nþ1;O nþ1;O (n + 1), Thf ;i and Thm;i are the corresponding measured values nþ1 nþ1 and Xf ;i and Xm;i are the weight given to the observations assumed to be equal to 1 if a measurement is taken at the edge (i) at the time step (n + 1) and 0 elsewhere.
pnew ¼ pold ðJf Xf Jf þ JmT Xm Jm þ gIÞ
0
þ
Ne h i XX 1 Nt1 nþ1 nþ1;O 2 nþ1 nþ1;O 2 Xnþ1 ðThf ;i Thf ;i Þ þ Xnþ1 Þ m;i ðThm;i Thm;i 2 n¼0 i¼1 f ;i
ð16Þ
" # E E0 E 0 E0 a þa 0 Dt 1 jEjSf þ jE jSf E;nþ1 0 E0 E SE þSE ðjEja þ jE ja Þe m m Thf ;i 3 Dt X E E;nþ1 X E0 E0 ;nþ1 Nij Thf ;j N ij Thf ;j j
JðpÞ ¼
nþ1 Thf ;i
where Dt is the time step size and index n refers to the time step level. Using Eq. (2), we can write
aE þaE00 Dt E E E;n E;nþ1 ¼ Thm;i Thf ;i e Sm þSm
In order to estimate the set of parameters fK f ; Sf ; ag, measurements of the state variables must be supplied. We assume that the observed hydraulic heads (in the matrix and in the fractures) are located at the element edges of the computational mesh. For modeling purposes, the objective is to identify K f , Sf and a from a limited number of observations scattered in the field. We use the classical least square error function J as the objective function. For a given set of Np parameters p ¼ ðp1 ; p2 ; . . . ; pNp ÞT , the objective function measures the errors between predicted and observed heads in the fractures and the matrix,
1
rJðpold Þ nþ1 f ;i
ð17Þ
where, Xf and Xm are the weight matrices having X and Xnþ1 m;i (defined in the previous sections) as coefficients, rJ is the gradient
182
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
of the objective function, I is the identity matrix, and Jf and Jm are the sensitivity matrices (also named Jacobian). The coefficients of these matrices are given by: n
n
Jf i;k ¼
@Thf ;i @pk
ð18Þ
n
Jmni;k ¼
@Thm;i @pk
ð19Þ
These derivatives are often called sensitivity coefficients. They represent respectively the sensitivity of the fracture and the matrix hydraulic heads to the parameters. The termination-criterion used for the optimization algorithm is:
krJk < eg
ð20Þ
krJk is the rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2
where
1 krJk ¼ Np
np k¼1
@JðpÞ @pk
norm and
of
the
gradient
defined
by
eg is the tolerance for the gradient.
3.3. Calculation of the sensitivity coefficients Sensitivity coefficients are calculated using the sensitivity method [47]. This method provides an exact calculation of the coefficients and leads to a linear system having the same matrix as the direct system (15). The idea of the sensitivity method is to derive the discrete system (15) with respect to the unknown parameters:
8 nþ1 n > < A @Thf ¼ @A Thnþ1 þ @Ff f @p @p @p k
k
nþ1 > : I @Thm ¼ @Fnm @p @p k
k
ð21Þ
k
The linear system (21) has the same matrix as the direct linear system (15). Hence the computational effort of the optimization procedure can be reduced by using a direct linear solver (with LU factorization), by factorizing the matrix only when solving the direct problem and by applying the solving step that performs only the forward elimination followed by back substitution for the sensitivity system. In this work, all linear systems are solved using the direct solver UMFPCAK [13–15]. This solver is first called to solve the direct system (15). At this level, the solver proceeds to the solution via three steps: (1) An ordering and analysis step to determine a pivot sequence and data structures to perform an efficient factorization, (2) a factorization step that uses the pivot sequence and (3) a solving step that performs forward elimination followed by back substitution. The first and the second steps are known to be largely more expensive in terms of CPU time than step (3). Then the solver is called again to solve the sensitivity system. Here, steps (1) and (2) are circumvented and only step (3) is executed.
of zone and allow merging and/or splitting of zones during the inverse procedure, they are known to be difficult to implement and very time consuming in CPU time since they require the resolution of a hyperbolic equation on the computational grid of the direct problem [28]. The refinement indicators, which are used in this work, might be more efficient in computational effort. They do not require any supplementary calculations to define the adaptive progressing parameterization but they require a specific algorithm to determine the shape of the zones [28]. 4.1. Refinement indicator The role of the refinement indicator is to determine the location of the cutting that allows the best decrease in the objective function. The refinement indicators were first developed based upon the optimal conditions [3,12,30]. Recently, Hayek et al. [28] give a new mathematical development to define the refinement indicator and extend it to the case of multi-dimensional parameters. In this paper, we adopt this formulation to define the global refinement indicator corresponding to the set of parameters {K f ; Sf ; a}. The formulation developed by Hayek et al. [28] consists in calculating the indicator of each parameter individually and then calculating a total indicator based on the summation of the individual parameter indicators. For simplicity, the method is explained by considering that one parameter (K f for example) is to be estimated. Generally, the refinement indicators are calculated after the optimization procedure. Let us note K f and J K f the optimal parameter and objective function corresponding to an initial parameterization which is a homogeneous domain (Fig. 1a). If J K f is small enough, the current parameterization is accepted. On the contrary, an additional degree of freedom (additional parameter) is necessary and is obtained by dividing the domain into two zones (Z 1 and Z 2 ). The main problem now is to find the optimal location of the boundary between these zones. This can be achieved by calculating the refinement indicators for a certain number of cutting configurations usually prescribed by the user. Let us consider the cutting configuration of Fig. 1b. The refinement indicator corresponding to this configuration is given by:
I‘K f
X @JðK Þ X @JðK Þ f f ¼ ¼ i2Z @K f ;i i2Z @K f ;i 1
ð22Þ
2
Superscript ð‘Þ denotes the cutting configuration of Fig. 1b. In Appendix A we gives a detailed development of the refinement indicator defined by Eq. (22). By definition, the refinement indicator does not require the resolution of the optimization problem corresponding to the current tentative parameterization. It can be computed from the partial derivatives of the objective function with respect to the local parameters (the parameter in each mesh element) in the zone. It depends only on the solution K f of the current parameterization.
4. Parameterization As discussed in the introduction, the adaptive multiscale zonation method is used in this work for the parameterization. With this method, the parameters are assumed to be constant per zone. The number of zones is increased during the calculation until the difference between the predicted and measured data becomes invariant. The zones are refined only in the specific regions depending on the sensitivity of the state variable. A specific procedure should be used to identify the zones splitting that allow the maximum decrease in the objective function. Two appropriate procedures are developed in the literature. They are based on the level-set methods [6,38,40] and the refinement indicators [3,12,27,28], respectively. Despite the fact that the level-set methods can handle any shape
a
b
Fig. 1. (a) Initial parameterization, (b) parameterization associated to a tentative parameter discontinuity.
183
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
In practice, and as defined in [28], a dimensionless refinement indicator is implemented in the numerical model. This dimensionless refinement indicator is defined as follows:
~I‘ ¼ Kf
I‘K f
ð23Þ
max8‘ ðI‘K f Þ
where max8‘ ðI‘K f Þ represents the largest value of the indicator in all zone configurations. Similar to K f we define the dimensionless refinement indicators ~I‘Sf and ~I‘a corresponding to parameters Sf and a, respectively. The zonation is assumed to be similar for the triplet of parameters (K f , Sf and a). Consequently, these parameters are discontinuous at the same locations of the domain. Therefore, it is coherent to use an overall refinement indicator in order to define the interfaces between zones. The global refinement indicator corresponding to the triplet of parameters fK f ; Sf ; ag is defined as the summation of ~I‘K f , ~I‘Sf and ~I‘a :
~k‘ ¼ ~I‘ þ ~I‘ þ ~I‘ Kf Sf a
ð24Þ
4.2. Refinement strategy Once the location of the best cut is defined, infinite possibilities to refine a zone can be deduced. Numerically, we are unable to test all these possibilities, and a strategy of refinement must be adopted. Ben Ameur et al. [3] and Hayek and Ackerer [27] used the so-called vertical and horizontal refinement indicators in order to define interfaces between zones. In their works, the mesh used for solving the direct problem is rectangular and refinement indicators are supported by the edges of the mesh elements themselves. Consequently, interfaces between zones are limited to be verticals or horizontals. In this work, we extend the study to the case of an irregular mesh (here triangular). Two different mesh grids are used for the spatial discretization and for the refinement indicators. We use vertical and horizontal indicators for simplicity but other geometries can be handled. Note that, the term vertical (resp. horizontal) indicator is used to denote the scalar value of an indicator corresponding to a vertical (resp. horizontal) cutting configuration of the domain. Vertical indicators are numbered starting from the bottom to the top and horizontal indicators are numbered starting from down to up. The distance between two vertical (resp. horizontal) indicators is defined by the user. Fig. 2a represents a grid of refinement indicators composed of four vertical and four horizontal indicators overlapped with the mesh used for the spatial discretization. Now, let us explain how the refinement indicators can be used in order to construct the interface between two zones. Let us take, for example the third vertical refinement indicator (see Fig. 2b). This indicator divides the domain into two zones left and right. It intersects six mesh elements. The method requires that each zone is composed of a certain number of mesh elements. Some parts of
a H3 H2 H1 V3
V4
The calculation of the refinement indicators requires the values of the gradients of the objective function with respect to the parameter values in each mesh element (see Eq. (24)). These gradients are often called local gradients. They can be computed using the sensitivity method. However, this method requires the resolution of Np Nm Nt linear systems which is very expensive in CPU time (Nm is the number of elements in the computational mesh). An alternative is to use the adjoint-state method [10,49] that calculates the local gradient by solving only one supplementary linear system (the adjoint system) per time step (i.e., Nt supplementary linear systems). In this work the adjoint-state method is modified by using two Lagrange multipliers in order to take into account the two state variables (hydraulic head in the fracture and the matrix). The main stages of the adjoint-state method are recalled here. Consider the following optimization problem:
8 > Find p 2 RNp that minimizes JðpÞ under constraints : > < nþ1 A Thf Fnf ¼ 0 > > : nþ1 I Thm Fnm ¼ 0
ð25Þ
The Lagrangian operator associated with this optimization problem is given by: 1
Nt
1
Nt
1 Nt Lðp; Thf ; . . . ; Thf ; Thm ; . . . ; Thm ; c1f ; . . . ; cNt f ; cm ; . . . ; cm Þ
¼ JðpÞ þ
Nt1 X
nþ1
hA Thf
Fnf ; cnþ1 iRNe þ f
n¼0
Nt1 X
nþ1
hI Thm Fnm ; cnþ1 m iRNe
n¼0
ð26Þ nþ1 f
nþ1 m
where c and c are two vectors representing respectively the two Lagrange multipliers at the time step (n + 1) and h; iRNe is the 1 Nt 1 Nt scalar product in RNe . In Eq. (26) p, Thf ; . . . ; Thf and Thm ; . . . ; Thm are considered as independent variables but are related by the constraints in Eq. (25). The Lagrange multipliers are solutions of the adjoint-system:
8 Nt1 X > nþ1 nþ1 > @L > ¼ 0 8dThf nþ1 ; dThf > @Thf < Ne R n¼0
Nt1 E > XD > nþ1 > @L > : nþ1 dThm @Th n¼0
m
ð27Þ
nþ1
RNe
¼ 0 8dThm
8 A cnf ¼ cnþ1 :cnþ1 þ pnþ1 cnþ1 > m f > > > > < Xnþ1 ðThnþ1 ThO;nþ1 Þ f f f for n ¼ Nt 1; . . . ; 0 > þ Knþ1 cnþ1 > I cnm ¼ Hnþ1 cnþ1 m > f > > : nþ1 O;nþ1 Xnþ1 Þ m ðThm Thm
H4
V2
5. Adjoint-state method
which is discretized at each time step in:
b
V1
the mesh elements which intersect the indicator must be gathered or removed from the left zone. The idea is as follows: if the center of mass of the mesh element is located at the left of the indicator then the entire element is gathered to the left zone. If not, then the element is removed from the left and gathered to the right zone. Consequently, for this example, the third indicator gives the interface plotted in Fig. 2b. The same strategy is used for the horizontal indicators.
V3
Fig. 2. (a) Vertical and horizontal indicators, (b) interface between zones corresponds to a vertical indicator.
ð28Þ
where Xnþ1 and Xnþ1 are diagonal matrices whose diagonal elef m nþ1 ments are Xnþ1 and Xnþ1 , pnþ1 , Hnþ1 and f ;i m;i and the matrices c nþ1 K have respectively the following coefficients:
Cnþ1 ¼ i;j
@F nþ1 f ;i nþ1
@Thf ;j
ð29Þ
184
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
Pnþ1 ¼ i;j
Hnþ1 ¼ i;j
@F nþ1 m;i
ð30Þ
nþ1 @Thf ;j
@F nþ1 m;i
ð31Þ
nþ1
@Thm;j @F nþ1 f ;i
Knþ1 ¼ i;j
ð32Þ
nþ1
@Thm;j
The system (28) is retrograde in time. It is solved from n ¼ Nt 1 to Nt n ¼ 0 with initial conditions data: cNt f ¼ cm ¼ 0. The adjoint system requires as much resolution of linear systems as the number of time steps. Note that, system (28) has the same matrix as the direct system (15). Hence, the use of a direct solver, allows saving the computational effort by avoiding matrix factorization as in the sensitivity system. The direct solver also allows avoiding problems that may arise when solving adjoint equations with iterative solvers due to the requirement of tight tolerances. Once the Lagrange multipliers are obtained, they can be used to calculate the local gradient of the objective function. Let us take as an example the local gradient of J with respect to pik00 which represents the value of the parameter ðpk0 Þ at the mesh element ði0 Þ. This local gradient can be obtained using the Lagrangian differential which writes:
dL ¼
Np Nm X X @L
dpik þ
@pik * N t1 X @L
i¼1 k¼1
þ
nþ1
n¼0
@Thm
Nt1 X
*
n¼0
@L
+ nþ1
; dThf nþ1 @Thf +
RNe
nþ1
ð33Þ
; dThm
RNe
where pik is the value of the parameter pk on the mesh element ðiÞ. Using Eq.(27) and taking dpik ¼ 0 for i–i0 and k–k0 , Eq. (33) becomes:
dL dpik00
¼
@L
ð34Þ
@pik00
dJ
¼
@L
ð35Þ
@pik00
Consequently, the gradient of the objective function with respect to a parameter value in certain mesh is equal to the partial derivative of the Lagrangian operator with respect to this value. This can be written as follows:
@JðpÞ @L ¼ @pk;i @pk;i
This section is devoted to the algorithm of the developed numerical model. This algorithm allows the triplet of parameters fK f ; Sf ; ag and the spatial distribution of their zones of heterogeneity to be determined by minimizing the difference between the predicted and calculated data. We assume that the total maximum number of zones is known and we denote it by nz. Using appropriate grids for the refinement indicators and the spatial discretization, the algorithm proceeds as follows: 1. Choose an initial parameterization P 1 (usually corresponds to the homogeneous case) and an initial vector of parameters p1 . 2. Minimize the objective function and compute the optimal solution p1; using the procedure described in Section 3.2. 3. For k ¼ 1; nz 1 do Denote by P k a parameterization containing k zones and by pk the corresponding parameters. For each zone of P k do a. Compute all global refinement indicators ~ k‘k using the k‘k repoptimal parameters ðpk; Þ as in Eq. (24). Note that ~ resents the restriction of the refinement indicator ~ k‘ for the zone (k). End do b. Let kmax ¼ maxk;‘ kk‘ (the largest value of the global indicators in all zones). c. Select a number of partitions corresponding to indicators which have values kk‘ P h kmax . Here h is a user defined value. It is fixed to be 0.6 in this study (see Remark below). For each selected partition d. Generate a new parameterization with k þ 1 zones. e. Minimize the objective function for this parameterization and calculate the optimal parameters pkþ1; . End do f. Retain the parameterization corresponding to the smallest objective function for the next iteration. End do
Since p, Thf and Thm are considered as independent variables and using Eq. (26) we obtain:
dpik00
6. Overall algorithm of the numerical model
8i ¼ 1; . . . ; Nm and k ¼ 1; . . . ; Np
4. The current parameterization now contains nz zones. It is denoted by P nz and the corresponding set of optimal parameters are pnz; . If the stopping criterion is reached then the solution has been found. 5. If not the global refinement indicators ð~ k‘k Þ are calculated for each zone of the Pnz parameterization using the optimal parameters pnz; .
Z1
ð36Þ
Z2 Z4
Using Eq. (26) we obtain:
Z1
* + t1 @Fnf nþ1 @JðpÞ NX @A nþ1 ¼ :Th ; c f @pik @pik @pik f n¼0 RNe * + n Nt1 X @F þ mi ; cnþ1 @pk m Ne n¼0
Z3
Z21
Z22
Z4
Z3
Z1
ð37Þ
R
Note that in Eq. (37) all derivatives can be calculated analytically since matrix A and vectors Fnf and Fnm are explicitly expressed in terms of the mesh parameters pik .
Z3
Z2 Z4
Fig. 3. Description of stage (5) of the algorithm for n ¼ 4.
185
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
9. Select the partition which corresponds to the smallest objective function. 10. If the stopping criterion is reached then the solution has been found. If not, we return to stage (4). Three stopping criteria are used in the above algorithm: (i) The objective function and the norm of the gradient are less than the prescribed tolerances.
krJk < eg and Jðpkþ1 Þ < eaJ where
ð38Þ
eaJ is the absolute tolerance for the objective function.
(ii) One of the already tested zone configurations is repeated again
jJðpkþ1 Þ Jðpk0 Þj < eaJ k0 ¼ 0; . . . ; k Fig. 4. Schematic representation of the computational mesh and boundary conditions of the 2D domain.
Fig. 5. The exact three-zone parameterization, solution of the inverse problem.
6. Calculate the maximum indicator for all cutting configurations kmax (as in steps (b) and (c)) and retain the partitions corresponding to indicators which have kk‘ P h kmax for the next step. 7. At this stage, the selected partitions now contain ðnz þ 1Þ zones (which is not accepted since the number of zones is imposed to be equal to nz). To avoid this problem two zones are gathered as in the remark below. 8. Minimize the objective function for all possible partitions (of nz zones) deriving from the partition of nz þ 1 zones (see remark below).
ð39Þ
In this case, parameterization Pkþ1 and Pk0 are similar, then they will be progressively repeated and the objective function cannot be decreased any more. In this case, the parameterization corresponding to the minimum value of the objective function is considered as the best possible solution. Note that, this criterion is fixed empirically after several numerical observations. (iii) The maximum number of iteration is reached. Remark. Note that, logically in steps (3.c) and (5) only the maximum indicator should be selected for the new parameterization. However here we select all indicators having a value larger than h kmax since refinement indicators are defined at the first order and do not provide exact and complete information. Based on trial and error numerical experiments the value of h is fixed to be 0.6. In stage (7), since we seek an nz-parameterization (parameterization contains nz zones), we suggest a way to gather three zones into two zones. Fig. 3 describes stage (7) of the algorithm for nz ¼ 4. The indicator (dashed line) divides zone Z 2 into two sub-zones Z 21 and Z 22 , then the current partition contains five zones. This partition leads to two partitions of four zones each. The first partition is obtained by gathering Z 21 to Z 4 and the second one is obtained by gathering Z 22 to Z 4 . We note that the two new partitions conserve the discontinuity of the parameters at the indicator location. 7. Numerical examples The applicability and the efficiency of the developed numerical model are evaluated using numerical experiments based upon a theoretical test case. The test case considered here is a square domain of side 200 m (Fig. 4). Dirichlet boundary conditions with constant head hf ¼ 1 m are imposed at the boundaries (Fig. 4). A 1 pumping well of flow rate Q ¼ 50 m3 h is located at the center. The duration of the simulation is 5000 s. The domain is assumed to be heterogeneous with a complicated zones distribution
Table 1 Results of the inverse problem with known parameterization. Parameters Exact solution
Initial solution
Optimal solution
Zone 1
Zone 2
Zone 3
JðpÞ
Kf
1:0 10
8:0 10
4:0 10
Sf
1:0 106
7:0 106
3:5 106
a
5:0 1010
1:0 1010
3:0 1010
Kf
5
5
5
1:0 10
6
6
Sf
1:0 10
5
1:0 106
5
5
1:0 10
1:0 10
a
1:0 1011
1:0 1011
1:0 1011
Kf
0:999 105
8:0 105
4:0 105
Sf
1:0 106
7:0 106
3:499 106
a
4:999 1010
0:999 1010
3:0 1010
1:0 10
0:0
8:24 108
3:07 1020
186
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
(Fig. 5). Parameters of the zones are given in Table 1 (exact solution). The parameter Sm is imposed to be 3 105 m1 in the whole domain. The test case is simulated by assuming that the distribution of the zones and the values of the corresponding parameters are unknowns. Measurement data are obtained by solving the direct problem using the exact parameters. The time step size is fixed to be 500 s. For the spatial discretization we use an unstructured triangular mesh with 388 elements (Fig. 5). Note that the magnitude of the hydraulic head for all studied problems is of the order of tens of meters and the used tolerances are as following:
Fig. 6. Evolution of the objective function as a function of iteration level in the developed inverse procedure.
(i) The tolerance for the gradients is eg ¼ 103 ; (ii) The absolute tolerance for the objective function is eaJ ¼ 1010 ; (iii) The maximum number of zonation attempts is fixed to 10.
Fig. 7. Refinement indicator values: (a) for the first iteration, (b) for the second iteration, (c) for the third iteration, (d) for the fourth iteration, (e) for the fifth iteration and (f) for the sixth iteration.
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
187
Fig. 8. The optimal refinement process required to recover the three-zone parameterization. (a) Initial parameterization and locations of refinement indicators selected for the first iteration. (b), (c), (d), (e), and (f) represent cuts retained at the first, second, third, fourth and fifth iteration, respectively, and locations of refinement indicators selected for the second, third, fourth, fifth and sixth iteration.
7.1. Parameters estimation At first, the convergence of the optimization procedure of the numerical model is studied by assuming that the distribution of
the zones is given and only the values of the parameters should be estimated. Full measurement data are considered in this case (overall the computational mesh). Arbitrarily, half of the observed points are used for the measurement in the fracture while the
188
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
Table 2 Results obtained by the algorithm at the first five iterations. Parameter
Zone 1
Zone 2
First iteration Kf
4:274 106
3:912 105
7
3:284 106
Sf
9:168 10
a
5:655 1010
2:665 1010 25:62 104
JðpÞ Second iteration Kf
Zone 3
1:086 105
4:563 105
1:091 105
7
6
4:080 10
3:736 107
3:353 1010
5:069 1010
Sf
4:362 10
a
5:024 1010
JðpÞ
3:85 10
4
Third iteration Kf
1:626 106
6:856 105
1:208 105
Sf
2:034 107
7:315 106
6:779 107
a
3:048 1010
3:238 1010
3:261 1010
JðpÞ
4:04 104
Fourth iteration Kf
9:409 106
7:178 105
3:635 105
Sf
9:372 107
6:853 106
1:018 106
a
4:953 10
10
3:548 10
JðpÞ
10
3:028 1010
3:68 104
Fifth iteration Kf
9:629 106
6:818 105
3:896 105
Sf
1:044 106
7:483 106
3:221 106
a
4:994 10
10
2:141 10
10
3:015 1010
0:51 104
JðpÞ
Fig. 9. Schematic representing the exact three-zone parameterization and the space-distribution of 50 measurement points; symbols ( ) indicate locations of measurements points.
other points are used for the measurement in the matrix. The pressure in the fracture and the matrix are measured for all time steps. The total number of parameters to be estimated is nine. The parameters used as initial values for the optimization procedure are given in Table 1 (initial solution). These values are selected to be as far as possible from the optimal solution. The behavior of the proposed inverse procedure is depicted in Fig. 6. This figure represents the evolution of the objective function versus the iteration level. It shows a relatively quick convergence of the inverse procedure. The obtained estimated parameters are close to the exact parameters Table 1 (optimal solution). Various initial parameter sets were tested and the results show a good convergence of the optimization procedure. 7.2. Parameters estimation and zones identification In this section, we discuss the accuracy and the efficiency of the developed numerical model. To do so we consider the previous test case and we assume that the distribution of the zones and their corresponding parameters are unknowns. Thus, we aim here to find the zonation and the parameters by simulating the test case with the developed numerical model. The measurement data are
the same as the previous section. For the simulation, we use a set of refinement indicators consisting of nine vertical and nine horizontal indicators (regular mesh). The initial parameterization consists in 1 zone with the following initial parameter set: K f ¼ 106 m s1 , Sf ¼ 105 m1 , a ¼ 1011 m s1 . After minimizing the objective function, the optimal solution for this parameterization is K f ¼ 2:985 105 m s1 , Sf ¼ 2:620 106 m1 , 11 1 a ¼ 4:113 10 m s and the value of the optimal objective function is 73:75 104 . Therefore, we have 18 possible locations of discontinuity. At this stage, we compute the global refinement indicators corresponding to these discontinuities by using Eq. (24). The global refinement indicator values are represented in Fig. 7a. The six selected discontinuities (V5, V6, V7, H3, H4, H5) having indicators greater than 0:6k1max are plotted in Fig. 8a (vertical and horizontal lines). The optimization is performed for these six discontinuities and the smallest value of the objective function is obtained for V7 (solid line in Fig. 8a). The parameterization corresponding to V7 is retained for the first iteration. It is composed of two zones (Fig. 8b). The parameter values in each zone and the value of the corresponding objective function are represented in Table 2. For the second iteration, the computation of the refinement indicators is based on the solution obtained previously at the first iteration. The global refinement indicator values in both zones are plotted in Fig. 7b. For this iteration, nine refinement indicators have values greater than 0:6k2max V8 in zone 1 and (V4, V5, V6, H3, H4, H5, H6 and H7) in zone 2 (see Figs. 7b and 8b). The parameterization corresponding to H3 provides the smallest objective function (J ¼ 3:85 104 ), and is retained for the second iteration (see Fig. 7c and Table 2 for the corresponding parameters). Since the optimal objective function does not reach the prescribed
Table 3 Comparison between the exact solution and the numerical solution obtained by the algorithm at the sixth iteration (the optimal solution). Parameter Exact solution
Optimal solution
Zone 1
Zone 2
Zone 3
JðpÞ
Kf
1:0 10
Sf
1:0 106
7:0 106
3:5 106
a
5:0 1010
1:0 1010
3:0 1010
Kf
0:999 105
8:0 105
4:0 105
Sf
1:0 106
7:0 106
3:499 106
a
4:999 1010
0:999 1010
3:0 1010
5
8:0 10
5
5
4:0 10
0:0
3:07 1020
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
189
Fig. 10. Cuts retained at: (a) the first iteration, (b) the second iteration, (c) the third iteration, (d) the fourth iteration, (e) the fifth iteration and (f) the sixth iteration, respectively, for the case of noisy data.
minimum, the next iteration is run by computing all global refinement indicators in each zone of the three-zone parameterization retained for the second iteration. For the third iteration, six refinement indicators are selected. All these indicators are in zone 2 (V3, V4, V5, H5, H6 and H7) (see Figs. 7c and 8c). Each indicator gives two parameterizations of three zones each. Then, the objective function is minimized for
the corresponding 12 parameterizations. The retained parameterization at the third iteration corresponding to V5 is represented in Fig. 8d. Again, the optimal objective function is greater than the fixed minimum (see Table 2), and the procedure is repeated three more times. Fig. 8d, e, and f represent the parameterization retained at the fourth, fifth, and sixth iterations, respectively. Fig. 7d, e, and f represent the values of the corresponding
190
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
Table 4 Comparison between the exact solution and the numerical solution obtained by the algorithm for the case of noisy data. Zone 2
Zone 3
JðpÞ
Exact solution K f
Parameter Zone 1 1:0 105
8:0 105
4:0 105
0:0
Sf
1:0 106
7:0 106
3:5 106
a
5:0 1010
1:0 1010
3:0 1010
Kf
0:999 105 8:0 105
Sf
0:999 106 6:998 106 3:5 106
Optimal solution
a
10
5:015 10
0:993 10
0:476
4:0 105 10
Table 5 Comparison between the exact solution and the numerical solution obtained by the algorithm for the case of lack of data.
Exact solution
Optimal solution
10
3:004 10
Parameter Zone 1
Zone 2
Zone 3
JðpÞ
Kf
1:0 105
8:0 105
4:0 105
0:0
Sf
1:0 106
7:0 106
3:5 106
a
5:0 1010
1:0 1010
3:0 1010
Kf
0:992 105 5:984 105 4:032 105 196:157
Sf
0:942 106 4:588 106 3:555 106
a
4:992 1010 3:013 1010 3:013 1010
based on the minimum value of the objective function is reached and the iterative procedure is stopped. The zonation is well identified and the parameter set for each zone is estimated with high accuracy. The solution obtained by the algorithm is the exact solution in terms of parameter values and zone distributions. Note that, in this case the refinement procedure is terminated by the first stopping criterion of Eq. (40). This test case is also used to study the efficiency of the numerical model in terms of CPU time. To do so, two runs have been performed. In the first one, the LU decomposition (in the linear system solver) is applied for the direct, the sensitivity and the adjoint system, respectively. In the second run, the LU decomposition is applied only for the direct system and it is retained for the resolution of the sensitivity and the adjoint systems. The CPU time of the two runs are estimated to be 6800 s and 4100 s, respectively. The direct solver implemented reduces the computational effort of the developed numerical model by approximately 40%.
Fig. 11. Hydraulic head profile in the fractures and distribution of 20 measurement points at locations of high head values.
7.3. Effect of number of data points The robustness of the algorithm was examined with respect to the number of data points. Four simulations were realized using respectively 300, 200, 100 and 50 data points. In all these simulations, the data-point locations were arbitrarily selected while respecting that each zone contains at least one data point (Fig. 9). Also arbitrarily, half of the observed points are used for the measurement in the fracture while the other points are used for measurement in the matrix. The pressures in the fracture and the matrix were measured every 1000 s. For the four simulations, the numerical solution converged to the exact parameterization using the same iterations as in the case of the total data. However, when the number of data points became smaller than 50, the algorithm became unable to detect the exact parameterization and therefore special attention should be focused here on the location of the data points. This problem is discussed further in Section 7.5. 7.4. Effect of measurement errors
Fig. 12. The optimal parameterization corresponding to the case of lack of data in zone 2.
refinement indicators, respectively. The corresponding parameters and optimal objective functions are presented in Table 2. At the sixth iteration, one of the parameterization corresponding to indicator V4 in zone 2 (see Fig. 8f) gives the smallest optimal objective function. The value of the optimal objective function is very close to 0 (see Table 3). The stopping criterion of the parameterization
We tested the behavior of the algorithm with respect to noise on the measured data. We considered the case depicted in Fig. 9 (50 data points). Measurements are taken every 1000 s as in the previous case. We assumed that data include uncorrelated (space and time) noise and we tried to reconstruct this parameterization. Measurement errors were simulated by adding a random numnþ1 nþ1 ber to the exact values Thf ;i ðpÞ and Thm;i ðpÞ: nþ1;O
¼ Thf ;i ðpÞ þ ef mnþ1 f ;i
nþ1;O
¼ Thm;i ðpÞ þ em mnþ1 m;i
Thf ;i
Thm;i
where
nþ1 nþ1
ð40Þ
nþ1 ef ¼ em ¼ 0:2 and mnþ1 f ;i ; mm;i 2 1; þ1½ are random numbers.
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
191
Fig. 13. Cuts retained at: (a) the first iteration, (b) the second iteration, (c) the third iteration, (d) the fourth iteration for the case of central inclusion.
In this case, at the optimum, the difference between calculated and measured data should satisfy the following relation:
JðpÞ <
1 f 2 2 ½N e þ Nm d em 2 d f
ð41Þ
where Nfd and N m d are respectively the number of measured data in fracture and matrix. Here, the total number of measured data is 250. Hence the corresponding criterion is fixed to be 10. Fig. 10 shows the successive parameterizations obtained by the algorithm. The exact parameterization was detected after six iterations as in the case of one data point per element and without measurement errors. However, the selected parameterizations are not similar for all iterations. Indeed, the parameterizations retained for the first and second iterations are identical in the two cases but those retained for the third, fourth and fifth iterations are different (see Figs. 8 and 10). The results concerning the parameters are given in Table 4. 7.5. Locations of data points As discussed previously, when the number of data points becomes less than 50 and if these points are distributed arbitrary in the domain, the algorithm does not converge and is unable to detect the exact parameterization. In this case more attention should be paid to the location of the measurement points. Indeed, measurement points located in the regions where there is a non-significant change of the head cannot provide sufficient
information for the parameterization. In our case this may be occur on the measurement points that are too far from the source. The objective of this part is to show the influence of the measurement point locations on the convergence of the algorithm. To do so we consider the same test case and we realize a simulation with 20 measurement points distributed arbitrarily in the domain. In this case, the algorithm appears to be unable to converge. We repeat the simulation with the same number of measurement points, but we distribute these points as near as possible to the source while considering that each zone contains at least one measurement point. The distribution of the measurement points is illustrated in Fig. 11. This figure shows also the hydraulic head distribution obtained with the exact parameters. The pressures in fracture and matrix are measured every 1000 s. Odd (resp. even) numbered points are used for the measurement in the fractures (resp. matrix). The simulation is performed here without measurement errors. The results show that, after six iterations, the exact parameterization is found and the estimated parameters are close to the exact ones. We note that the parameterizations selected during the iterations are exactly the same as in the case of full data points. 7.6. Lack of data In this section we study the behavior of the suggested numerical model in the case of lack of data. Indeed, in all the previous examples, measurement points were distributed in such a manner that each zone contains at least one data point. Here, we consider
192
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193
the case where, in a certain zone, there are no measurement points. To do so we modified the case depicted in Fig. 11 by removing the measurement points from the second zone and adding them into zones 1 and 3. We realize a simulation without measurement errors. The algorithm performed eight iterations and terminated with the second stopping criterion (41) since the eighth and the sixth iterations are similar. The parameterization which corresponds to the smallest objective function is presented in Fig. 12. This parameterization is obtained after six iterations and the corresponding objective function value is about 200. The obtained parameterization is different from the exact one. Indeed, zone 1 is the same, however a part of zone 3 (the dashed boundary) is gathered to zone 2. The lack of data affects locally the parameterization and consequently the parameter values (see Table 5). Only the interface between zone 2 (which does not contain measurement points) and zone 3 has been changed. The parameter values in zones 1 and 3 are well estimated compared to the values in zone 2 (Table 5). Finally it is relevant to recall that the numerical model has been tested for several configurations of heterogeneity distribution. Good convergence has been obtained even for complicated cases. As an example, we present the required iterations for the case of central inclusion (Fig. 13).
1 2 J K f ¼ J K f ;1 ; . . . ; K f ;Nm and J 12 K f ¼ J K f ;ii2Z ; K f ;ii2Z 1
2
dK f ;i the perturbation of the parameter K f ;i is defined by:
( dK f ;i
¼
K 1 f ;i K f ;i
if i 2 Z 1
K 2 f ;i
if i 2 Z 2
K f ;i
ðA:5Þ
Then, the variation of the objective function corresponding to the perturbation dK f ;i (i.e., corresponding to the addition of the new degree of freedom) is: 2 12 dJ K f ¼ JðK 1 f ;ii2Z ; K f ;ii2Z Þ JðK f ;1 ; . . . ; K f ;Nm Þ ¼ J K f J K f 1
2
This paper describes a numerical model for efficiently identifying the hydrodynamic parameters in 2D dual-porosity media. An adaptive multiscale algorithm, based on the concept of global refinement indicators is proposed for the parameterization. The principle is as follows: We search the parameters vector as a piecewise constant and unknowns are both the parameters values and the geometry of zones. At each step, the algorithm ranks new possible parameterizations in the domain by using the global refinement indicators and selects from the remaining nominated parameterizations, the one which gives the best actual decrease in the objective function. The algorithm is tested on synthetic cases. In ideal cases (‘‘enough’’ measurements, spread over a wide range of values), the suggested algorithm allows to identify the proper parameterization and the hydraulic parameters for each zone. The measurement errors do not seriously affect the algorithm behavior. The algorithm provides a good parameterization even in case of lack of data. Appendix A. The concept of refinement indicator K f
Using a first order Taylor development, we obtain: Nm X @J dK f ;i @K f ;i i¼1 X @J X @J ¼ dK f ;i þ dK f ;i @K @K f ;i f ;i i2Z i2Z
J 12 K f J K f ¼ dJ K f
1
2
X @J X @J ðK 1 ðK 2 K f ;i Þ f ;i K f ;i Þ þ @K f ;i @K f ;i f ;i i2Z i2Z 1
Ne X @JðK f Þ i¼1
¼0
ðA:8Þ
Eq. (A.8) gives
X @JðK f Þ i2Z 1
@K f ;i
¼
X @JðK f Þ i2Z 2
ðA:9Þ
@K f ;i
Inserting (A.2) and (A.3) in (A.7) and using (A.9) leads to: J 12 K f J K f ¼ rK f
X @JðK f Þ i2Z 1
@K f ;i
¼ rK f
X @JðK f Þ i2Z 2
@K f ;i
X @J K X @JðK f Þ f ¼ ¼ i2Z1 @K f ;i i2Z2 @K f ;i
K f ;i ¼ K f
References
8i ¼ 1; . . . ; Nm
ðA:1Þ
where Nm is the number of elements in the computational mesh. To develop the calculation of the refinement indicator, let us consider the cutting configuration of Fig. 1b. And let us denote 2 12 by K 1 f , K f and J K f the optimal parameters and objective function corresponding to this tentative configuration. The value of K f in each mesh element is given by: 1 K 1 f ;i ¼ K f
8i 2 Z 1
K 2 f ;i
8i 2 Z 2
ðA:2Þ
2 rK f ¼ K 1 f K f , then we have
2 K 1 f ;i K f ;j ¼ rK f
8ði; jÞ 2 Z 1 Z 2 JK f
ðA:3Þ J 12 Kf
The optimal objective functions and can be written as functions of the parameter values in the mesh elements as follows:
ðA:10Þ
Eq. (A.10) shows clearly that, the largest decrease of the objective function (descending from J K f ) is measured by the absolute value P @JðK f Þ of the quantity i2Z1 @K . This absolute value is called the refinef ;i ment indicator associated with the splitting of the whole domain into two zones Z 1 and Z 2 , and it is defined by:
I‘K f
J K f
@K f ;i
Superscript ð‘Þ denotes the cutting configuration of Fig. 1b.
Let us denote
ðA:7Þ
2
Let us denote and the optimal parameter and objective function corresponding to the homogeneous case. Hence, we have the value of K f in each mesh element
¼
ðA:6Þ
Since the solution K f of the initial problem corresponds to a minimum, we have:
8. Conclusion
K 2 f
ðA:4Þ
ðA:11Þ
[1] Ackerer P, Younes A. Efficient approximations for the simulation of density driven flow in porous media. Adv Water Res 2006;31:15–27. http://dx.doi.org/ 10.1016/j.advwatres.2007.06.001. [2] Barenblatt GI, Zheltov YP, Kochina IN. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J Appl Math Mech 1960;24:1286–303. http://dx.doi.org/10.1016/0021-8928(60)90107-6. [3] Ben Ameur H, Chavent G, Jaffré J. Refinement and coarsening indicators for adaptive parameterization: application of the estimation of hydraulic transmissivities. Inverse Prob 2002;18(3):775–94. http://dx.doi.org/10.1088/ 0266-5611/18/3/317. [4] Boudreau BP. Diagenetic models and their implementation: modeling transport and reactions in aquatic sediments. Springer-Verlag; 1997. p. 414. [5] Brezzi F, Fortin M. Mixed and hybrid finite element methods. Berlin: Springer; 1991. [6] Burger M. A level set method for inverse problems. Inverse prob 2001;17:1327–56. http://dx.doi.org/10.1088/0266-5611/17/5/307. [7] Cacas MC, Ledoux E, de Marsily G, Tillie B, Barbreau A, Durand E, Feuga B, Peaudecerf P. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation 1. The flow model. Water Resour Res 1990;26(3):479–89. http://dx.doi.org/10.1029/WR026i003p00479.
H. Fahs et al. / Advances in Water Resources 63 (2014) 179–193 [8] Carrera J, Neuman SP. Estimation of aquifer parameters under transient and steady state conditions: 3. Application to synthetic and field data. Water Resour Res 1986;22(2):228–42. http://dx.doi.org/10.1029/WR022i002p00228. [9] Carrera J, Alcolea A, Medina A, Hidalgo J, Slooten J. Inverse problem in hydrogeology. Hydrogeol J 2005;13:206–22. http://dx.doi.org/10.1007/ s10040-004-0404-7. [10] Chavent G. Identification of functional parameter in partial differential equations. In: Goodson RE, Polis M, editors. Identification of parameters in distributed systems. New York: Am. Soc. Mech. Eng; 1974. p. 31–48. [11] Chavent G, Roberts JE. A unified physical presentation of mixed, mixed hybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems. Adv Water Res 1991;14:329–48. http://dx.doi.org/10.1016/0309-1708(91)90020-O. [12] Chavent G, Bissel R. Indicator for the refinement of parameterization. In: Tanaka M, Dulikravich GS, editors. Inverse problems in engineering mechanics. Elsevier; 1998. p. 309–14. [13] Davis TA. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Trans Math Softw 2004;30(2):165–95. http:// dx.doi.org/10.1145/992200.992205. [14] Davis TA. Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method. ACM Trans Math Softw 2004;30(2):196–9. http://dx.doi.org/10.1145/ 992200.992206. [15] Davis TA, Duff IS. A combined unifrontal/multifrontal method for unsymmetric sparse matrices. Technical Report, TR-97-016. Gainesville, FL: Computer and Information Science and Engineering Department, University of Florida; 1999. [16] Delay F, Kaczmaryk A, Ackerer P. Inversion of interference hydraulic pumping tests in both homogeneous and fractal dual media. Adv Water Res 2007;30:314–34. http://dx.doi.org/10.1016/j.advwatres.2006.06.008. [17] Dershowitz WS, Einstein HH. Characterizing rock joint geometry with joint system models. Rock Mech Rock Eng 1988;1(1):21–51. http://dx.doi.org/ 10.1007/BF01019674. [18] Durlofsky LJ. Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities. Water Resour Res 1994;30(4):965–73. http://dx.doi.org/10.1029/94WR00061. [19] Dverstorp B, Andersson J. Application of the discrete fracture network concept with field data: possibilities of model calibration and validation. Water Resour Res 1989;25(3):540–50. http://dx.doi.org/10.1029/WR025i003p00540. [20] Eppstein MJ, Dougherty DE. Simultaneous estimation of transmissivity values and zonation. Water Resour Res 1996;32(11):3321–36. http://dx.doi.org/ 10.1029/96WR02283. [21] Gerke HH, van Genuchten MT. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour Res 1993;29(2):305–19. http://dx.doi.org/10.1029/92WR02339. [22] Gerke HH, van Genuchten MT. Evaluation of a first order water transfer term for variably saturated dual-porosity flow models. Water Resour. Res. 1993;29(4):1225–38. http://dx.doi.org/10.1029/92WR02467. [23] Gerke HH, van Genuchten MT. Macroscopic representation of structural geometry of simulating water and solute movement in dual-porosity media. Adv Water Res 1996;19:343–57. http://dx.doi.org/10.1016/03091708(96)00012-7. [24] Hantush MM, Mariño MA, Islam MR. Models for leaching of pesticides in soils and groundwater. J Hydrol 2000;227:66–83. http://dx.doi.org/10.1016/S00221694(99)00166-3. [25] Hantush MM, Govindaraju RS, Mariño MA, Zhang Z. Screening mode for volatile pollutants in dual-porosity soils. J Hydrol 2002;260:58–74. http:// dx.doi.org/10.1016/S0022-1694(01)00597-2. [26] Hantush MM, Govindaraju RS. Theoretical development and analytical solutions for transport of volatile organic compounds in dual-porosity soils. J Hydrol 2003;279:18–42. http://dx.doi.org/10.1016/S0022-1694(03)00157-4. [27] Hayek M, Ackerer P. An adaptive subdivision algorithm for the identification of the diffusion coefficient in two-dimensional elliptic problems. J Math Model Algorithms 2007;6:529–45. http://dx.doi.org/10.1007/s10852-006-9046-1. [28] Hayek M, Lehmann F, Ackerer P. Adaptive multi-scale parameterization for one-dimensional flow in unsaturated porous media. Adv Water Res 2008;31:28–43. http://dx.doi.org/10.1016/j.advwatres.2007.06.009. [29] Herbert AW, Lanyon GW. Discrete fracture network modelling of flow and transport within a fracture zone at Stripa. In: Proceedings of the ISRM Regional Conference on Fractured and Jointed Rock Masses, Lake Tahoe, vol. 3. Berkeley, California: Lawrence Berkeley National Laboratory; 1992. p. 699–708. [30] Jarvis NJ, Jansson PE, Dik PE. Modelling water and solute transport in macroporous soil. 1. Model description and sensitivity analysis. J Soil Soc 1991;42:59–70. http://dx.doi.org/10.1111/j.1365-2389.1991.tb00091.x. [31] Jarvis NJ, Jansson PE, Dik PE. Modelling water and solute transport in macroporous soil. 2. Chloride breakthrough under non-steady flow. J Soil Soc 1991;42:71–81. http://dx.doi.org/10.1111/j.1365-2389.1991.tb00092.x.
193
[32] Kätterer TB, Schmied KC, Abbaspour KC, Schulin R. Single and dual-porosity modelling of multiple tracer transport through soil columns: effects of initial moisture and mode of application. Eur J Soil Soc 2001;52(1):25–36. http:// dx.doi.org/10.1046/j.1365-2389.2001.00355.x. [33] Khaleel R. Scale dependence of continuum models for fractured basalts. Water Resour Res 1989;25(8):1847–55. http://dx.doi.org/10.1029/WR025i008 p01847. [34] Larsbo M, Jarvis N. Simulating solute transport in a structured field soil: uncertainty in parameter identification and predictions. J Environ Qual 2005;34:621–34. http://dx.doi.org/10.2134/jeq2005.0621. [35] Larsson MH, Jarvis NJ. Evaluation of a dual-porosity model to predict field scale solute transport in a macroporous soil. J Hydrol 1999;215:153–71. http:// dx.doi.org/10.1016/S0022-1694(98)00267-4. [36] Levenberg K. A method for the solution of certain nonlinear problems in least squares. Quart Appl Math. 1944;2:164–8. MR0010666 (6,52a). [37] Lichtner PC. Critique of dual continuum formulation of multi-component reactive transport in fractured porous media. Geophys Monogr 2000;122:281–98. http://dx.doi.org/10.1029/GM122p0281. [38] Lien M, Berre I, Mannseth T. Combined adaptative and level-set parameter estimation. Multiscale Model Simul 2007;4(4):1349–72. http://dx.doi.org/ 10.1137/050623152. [39] Long JCS, Remer JS, Wilson CR, Witherspoon PA. Porous media equivalents for networks of discontinuous fractures. Water Resour Res 1982;18(3):645–58. http://dx.doi.org/10.1029/WR018i003p00645. [40] Lu Z, Robinson BA. Parameter identification using the level set method. Geophys Res Lett 2006;33:L06404. http://dx.doi.org/10.1029/2005GL025541. [41] Marquardt DW. An algorithm for least squares estimation of non-linear parameters. J Soc Ind Appl Math 1963;11(431):441. http://dx.doi.org/10.1137/ 0111030. [42] Mosé R, Siegel P, Ackerer P, Chavent G. Application of the mixed hybrid finite element approximation in a groundwater flow model: luxury or necessity? Water Resour Res 1994;30(11):3001–12. http://dx.doi.org/10.1029/ 94WR01786. [43] Oda M. Permeability tensor for discontinuous rock masses. Geotechnique 1985;35:483. http://dx.doi.org/10.1680/geot.1985.35.4.483. [44] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes: the art of scientific computing. 2nd ed. Cambridge: Cambridge University Press; 1992. [45] Ray C, Ellsworth TR, Valocchi AJ, Boast CW. An improved dual-porosity model for chemical transport in macroporous soils. J Hydrol 1997;193:270–92. http://dx.doi.org/10.1016/S0022-1694(96)03141-1. [46] Schwartz RC, Juo ASR, McInnes KJ. Estimating parameters for a dual-porosity model to describe non-equilibrium, reactive transport in an fine-textured soil. J Hydrol 2000;229:149–67. http://dx.doi.org/10.1016/S0022-1694(00)00164-5. [47] Šimunek J, Jarvis NJ, van Genuchten MT, Gårdenäs A. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J Hydrol 2003;272:14–35. http://dx.doi.org/10.1016/ S0022-1694(02)00252-4. [48] Sun NZ. Inverse problems in groundwater modelling. Dordrecht, The Netherlands: Kluwer Academic; 1994. p. 337. [49] Sun NZ, Yeh WG. Identification of parameter structure in groundwater inverse problem. Water Resour Res 1985;21(6):869–83. http://dx.doi.org/10.1029/ WR021i006p00869. [50] Svensson U. A continuum representation of fracture networks. Part I: method and basic test cases. J Hydrol 2001;250:170–86. http://dx.doi.org/10.1016/ S0022-1694(01)00435-8. [51] Tsai FTC, Sun NZ, Yeh WG. Global–local optimization for parameter structure identification in three-dimensional groundwater modeling. Water Resour Res 2003;39(2):1043. http://dx.doi.org/10.1029/2001WR001135. [52] Wilson CR, Witherspoon PA, Long JCS, Galbraith RM, Dubois AO, Macpherson MJ. Large-scale hydraulic conductivity measurements in fractured granite. Int J Rock Mech Min Sci Geomech Abstr 1983;20(6):71–83. http://dx.doi.org/ 10.1016/0148-9062(83)90596-X. [53] Younes A, Ackerer P, Lehmann F. A new mass lumping scheme for the mixed hybrid finite element method. Int J Numer Methods Eng 2006;67:89–107. http://dx.doi.org/10.1002/nme.1628. [54] Younes A, Ackerer P, Delay F. Mixed finite elements for solving 2-D diffusiontype equations. Rev Geophys 2010;48:RG1004. http://dx.doi.org/10.1029/ 2008RG000277. [55] Zhang X, Sandersson DJ, Harkness RM, Last NC. Evaluation of the 2-D permeability tensor for fractured rock mass. Int J Rock Mech Min Sci Geomech Abstr 1996;33(1):17–37. http://dx.doi.org/10.1016/01489062(95)00042-9.
