Comput Geosci (2008) 12:91104

DOI 10.1007/s10596-007-9062-x
ORIGINAL PAPER
Modeling fractures as interfaces: a model
for Forchheimer fractures
Najla Frih Jean E. Roberts Ali Saada
Received: 18 May 2007 / Accepted: 6 November 2007 / Published online: 17 January 2008
Springer Science + Business Media B.V. 2007
Abstract In this paper we are concerned with modeling
single phase ow in a medium with known fractures.
In particular we are interested in the case in which the
ow rate in the fractures is large enough to make it
appropriate to use Forchheimers law for modeling the
ow in the fractures even though the ow in the sur-
rounding domain is such that Darcys law is adequate.
We describe a model in which the fractures are treated
as interfaces. We also consider the case of intersecting
fractures and the case of nonconforming meshes.
Keywords Porous media Fractures
Forchheimers law Mixed nite elements
Domain decomposition Nonconforming meshes
1 Introduction
Modeling ow in porous media is made difcult by
the presence of heterogeneities in the characteristics
of the medium. In particular, it is difcult to take into
account fractures in a numerical model, because of their
width which is very small in comparison to the size of
the domain. Yet because their permeability may also
N. Frih (B) A. Saada
ENIT-LAMSIN, BP 37, 1002 Tunis-le Belvdre, Tunisia
e-mail: najla.frih@lamsin.rnu.tn, najla.frih@inria.fr
A. Saada
e-mail: ali.saada@ipein.rnu.tn
J. E. Roberts
INRIA-Rocquencourt, BP 105,
78153 Le Chesnay Cedex, France
e-mail: jean.roberts@inria.fr
differ greatly from that of the surrounding medium,
they may have a very important inuence on the ow
in the medium. They may act as privileged channels for
ow or they may act as barriers.
Fractures may occur as networks of very ne frac-
tures, impossible to localize individually, [9]. These are
taken into account using homogenization techniques
or double porosity models, [7, 8]. Here however, we
are concerned with larger known fractures. In [1], a
model was introduced in which such larger fractures are
modeled as (n 1) dimensional interfaces in an n di-
mensional medium. The fractures were supposed to be
of high permeability so the pressure was assumed to be
continuous across the fractures. However the ux was
not supposed to be continuous as the uid could ow
into and out of as well as along the fractures. A model
introduced in [14] generalizes the earlier model so that
it can handle both large and small permeability in the
fracture. For this model it is no longer assumed that the
pressure is continuous across the fracture. In [2], a 3-D
model with intersecting fractures was introduced. For
all of these models, ow in the fractures as well as in
the surrounding medium is governed by Darcys law.
In an earlier article, [12], a model was introduced in
which the ow in the fracture is sufciently rapid that it
is governed by a nonlinear law, Forchheimers law, see
[6, 11, 17] and [20]. In the present article we develop
the model of [12] in more detail, extending it to the
case of intersecting fractures and giving more numerical
results. We also show that nonmatching grids may be
used for the subdomains and for the fractures.
We mention that others have also treated fractures
as interfaces: for Darcy ow, a model introduced by
Angot et al. in [3] is based on Robin boundary condi-
tions at the interface and assumes the continuity of the
92 Comput Geosci (2008) 12:91104
ux across the fracture. A model similar to that of [14]
was studied by Flauraud et al. in [10]. For a problem
of multiphase ow in a fractured medium, Helmig et al.
in [15], represented the fracture by lower dimensional
nite elements.
In Section 2, we describe the different equations
governing single phase, incompressible ow in a porous
media. In Section 3, we describe a simple model prob-
lem with a domain containing a single fracture. Then,
in Section 4, the interface model of [12] is derived
for this simple problem and in Section 5 we give the
discrete formulation along with some numerical results.
In Section 6, we extend the simple model problem to a
problemwith intersecting fractures. In Section 7, we use
nonconforming meshes and we show some numerical
results with several intersecting fractures.
2 Equations governing ow in a porous medium
Flow of a single-phase, incompressible uid in a porous
medium is governed by the law of mass conservation:
div (u) = q in , (1)
where q is a source term and u is a volumetric ow
rate or velocity. The velocity u is related to the gra-
dient of the pressure p (in the absence of gravity) by
Darcys law:
u + K p = 0, in , (2)
where K is the permeability of the medium.
Darcys law is valid for low ow rates for which
inertial effects are negligible. For higher ow rates
however, the results obtained by Darcys law do not
coincide with experimental results; the ow rate pre-
dicted by Darcys law is too high. Forchheimers law
introduces an inertial term, a quadratic term, which
slows down the ow. So, when the ow is sufciently
rapid, Forchheimers law gives a more accurate relation
between the gradient of the pressure and the ow rate.
This law is a nonlinear law given by
(1 +b |u|) u + K p = 0 in , (3)
where b denotes the Forchheimer constant also called
the inertial constant. In some works b is taken to be
a tensor, but here we have restricted our attention
to the case b is a scalar. For both experimental and
theoretical derivations of Forchheimers law see the list
of references in the introduction of [13].
3 Description of a simple model problem
Suppose that is a convex domain in R
n
, n = 2 or 3,
and denote by = the boundary of . We sup-
pose that the ow in is governed by a conservation
equation together with Forchheimers law relating the
gradient of the pressure p to the ow velocity u:
div u = q in
(1 +b|u|)u = K p in
p = p
D
on , (4)
where q is a source term, b the Forchheimer coefcient,
p
D
the given pressure on the boundary and K a di-
agonal permeability (or hydraulic conductivity) tensor
which is supposed to satisfy
0 < K
min
x
2
0
(Kx, x) K
max
x
2
0
< , x = 0.
For this simple model problem, we suppose (see
Fig. 1) that the fracture
f
, a subdomain of , is such
that there is a hyperplane with unit normal vector n
so that

