375
Views
1
CrossRef citations to date
0
Altmetric
Articles

First-order indicators for the estimation of discrete fractures in porous media

, , , , &
Pages 1-32 | Received 23 Feb 2016, Accepted 29 Jan 2017, Published online: 06 Mar 2017

Abstract

Faults and geological barriers can drastically affect the flow patterns in porous media. Such fractures can be modelled as interfaces that interact with the surrounding matrix. We propose a new technique for the estimation of the location and hydrogeological properties of a small number of large fractures in a porous medium from given distributed pressure or flow data. At each iteration, the algorithm builds a short list of candidates by comparing fracture indicators. These indicators quantify at the first order the decrease of a data misfit function; they are cheap to compute. Then, the best candidate is picked up by minimization of the objective function for each candidate. Optimally driven by the fit to the data, the approach has the great advantage of not requiring remeshing, nor shape derivation. The stability of the algorithm is shown on a series of numerical examples representative of typical situations.

AMS Subject Classifications:

1. Introduction

The accurate simulation of flow in porous media is important for many applications such as petroleum reservoir management, monitoring and clean-up of underground pollutants, and planning for underground nuclear-waste disposal. In porous media, large fractures are often present [Citation1,Citation2], and they can modify drastically the flow pattern: the fracture, sometimes called a fault in this case, can be very permeable and thus channel rapidly the fluid from and to the surrounding rock. Another typical regime arises when the fracture is almost impervious to the flow: it acts as a geological barrier. Accurately simulating the interaction of flow between porous media and fractures is therefore important. For this reason, it is desirable to determine as precisely as possible the location of any important fracture and its hydrogeological signature, i.e. its impact on the flow.

The goal of this article is to present a new methodology for estimating the location and the hydrogeological properties of a small number of large fractures in a porous medium. Geological experiments [Citation3] indeed show that the fluid tends to choose its pathway along only some few of the many existing fractures. So we do not address the case of a large number of small fractures that could be treated using a double continuum model [Citation4,Citation5], but search only for a limited number of fractures, the ones having the major impact on the flow.

Imaging techniques for assessing underground media such as seismic imaging [Citation6,Citation7] can give important information about the location of fractures, but they are not able to tell to what extent if any of the fracture influences the flow. These techniques should be seen as complementary to the present work: they could for instance be used as a priori information for the inversion algorithm designed here.

The numerical flow model for the fractured medium must take into account the strong heterogeneity and the very different spatial scales involved. Reduced models or co-dimension 1 models treat the fractures as interfaces in the porous medium, and allow flow in the fracture as well as fluid exchange between the fracture and the rock matrix. Such interface models have been studied extensively in the engineering literature, e.g. see [Citation8Citation13], as well as in the mathematical literature, e.g. see [Citation14Citation22], to name just a few. We use here the model developed in [Citation14], where flow in the fracture as well as in the rock is governed by Darcy’s law. The model was first presented only for a fault regime [Citation23], but was extended to treat both fault and barrier regimes in [Citation14,Citation24,Citation25]. In this model, fractures lie along the edges of the mesh and can be easily opened or closed by adjusting the fracture parameters on edges, which makes it convenient for the purpose of fracture determination. This model has been extended so that one may use non-matching grids [Citation18,Citation26], or disconnect the fracture mesh from that of the domain [Citation17,Citation27]. It has also been extended to treat Forchheimer flow in the fracture [Citation20,Citation28], and multiphase flow [Citation11], but these extensions are not considered in this paper.

We suppose that the permeability is known outside the fracture, and concentrate on the determination of the position and intensities of a few fractures that explain pressure and/or flow data. We have considered both distributed pressure and distributed velocity measurements throughout the domain. These data are used to build a least-squares objective function which measures the misfit between data and the corresponding computed values. Depending on the available data, the determination of the position and intensity of fractures which minimize this objective function can be poorly conditioned or under determined. We count on the parsimonious introduction of fractures, which is inherent to our inversion algorithm, to have a beneficial regularizing effect, see [Citation29,Citation30] for instance.

The key ingredient in our approach is the notion of indicator [Citation30, section 3.7]. This idea was originally developed in [Citation31], for the estimation of hydraulic transmissivities, and was then extended for the estimation of vector valued distributed parameters in [Citation32,Citation33]. In this paper, we adapt this notion to construct fracture indicators, which are used to locate fractures. Given the current model of the fractured porous medium one associates any additional candidate fracture, i.e. any set of new contiguous edges not belonging to any existing fracture, with a fracture indicator. This indicator measures, up to the first order, the rate of decrease of the objective function which would be achieved by adding the candidate fracture to the current model. But, the important feature of the method is that, once both direct and adjoint systems have been solved for the current model, the fracture indicator associated with any candidate fracture can be computed at a very small cost.

Ideally, one should compute fracture indicators for the extremely long list of all possible candidate fractures, and retain only a short list of those with the largest indicators. However, even with this cheap indicator, such an exhaustive search is way out of reach. Hence, one has to resort to heuristic exploration strategies to reduce the length of this long seminal list. One possible strategy [Citation31Citation33] is to select candidate parameterizations in a long list made of predefined families. Here, this would have required to consider predefined families of candidate fractures. Due to the lack of a priori information for the choice of such families, we propose instead a new constructive strategy based on elementary candidate fractures. In this approach, the fracture indicators are first computed for all items of the long list of elementary candidate fractures. Then, proceeding by aggregation and extension, a short list of candidate fractures with large indicators is produced.

The algorithm is iterative. Each iteration provides an additional fracture that produces a better fit to the data. The sketch of the algorithm goes as follows. For a given set of data measures, an initial model with no fracture is chosen. Then, each iteration is made of two steps: an indicator step, which provides a short list of candidate fractures as explained above, and an optimization step, which selects the winning fracture in the short list of candidate fractures and determines the associated intensities. In this optimization step, pre-optimal fracture intensities are computed for all models obtained by adding in turn each candidate fracture of the short list to the current model. Note that only a limited number of such minimizations are necessary, as only a small number of candidate fractures have been retained in the short list during the indicator step. Note also that these optimizations are performed on a space of small dimension, that is equal to the number of fractures (say less than 10) times the number of physical parameters (1 or 2). At the end of the iteration, the winning candidate fracture is added to the collection of already retained fractures. This new fracture may be either an extension of a previously found fracture or an entirely new fracture. Then, the iterative process may resume with the new estimated model.

With this approach, the underlying mesh is fixed during the procedure, and the “opening” of a fracture requires only changing its intensity parameter to a nonzero value. Hence, one avoids remeshing, most of the assembly process, and the computation of shape derivatives and topological gradients, or similar techniques to “open the fracture” (in other words, to compute a gradient with respect to a fracture that does not exist yet). The one drawback of the method is that the search is limited to fractures located on the edges of a fixed mesh. The fracture indicator algorithm is presented and tested numerically for a two-dimensional model; however, it remains valid for a three-dimensional model.

The paper is organized as follows: the direct model for the numerical simulation of interactions between a porous medium and a discrete fracture is given in detail in Section 2. The inverse problem for the estimation of fracture parameters is presented in Section 3. The fracture indicators are introduced in Section 4, and the fracture indicator algorithm is detailed in Section 5. Finally, Section 6 is devoted to numerical results.

2. The direct model

We consider the simulation of single-phase flow in a porous medium Ω containing a fracture. For simplicity of exposition we describe both the continuous and discrete flow models in a two-dimensional setting, though we stress that the methodology is equally valid in a three-dimensional setting. Also for simplicity, in the description of the models only, we suppose that Ω contains one fracture γ.

Figure 1. Notations associated with a fracture γ in the domain Ω.

Figure 1. Notations associated with a fracture γ in the domain Ω.

2.1. The model problem

Let the porous medium Ω be a bounded domain in R2. We suppose that the fracture γΩ is a regular, non-self-intersecting, curve of finite length and of bounded curvature. In this two-dimensional setting γ, the boundary of γ, consists of the endpoints of γ which we assume are distinct and do not lie on the boundary of Ω. Let n+ be a unit normal vector field on γ, and let n-=-n+. We distinguish the two sides of the fracture, γ+ and γ-, such that n+ (resp. n-) can be seen as the exterior normal to Ω\γ on γ+ (resp. γ-), see Figure . For a sufficiently regular function φ in Ω\γ, we consider two traces φ+ and φ- on γ, one for each side; see [Citation15] for a more rigorous definition. We also introduce the tangential gradient τ and the tangential divergence divτ operators along the fracture.

The fracture interface model [Citation14] may be written:(1) u=-KpinΩ\γ,divu=finΩ\γ,uγ=-ατpγinγ,divτuγ=fγ+[[u·n]]inγ,-βu+·n++p+=pγonγ+,-βu-·n-+p-=pγonγ-.(1) (With respect to the model (4.1) in [Citation14], the numerical parameter ξ is taken equal to 1, since it has very little influence on the flow in most cases, see [Citation14].)

The first two equations in (Equation1) are Darcy’s law and the law of mass conservation with a source term f in the domain Ω\γ. They relate the Darcy velocity u with the pressure p through the permeability tensor K. The next equation is Darcy’s law, relating the tangential Darcy velocity uγ to the fluid pressure in the fracture pγ. The fourth equation is the law of mass conservation along the fracture γ, with a possible fracture source term fγ, supplemented by an exchange term [[u·n]]=u+·n++u-·n- accounting for fluid exchange between γ and Ω\γ. The last two equations that close the system represent an averaged Darcy law across normal cross sections of the physical fracture represented by γ.