f
=
_
x : x = s +rn for some s
and some r in the interval
_

d(s)
2
,
d(s)
2
__
,
Fig. 1 Left the domain with the fracture
f
. Right the subdomains
1
and
2
separated by the fracture considered as an interface
Comput Geosci (2008) 12:91104 93
where d(s) denotes the thickness of the fracture at
s . We also suppose that
f
separates into two
disjoint, connected subdomains:
\
f
=
1

2
,
1

2
= .
We will use the notation
i
for the part of the boundary
of
i
which lies on , i = 1, 2, f

i
=
i
, i = 1, 2, f
and we denote by
i
the part of the boundary of
i
,
i = 1, 2 which lies on
f

i
=
i

f
, i = 1, 2.
The case of interest for us is that the velocity in
the subdomains is small enough for Darcys law to be
sufcient while that in the fracture requires the use of
Forchheimers law. In this case in order to avoid having
to solve a nonlinear problem in all of it would seem
appropriate to formulate the problem as a transmission
problem using the nonlinear law only in the fracture. If
we assume that the ow in the subdomains is governed
by Darcys law and if we denote by p
i
, u
i
, K
i
, and q
i
the restrictions of p, u, K, and q respectively to
i
, i =
1, 2, f , and by p
D
i
the restriction of p
D
to
i
, i = 1, 2, f
we can rewrite the above problem (4) as a transmission
problem:
div u
i
= q
i
in
i
, i = 1, 2, f,
u
i
= K
i
p
i
in
i
, i = 1, 2,
(1 +b|u
f
|)u
f
= K
f
p
f
in
f
,
p
i
= p
D
i
on
i
, i = 1, 2, f,
p
i
= p
f
on
i
, i = 1, 2,
u
i
n = u
f
n on
i
, i = 1, 2. (5)
The system (5) could then be solved by using a con-
ventional nonoverlapping domain decomposition tech-
nique in which
f
would be considered simply as a
third subdomain. Following [1, 14] and [2], we pro-
pose here as in [12] the alternative technique of treat-
ing
f
not as a subdomain, but as an interface
between the subdomains
1
and
2
(see Fig. 1) on
which nonlocal transmission conditions are imposed.
The advantage of using domain decomposition for
modeling the fractured medium is that one no longer
needs to have local renement around the fracture.
However if the fracture is taken to be a subdomain
itself this introduces a subdomain much thinner than
other subdomains which is not usually desirable for
domain decomposition algorithms. Further with a stan-
dard domain decomposition approach the elements in
the fractures would either have one dimension much
smaller than the other dimensions or the grids would
be nonconforming requiring the use of mortar elements
or other special techniques. By treating the fracture
not as a thin subdomain but as an interface this prob-
lem is eliminated. An added advantage of the inter-
face approach in the present case with Forchheimer
ow in the fracture is that one has only an (n 1)
dimensional nonlinear problem to solve instead of an n
dimensional one.
4 Derivation of the interface model
The model in which the subdomain
f
is replaced
by the interface , is obtained by using the technique
of averaging across the fracture. For this however, a
simplifying hypothesis, that ow in the fracture in the
direction normal to the central hyperplane is much
smaller than that in the tangential direction, is used.
This hypothesis seems reasonable in light of the very
small ratio, width to length of the fracture. Thus we sup-
pose that the owin the direction normal to the fracture
is adequately described by Darcys law and that in the
tangential direction, it is described by Forchheimers
law.
The rst step toward deriving the model is to decom-
pose u
f
as u
f
= u
f,n
+u
f,
with u
f,n
= (u
f
n) n (recall
that n = n
1
= n
2
), and to introduce the notation

and div

for the tangential gradient and divergence

operators and
n
and div
n
for the normal gradient and
divergence operators.
4.1 Averaging the conservation equation
With the above notation, the rst equation of Eq. 5 for
i = f , may be rewritten as
div
n
u
f
+div

u
f
= q
f
in
f
. (6)
Integrating in the direction normal to the fracture, one
obtains
u
f
n
|
2
u
f
n
|
1
+div

= Q

on , (7)
where U

=
_ d
2

d
2
u
f,
dn is the ow rate through a nor-
mal cross section of the fracture and Q

=
_ d
2

d
2
q
f
dn.
Then using the continuity of the uxes across
1
and

2
, the last equation of Eq. 5 for i = 1 and 2, we may
write
div

= Q

+(u
1
n
1|
1
+u
2
n
2|
2
) on . (8)
94 Comput Geosci (2008) 12:91104
This is the conservation equation on with the addi-
tional source term u
1
n
1|
1
+u
2
n
2|
2
representing the
difference between what ows into the fracture and
what ows out of the fracture at a given point of .
4.2 Averaging Forchheimers law
With the hypothesis that ow in the fracture in the
direction normal to the fracture is described by Darcys
law, the decomposition of the second equation of the
system (5) into its tangential direction and the normal
direction parts is given by
(1 +b |u
f
|) u
f,
= K
f,

p
f
(a)
u
f,n
= K
f,n

n
p
f
. (b) (9)
Integrating the rst equation of the system (Eq. 9, a)
along the line segments
_

d
2
,
d
2
_
, yields:
_ d
2

d
2
(1 +b |u
f
|) u
f,
dn = K
f,