The model (Equation1) is supplemented with the boundary conditions(2) p=p¯onΩ,uγ·nγ=0onγ,(2)

where p¯ is a given function, and nγ is the unit normal vector field on γ pointing outwards from γ (and tangent to γ at its endpoints), see Figure . The condition uγ·nγ=0 is coherent with the hypothesis that the endpoints γ are inside the porous medium.

The coefficients α and β are hydrogeological parameters characterizing the fracture: α represents an effective tangential permeability. In the case of (conducting) geological faults, α is positive and the Darcy velocity is generally discontinuous across the fracture ([[u·n]]0). In a two-dimension setting, α is the product of a thickness and a permeability divided by a length. The coefficient β represents the inverse of an effective permeability normal to the fracture. So β models the resistivity of the fracture to flow across it. In the case of geological barriers, β is positive and the pressure is generally discontinuous across the fracture (p+p-). In two dimensions, β is the product of a thickness, a permeability and a length.

Our objective will be to determine both the location of the fracture γ and the physical properties α and β of the fracture.

We assume that the permeability tensor K is symmetric and positive definite, and bounded away from zero, and that the parameters α and β satisfy the hypotheses(3) ααmin>0andβ0.(3)

Under these hypotheses, the model (Equation1) and (Equation2) is well-posed, see [Citation14].

2.2. Discretization of the model problem

Let Th be a finite element mesh of Ω, conforming with the fracture γ in the sense that γ is a subset of the union of the closures of the edges of the elements of Th. We denote by TEN, respectively, the elements (e.g. triangles), edges and nodes of Th, by ET the boundary edges of an element T, by TE the elements having E as an edge, by NE the endpoints of an edge E, and by EN the edges having N as an endpoint. Let Ehi be the set of internal edges of Th and Ehe the set of boundary edges of Th, so that Ω=EEheE¯ and EheEhi=. We denote by Eh=defEhiEhe the set of all edges of Th.

We say that the fracture γ is supported by the internal mesh edges to mean that γEEhiE¯, and we assume that the parameters α and β are discretized by piecewise constant functions on these fracture edges. To simplify the search for independent faults and barriers, we discretize the fracture γ as(4) Fh=def(Fα,Fβ),(4)

where we have separated fault edges in Fα from barrier edges in Fβ,(5) Fα=def{EEhi|αE>0}andFβ=def{EEhi|βE>0}.(5)

We denote the collections of fault (or tangential) and barrier (or normal) parameters by(6) α=def(αE)EFαandβ=def(βE)EFβ.(6)

The sets of fault and barrier edges need not be disjoint. Let Nα be the set of all nodes of Fα.

In order to simplify the presentation, we assume that the mesh is made up of square elements of size h, that the permeability tensor in Ω is constant and scalar K=KI (with K>0), and that the model (Equation1) and (Equation2) is discretized using a cell-centered finite volume method, see [Citation25]. The discretization in the more general case of lowest order Raviart–Thomas mixed hybrid finite elements with general meshes is presented in [Citation34]. These two discretizations in rectangular meshes are equivalent up to an approximate quadrature formula, see [Citation35].

Let PT, PE and PN be the pressure unknowns on an element TTh, an edge EEh and a node NNα. Let UT,E be the flow rate of fluid leaving an element TTh through an edge EET, let UE,N be the flow rate of fluid leaving an edge EFα through a node NNE. Using the notations κ=defK-1h, fT=defTf, fE=defEfγ, and P¯E=def1|E|Ep¯, the finite volume formulation of (Equation1) and (Equation2) reads(7) EETUT,E=fT,for allTTh,κ2UT,E=PT-PE,for allEFβ,TTE,κ2+βEUT,E=PT-PE,for allEFβ,TTE,PE=P¯E,for allEEhe,-TTEUT,E=0,for allEFα,EEhi,(7)

and(8) -TTEUT,E+NNEUE,N=fE,for allEFα,h2αEUE,N=PE-PN,for allEFα,NNE,EFαEENUE,N=0,for allNNα,(8)

where the last equation already contains, for endpoints of the fault Fα, the no-flow condition uγ·nγ=0 of (Equation2).

In Section 4, the construction of fracture indicators will require to pass to the limit in the above equations when αE0 on some edges of Fα. But a quick look at the second equation in (Equation8) let us foresee some problem in doing so. So we replace on all fault edges the flow rate unknown UE,N by a new unknown VE,N defined by (remember that αE>0 in this section)(9) αEVE,N=UE,N,for allEFα,NNE.(9)

With these unknowns, (Equation8) is replaced by(10) -TTEUT,E+αENNEVE,N=fE,for allEFα,h2VE,N=PE-PN,for allEFα,NNE,EFαEENαEVE,N=0,for allNNα.(10)

The second equation in (Equation10) shows that the new unknown VE,N is the (discrete) pressure gradient between the middle of edge E and its end node N.

For the sake of simplicity, we suppose, in the sequel, the fracture source term to be zero (fE=0 in (Equation10)). To sum up, the (discrete) direct model for the simulation of flow in a porous domain Ω containing a discrete fracture (Fα,Fβ) with fault and barrier parameters (α,β) is made up of equations (Equation7) and (Equation10). Thus the direct problem associates to the fracture, i.e. the pair Fh=(Fα,Fβ) together with the fracture parameters, the solution of the direct model.

3. The inverse problem

Focusing now on the inverse problem, we assume that pressure and/or flow measurements are available at some points of the domain Ω, and consider the problem of estimating the location and parameters of the corresponding fracture. For the sake of simplicity, we limit ourselves in this section to the simple case where a pressure measurement P¯T is available for each element T in Th. More realistic measurements will be considered for numerical experiments in Section 6.

3.1. The least-squares formulation

In the framework of an iterative process for solving the inverse problem, let the current fracture be described by its edges Fh=(Fα,Fβ) and parameters (α,β). We associate the current fracture with the least-square data misfit(11) JFh(α,β)=def12TTh(PT-P¯T)2,(11)

where PT is given by the current model (Equation7) and (Equation10) set on the current fracture Fh.

If this misfit is above the uncertainty level of the data measurements, one is led to consider adding new candidate edges ΔFh=(ΔFα,ΔFβ) to the current fracture Fh. We suppose that current and candidate fractures have no common edge,(12) FαΔFα=andFβΔFβ=.(12)

Let ΔNα be the set of new nodes of the candidate fault ΔFα, so that(13) NαΔNα=.(13)

We use nominal parameters α¯ and β¯ to describe typical parameter profiles expected along the new candidate fault and barrier edges. In the absence of such information, nominal parameters are chosen to be constant. We model the introduction of the candidate fracture by new intensity parameters εα and εβ, which let the fracture parameters α and β grow from zero to their nominal value and beyond,(14) αE=εαα¯E,for allEΔFα,εα0,(14) (15) βE=εββ¯E,for allEΔFβ,εβ0.(15)

Because of (Equation14), the last equation of (Equation10) can be divided by εα on ΔNα, thus the model (Equation7) and (Equation10) set on both current and candidate fractures Fh and ΔFh simplifies to the tentative model (16) EETUT,E=fT,for allTTh,κ2UT,E=PT-PE,for allEFβΔFβ,TTE,κ2+βEUT,E=PT-PE,for allEFβΔFβ,TTE,PE=P¯E,for allEEhe,-TTEUT,E=0,for allEFαΔFα,EEhi,-TTEUT,E+αENNEVE,N=0,for allEFαΔFα,h2VE,N=PE-PN,for allEFαΔFα,NNE,EFαΔFαEENαEVE,N=0,for allNNα,EΔFαEENα¯EVE,N=0,for allNΔNα.(16)

When εα=εβ=0, one easily checks that the tentative model (Equation16) reduces to the current model (Equation7) and (Equation10) complemented with the limit fault model set on the candidate fault,(17) h2VE,N=PE-PNfor allEΔFα,NNE,EΔFαEENα¯EVE,N=0for allNΔNα.(17)

Hence, we see that, for any εα,εβ0, the tentative model (Equation16) determines uniquely PT, PE, and UT,E on all elements and edges, and PN and VE,N on all nodes and edges of the current and candidate faults. So, we can extend the definition (Equation11) of the least-squares data misfit to a function of the intensity parameters εα,εβ0 on the candidate fracture by(18) JΔFh(εα,εβ)=def12TTh(PT-P¯T)2,(18)

where PT is now given by the tentative model (Equation14)–(Equation16). The purpose of the indicator algorithm developed in the next sections will be to select, in a computationally efficient way, the candidate fractures ΔFh for which there exists εα,εβ>0 such that(19) JΔFh(εα,εβ)<JΔFh(0,0)=JFh(α,β),(19)

with a preferably large decrease of the data misfit.

Remark 3.1:

Replacing the hypothesis εα=εβ=0 in (Equation14) and (Equation15) by the weaker assumption that α,β0 for all candidate edges of ΔFh, makes the last equation in (Equation10) disappear on nodes of ΔNα: unknowns (PN)NΔNα and (VE,N)EΔFα,NNE become undetermined for the limit equations, but all other unknowns remain uniquely determined.

Actually, in Section 4, the key ingredient for the definition of fracture indicators is the fading of fracture parameters on the candidate fracture. And, choosing a specific direction in the space of fracture edge parameters through (Equation14) and (Equation15), makes the solution of the tentative model (Equation16) behave continuously for vanishing intensity parameters.

3.2. Gradient with respect to the intensity parameters

We compute now the gradient of the objective function with respect to εα and εβ by the adjoint state method, for instance see [Citation30, Section 2.4]. For the sake of simplicity, the Dirichlet boundary condition (fourth equation in (Equation16)) is not written explicitly, it is understood that whenever E is a boundary edge, PE means P¯E (and hence its differential vanishes). With this convention, the vector of edge pressure unknowns becomes (PE)EEhi.

The Lagrangian associated with the minimization of the objective function (Equation18) under the constraint of the state equations (Equation16) for the current and candidate fractures FhΔFh is then(20) L(X,Λ;εα,εβ)=def12TTh(PT-P¯T)2+TTh-EETUT,E+fTλT+EFβΔFβTTEκ2UT,E-PT+PEμT,E+EFβΔFβTTEκ2+βEUT,E-PT+PEμT,E+EFαΔFαEEhiTTEUT,EλE+EFαΔFαTTEUT,E-αENNEVE,NλE+EFαΔFαNNEh2VE,N-PE+PNμE,N+NNαEFαΔFαEENαEVE,NλN+NΔNαEΔFαEENα¯EVE,Nλ¯N,(20)

where εα and εβ are the intensity parameters, X=def(Pt,Pe,Pn,Ute,Ven) is the state vector, and Λ=def(λt,λe,λn,λ¯n,μte,μen) is the adjoint state vector, with the notations(21) Pt=def(PT)TTh,Ute=def(UT,E)TTh,EET,Pe=def(PE)EEhi,Ven=def(VE,N)EFαΔFα,NNE,Pn=def(PN)NNαΔNα,λt=def(λT)TTh,μte=def(μT,E)TTh,EET,λe=def(λE)EEhi,μen=def(μE,N)EFαΔFα,NNE,λn=def(λN)NNα,λ¯n=def(λ¯N)NΔNα.(21)

Let X(εα,εβ) denote the state vector corresponding to the solution of the tentative model (Equation16) with the fracture parameters α and β defined by (Equation14) and (Equation15) on the candidate fracture ΔFα. Let Λ be a given adjoint state vector. Then, for all εα,εβ0, we have the identity(22) JΔFh(εα,εβ)=L(X(εα,εβ),Λ;εα,εβ).(22)

Differentiation with respect to parameters εα and εβ gives (remember that Λ is given)(23) δJΔFh(εα,εβ)=δL(X(εα,εβ),Λ;εα,εβ),(23)

hence, for all εα,εβ0,(24) δJΔFh(εα,εβ)=LX(X(εα,εβ),Λ;εα,εβ)δX+L(εα,εβ)(X(εα,εβ),Λ;εα,εβ)δ(εα,εβ).(24)

This identity holds for any adjoint state vector Λ. Thus, nothing can prevent us from choosing for Λ precisely the unique solution Λ(εα,εβ) of the adjoint equations such that(25) LX(X(εα,εβ),Λ;εα,εβ)δX=0,δX.(25)

Then, (Equation24) simplifies to the desired gradient equations, for all εα,εβ0,(26) δJΔFh(εα,εβ)=L(εα,εβ)(X(εα,εβ),Λ(εα,εβ);εα,εβ)δ(εα,εβ),(26)

and the gradient of JΔFh is simply obtained by differentiation of the Lagrangian (Equation20) with respect to the intensity parameters εα and εβ. This gives, for all εα,εβ0,(27) JΔFhεα(εα,εβ)=-EΔFαα¯ENNEVE,NλE+EΔFαα¯ENNENNαVE,NλN,JΔFhεβ(εα,εβ)=EΔFββ¯ETTEUT,EμT,E,(27)

where VE,N and UT,E are solutions of the state equations (Equation16), and λE, λN, and μT,E are solution of the adjoint equations (Equation25), which we explicit now.

Differentiation of the Lagrangian (Equation20) with respect to the state vector X gives(28) LX(X,Λ;εα,εβ)δX=TTh(PT-P¯T)δPT-TThEETδUT,EλT+EFβΔFβTTEκ2δUT,E-δPT+δPEμT,E+EFβΔFβTTEκ2+βEδUT,E-δPT+δPEμT,E+EFαΔFαEEhiTTEδUT,EλE+EFαΔFαTTEδUT,E-αENNEδVE,NλE+EFαΔFαNNEh2δVE,N-δPE+δPNμE,N+NNαEFαΔFαEENαEδVE,NλN+NΔNαEΔFαEENα¯EδVE,Nλ¯N,(28)

where, according to our convention, δPE=0 on boundary edges. The adjoint state equation (Equation25) is satisfied if and only if all coefficients of δPT, δPE, δPN, δUT,E and δVE,N are equal to zero. This gives the following tentative adjoint model for the determination of the adjoint variables λT, λE, λN, λ¯N, μT,E and μE,N,(29) EETμT,E=PT-P¯T,for allTTh,κ2μT,E=λT-λE,for allEFβΔFβ,TTE,κ2+βEμT,E=λT-λE,for allEFβΔFβ,TTE,λE=0,for allEEhe,-TTEμT,E=0,for allEFαΔFα,EEhi,(29) (30) -TTEμT,E+NNEμE,N=0,for allEFαΔFα,h2μE,N=αE(λE-λN),for allEFαΔFα,NNαNE,h2μE,N=αEλE-α¯Eλ¯N,for allEΔFα,NΔNαNE,EFαΔFαEENμE,N=0,for allNNαΔNα.(30)

We rewrite now these adjoint equations in a form similar to that of the tentative model (Equation16). First, define λN for NΔNα and νE,N for all current and candidate fault edges by (compare with the seventh equation in (Equation16))(31) λN=EΔFαEENα¯EλEEΔFαEENα¯E,for allNΔNα,h2νE,N=λE-λN,for allEFαΔFα,NNE.(31)

Then, from the last two equations in (Equation30), one checks that the adjoint variables λ¯N and μE,N satisfy (compare with (Equation9))(32) λ¯N=εαλNfor allNΔNα,μE,N=αEνE,Nfor allEFαΔFα,NNE.(32)

With the definition (Equation31), (Equation30) becomes (compare with the last four equations in (Equation16))(33) -TTEμT,E+αENNEνE,N=0,for allEFαΔFα,h2νE,N=λE-λN,for allEFαΔFα,NNE,EFαΔFαEENαEνE,N=0,for allNNα,EΔFαEENα¯EνE,N=0,for allNΔNα.(33)

Remark 3.2:

Unsurprisingly, the tentative adjoint model (Equation29) and (Equation30), or (Equation29) and (Equation33), makes sense for εα,εβ0, as the tentative direct model (Equation16) does. When εα=εβ=0, the adjoint model reduces to the current adjoint model, that will be defined in (Equation37), complemented with a limit fault adjoint model similar to (Equation17) and whose solution is given by (Equation31) on the candidate fault.

The introduction of λN and νE,N gives a new formula for JΔFh/εα. Indeed, the multiplication of the last equation in (Equation16) by λN and the summation over all nodes of ΔNα give(34) NΔNαEΔFαEENα¯EVE,NλN=0.(34)

Then, by adding this zero term to the right-hand side of the first gradient equation in (Equation27), one obtains a new formula for the derivative of JΔFh with respect to εα, for all εα,εβ0,(35) JΔFhεα(εα,εβ)=-h2EΔFαα¯ENNEVE,NνE,N.(35)

This formula is more elegant than the one in (Equation27). It requires to solve the adjoint limit fault model (Equation31) for λN and νE,N, but only on the candidate fault.

Remark 3.3:

When we do not have a pressure measurement for all elements of the mesh, the right-hand side of the first equation in (Equation29) vanishes for elements with no measurement. In the same way, when a Darcy velocity measurement U¯T,E is provided for element TTh and for edge EET, the additional source term UT,E-U¯T,E appears on the right-hand side in either the second or the third equation in (Equation29), following that E is a barrier edge or not.

To conclude this section, we recall that the tentative direct and adjoint systems (Equation16), (Equation29), and (Equation33) define uniquely the direct and adjoint variables X and Λ for all εα,εβ0, including PN, λN, VE,N and νE,N on the candidate fault. Hence, the gradient (Equation27) or (Equation35) of the objective function JΔFh(εα,εβ) is also defined when the intensity parameters εα and εβ vanish on the candidate fracture.

4. Fracture indicators

We address in this section the definition of fault and barrier indicators for the selection of additional candidate faults and barriers which are likely to produce a better fit to the data.

Let the current fracture, described by its locations Fh=(Fα,Fβ) and parameters α and β, be given by a priori knowledge or from previous computations. Fh may be void, for instance at the beginning of the fracture determination procedure when no a priori information is available. Let X=(Pt,Pe,Pn,Ute,Ven) be the solution of the current model (Equation7) and (Equation10), which reads(36) EETUT,E=fT,for allTTh,κ2UT,E=PT-PE,for allEFβ,TTE,κ2+βEUT,E=PT-PE,for allEFβ,TTE,PE=P¯E,for allEEhe,-TTEUT,E=0,for allEFα,EEhi,-TTEUT,E+αENNEVE,N=0,for allEFα,h2VE,N=PE-PN,for allEFα,NNE,EFαEENαEVE,N=0,for allNNα,(36)

and let Λ=def(λt,λe,λn,μte,νen) be the solution of the current adjoint model obtained from (Equation29) and (Equation33) for a vanishing candidate fracture with zero intensity parameters,(37) EETμT,E=PT-P¯T,for allTTh,κ2μT,E=λT-λE,for allEFβ,TTE,κ2+βEμT,E=λT-λE,for allEFβ,TTE,λE=0,for allEEhe,-TTEμT,E=0,for allEFα,EEhi,-TTEμT,E+αENNEνE,N=0,for allEFα,h2νE,N=λE-λN,for allEFα,EEN,EFαEENαEνE,N=0,for allNNα.(37)