_ d
2

d
2
p
f
dn, (10)
where we have assumed that K
f,
is constant along
the segments
_

d
2
,
d
2
_
. We approximate |u
f
| = |u
f,
+
u
f,n
|by|u
f,
|
_
1 +
1
2
u
2
f,n
u
2
f,
_
.
Then using the hypothesis that u
f,n
is very small in
comparison to u
f,
, we obtain the approximation
|u
f
| |u
f,
|
|U

|
d
.
Replacing |u
f
| by
|U

|
d
in Eq. 10 and integrating one
obtains
_
1 +
b
d
|U

|
_
U

= K
f,
d

, (11)
where P

=
1
d
_ d
2

d
2
p
f
dn is the average pressure along
a normal cross section of the fracture.
Equation 11 is Forchheimers law in the (n 1) di-
mensional domain . Together Eqs. 8 and 11 give a
ow equation in with a source term representing the
ow from the subdomains
1
and
2
into the fracture.
The second equation of Eq. 9 can now be used to give
boundary conditions along for the subdomains
1
and
2
.
Integrating the second equation of the system
(Eq. 9, b) along the line segments
_

d
2
,
d
2
_
, yields:
_ d
2

d
2
u
f,n
= K
f,n
_ d
2

d
2

n
p
f
, (12)
where we have assumed that K
f,n
is constant along the
segments
_

d
2
,
d
2
_
.
Equation 12 is equivalent to
p
2|
p
1|
=
1
K
f,n
_ d
2

d
2
u
f,n
dn (13)
With the hypotheses that the fracture width d is much
smaller than the fracture length, and the high perme-
ability in the fracture, i.e K
f,n
is large, then we can
suppose that the pressure p is continuous through
so that
p
1|
= P

= p
2|
.
4.3 The strong formulation of the interface model
problem
The interface model for the fracture problem can now
be written
div u
i
= q
i
in
i
, i =1, 2
u
i
= K
i
p
i
in
i
, i =1, 2
div

= Q

+(u
1
nu
2
n) on
_
1+
b
d
|U

|
_
U

=K
f,
d

on
p
i
= p
Di
on
i
, i =1, 2
p
i
= P

on , i =1, 2
P

= P
D
on ,
(14)
where P
D
is the average value of P
Df
along the seg-
ment
_

d
2
,
d
2
_
in
f
.
5 Numerical discretization for the interface model
5.1 The discrete formulation of the model problem
Let T
h,i
be a conforming nite element partition of

i
, i = 1, 2, and suppose that the meshes T
h,i
, i = 1, 2
match at the interface between the subdomains
i
,
i.e. suppose that these partitions are such that the
two partitions induced on by these partitions T
h,i
coincide. In other words T
h,1
T
h,2
forms a standard
nite element partition of . Then denote by T
h,
the
Comput Geosci (2008) 12:91104 95
induced mesh on consisting of all edges of elements
of T
h,i
that lie on . Then let T
h
= T
h,i
, i = 1, 2, . Thus
T
h
consists of both n-dimensional elements in
i
and
(n 1)-dimensional elements on .
To dene an approximation space for the pressure,
let M
h,i
, i = 1, 2 be the space of piecewise constant
functions constant on each element T of T
h,i
and let
M
h,
be the approximation space of piecewise constant
functions constant on each element E of T
h,
. Then let
M
h
= M
h,1
M
h,2
M
h,
.
Similarly, to dene an approximation space for the
velocity, let W
h,i
, i = 1, 2 be the RaviartThomas
Nedelec space of lowest order dened on the n-
dimensional space
i
subordinate to the mesh T
h,i
,
and let W
h,
be the RaviartThomas space of lowest
order dened on the (n 1)-dimensional interface
subordinate to the mesh T
h,
. Then let
W
h
= W
h,1
W
h,2
W
h,
.
The discrete mixed nite element approximation
for the interface problem (14) can now be written as
follows:
Find (u
h,1
, u
h,2
, U
h,
)W
h
and ( p
h,1
, p
h,2
, P
h,
)M
h
such that for all (v
h,1
,v
h,2
,V
h,
)W
h
and (r
h,1
, r
h,2
,
r
h,
) M
h
_

i
K
1
i
u
h,i
v
h,i

_

i
p
h,i
divv
h,i
=
_

P
h,
v
h,i
n
i

_

i
p
D
i
v
h,i
n
i
, i = 1, 2
(P
h
)
_

i
divu
h,i
r
h,i
=
_

i
q
i
r
h,i
, i = 1, 2
_

(d K
f,
)
1
_
1 +
b
d
|U
h,
|
_
U
h,
V
h,

P
h,
div

V
h,
=
_

P
D
V
h,
n

div

U
h,
r
h,
=
_

r
h,
+
_

(u
h,1
n
1
+u
h,2
n
2
)r
h,
.
(15)
Following [1, 2, 14], we use a domain decomposition
type method [18] to solve the discrete problem (P
h
).
In this way, the problem (15) is reduced to a problem
only on the interface . For i = 1, 2, u
h,i
and p
h,i
are
decomposed into two parts,
u
h,i
= u
0
h,i
+u

h,i
, p
h,i
= p
0
h,i
+ p

h,i
.
For each i, we use the SteklovPoincar operator
S
i
: M
h,
M
h,
r
h,
u
0
h,i
n
i
which associates to an element r
h,
M
h,
the element
u
0
h,i
n
i
M
h,
where (u
0
h,i
, i, p
0
h,i
) W
h,i
M
h,i
is the
solution of the problem
Find (u
0
h,i
, p
0
h,i
) W
h,i
M
h,i
such that
_
P
0
i
_
_