If the current fracture Fh does not give a small enough value to the objective function JFh(α,β) defined in (Equation11), one considers testing the effect of adding a set ΔFαEhi of new fault edges and/or a set ΔFβEhi of new barrier edges. As in Section 3, we suppose that these candidate fractures have no common edge with the current fracture, i.e. satisfy (Equation12) and (Equation13).

In order to select candidate fractures ΔFh=(ΔFα,ΔFβ) which are likely to decrease the objective function in the sense of (Equation19), we make a first-order development of JΔFh defined in (Equation18),(38) JΔFh(εα,εβ)-JΔFh(0,0)JFh(α,β)JΔFhεα(0,0)εα+JΔFhεβ(0,0)εβ.(38)

Then, following [Citation30, Section 3.7], we define the fault indicator Iα and the barrier indicator Iβ associated, respectively, with the candidate fault ΔFα and the candidate barrier ΔFβ by(39) Iα(ΔFα)=defJΔFhεα(0,0)andIβ(ΔFβ)=defJΔFhεβ(0,0).(39)

When these indicators are negative and of large absolute value, ΔFα and/or ΔFβ will be good candidate fault or barrier, as they should provide a strong decrease of the cost function (at least at first order) when the intensity parameters εα or εβ increase. Of course, these fracture indicators only give first-order information, but they are inexpensive to compute, as we see now.

Equation (Equation35) gives for the fault indicator the formula(40) Iα(ΔFα)=-h2EΔFαα¯ENNEVE,NνE,N,(40)

where (VE,N,νE,N)EΔFα,NNE are given by the direct limit fault model (Equation17) and its adjoint counterpart set on the candidate fault edges of ΔFα only, whose solutions are given by(41) PN=EΔFαEENα¯EPEEΔFαEENα¯E,λN=EΔFαEENα¯EλEEΔFαEENα¯E,for allNΔNα,h2VE,N=PE-PN,h2νE,N=λE-λN,for allEΔFα,NNE.(41)

Here, (PE,λE)EΔFα and (PN)NNαNE,EΔFα are simply extracted from the solutions X and Λ of the current direct and adjoint models (Equation36) and (Equation37).

Similarly, the second equation in (Equation27) gives for the barrier indicator the formula(42) Iβ(ΔFβ)=EΔFββ¯ETTEUT,EμT,E,(42)

where (UT,E,μT,E)EΔFβ,TTE are simply extracted from the solutions X and Λ of the current direct and adjoint models (Equation36) and (Equation37).

Remark 4.1:

The fault indicator formula (Equation40) depends on the direct and adjoint tangential Darcy velocities along the candidate fault. From Section 2.1, this is physically sound, since the fault parameter α is linked to the flow along the candidate fault. In the same way, the barrier indicator formula (Equation42) depends on the direct and adjoint Darcy velocities normal to the candidate barrier, which is natural, since β is linked to the flow across to the candidate barrier.

Remark 4.2:

If ΔFα is made of a single edge E with no endpoint on the current fault Fα, one checks easily from the last equation of (Equation16) that JΔFh(εα,0)=JFh(α,β) for all εα0, i.e. opening a fault on a single edge does not change the pressure and flow patterns. Notice also that in this case, (Equation41) implies immediately that VE,N=0 at the endpoints of E, and the corresponding fault indicator is zero (which is coherent). Thus, it will be necessary to use only candidate faults constituted of at least two connected edges, or having a common node with the current fault.

To conclude this section, we remark that once the current direct and adjoint models (Equation36) and (Equation37) have been solved, the computation of the indicators Iα(ΔFα) or Iβ(ΔFβ) only requires the summation of known quantities along the candidate fractures ΔFα or ΔFβ (plus the resolution of (Equation41) in the case of faults). Hence, the computation of the indicators for many candidate faults or barriers is inexpensive since it requires the resolution of only one set of direct and adjoint model with the current fracture Fh. This is the basis of the indicator algorithm developed in Section 5.

5. Algorithm for the estimation of fractures

We have now all the ingredients for the presentation of an iterative algorithm, based on fracture indicators, designed for the estimation of discrete fractures in a porous medium from given pressure and/or flow measurements.

We suppose in this section that pressure data are available in all elements of the computational mesh, but the case of coarser measurements, or of flow data can be treated as well by modifying the objective function (Equation18) and the adjoint state model (Equation37) according to Remark 3.3.

The simultaneous determination of faults and barriers is a delicate matter: the number of parameters is multiplied by two, and this can lead to underdetermination unless a large number of measurements are available. Hence, we limit ourselves to an algorithm for estimating either faults or barriers, but not both at the same time. We present the algorithm under a generic form, hence throughout this section, we drop the subscripts α and β: F, ΔF, and ε will refer to either Fα, ΔFα and εα, or Fβ, ΔFβ and εβ, depending on the nature of the sought fracture. When necessary, the fracture parameter of interest α or β will be denoted ξ.

Fractures have been assumed so far to be any set of internal edges. Again, this can lead to underdetermination, which is usually associated with a poor conditioning of the inverse problem and with the presence of many meaningless local minima for the objective function. Moreover, we are only interested in locating fractures that have a macroscopic impact on the flow. Therefore, we make the assumption that fractures are defined on a coarser mesh than the computational mesh. Nevertheless, the amount of possible fractures for such a coarse representation can still be quite large, so in order to avoid overparameterization, we use the fracture indicators defined in Section 4 to introduce fractures one at a time, and we use (Equation14) and (Equation15) to parameterize each fracture with only one parameter ε.

5.1. The fracture mesh

To exclude very small fractures, typically made of only a few edges, we decouple the computational mesh Th from the fracture research mesh Thc. The computational mesh should be fine enough to perform accurate calculations, and the fracture mesh should correspond to the expected size of the fracture. To avoid interpolation issues, we simply choose the fine mesh Th to be a refinement of the coarse mesh Thc. Hence, coarse edges are made up of several computational edges.

From now on, we suppose that both current and candidate fractures are made of internal coarse edges, and we denote by Fc (resp. ΔFc) the set of coarse edges belonging to the current fractures (resp. a candidate fracture). In terms of computational edges, these fractures are given by(43) F={EEhi|EcFc,EEc},ΔF={EEhi|EcΔFc,EEc}.(43)

It is understood that in the sequel all current and candidate fractures F and ΔF are of the form (Equation43) for some Fc and ΔFc, i.e. actually made of coarse edges.

5.2. A collection of fractures

At each iteration 1, a new fracture ΔF is added to the current collection of fractures F-1. Hence, after k-1 iterations, the current collection of fractures is of the form(44) Fk-1=F0ΔF1ΔFk-1,(44)

where F0 is the (usually void) set of given a priori fractures.

The setting of Sections 3 and 4 is replicated for each fracture. Let εk-1 be the intensity parameter associated with the fracture ΔF. Equations (Equation14) and (Equation15) become (different fractures have no common edges)(45) ξE=εk-1ξ¯E,EΔF,=1,,k-1.(45)

When a priori information is available, values ξ0 of the fracture parameter on the initial set of fractures F0 are also provided; they are kept fixed. Then, the data misfit of equation (Equation18) becomes a function of the vector of intensity parameters,(46) JFk-1(ε1k-1,,εk-1k-1)=def12TTh(PT-P¯T)2,(46)

where PT is given by the current model (Equation36) with the current collection of fractures Fk-1. The optimal intensities εk-1 are the minimizer of the objective function satisfying(47) εk-1=(ε1k-1,,εk-1k-1)=defargminεk-1>0JFk-1(ε1k-1,,εk-1k-1).(47)

As the current collection of fractures is known by its locations Fk-1 and intensities εk-1, one can solve the associated current direct and adjoint models (Equation36) and (Equation37), complemented with (Equation45) for the optimal intensities εk-1. This determines the current direct and adjoint variables Xk-1=(Pt,Pe,Pn,Ute,Ven) and Λk-1=(λt,λe,λn,μte,νen).

To sum up, after k-1 iterations, the following quantities are available: the collection of k-1 estimated fracture locations ΔF1, ..., ΔFk-1, the associated optimal intensities εk-1=(ε1k-1,,εk-1k-1), and the corresponding vectors of current direct and adjoint variables Xk-1 and Λk-1.

The case k=1 corresponds to the initialization of the algorithm. It is specific: there is no estimated fracture (i.e. no intensity), X0 and Λ0 are obtained as solution of (Equation36) and (Equation37) with the set of given a priori fractures F0.

5.3. The next iteration

When the fit to the data obtained in (Equation47) after k-1 iterations is not satisfactory, one needs to perform an additional iteration. This means selecting a new fracture ΔFk enabling the data misfit to decrease, and determining the resulting optimal fracture intensities εk=(ε1k,,εk-1k,εkk) for the new collection of fractures Fk=Fk-1ΔFk. This is achieved in three steps.

Indicator step. A long list of candidate fractures is chosen. As in Sections 3 and 4, the candidate fractures must have no common edge with the current collection of fractures Fk-1. The fracture indicators I(ΔF) for all candidate fractures ΔF of the long list are computed by (Equation40) and (Equation41) for faults, or (Equation42) for barriers. These computations are very fast as all terms are already available from the current direct and adjoint variables Xk-1 and Λk-1 provided by the previous iteration.

Then, a short list of candidate fractures is built. The candidate fractures of the short list are associated with (negative) fracture indicators of large magnitude. Thus, according to the first-order information carried by the indicators, they are the most likely to produce a significant enhancement of the fit to the data.