i
K
1
i
u
0
h,i
v
h,i

_

i
p
0
h,i
divv
h,i
=
_

r
h,
v
h,i
n
i
v
h,i
W
h,i
_

i
divu
0
h,i
r
h,i
= 0 r
h,i
M
h,i
. (16)
We also use the notation
i
= u

h,i
n
i
i = 1, 2, where
(u

h,i
, p

h,i
) W
h,i
M
h,i
is the solution of the problem
Find (u

h,i
, p

h,i
) W
h,i
M
h,i
such that
_
P

i
_
_

i
K
1
i
u

h,i
v
h,i

_

i
p

h,i
divv
h,i
=
_

i
p
D
v
h,i
n
i
v
h,i
W
h,i
_

i
divu

h,i
r
h,i
=
_

i
q
i
r
h,i
r
h,i
M
h,i
. (17)
Then if ( p
h,1
, p
h,2
, P
h,
) with (u
h,1
, u
h,2
, U
h,
) is the
solution of (P
h
); (u
0
h,i
, p
0
h,i
) is the solution of (P
0
i
) with
r
h,
= P
h,
and (u

h,i
, p

h,i
) is the solution of (P

i
) one has
p
h,i
= p
0
h,i
+ p

h,i
and u
h,i
= u
0
h,i
+u

h,i
. (18)
Thus u
h,i
n
i
= S
i
(P
h,
) +
i
, and (U
h,
, P
h,
) is the
solution of the following problem:
Find (U
h,
, P
h,
) W
h,
M
h,
such that for all
r
h,
M
h,
and V
h,
W
h,
(P

)
_

(d K
f,
)
1
_
1 +
b
d
|U
h,
|
_
U
h,
V
h,

P
h,
div

V
h,
=
_

p
D
V
h,
n

div

U
h,
r
h,
+
_

(S
1
(P
h,
) +S
2
(P
h,
))r
h,
=
_

r
h,
+
_

(
1
+
2
)r
h,
. (19)
This problem is nonlinear but it is only an (n-1) di-
mensional nonlinear problem, and it can be solved by
96 Comput Geosci (2008) 12:91104
Fig. 2 Test-case
using a xed point iteration method or a quasi-Newton
method. In either case the iteration can be initialized
with the solution for Darcy ow in the fracture. The
solutions (u
h,i
, p
h,i
) in the subdomains can then be
calculated using Eq. 18.
5.2 Some 2-D numerical results
A rst numerical result for this type of problem is to
obtain a good correspondence between the solution
obtained by the interface problem and a reference solu-
tion. We obtain a reference solution by using the model
problem where the fracture is considered as a third
subdomain and using a mixed nite element method.
5.2.1 A simple test case
In this section numerical results obtained with the inter-
face model will be compared with those obtained using
a standard model with 2-D fracture. The results are
obtained for a simple model problem: the subdomains

1
and
2
are both squares of unit length and width,
and the fracture
f
separating the subdomains is of
unit length and of width d = 0.01. The permeability
in each of the subdomains is assumed to be constant
and equal to 10
9
, and the fracture is assumed to be of
much higher (also constant) permeability, K
f
= 10
6
.
The Forchheimer coefcient b is taken to be 10. The
upper and lower boundaries of the two subdomains are
assumed to be impermeable, and there is a pressure
drop from right to left of 10
6
. The same pressure drop
from top to bottom of the fracture is imposed. See
Fig. 2.
To compare the different results, one can see in Fig. 3
the reference pressure and also the reference velocity,
both obtained by using the standard model (model with
2-D fracture). In Fig. 4, those obtained by using the
interface model are given. In both cases a 15 by 15
grid was used for each of the subdomains
1
and
2
,
and in the case of the reference model a 15 by 15 grid
was also used for the fracture domain
f
. As we can
see, these results show good agreement between the
solutions obtained by the two models.
5.2.2 Comparison between Darcy and Forchheimer
velocity in the fracture
In the experiment above the data were chosen so that
the Forchheimer term would have a non negligible ef-
fect. Here a comparison between the Darcy and Forch-
heimer velocities in the fracture is given. The same data
as in the above experiment are used only in one case the
Forchheimer term in the fracture is omitted.
In Fig. 5 the tangential Darcy velocity (which for
this simple 2-D model is just a scalar) is represented
in blue while the Forchheimer velocity is given by the
0
0.5
1
1.5
2
2.5
0
0.2
0.4
0.6
0.8
1
0
2
4
6
8
10
x 10
5
x
Standard model
y
P
r
e
s
s
u
r
e
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
V
e
l
o
c
i
t
y

i
n

t
h
e

t
o
t
a
l

d
o
m
a
i
n
Standard Model
Fig. 3 Reference pressure (left) and reference velocity (right) given by the standard model
Comput Geosci (2008) 12:91104 97
0
0.5
1
1.5
2
2.5
0
0.2
0.4
0.6
0.8
1
0
2
4
6
8
10
x 10
5
x
Interface Model
y
P
r
e
s
s
u
r
e
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
V
e
l
o
c
i
t
y

i
n

t
h
e

t
o
t
a
l