Strategies to choose the long list and to build the short list are discussed in Section 5.6.

Optimization step. For each candidate fracture ΔF of the short list, solve the minimization problem(48) ε~k(ΔF)=(ε~1k,,ε~kk)=argminεk>0JFk-1ΔF(ε1k,,εkk),(48)

where the objective function JFk-1ΔF is similar to (Equation46), here PT is still given by the current model (Equation36), but set on Fk-1ΔF with the k fracture intensities (ε1k,,εkk). Thanks to the use of a nominal parameter ξ¯, the natural initial guess for the minimization is (ε1k-1,,εk-1k-1,1).

The computation of the gradient of JFk-1ΔF is required at each iteration of the minimization algorithm. Hence, the need for the resolution of both direct and adjoint models (Equation36) and (Equation37) set on Fk-1ΔF with the current value of (ε1k,,εkk) at each minimization iteration. Thus, each minimization is much more computationally intensive than the calculation of the fracture indicators by (Equation40) and (Equation41), or (Equation42). This is why a short list of a small number of candidate fractures is built. The optimization step is usually more expensive than the indicator step, by several orders of magnitude.

Update step. The winner ΔFk is the candidate fracture that gives the smallest minimum value to the objective function. The new current collection of fractures is Fk=Fk-1ΔFk, and the new optimal fracture intensity vector is εk=ε~k(ΔFk).

Finally, the new vectors of direct and adjoint variables Xk and Λk are determined by solving direct and adjoint models (Equation36) and (Equation37) set on Fk with intensities εk.

5.4. Stopping the algorithm

The algorithm stops either when the maximal number of fractures is reached, when adding a new fracture does not significantly improve the objective function, or when the data misfit is of the same magnitude as the uncertainty on the given data measures. The constants driving the stopping criteria are: the absolute uncertainty on the data measures Δmeas0, the relative tolerance for convergence ηconv>0, the relative tolerance for stationary sequence ηstat>0, and the maximum number of admissible fractures kmax.

When the algorithm stops, the result is given in terms of the fracture parameter: ξk is restored from the a priori given values ξ0 on the initial set of fractures F0, and from the last estimated intensities εk through (Equation45) (expressed for k instead of k-1).

5.5. Algorithm

Here comes the precise description of the algorithm. Remember that the fault and barrier indices α and β are omitted: F, ΔF and ε stand either for Fα, ΔFα and εα, or for Fβ, ΔFβ, and εβ. Moreover, ξ stands either for α or β.

A priori information, or maybe a previous computation, provides the initial guess for fractures (location F0 and fracture parameter ξ0, which may be void), and a value ξ¯E of the nominal parameter in each edge E of the mesh.

Initialization.

(1)

Compute the initial solutions X0 and Λ0 of the direct and adjoint models (Equation36) and (Equation37) set on the initial fractures F0 with fracture parameter ξ0.

(2)

Compute the initial objective function J0=JF0(ξ0) with (Equation11).

Iterations. For k1, do:
(1)

Indicator step. Build the short list of candidate fractures according to the chosen strategy (see Section 5.6). This uses fracture indicators (Equation40) and (Equation41), or (Equation42), that depends on Xk-1 and Λk-1.

(2)

If the short list is empty, then the algorithm stops, the result is Fk-1 (with k-1 estimated fractures), and ξk-1 restored from ξ0 and intensities εk-1 through (Equation45).

(3)

Optimization step. For each candidate fracture ΔF in the short list, solve (Equation48):

  • minimize the objective function JFk-1ΔF with respect to the k fracture intensity parameters εk=(ε1k,,εkk);

  • call ε~k(ΔF) the minimizer.

(4)

Update step. Retain in the short list the candidate fracture ΔFk that gives the smallest value to JFk-1ΔF(ε~k(ΔF)), and set Fk=Fk-1ΔFk,εk=ε~k(ΔFk),Jk=JFk(εk). Then, compute the solutions Xk and Λk of the current direct and adjoint models (Equation36) and (Equation37) set on the current collection of fractures Fk with intensities εk.

(5)

If Jk-1-JkηstatJ0, then the algorithm is stationary, the result is Fk-1 (with k-1 estimated fractures), and ξk-1 restored from ξ0 and intensities εk-1.

(6)

If Jk-ΔmeasηconvJ0, then the algorithm has converged, the result is Fk (with k estimated fractures), and ξk restored from ξ0 and intensities εk.

(7)

If kkmax, then the algorithm has reached the maximum number of expected fractures, the result is Fk (with k estimated fractures), and ξk restored from ξ0 and intensities εk.

 

5.6. Strategy to build a short list of candidate fractures

The general idea is to choose a long list of candidate fractures, trying to be as exhaustive as possible, or to incorporate a priori knowledge, and then to use the inexpensive first-order indicators to select a short list of candidate fractures that are likely to provide a large decrease of the objective function.

For the design of refinement indicators in [Citation31Citation33], the unknowns were of a completely different nature: the parameters were distributed all over the domain, with a value in each element of the mesh. The algorithm sought for the location of parameter discontinuities represented by cutting curves. The long list of candidate cuttings was chosen from a small number of predefined families of curves, such as parallel lines, or circles. Here, some fractures have to be created at specific locations on the edges of the mesh. It would have been possible to use the same kind of choices for the long list of candidate fractures. Instead, we propose a constructive strategy.

The key idea is to choose the initial long list as a collection of elementary candidate fractures that allow for a cartography of the values of the fracture indicators all over the domain. Elementary fractures are made of a few edges to allow for exhaustiveness, and they are meant to grow into larger fractures having a macroscopic impact on the flow and associated with (negative) fracture indicators of larger magnitude.

The strategy depends on the choice of a coarse fracture research mesh Thc of which the computational mesh Th is a refinement. It is also parameterized by constants driving the selection criteria: two ratios θelem and θext (between 0 and 1) for indicator-based selection of candidate fractures, and the maximum admissible number of candidate fractures nmax.

The three stages of the strategy go as follows. First, according to Remark 4.2, the fracture indicators are computed for a long list of elementary fractures made of all pairs of contiguous internal coarse edges with at most one node in common with the current collection of fractures, complemented with all single internal coarse edges with one node on the current collection of fractures. Let I,elem be the best (i.e. minimum) indicator for all elementary fractures (remember that good indicators are very negative). If this minimum is positive, then the short list of candidate fractures is empty (and the fracture indicator algorithm stops). Otherwise, we select all elementary fractures whose indicator is lower or equal to θelemI,elem.

Then, we aggregate together all selected elementary fractures that share a common coarse edge to form the longest possible candidate fractures. The minimum indicator for all aggregates of elementary fractures I,aggr is usually lower than I,elem, and there is no need for a selection here since the number of aggregates is smaller than the number of selected pairs.

Finally, we allow for an extension of all aggregates by one coarse edge at one or two of its endpoints. Note that there could be three endpoints, or more. Note also that it can increase considerably the number of candidate fractures. Let I,ext be the best (i.e. minimum) indicator for all extended aggregates. Again, this minimum is lower or equal to I,aggr. We select all extended aggregates whose indicator is lower or equal to θextI,ext. At the end, in order to control the cost of the optimization step of the algorithm, we truncate the short list of candidate fractures to the nmax selected extended aggregates associated with the lowest indicators.

6. Numerical results

We present now numerical experiments conducted with a Matlab® implementation of the fracture indicator-based algorithm described in Section 5 for the estimation of discrete fractures in a porous medium from given pressure measurements.

First of all, as it must always be the case, the implementation of the closed-form formulas (Equation27), (Equation35), (Equation40) and (Equation42) for the gradient of the objective function with respect to the intensity parameters, and for the fracture indicators, have been successfully checked by comparison with a finite-difference approximation based on the sole computation of the objective function.

The test-cases are simple, but represent typical situations of interest: faults parallel to the main direction of the flow, and barriers normal to the flow. The synthetic pressure and flow data are obtained by using the same simulation programme: the tests range from the so-called “inverse crime” for which the very same model is used for the inversion and simulation of data, to situations where random noise of increasing level is added, and also to the much more difficult but interesting situation where the target fractures are not carried by the fracture research mesh.

We do not use any a priori information on the target fractures: the initial estimation is no fracture (F0=), and nominal parameters are taken constant α¯=β¯=1; thus εα=α and εβ=β. The constants for tuning the behaviour of the algorithm are chosen as follows: Δmeas is proportional to the noise level, ηconv=ηstat=10-2, and kmax=8 for the stopping criteria, and θelem=0.8, θext=0.9, and nmax=10 for the building of the short list of candidate fractures.

Three levels of mesh. To allow for the decimation of the measurements, in addition to the (fine) computational mesh Th and to the (coarse) fracture research mesh Thc, we also use a specific mesh for the synthetic measurements, denoted by Thm (the sum in (Equation46) is now meant for TThm).

Let Ω be a rectangular domain, and QNx,Ny denote the regular rectangular grid of Nx×Ny cells for positive integers Nx and Ny. In the sequel, the domain Ω is the unit square, the computational mesh is Th=Q72,72, the measurement mesh is Thm=QNm,Nm with Nm ranging from 72 down to 8, and the fracture mesh defined in Section 5.1 is Thc=QNf,Nf with Nf=12 or 9. Note that Th is a refinement of Thc since Nf is a divisor of 72.

Test-cases. The porous medium is homogeneous with a permeability K=1. Impervious Neumann conditions are imposed on top and bottom, and Dirichlet conditions impose a pressure drop from right (p¯=1) to left (p¯=0). Without any fracture, the solution of the direct model is a linear pressure (see Figure (b)) and a uniform Darcy velocity.

A fault normal to the flow, or a barrier tangential to the flow, has a null hydrogeological signature: the solution of the direct model is still the same uniform flow from right to left. Therefore, we consider either tangential faults, or normal barriers, for a maximum effect on the solution. In the following examples, the fractures have length 0.5 (half the size of the domain), and they are located either in the middle of the domain (single fracture), or centred at one quarter and at three quarters of the domain (two fractures). The fault parameter α (conductivity along the fracture) ranges from 0.2 up to 200. The barrier parameter β (resistivity across to the fracture) ranges from 0.02 up to 200.

Figure 2. Case of a single fault tangential to the flow. (a) Data pressure, solution of the direct model (Equation36) with a fault at y=0.5 (α=2). (b) Initial pressure, solution of the direct model with no fracture. (c) Distribution of indicators for the long list of elementary candidate faults: the lowest indicator for all six pairs of edges centred in each interior coarse node is represented. (d) In dark blue, the target fault; above in yellow, the selected pairs of edges (three of them are superimposed above the target fault). (e) In light blue, the candidate faults of the short list after aggregation and extension; above in green, the best candidate fault. (f) Best result after minimization for all candidate faults of the short list (the lowest permeabilities are in light blue, and the highest in pink).

Figure 2. Case of a single fault tangential to the flow. (a) Data pressure, solution of the direct model (Equation36(36) ∑E∈ETUT,E=fT,for allT∈Th,κ2UT,E=PT-PE,for allE∉Fβ,T∈TE,κ2+βEUT,E=PT-PE,for allE∈Fβ,T∈TE,PE=P¯E,for allE∈Ehe,-∑T∈TEUT,E=0,for allE∉Fα,E∈Ehi,-∑T∈TEUT,E+αE∑N∈NEVE,N=0,for allE∈Fα,h2VE,N=PE-PN,for allE∈Fα,N∈NE,∑E∈FαE∈ENαEVE,N=0,for allN∈Nα,(36) ) with a fault at y=0.5 (α=2). (b) Initial pressure, solution of the direct model with no fracture. (c) Distribution of indicators for the long list of elementary candidate faults: the lowest indicator for all six pairs of edges centred in each interior coarse node is represented. (d) In dark blue, the target fault; above in yellow, the selected pairs of edges (three of them are superimposed above the target fault). (e) In light blue, the candidate faults of the short list after aggregation and extension; above in green, the best candidate fault. (f) Best result after minimization for all candidate faults of the short list (the lowest permeabilities are in light blue, and the highest in pink).

6.1. Illustration of the algorithm for the estimation of faults

Let us first consider the simple case of a fault tangential to the flow located in the middle of the domain with α=2. The synthetic data are represented in Figure (a). The pressure is measured in all elements of the computational mesh (i.e. Nm=72). The fracture research mesh corresponds to Nf=12, i.e. each coarse edge is made of six fine edges.

The algorithm is initialized with no fracture (F0=). The initial direct and adjoint variables X0 and Λ0 are computed as solution to (Equation36) and (Equation37) with no fracture. The initial pressure is linear (Figure (b)).

The first iteration starts with the computation of the fault indicators through (Equation41) and (Equation40) from X0 and Λ0 for all elementary candidate faults of the long list made of the 726 pairs of contiguous interior coarse edges of Thc=Q12,12. The cartography of these indicators (Figure (c)) shows a nice blue spot in the middle of the domain, right where the target fault is located (remember that good indicators are large negative numbers). Only five pairs of coarse edges have an indicator lower than θelemIα,elem; they are parallel and close to the target fault (Figure (d)).

The aggregation stage collapses the three pairs of edges superimposed with the target fault, reducing the number of candidate faults down to three (Figure (d)). There is a significant improvement of the indicator as Iα,aggr1.8Iα,elem.

The extension stage adds one coarse edge at zero, one or two of the endpoints of the aggregated pairs of edges. This multiplies the number of candidate faults by 16 (6 ways to add one edge, 9 ways to add two edges, plus the initial one). The best indicator increases again, by 14%; it corresponds to the target fault (in green in Figure (e)). The short list of candidate faults is made of the 7 extended aggregated pairs having an indicator lower than θextIα,ext (in light blue in Figure (e)).

After the minimization of the objective function with respect to the fault parameter (remember that εα=α) for all 7 candidate faults of the short list, the best data misfit is divided by a huge factor of 1016 from its initial value. This corresponds to the candidate fault with the best indicator, which happens to be the target fault. And at the minimum, the target value α=2 is recovered (Figure (f)). The algorithm stops on the convergence criterion.

Figure 3. Case of two faults tangential to the flow. (a) Data pressure, solution of the direct model (Equation36) with two faults at y=0.25 (α=2) and y=0.75 (α=20). (b) Distribution of indicators for the long list of elementary candidate faults at the first iteration (see Figure (c)). (c) Best result after minimization at the first iteration for all candidate faults of the short list. (d) Pressure, solution of the direct model with the estimated fault of (c). (e) Distribution of indicators for the long list of elementary candidate faults at the second iteration. (f) Best result after minimization at the second iteration for all candidate faults of the short list.

Figure 3. Case of two faults tangential to the flow. (a) Data pressure, solution of the direct model (Equation36(36) ∑E∈ETUT,E=fT,for allT∈Th,κ2UT,E=PT-PE,for allE∉Fβ,T∈TE,κ2+βEUT,E=PT-PE,for allE∈Fβ,T∈TE,PE=P¯E,for allE∈Ehe,-∑T∈TEUT,E=0,for allE∉Fα,E∈Ehi,-∑T∈TEUT,E+αE∑N∈NEVE,N=0,for allE∈Fα,h2VE,N=PE-PN,for allE∈Fα,N∈NE,∑E∈FαE∈ENαEVE,N=0,for allN∈Nα,(36) ) with two faults at y=0.25 (α=2) and y=0.75 (α=20). (b) Distribution of indicators for the long list of elementary candidate faults at the first iteration (see Figure 2(c)). (c) Best result after minimization at the first iteration for all candidate faults of the short list. (d) Pressure, solution of the direct model with the estimated fault of (c). (e) Distribution of indicators for the long list of elementary candidate faults at the second iteration. (f) Best result after minimization at the second iteration for all candidate faults of the short list.

The picture is even more interesting with two faults. Let us now consider the case of two faults tangential to the flow centred at one quarter and at three quarters of the domain with α=2 (bottom) and 20 (top). The synthetic data are represented in Figure (a). The measurement and fracture meshes remain unchanged (Nm=72 and Nf=12).

At the first iteration, the cartography of indicators for all elementary candidate faults of the long list shows in Figure (b) two blue spots of highly negative values centred at the correct location of the target faults. Moreover, the best indicators are located around the upper fault having the highest permeability (in dark blue). As in the case of a single fault, the aggregation and extension stages multiply the best indicator by a factor of 2, and the best indicator corresponds to the most permeable target fault (the upper one). Again, the best minimizer among the 9 candidate faults of the short list corresponds to the best indicator (the upper target fault). The optimization step decreases the objective function by a factor of 3 only, and at the minimum, a highly overestimated value of about 3300 is recovered for α (Figure (c)). Yet, the pressure of Figure (d) recovers the main perturbation due to the upper target fault (compare with Figure (a)). As the stopping criteria are not satisfied, another iteration is performed.

At the second iteration, the cartography of indicators presents now in Figure (e) a strong dark blue spot located at the centre of the least permeable target fault. Clearly, there is nothing left to do around the upper target fault (the area is mostly dark red). The indicators are smaller in magnitude than in the first iteration by a factor of about one third. The optimization step exhibits the best candidate fault of the short list (among 10 of them) as the best minimizer. Again, this best minimizer is exactly the missing lower target fault, the data misfit is divided by a huge factor of 21013 from its initial value, and the target values α=20 and α=2 are perfectly recovered (Figure (f)). The algorithm stops again on the convergence criterion.

6.2. Beyond the inverse crime

We investigate now less favourable situations with less measurements, and where the synthetic data are no longer in the range of the model used for inversion. The comprehensive campaign of numerical tests is detailed in [Citation34].

Figure 4. Estimation of faults with less measurement data. The inverted faults are represented with colours ranging from light blue (low permeabilities) to pink (high permeabilities). The target fault is represented in dark blue. (a) Single target fault, no noise, Nm=6. (b) Two target faults, no noise, Nm=8. (c) Single target fault, 4% noise, Nm=8. (d) Two target faults, 4% noise, Nm=12. (e) Single target fault, 6% noise, Nm=8. (f) Two target faults, 2% noise, Nm=8.

Figure 4. Estimation of faults with less measurement data. The inverted faults are represented with colours ranging from light blue (low permeabilities) to pink (high permeabilities). The target fault is represented in dark blue. (a) Single target fault, no noise, Nm=6. (b) Two target faults, no noise, Nm=8. (c) Single target fault, 4% noise, Nm=8. (d) Two target faults, 4% noise, Nm=12. (e) Single target fault, 6% noise, Nm=8. (f) Two target faults, 2% noise, Nm=8.

Influence of decimation of measurements. When the measurement mesh Thm=QNm,Nm is reduced from Nm=72 down to Nm=8, numerical tests shown in [Citation34] indicate that the algorithm still provides the target fault as the best candidate fault, and the minimization step still yields the target value of α. Moreover, this remains true for target values of α ranging from 0.2 up to 200. For values above 200, the direct model outputs almost the same pressure field: the cost function becomes insensitive to the fault parameter α. It is the same for the case of two target faults: down to Nm=12, the algorithm still provides the most permeable target fault as the best candidate fault at the first iteration, the least permeable target fault at the second iteration, and then, the minimization step yields the target values of α for both faults.