d
o
m
a
i
n
Interface Model
Fig. 4 Numerical pressure (left) and numerical velocity (right) given by the interface model
red curve. On the left, the velocity calculated using
the standard model along a central vertical section is
shown, and on the right is the velocity obtained using
the interface model.
A good agreement between the results obtained by
the standard model and the interface model can be
observed and in addition, the difference between the
Forchheimer velocity and the Darcy velocity can be
seen. The Forchheimer velocity is of much smaller
magnitude than the Darcy velocity because of the large
inertial coefcient.
Indeed, in this example the Forchheimer coefcient
b was taken to be extremely large to show that the
model is valid even when the nonlinear term is very
important. For more physically correct values of b, the
percent decrease in the magnitude of the velocity when
the Forchheimer lawis used instead of simply the Darcy
law is much smaller. To illustrate this, we have shown
in Fig. 6, the average percent decrease in the velocity
due to the Forchheimer term for different values of
the Forchheimer coefcient b. The percentage calcu-
lated using the interface model is represented in blue
and that obtained using the standard model is given
in red.
5.2.3 Numerical study of the error
The object of this paragraph is to obtain a quantitative
appreciation of the error made when using the inter-
face model. A relative L
2
error in the pressure in the
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
x
V
e
l
o
c
i
t
y

i
n

t
h
e

i
n
t
e
r
f
a
c
e
Interface model
Forchheimer velocity
Darcy velocity
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
x
V
e
l
o
c
i
t
y

i
n

t
h
e

f
r
a
c
t
u
r
e
Standard model
Forchheimer velocity
Darcy velocity
Fig. 5 Darcy and Forchheimer velocity in the fracture given both by the interface model (left) and the standard model (right)
98 Comput Geosci (2008) 12:91104
0 1 2 3 4 5 6 7 8 9 10
0
10
20
30
40
50
60
70
80
90
100
b
a
v
e
r
a
g
e

%

v
e
l
o
c
i
t
y

d
e
c
r
e
a
s
e
Average percent velocity decrease
due to Forchheimer term for different b values
interface model
standard model
Fig. 6 Average percent velocity decrease due to Forchheimer
term for different values of the Forchheimer coefcient b
subdomains is calculated. A reference solution P

is
obtained using the standard model with 2-D fracture,
using mixed nite elements on a regular ne mesh T

of
square elements of side length , with = 1/192 in the
subdomains
1
and
2
. Approximate solutions P
h
for
different values of h are calculated using the interface
model with a uniform mesh of square elements of side
length h in each subdomain and a 1-D mesh in the
fracture with elements of length h.
The discretization parameter h for each of the cal-
culations below is chosen such that the ne mesh in the
subdomains is a renement of the mesh of size h. Then a
projection

P
h
of each approximate solution P
h
onto
the ne mesh is obtained in the obvious way. Then the
relative L
2
error in the pressure is given by
P
h
P

2
L
2
rel
()
=

P
h
P

)
2
|C

(P

)
2
|C

|
(20)
where |C

| is the measure of the cell C

.
In these experiments there are two sources of error,
the usual numerical error dependent upon the discreti-
sation parameter h and an error due to the interface
model dependent upon the width d of the fracture.
In Fig. 7, on the left, the log of the relative L
2
pres-
sure error in the total domain is plotted as a function of
the log of the fracture width d. Each curve corresponds
to a different value of the mesh size h: h = 1/6, h =
1/12, h = 1/24, h = 1/48, h = 1/96, h = 1/192. On the
right, the log of the relative L
2
pressure error in the
total domain is plotted as a function of the log of 1/h,
the inverse of the mesh size, for a constant fracture
width d = 0.01.
In the rst ve curves on the left there is no de-
crease in the error when the fracture width diminishes
because the error due to the model is dominated by the
error due to the numerical discretization. Only in the
sixth curve with h = 1/192 can we see that the error
diminishes with decreasing d just until the numerical
discretization error again becomes dominant. In the
curve on the right, as expected, the error decreases
as 1/h increases. The large drop in the error between
1/h = 96 and 1/h = 192 is due to the fact that the
reference solution is not an analytic solution but a
solution obtained with a mesh for which h = 1/192.
10
3
10
2
10
1
10
5
10
4
10
3
10
2
10
1
10
0
Log (width of the fracture)
L
o
g

(
L
2

E
r
r
o
r

p
r
e
s
s
u
r
e
)
Mesh 192X192
Mesh 96X96
Mesh 48X48
Mesh 24X24
Mesh 12X12
Mesh 6X6
10
0
10
1
10
2
10
3
10
5
10
4
10
3
10
2
10
1
10
0
Log (1/ h)
L
o
g

(
L
2

E
r
r
o
r

p
r
e
s
s
u
r
e
)
error= f(1/ h)
Fig. 7 Left the relative L
2
error in the pressure in the total
domain as a function of the fracture width d for different values
of the mesh size: h = 1/6, h = 1/12, h = 1/24, h = 1/48, h =
1/96, h = 1/192. Right the relative L
2
error variation in the total
domain as a function of 1/h for a constant fracture width d = 0.01
Comput Geosci (2008) 12:91104 99
6 A model with intersecting fractures
6.1 Description of the problem
For this model problem, as earlier, is supposed to
be a convex domain in R
n
, n = 2 or 3, with boundary
= decomposed into a Dirichlet part
D
and a
Neumann part
N
. For simplicity, we suppose that there
are three fractures
k
, k = 1, 2, 3, dividing

into three
subdomains
i
, i = 1, 2, 3,

=
3
_
i=1

i
,
i

j
= and

j
=
k
i, j, k
with i = j = k = i,
in such a way that the intersection of each pair of
distinct fractures is the same set T:

1

2
=
1

3
=
2