Problems appear when the distribution of measurement points becomes too loose. At best, the target faults are still located quite accurately, but the minimization step is unable to recover the target values of α, and the results may be misevaluated by several orders of magnitude as in the case of two target faults with Nm=8 (Figure (b)). Or the situation may degenerate progressively: the recovered locations become quite diffuse, but the effective permeabilities can still be quite accurate, as in the case of a single target fault with Nm=6 (Figure (a)).

Influence of random noise. When adding white Gaussian noise up to a relative noise level of 6% (the perturbation due to the presence of the fault becomes hardly visible by the eye), the algorithm still exhibits the target fault as the best candidate fault, and the minimization step recovers a value of α=2.01, almost identical to the target value α=2, see [Citation34]. When also decimating the measurements down to Nm=8, the algorithm still recovers the target fault as the best candidate fault up to a relative noise level of 4%, and the minimization step yields promising values of 1.5 (with 2% noise) and 1.2 (with 4% noise) for the fault parameter α (Figure (c)). The algorithm is not as successful with 6% noise and Nm=8: a reasonable estimated fault in the middle (in medium blue) is polluted by a very strong fault in the lower right corner of the domain (Figure (e)).

When considering the case of two faults, the conditions are still favourable up to a noise level of 4%: the location of the most permeable fault is recovered at the first iteration and that of the other at the second iteration. Moreover, the minimization step recovers the value α=17 for the upper fault and α=1.99 for the lower fault, which is still close to the target values 20 and 2. When also decimating the measurements down to Nm=12, the algorithm still recovers correctly the least permeable (lower) fault with an excellent inverted value of α=1.98, but the most permeable (upper) fault is highly overestimated (values over half a million, Figure (d)) and with 4% noise, the location of this fault exhibits artifacts associated with negligible permeabilities under 0.1. For Nm=8, only the strongest (upper) fault is estimated with 2% noise, it is located one cell under the target location and the permeability is highly overvalued to about 8700 (Figure (f)).

Figure 5. Estimation of faults with coherent noise: the target faults (dark blue) are not carried by the fracture research mesh (Nf=9). The inverted faults are represented with colours ranging from light blue (low permeabilities) to pink (high permeabilities). (a) Single target fault, Nm=72. (b) Single target fault, Nm=8 (c) Two target faults, Nm=72. (d) Two target faults, Nm=8.

Figure 5. Estimation of faults with coherent noise: the target faults (dark blue) are not carried by the fracture research mesh (Nf=9). The inverted faults are represented with colours ranging from light blue (low permeabilities) to pink (high permeabilities). (a) Single target fault, Nm=72. (b) Single target fault, Nm=8 (c) Two target faults, Nm=72. (d) Two target faults, Nm=8.

Influence of coherent noise. Finally, we investigate the much more difficult situation where the synthetic pressure data do not belong to the range of the direct model used for inversion, i.e. the global minimum of the objective function is not 0. Of course, we can no longer expect to perfectly retrieve the location and permeability of the target faults, but at least we can hope to recover groups of faults with similar hydrogeological signatures; in fact, values of α are adding up for nearby parallel faults tangential to the flow.

We consider Nf=9, hence the target faults located in the middle, or at one quarter and three quarters of the domain, are no longer carried by the fracture research mesh anymore. We also decimate the pressure data from Nm=72 down to Nm=8. The estimated locations of the faults contain most of the coarse edges of the fracture research mesh that are closest to the target faults. Moreover, with the single target fault, the recovered effective permeabilities range from 1.8 to 3.4, for a target value of 2 (Figure (a) and (b)). And with the two target faults, the recovered effective permeabilities are higher for the upper fault than for the lower, but they are underestimated up to an order of magnitude (Figures (c) and (d)): from 1.6 to 2.5 for the most permeable fault (the target value is 20), and from 1.0 to 1.6 for the least one (the target value is 2).

Figure 6. Cases of a single barrier (a)–(c), and of two barriers (d)–(g), normal to the flow. (a) Data pressure, solution of the direct model (Equation36) with a barrier at x=0.5 (β=2). (b) Distribution of indicators for the long list of elementary candidate barriers (see Figure (c)). (c) Best result after minimization for all candidate barriers of the short list. (d) Data pressure, solution of the direct model (Equation36) with two barriers at x=0.25 (β=2) and x=0.75 (β=20). (e) Distribution of indicators at the first iteration. (f) Distribution of indicators at the second iteration. (g) Best result after minimization for all candidate barriers of the short list.

Figure 6. Cases of a single barrier (a)–(c), and of two barriers (d)–(g), normal to the flow. (a) Data pressure, solution of the direct model (Equation36(36) ∑E∈ETUT,E=fT,for allT∈Th,κ2UT,E=PT-PE,for allE∉Fβ,T∈TE,κ2+βEUT,E=PT-PE,for allE∈Fβ,T∈TE,PE=P¯E,for allE∈Ehe,-∑T∈TEUT,E=0,for allE∉Fα,E∈Ehi,-∑T∈TEUT,E+αE∑N∈NEVE,N=0,for allE∈Fα,h2VE,N=PE-PN,for allE∈Fα,N∈NE,∑E∈FαE∈ENαEVE,N=0,for allN∈Nα,(36) ) with a barrier at x=0.5 (β=2). (b) Distribution of indicators for the long list of elementary candidate barriers (see Figure 2(c)). (c) Best result after minimization for all candidate barriers of the short list. (d) Data pressure, solution of the direct model (Equation36(36) ∑E∈ETUT,E=fT,for allT∈Th,κ2UT,E=PT-PE,for allE∉Fβ,T∈TE,κ2+βEUT,E=PT-PE,for allE∈Fβ,T∈TE,PE=P¯E,for allE∈Ehe,-∑T∈TEUT,E=0,for allE∉Fα,E∈Ehi,-∑T∈TEUT,E+αE∑N∈NEVE,N=0,for allE∈Fα,h2VE,N=PE-PN,for allE∈Fα,N∈NE,∑E∈FαE∈ENαEVE,N=0,for allN∈Nα,(36) ) with two barriers at x=0.25 (β=2) and x=0.75 (β=20). (e) Distribution of indicators at the first iteration. (f) Distribution of indicators at the second iteration. (g) Best result after minimization for all candidate barriers of the short list.

6.3. Estimation of barriers

The case of barriers is quite similar. We present in Figure some tests for one or two barriers normal to the flow. The barrier indicators for the long list of elementary candidate barriers are even more precise: the target barriers are almost already drawn at the beginning of the first iteration of the algorithm in the cartographies of indicators in Figure (b) and (e). Indeed, the location of the most resistive target barrier is associated with the strongest indicator, and the optimization step picks it up for it produces a perfect fit to the data (Figure (c) and (g)). Nevertheless, in the case of two barriers, during the extension stage of the second iteration, the location of the weakest target barrier is not associated with the strongest indicator, and it is not even picked up by the optimization step. Consequently, although the best location for the second barrier contains all edges of the least resistive target barrier, the optimization results are not perfect. The recovered values for the resistivity are β=82 and β=0.41, instead of the target values 20 and 2; they are quite far, but the hierarchy is correct.

Figure 7. Estimation of barriers with coherent noise: the target barriers (dark blue) are not carried by the fracture research mesh (Nf=9). The inverted barriers are represented with colours ranging from light blue (low resistivities) to pink (high resistivities). (a) Single target barrier, Nm=72. (b) Single target barrier, Nm=8. (c) Two target barriers, Nm=72. (d) Two target barriers, Nm=8.

Figure 7. Estimation of barriers with coherent noise: the target barriers (dark blue) are not carried by the fracture research mesh (Nf=9). The inverted barriers are represented with colours ranging from light blue (low resistivities) to pink (high resistivities). (a) Single target barrier, Nm=72. (b) Single target barrier, Nm=8. (c) Two target barriers, Nm=72. (d) Two target barriers, Nm=8.

When using a fracture research mesh that no longer carries the target barriers (Nf=9), the inversion results are still quite reasonable, see Figure . We also decimate the pressure data measures from Nm=72 down to Nm=8. The estimated locations of the barriers contain most of the coarse edges of the fracture research mesh that are closest to the target barriers. After minimization, the recovered effective resistivities (values of β are also adding up for nearby parallel barriers normal to the flow) are underestimated by one or two orders of magnitude, but the hierarchy is still preserved.

6.4. Discussion

The first lesson learned from the numerical experiments is that, as hoped, the inversion algorithm presented in Section 5 performs accurately in favourable situations. Indeed, exact fit to synthetic data that are in the range of the inversion model can be reached in a wide variety of situations: for one or two faults tangential to the flow, for a barrier normal to the flow, for a wide range of target parameter values (up to four orders of magnitude), and when loosening the number of distributed measurements (up to a factor of about 80).

Much more interesting is the stability property shown by the algorithm. The recovered fractures remain perfectly located, or at least superimposed with the target location, when adding random noise up to a level of 6% for faults, and of 8% for barriers, see [Citation34]. And this remains true when loosening the number of measurements by a factor of about 80 for up to 4% noise (for faults) and up to 5% noise for barriers. Moreover, when the situation deteriorates a bit more, for instance when the measurement mesh becomes just too loose, when the noise level becomes just too high, or when the target fractures no longer belong to the fracture research mesh, then the location of recovered fractures remains within a distance of no more than a coarse cell from the target location. Furthermore, even when fractures are estimated by several smaller fractures of close location, the value of the recovered effective parameter is of the same order than the target value for single fractures, and the hierarchy in values is kept in the case of two target fractures.