3
= T.
The intersection T is thus a point in the 2-D case and a
segment in the 3-D case. See Fig. 8.
As before, in each subdomain
i
, i = 1, 2, 3, we
have the law of mass conservation together with
Darcys law and the boundary conditions on the exte-
rior boundary:
div u
i
= q
i
in
i
u
i
= K
i
p
i
in
i
p
i
= p
Di
on
i

D
u
i
n
i
= 0 on
i

N
. (21)
Also as before, we assume that the pressure is con-
tinuous across the fracture interfaces
k
:
p
i
= P

k
= p
j
on
k
, i = j = k = i, (22)
Fig. 8 A domain with intersecting fractures
where P

k
is the pressure on the fracture
k
. In each
fracture we have the law of mass conservation together
with Forchheimers lawand the boundary conditions on
the exterior boundary:
div

k
=Q

k
+
_
u
i
n
i
+u
j
n
j
_
on
k
_
1+
b

k
d

k
|U

k
|
_
U

k
= d

k
K

k
on
k
P

k
= P
D
k
on
k
,
(23)
where U

k
, Q

k
, K

k
, b

k
and d

k
are respectively the
pressure, the ow velocity, the source term, the per-
meability tensor, the Forchheimer coefcient and the
width of the fracture for
k
, and P
D
k
is the average
value of p
D
on
k
. On the intersection T of the
fractures, following [1] and [2], we impose the continu-
ity of the pressure and also the continuity of the ux:
P

1
= P

2
= P

3
on T
U

1
+U

2
+U

3
= 0 on T, (24)
where

k
is the exterior normal to
k
on
k
, k =
1, 2, 3.
6.2 Numerical discretization
As in the case of a single simple fracture, to obtain
a numerical solution of the system (Eqs. 21, ..., 24),
we use mixed nite elements of dimension n in the
subdomains and mixed nite elements of dimension
(n 1) in the fractures. Thus as before the pressure
is approximated in each element (in dimension n and
in dimension n 1) by a constant, and the velocity is
approximated in the space of lowest order Raviart
ThomasNedelec elements. However, in addition to the
unknowns associated with these spaces there are scalar
unknowns associated with the intersection T. In the
case R
2
so that T is simply a point this amounts
to one unknown which is a multiplier for enforcing
continuity of the pressure at the fracture intersection.
In the case R
3
so that T is a segment, there is a
naturally induced mesh on T from each fracture having
T as part of its boundary, and as we have assumed here
that the meshes are compatible, i.e. the union of the
meshes for the subdomains forms a conforming mesh
on all , the meshes induced on T from the different
fractures having T as part of their boundaries are all the
same. Thus there is one additional unknown for each
element of this mesh.
Again as the nonlinearity of the problem is asso-
ciated only with the fractures, domain decomposition
techniques can be used to reduced the problem to a
100 Comput Geosci (2008) 12:91104
problem posed only on the fractures. The resulting
nonlinear system can now be written as follows:
_
_
_
_
_
_
_
_
_
_
A
1
(U

1
) B

1
0 0 0 0 M

1
B
1
S
1,1
0 S
1,2
0 S
1,3
0
0 0 A
2
(U

2
) B

2
0 0 M

2
0 S
2,1
B
2
S
2,2
0 S
2,3
0
0 0 0 0 A
3
(U

3
) B

3
M

3
0 S
3,1
0 S
3,2
B
3
S
3,3
0
M
1
0 M
2
0 M
3
0 M
T
_

_
_
_
_
_
_
_
_
_
_
_
U

1
P

1
U

2
P

2
U

3
P

3
P
T
_

_
=
_
_
_
_
_
_
_
_
_
_

0
_

_
(25)
where U

k
and P

k
for k = 1, 2, 3, are the velocity and
the pressure respectively on the interface
k
and P
T
the trace of the pressure on the intersection T of the
fractures. The matricies A
k
(U

k
) and B
k
, represent the
mixed nite element matricies for the ow equations
in the fracture
k
, k = 1, 2, 3. The matrix S
i, j
, i, j =
1, 2, 3, represents the domain decomposition matrix
associated with the SteklovPoincar operator taking
into account the values of the pressure on the fracture

i
and returning the values of the ux on the fracture
j
.
The matricies M
i
, i = 1, 2, 3, T, are the matrices associ-
ated with the equations of continuity of the pressure
and continuity of the ux on the intersection of the
fractures.
6.3 A numerical result for the case of intersecting
fractures
In this experiment the domain R
2
is divided into
4 subdomains
i
, i = 1, . . . , 4 which are separated by
4 fractures
k
, k = 1, . . . , 4. These fractures intersect at
the point T; see Fig. 9. The width d

k
is the same in
all fractures and is equal to 0.01. The permeability in
each of the subdomains is assumed to be constant and
equal to 10
9
, and the permeability in each fracture
is assumed to be much higher (also constant) K

k
=
10
6
. The Forchheimer coefcient b

k
is the same in all
Fig. 9 A test-case
fractures and is taken to be 10. The upper and lower
exterior boundaries of the subdomains are assumed to
be impermeable, and there is a pressure drop from right
to left of 10
6
. A pressure of 10
6
is imposed at the ends of
the fractures lying on the upper and right-hand sides of
the domain and a pressure of 0 is imposed on the ends
of the fractures on the left-hand and lower sides of the
domain.
In Fig. 10 the conforming mesh used for this test case
is shown. In Fig. 11 the calculated pressure is shown on
the left and the calculated velocity on the right.
7 Nonconforming meshes
In many applications one may want to use nonconform-
ing meshes. For simplicity, in this section, we suppose
that the domain is divided into only two subdomains
separated by a single fracture. In the numerical discre-
tization, the domain will be decomposed into nonover-
lapping subdomains with the fracture as the interface
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 10 The mesh
Comput Geosci (2008) 12:91104 101
0 0.5 1 1.5 2
0
0.2
0.4
0.6
0.8
1
Fig. 11 The numerical pressure (left) and the numerical velocity (right)
between the two subdomains. The grids will be dened
independently in each subdomain and in the fracture,
see Fig. 12.
Let T
h,i
be a conforming nite element partition of