However, the heuristic strategy detailed in Section 5.6 was designed for the case of fault detection, and it may be suboptimal for the estimation of barriers. Indeed, even in the most advantageous case of the inverse crime with full measurements, the indicator step of the algorithm is unable to nominate the target location of the least resistive barrier in the case of two target barriers normal to the flow.

Finally, the amount of measurements necessary to accurately recover the location of fractures may seem prohibitive (64 in our tests). This indicates that one needs a fine enough distribution of measurement to hope solving this difficult inverse problem. Some further experiments could deal with a study of the influence of the location of measurements: a priori information may help driving the measurements in the vicinity of the target fracture in order to reduce their number.

7. Conclusion

The estimation of fractures in porous media from pressure and/or flow data is known to be an ill-posed inverse problem. We propose a new approach for the recovering of both location and hydrogeological properties of a small number of large fractures. The method does not need any remeshing, nor shape derivation.

We have chosen a numerical flow model that treats fractures as interfaces interacting with the surrounding porous matrix, and we have derived a forward model for the simulation of both faults and barriers with vanishing intensities, provided they are located on edges of the simulation mesh. Then, the approach is based on the minimization of an error function that measures the distance between data and simulated measurements.

We have defined specific fracture indicators, whose sign and modulus give the effect on the data misfit of opening any candidate fault or barrier. These indicators give only first-order information, but are computationally inexpensive.

Overparameterization is avoided by restricting the search for fractures to edges of a coarser mesh than the computational grid, and by introducing fractures one at a time. When the data misfit is too high, a short list of new candidate fracture locations is built from a large number of elementary fracture locations that span the entire domain of interest, with the assistance of fracture indicators. The actual enhancement brought by each candidate is computed by optimization, and the best performing fracture is retained.

Numerical tests have been performed on synthetic data corresponding to simple typical situations. They demonstrate the ability of the proposed algorithm to automatically retrieve one or two faults parallel to the flow, or one or two barriers perpendicular to the flow. The algorithm is shown to be fairly stable when noise is added to the data, and when the fracture to be detected is not located on the fracture search mesh.

Notes

No potential conflict of interest was reported by the authors.

References

  • Dietrich P, Helmig R, Sauter M, et al. Flow and transport in fractured porous media. Berlin: Springer; 2005.
  • Adler PM, Thovert JF, Mourzenko VV. Fractured porous media. Oxford: Oxford University Press; 2013.
  • Le Goc R, de Dreuzy JR, Davy P. Statistical characteristics of flow as indicators of channeling in heterogeneous porous and fractured media. Adv Water Res. 2010;33(3):257–269. Available from: http://www.sciencedirect.com/science/article/pii/S0309170809001936.
  • Barenblatt GI, Aheltov IP, Kochina IN. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J Appl Math Mech. 1960;24:1286–1303.
  • Warren JE, Root PJ. The behavior of naturally fractured reservoirs. Soc Pet Eng J. 1963;245–255.
  • Biondi BL. 3D seismic imaging, Investigations in geophysics. Vol. 14. Tulsa (OK): Society of exploration geophysicists; 2006.
  • Bimpas M, Anditis A, Uzunoglu N. Integrated high resolution imaging radar and decision support system for the rehabilitation of water pipelines. London: IWA publishing; 2010.
  • Baca RG, Arnett RC, Langford DW. Modelling fluid flow in fracture-porous rock masses by finite-element techniques. Int J Num Methods Fluids. 1984;4:337–348.
  • Karimi-Fard M, Firoozabadi A. Numerical simulation of water injection in fractured media using the discrete-fracture model and the galerkin method. SPE REE. 2003;117–126.
  • Karimi-Fard M, Durlofsky LJ, Aziz K. An efficient discrete fracture model applicable for general purpose reservoir simulators. SPE J. 2004;227–236.
  • Reichenberger V, Jakobs H, Bastian P, Helmig R. A mixed-dimensional finite volume method for multiphase flow in fractured porous media. Adv Water Res. 2006;29(7):1020–1036.
  • Kim Deo. Finite element discrete-fracture model for multiphase flow in porous media. AIChE J. 2006;463(6):1120–1130.
  • Hoteit H, Firoozabadi A. An efficient numerical model for incompressible twophase flow in fractured media. Adv Water Res. 2008;31:891–905.
  • Martin V, Jaffré J, Roberts JE. Modeling fractures and barriers as interfaces for flow in porous media. SIAM J Sci Comput. 2005;26(5):1667–1691.
  • Angot P, Boyer F, Hubert F. Asymptotic and numerical modelling of flows in fractured porous media. M2AN. 2009;43(2):239–275.
  • Lesinigo M, D’Angelo C, Quarteroni A. A multiscale darcy-brinkman model for fluid flow in fractured porous media. Num Math. 2011;117(4):717–752.
  • D’Angelo C, Scotti A. A mixed finite element method for Darcy flow in fractured porous media with non-matching grids. M2AN Math Model Numer Anal. 2011;46(2):465–489.
  • Frih N, Martin V, Roberts JE, et al. Modeling fractures as interfaces with non-matching grids. Computational Geosciences. 2012;16:1043–1060. DOI:10.1007/s10596-012-9302-6.
  • Tunc X, Faille I, Gallouët T, et al. A model for conductive faults with non-matching grids. Computational Geosciences. 2012;16(2):277–296. DOI:10.1007/s10596-011-9267-x.
  • Knabner P, Roberts JE. Mathematical analysis of a discrete fracture model coupling darcy flow in the matrix with darcy-forchheimer flow in the fracture. ESAIM M2AN. 2014;48:1451–1472.
  • Brenner K, Groza M, Guichard C, et al. Gradient discretization of hybrid dimensional Darcy flows in fractured porous media. In: Fuhrmann J, Ohlberger M, Rohde C, editors. Finite volumes for complex applications VII Elliptic parabolic and hyperbolic problems. Springer Proceedings in Mathematics & Statistics. Vol. 78. Cham: Springer; 2014. p. 527–535. DOI:10.1007/978-3-319-05591-6\_52.
  • Brenner K, Groza M, Guichard C, et al. Vertex approximate gradient scheme for hybrid dimensional two-phase Darcy flows in fractured porous media. ESAIM Math Model Numer Anal. 2015;49(2):303–330. DOI:10.1051/m2an/2014034.
  • Alboin C, Jaffré J, Roberts JE, et al. Domain decomposition for flow in porous fractured media. In: Lai CH, Bjorstad PE, Cross M, et al, editors. Domain decomposition methods in sciences and engineering. Bergen: Domain Decomposition Press; 1999. p. 365–373.
  • Jaffré J, Martin V, Roberts JE. Generalized cell-centered finite volume methods for flow in porous media with faults. In: Herbin R, Kr{\"o}ner D, editors. Finite volumes for complex applications, III. Paris: Hermes Science Publication; 2002. p. 367–364.
  • Faille I, Flauraud E, Nataf F, et al. A new fault model in geological basin modelling, application to finite volume scheme and domain decomposition methods. In: Herbin R, Kroner D, editors. Finite volumes for complex applications III. Paris: Hermés Penton Sci; 2002. p. 543–550.
  • Tunc X, Faille I, Gallouët T, et al. A model for conductive faults with non-matching grids. Comput Geosci. 2011;16(2):277–296. DOI:10.1007/s10596-011-9267-x.
  • Moreles F, Showalter RE. The narrow fracture approximation by channeled flow. JMAA. 2010;365(1):320–331.
  • Frih N, Roberts JE, Saâda A. Modeling fractures as interfaces: a model for Forchheimer fractures. Comput Geosci. 2008;12:91–104.
  • Engl HW, Hanke M, Neubauer A. Regularization of inverse problems. Mathematics and its applications. Vol. 375. Dordrecht: Kluwer Academic; 1996. DOI:10.1007/978-94-009-1740-8
  • Chavent G. Nonlinear least squares for inverse problems. Scientific computation. New York (NY): Springer; 2009; theoretical foundations and step-by-step guide for applications.
  • Ameur H, Chavent G, Jaffré J. Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities. Inverse Prob. 2002;18(3):775–794. DOI:10.1088/0266-5611/18/3/317.
  • Ameur H, Clément F, Weis P, et al. The multidimensional refinement indicators algorithm for optimal parameterization. J Inverse Ill-Posed Prob. 2008;16(2):107–126. DOI:10.1515/JIIP.2008.008.
  • Ameur H, Chavent G, Clément F, et al. Image segmentation with multidimensional refinement indicators. Inverse Prob Sci Eng. 2011;19(5):577–597; special Issue: Proceedings of the 5th International Conference on Inverse Problems: Modeling and Simulation; 2010 May 24--29; Antalya, Turkey. Available from: http://dx.doi.org/10.1080/17415977.2011.579609.
  • Cheikh F. Identification des fractures dans un milieu poreux [dissertation]. ENIT Tunis (Tunisia), U. Paris 6 (France); Forthcoming. French.
  • 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 aterflow problems. Adv Water Res. 1991;14(6):329–348.

Reprints and Corporate Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

To request a reprint or corporate permissions for this article, please click on the relevant link below:

Academic Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

Obtain permissions instantly via Rightslink by clicking on the button below:

If you are unable to obtain permissions via Rightslink, please complete and submit this Permissions form. For more information, please visit our Permissions help page.