i
, i = 1, 2, and let T
h,
be a conforming nite element
partition of the interface in (n 1)-dimensions. The
meshes T
h,1
and T
h,2
need not form a conforming mesh
on and T
h,
is independent of T
h,1
and of T
h,2
. The
notation T
h,,i
, i = 1, 2, will be used for the partition
induced on by the partition T
h,i
on
i
; i.e. T
h,,i
,
i = 1, 2 consists of all of the (n 1) dimensional faces
of elements of T
h,i
which lie on . See Fig. 13.
As above, we use the lowest order RaviartThomas
Nedelec mixed nite element spaces associated with
T
h,i
, i = 1, 2, and T
h,
so that the vector variable is
approximated in a space W
h
= W
h,1
W
h,2
W
h,
and the spaces of piecewise constants associated with
T
h,i
, i = 1, 2, and T
h,
so that the scalar variable is ap-
proximated in a space M
h
= M
h,1
M
h,2
M
h,
.
Also as above, domain decomposition techniques
are used to reduce the discrete problem to a problem
posed only on the interface , but here we can not
use directly the operator S
i
: M
h,
M
h,
as dened
in Section 7 because of the incompatibility between
the grids. Letting M
h,,i
denote the space of piecewise
constants associated with T
h,,i
, i = 1, 2, we will, with
some slight abuse of notation, now denote by S
i
the
mapping S
i
: M
h,,i
M
h,,i
dened in the obvious
manner. Then, let R
i
: M
h,
M
h,,i
be the projec-
tion operator dened by
R
i
: M
h,
M
h,,i
r
h,
r
h,,i
where
r
h,,i
|E
=
1
|E|
_
E
r
h,
E T
h,,i
.
Fig. 12 An example of
nonconforming grids for the
fracture and the subdomains
102 Comput Geosci (2008) 12:91104
Fig. 13 An example of a
numerical discretization with
nonconforming grids for the
fracture and the subdomains
So, for the case of nonmatching grids the Steklov
Poincar operator S
i
: M
h,
M
h,
of Section is re-
placed by the operator

S
i
: M
h,
M
h,
with

S
i
= R
t
i
S
i
R
i
,
where R
t
i
is just the transpose of the projection oper-
ator R
i
.
In the numerical experiments below, no mortar ele-
ments were introduced as there are already (n 1)-D
elements on the fracture interfaces, nor were the usual
size restrictions associated with mortar elements, c.f.
[4, 5] and [19], respected. We have shown (Frih et al.,
in preparation) that this is sufcient in the linear case,
i.e. in the case of Darcy ow in the fractures, for the
model of [14].
Remark In the case of intersecting fractures in a
3-D domain, one would need either to have compatible
meshes at the intersection T or to use mortar elements
on T.
7.1 A numerical experiment for nonconforming
meshes
In this experiment the domain R
2
is divided into
4 subdomains
i
, i = 1, ..., 4, by 5 fracture interfaces

k
, k = 1, ...5. The fractures
1
,
2
, and
3
intersect at
Fig. 14 Test-case
a point T
1
and the fractures
3
,
4
, and
5
intersect at a
point T
2
, see Fig. 14.
The width d

k
is the same for all the fractures and is
equal to 0.01. The permeability in each of the subdo-
mains is assumed to be constant and equal to 10
9
, and
the permeability in each fracture is assumed to be much
higher and also constant: K

k
= 10
6
for each k. The
Forchheimer coefcient b

k
is taken to be 10 in all of
the fractures. The upper and lower exterior boundaries
of the subdomains are assumed to be impermeable, and
there is a pressure drop from right to left of 10
6
. A
pressure of 10
6
is imposed at the exterior ends of
4
and

5
and a pressure of 0 is imposed on the exterior ends
of
1
and
2
.
The nonconforming mesh used for this test case is
shown in Fig. 15. It is made up of independently chosen
grids for the subdomains and grids for the fracture
interfaces. Here for simplicity each of the one dimen-
sional grids is a uniform grid. The mesh parameters
for the two dimensional subdomain grids and for the
one dimensional fracture interface grids are given in
the tables below, where by mesh parameter of a two or
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 15 The nonconforming meshes
Comput Geosci (2008) 12:91104 103
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 16 The pressure (left) and the velocity (right) computed solutions
three dimensional grid is meant the maximum length
of an edge of an element of the mesh and of a one
dimensional grid the maximum length of a segment of
the mesh.
Mesh parameters for the subdomain girds
Subdomain
1

2

3

4
Mesh parameter 0.20 0.10 0.08 0.05
Mesh parameters for the fracture interface grids
Interface
1

2

3

4

5
Mesh 0.035 0.029 0.07 0.048 0.11
parameter
The numerical solution obtained with the mesh
shown in Fig. 15 is given in Fig. 16, where the approxi-
mate pressure is shown on the left and the approximate
velocity eld on the right. Agood solution was obtained
in spite of the nonconforming gird, even though no
additional mortar elements were introduced.
8 Conclusion
This article is concerned with a numerical model, in-
troduced in the short article [12], for ow in a porous
medium with fractures. For this model it was assumed
that the ow rate in the fractures is large enough to
make it appropriate to use Forchheimers law for mod-
eling the ow in the fractures even though the ow
in the surrounding domain is such that Darcys law
is adequate. The fractures were treated as interfaces
between subdomains and nonlocal, nonlinear transmis-
sion conditions were imposed on the interfaces. Then
domain decomposition techniques were used to reduce
the problem to a problem posed on the interfaces.
In this article this model is described in more detail
and numerical studies are given. The model is extended
to the case of intersecting fractures and to the case of
nonconforming meshes.
Acknowledgements The authors would like to thank Rachida
Bouhlila and Jrome Jaffr for helpful conversations concerning
this article.
References
1. Alboin, C., Jaffr, J., Roberts, J.E.: Domain decomposition
for ow in fractured porous media. In: Domain Decompo-
sition Methods in Sciences and Engineering, pp. 365373.
Domain Decomposition Press, Bergen (1999)
2. Amir, L., Kern, M., Martin, V., Roberts, J.E.: Dcomposition
de domaine pour un milieu poreux fractur: un modle en 3D
avec fractures qui sintersectent. ARIMA 5, 1125 (2006)
3. Angot, Ph., Gallouet, T., Herbin, R.: Convergence of -
nite volume method on general meshes for non smooth
solution of elliptic problems with cracks. In: Vilsmeier, R.,
Benkhaldoun, F., Hanel, D. (eds.) In: Proceedings of the
2nd International Symposium Finite Volumes for Complex
Applications, pp. 215222. Herms Duisberg (1999)
4. Arbogast, T., Cowsar, L.C., Wheeler, M.F., Yotov, I.: Mixed
nite element methods on nonmatching multiblock grids.
SIAM J. Numer. Anal. 37, 12951315 (2000)
5. Bernardi, C., Maday, Y., Patera, A.T.: A new nonconform-
ing approach to domain decomposition: the mortar element
method. In: Brezis, H., Lions, J.-L. (eds.) In: Nonlinear Par-
tial Differential Equations and Their Applications, pp. 1351.
Pitman, NewYork (1994)
6. Douglas, J., Paes Leme, P.J.S., Giorgi, T.: Generalized
Forchheimer ow in porous media. In: Lions, J.-L.,
Baiocchi, C. (eds.) Boundary Value Problems for Partial
Differential Equations and Applications. Research Notes in
104 Comput Geosci (2008) 12:91104
Applied Mathematics, vol. 29, pp. 99113. Masson, Paris
(1993)
7. Douglas, Jr. J., Arbogast, T.: Dual porosity models for
ow in naturally fractured reservoirs. In: Cushman, J.H.
(ed.) Dynamics of Fluids in Hierarchial Porous Formations,
pp. 177221. Academic Press (1990)
8. Douglas, J., Arbogast, T., Hornung, U.: Derivation of the
double porosity model of single phase ow via homogeniza-
tion theory. SIAM J. Math. Anal. 21, 823836 (1990)
9. Erhel, J., Dreuzy, J.-R.: Efcient algorithms for the determi-
nation of the connected fracture network and the solution of
the steady-state ow equation in fracture networks. Comput.
Geosci. 29(1), 107111 (2003)
10. Faille, I., Flauraud, E., Nataf, F., Pegaz-Fiornet, S.,
Schneider, F., Willien, F.: A new fault model in geologi-
cal basin modelling, application to nite volume scheme and
domain decomposition methods. In: Herbin, R., Kroner, D.
(eds.) In: Finite Volumes for Complex Applications III, pp.
543550. Herms Penton Sci. (2002)
11. Forchheimer, P.: Wasserbewegung durch Boden. Z. Ver.
Deutsch. Ing. 45, 17821788 (1901)
12. Frih, N., Roberts, J.E., Sada, A.: Un modle Darcy
Forchheimer pour un coulement dans un milieu poreux
fractur. ARIMA 5, 129143 (2006)
13. Giorgi, T.: Derivation of the Forchheimer lawvia matched as-
ymptotic expansions. Trans. Porous Media 29, 191206 (1997)
14. Martin, V., Jaffr, J., Roberts, J.E.: Modeling fractures and
barriers as interfaces for ow in porous media. SIAM J. Sci.
Comput. 26, 16671691 (2005)
15. Reichenberger, V., Jakobs, H., Bastian, P., Helmig, R.:
A mixed-dimensional nite volume method for multiphase
ow in fractured porous media. Adv. Water Resour. 29(7),
10201036 (2006)
16. Roberts, J.E., Thomas, J.-M.: Mixed and hybrid methods. In:
Ciarlet et, P.G., Lions, J.-L. (eds.) Handbook of Numerical
Analysis 2. Finite Element Methods part 1, pp. 523639.
Elsevier, North Holland (1991)
17. Panlov, M., Fourar, M.: Physical splitting of nonlinear
effects in high-velocity stable ow through porous media.
Adv. Water Resour. 29, 3041 (2006)
18. Quateroni, A., Valli, A.: Domain Decomposition Methods
for Partial Differential Equations. Claderon Press, Oxford
(1999)
19. Wheeler, M.F., Yotov, I.: Multigrid on the interface for
mortar mixed nite element methods for elliptic problems.
Comput. Methods Appl. Mech. Eng. 184, 287302 (2000)
20. Witaker, S.: The Forchheimer equation: theoretical develop-
ment. Trans. Porous Media 25, 2761 (1996)