Benchmarks for single-phase flow in
fractured porous media
Bernd Flemisch1, Inga Berre2, Wietse Boon2, Alessio Fumagalli2,
Nicolas Schwenck1, Anna Scotti3, Ivar Stefansson2,
and Alexandru Tatomir4

arXiv:1701.01496v1 [math.NA] 5 Jan 2017

1

Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart,
Pfaffenwaldring 61, 70569 Stuttgart, Germany, bernd@iws.uni-stuttgart.de
2

Department of Mathematics, University of Bergen, Allégaten 41, 5007 Bergen, Norway
3

Laboratory for Modeling and Scientific Computing MOX, Politecnico di Milano,
p.za Leonardo da Vinci 32, 20133 Milano, Italy
4

Department of Applied Geology, Geosciences Center, University of Göttingen,
Goldschmidtstrasse 3, 37077 Göttingen, Germany

January 9, 2017
This paper presents several test cases intended to be benchmarks for numerical schemes for single-phase fluid flow in fractured porous media. A number
of solution strategies are compared, including a vertex and a cell-centered finite volume method, a non-conforming embedded discrete fracture model, a
primal and a dual extended finite element formulation, and a mortar discrete
fracture model. The proposed benchmarks test the schemes by increasing the
difficulties in terms of network geometry, e.g. intersecting fractures, and physical parameters, e.g. low and high fracture-matrix permeability ratio as well as
heterogeneous fracture permeabilities. For each problem, the results presented
by the participants are the number of unknowns, the approximation errors in
the porous matrix and in the fractures with respect to a reference solution, and
the sparsity and condition number of the discretized linear system. All data
and meshes used in this study are publicly available for further comparisons.

1 Introduction
In porous-media flow applications, the domains of interest often contain geometrically
anisotropic inclusions and strongly discontinuous material coefficients that can span several orders of magnitude. If the size of these heterogeneities is small in normal direction
compared to the tangential directions, these features are called fractures. Fractures can

1

act both as conduits and barriers and affect flow patterns severely. Target applications
concerning fractured porous-media systems in earth sciences include groundwater resource
management, renewable energy storage, recovery of petroleum resources, radioactive waste
reposition, coal bed methane migration in mines, and geothermal energy production.
The analysis and prediction of flow in fractured porous media systems are important
for all the aforementioned applications. Many different conceptual and numerical models
of flow in fractured porous-media systems can be found in the literature. Even though
fractured porous-media systems have been of interest to modelers for a long time, they still
present challenges for simulators. During the last 70 years, different modeling approaches
have been developed and gradually improved. Comprehensive reviews can be found in
Berkowitz (2002); Dietrich et al. (2005); Hoteit and Firoozabadi (2008); Neumann (2005);
Sahimi (2011); Singhal and Gupta (2010). Roughly, the fractured porous media systems are
classified in two broad categories: discrete fracture-matrix (DFM) models and continuum
fracture models. Within this paper, we will only consider DFM models.
The DFM models consider flow occurring in both the fracture network and the surrounding rock matrix. They account explicitly for the effects of individual fractures on the fluid
flow. An efficient way to represent fractures in DFMs is the hybrid-dimensional approach,
e.g. Helmig (1997); Firoozabadi and Monteagudo (2004); Karimi-Fard et al. (2004); Martin et al. (2005); Reichenberger et al. (2006), where fractures in the geometrical domain
are discretized with elements of co-dimension one with respect to the dimension of the
surrounding matrix, such as one-dimensional elements in two-dimensional settings. The
aforementioned classical DFM approaches all rely on matching fracture and matrix grids in
the sense that a fracture element coincides geometrically with co-dimension-one mesh entities, i.e. faces of matrix grid elements. In addition to the classical models, several so-called
non-conforming DFM models have been developed in recent years, such as EDFM (Moinfar
et al., 2014a; Hajibeygi et al., 2011), XFEM-based approaches (D’Angelo and Scotti, 2012;
Schwenck et al., 2015; Huang et al., 2011), or mortar-type methods (Frih et al., 2012).
Benchmarking represents a methodology for verifying, testing and comparing the modeling tools. Various codes have been developed by academic institutions or companies based
on different conceptual, mathematical, and numerical models. Even though benchmarking
studies are increasing in all fields of engineering and workshops have been organized around
specific problems (e.g. Class et al. (2009)), there are still only a limited number of studies
of this type in the field of geoscience. Some are related to a specific application and are
flexible as to how the problem is modeled in terms of assumptions regarding the physics and
the selection of the domain, see Dahle et al. (2010); Nordbotten et al. (2012); Caers (2013);
Kolditz et al. (2015). Others (De Dreuzy et al., 2013; Caers, 2013), like ours, focus on the
comparison of numerical schemes. One of the common requirements when selecting the
test problems for comparing numerical schemes is that they allow the examination of the
capabilities of each of the compared methods. Therefore, our benchmark study proposes
a set of problems starting from simple geometries and then gradually increasing the geometrical complexity. The test problems are specifically selected to make clear distinctions
between the different methods.
The main focus of this work is to use existing and new computational benchmarks for
fluid flow in fractured porous media to compare several DFM-based numerical schemes in
a systematic way. We would also like to invite the scientific community to follow up on
this study and evaluate further methods by means of the proposed benchmarks. In order
to facilitate this, the paper is accompanied by grid and result files in the form of a Git
repository at https://git.iws.uni-stuttgart.de/benchmarks/fracture-flow.

2

The remainder of this paper is organized as follows. In Section 2, we formulate the
model problem in terms of the partial differential equation to be solved. The participating
DFM models are described in Section 3. The central Section 4 proposes the benchmarks
and compares the results of the different methods. Finally, Section 5 concludes with a
summary and outlook.

2 The model problem
We are considering an incompressible single-phase flow through a porous medium, assumed
to be described by Darcy’s law, resulting in the governing sytem of equations
u = −K grad p,
div u = q,

(1a)
(1b)

in an open bounded domain D ⊂ RN , subject to boundary conditions
p = pD on ∂DD ,
u · n = qN on ∂DN ,

(1c)
(1d)

with ∂D = ∂DD ∪ ∂DN and ∂DD ∩ ∂DN = ∅. In equations (1), u denotes the macroscopic
fluid velocity whereas K and p stand for absolute permeability and pressure.
Let us assume that D contains several fractures, that all together constitute a single
domain Γ of spatial dimension N such that Γ ⊂ D, which is a possibly unconnected, open
subset of D. The surrounding porous rock, namely, the remaining part of D, is called
Ω = D \ Γ. Assuming that the fracture aperture ε at each point of Γ is small compared
to other characteristic dimensions of the fractures, the full-dimensional domain Γ can be
reduced to the (N −1)-dimensional fracture network γ. This reduction involves modeling
choices resulting in different hybrid-dimensional problem formulations that form the basis
for the methods presented in the following section.

3 Participating discretization methods
Within this section, the discretization methods participating in this benchmark study are
described. The purpose of this article is the comparison of well-known, established and/or
at least published methods. Therefore, only the most significant aspects of each method are
summarized. We do not show a comparison against analytical solutions here. The analysis
of the methods and theoretical results such as proofs of optimal convergence can be found
in the corresponding references. A summary of all participating methods is provided in
Table 1. In the sequel, we will denote with d.o.f. the degrees of freedom associated to a
specific method. We indicate also the type of conformity required to the computational grid
with respect to the fractures and the assumption that the pressure is considered continuous
across the fractures.
3.1 Vertex-centered, continuous-pressure, conforming lower-dimensional DFM
(Box-DFM)
The lower-dimensional representation of fractures allows easier mesh generation in comparison to the equi-dimensional approach, as it circumvents the appearance of very small

3

method
d.o.f.
frac-dim conforming p-cont.
Box-DFM
p (vert)
dim-1
yes
yes
CC-DFM
p (elem)
dim-1
yes
no
EDFM
p (elem)
dim-1
no
yes
mortar-DFM p (elem), u (faces)
dim-1
geometrically
no
P-XFEM
p (vert)
dim-1
no
no
D-XFEM
p (elem), u (faces)
dim-1
no
no
MFD
p (faces)
dim
geometrically
no
Table 1: Participating discretization methods.
elements when discretizing the interior of the fracture (i.e. , within the fracture width).
The conforming mesh generation algorithm honors the geometrical characteristics of the
fracture system. Conform meshing implies that the fractures are discretized with a set of
line elements (in a 2D domain) that are also the edges of the triangular finite elements.
The spatial discretization in Box-DFM is performed with the Box method, a vertexcentered finite-volume method proposed in, e.g. Helmig (1997) which combines the advantages of finite element and finite volume grids, allowing unstructured grids and guaranteeing a locally conservative scheme (Reichenberger et al. (2006)). Figure 1 illustrates a
two-dimensional representation of the dual-grid with two finite elements E1 and E5 sharing
the same edge (ij1 ) that represents a lower-dimensional fracture with the aperture εij . The
main characteristic in terms of the fractured system is that the pressure is required to be
continuous, in particular in those vertices whose control volumes overlap both fracture and
matrix regions.
The Box-DFM method used for this paper is implemented in the open-source numerical
simulator DuMux . A detailed description of the conceptual, mathematical and numerical
model and code implementation is published in Tatomir (2012). The Box-DFM simulation
code used for the benchmark studies is publicly available under https://git.iws.unistuttgart.de/dumux-pub/Flemisch2016a.git.
3.2 Cell-centered, discontinuous-pressure, conforming DFM (CC-DFM)
The control volume finite difference method uses a two-point flux approximation (TPFA)
based on the cell-center pressure values for the evaluation of the face fluxes, and is a
widely applied and standard method for simulation of flow in porous media. The domain
is partitioned with fractures coinciding with the interior faces between matrix cells just as
described in Section 3.1. The flux over the face between matrix cells i and j is approximated
by
uij = Tij (pi − pj ),
(2)
where pi and pj are the pressures in the neighboring cells and Tij is the face transmissibility,
computed as the harmonic average of the two half transmissibilities corresponding to the
face and the two cells. The half transmissibility of cell-face pair i is in turn given as
Ai n>
i Ki
· di ,
αi =
>
di di

(3)

where Ai and ni are the area and unit normal vector of the face, Ki is the permeability
assigned to the cell and di is the distance vector from cell center to face centroid.

4

Figure 1: Conceptual representation of the Box-DFM method: (left-hand side) The dual
finite element and finite volume mesh from which the control volume Bi around
node i is created. Node i is surrounded by nodes {j1 , j2 , j3 , j4 , j5 }, where segment
ij1 represents both a fracture and a shared FE edge; (right-hand side) Sub-control
1
volume (SCV) bE
in element E1 has barycenter G1 and the mid-points of the
i
E1
edges ij1 and ij2 are Mij1 , respectively Mij2 . The SCV face fij1
is the segment
E1
G1 Mij1 which contains the integration point xij1 where the normal vector nE1
ij1 is
applied.
In addition to the unknowns given at the centroids of the matrix cells, unknowns are
associated to the centroids of the fracture cells. The fracture cells are associated with
apertures, which multiplied with the length give the volume of these cells. The aperture is
also used to construct hybrid faces for the matrix-fracture interfaces. These faces, parallel
to the fracture but displaced half an aperture to either side, enable us to compute the half
transmissibilities between the fracture cell and the matrix cells on the two sides. These
faces are indicated by the dashed blue lines in Figure 2, where the computational domain
is superimposed on the geometrical grid. The result is a hybrid grid with fractures which
are lower dimensional in the grid, but equidimensional in the computational domain at the
cost of a small matrix volume error corresponding to the overlap of the matrix cells with
the fracture cells.
Following the method proposed by Karimi-Fard et al. (2004), the intermediate fracture
intersection cell drawn with dashed red lines in Figure 2 is removed, leading to direct
coupling of the fracture cells neighbor to the intersection. The purpose of this is both to
obtain a smaller condition number and to avoid severe time-step restrictions associated
with small cells in transport simulations. To each new face between cell i and j, face
transmissibilities are assigned, calculated using the star delta transformation as described
in Karimi-Fard et al. (2004):
αi αj
Tij = P
,
(4)
n
αk
k=1

with n denoting the number of fracture cells meeting at the intersection. As this elimination
disregards all information on the permeability of the intersection, it should be used with
caution in cases of crossing fractures of different permeability. We encounter this feature
in 4.3, and include results both with and without the elimination for one of the test cases
presented in that section.
Inspired by the CC-DFM method by Karimi-Fard et al. (2004) presented above, a method

5

(a)

(b)

Figure 2: (a) Conceptual decomposition of the domain according to element dimension with
the matrix depicted in black, fractures in blue and their intersections in red. (b)
The computational domain of the CC-DFM. Dashed lines are faces of the fracture
cells.
based on the multi-point flux approximation has also been developed Sandve et al. (2012).
The MPFA variant of the method reduces errors associated with the TPFA approach for
grids that are not close to K-orthogonal, and avoids errors related to the splitting of the
fluxes in the star-delta transformation. We refer to Sandve et al. (2012) for a thorough
comparison of the TPFA and MPFA CC-DFM approaches. The implementation of both
methods is available in the open-source Matlab Reservoir Simulation Toolbox http://www.
sintef.no/projectweb/mrst/.
3.3 Continuous-pressure, non-conforming embedded DFM (EDFM)
Recently, non-conforming methods for the treatment of lower-dimensional fractures have
been developed, for example in Moinfar et al. (2014a, 2011); Hajibeygi et al. (2011), to
avoid the time-consuming construction of complex matrix grids which explicitly represent
the fractures. They are mostly used in the context of single and multi-phase flow simulations for petroleum engineering applications and require the normal fracture permeability
to be orders of magnitude higher than the matrix permeability, as in the case of enhanced
reservoir exploitation and fractures stimulation. In this field of applications corner-point
grids are normally employed to describe the geological layers, e.g. different rock type, of the
reservoir. An adaptation of such computational grids to the fractures could be unaffordable
for real cases. The numerical method belongs to the family of two-point schemes, where
a one-to-one connection between the degrees of freedom is considered through the transmissibility concept (Eymard et al. (2000)). References on the embedded discrete fracture
method (EDFM) can be found, for example, in Li and Lee (2008); Panfili et al. (2013);
Moinfar et al. (2014b); Panfili and Cominelli (2014); de Araujo Cavalcante Filho et al.
(2015); Fumagalli et al. (2016).
In practice, the meshes of the fractures are generated on top of the rock grid so that
each rock cell cut by fractures contains exactly one fracture cell per fracture. Intersections
between fractures are computed without affecting the creation of the grids of fractures
and rock and used to compute approximate transmissibilities between different fracture

6

cells. See Figure 3 as an example. A degree of freedom that represents a pressure or

Figure 3: Example of meshes, for both fractures and rock matrix, suited for EDFM. The
rock matrix is considered as a background mesh. Each fracture cell is represented
by two blue dots and the green dots are the non-matching intersection among
fractures.
a saturation value is assigned to each matrix cell and to each fracture cell. This means
that transmissibilities between matrix and fracture cells, as well as those between different
fracture cells, need to be computed. We compute the transmissibility between a fracture
cell and a matrix cell Tf m and the half-transmissibility Ti between two intersecting fracture
cells (related to the fracture i) through the following approximate expressions
Tf m = A

n>
f K · nf
df,m

and Ti = s

ki εi
.
di,s

Here A is the measure of the fracture cell in the current rock cell, nf is the normal of the
fracture cell and df,m is an average distance between the fracture cell and the matrix cell,
see Li and Lee (2008). For the fracture-fracture transmissibility, s indicates the measure of
the intersecting segment, ki the scalar permeability of the fracture, εi the aperture and di,s
is the average distance between the fracture cell and the intersecting segment. The standard
harmonic average is considered to compute the transmissibility between the two fracture
cells. Standard formulae for fracture-fracture as well as matrix-matrix transmissibilities
are computed by means of a two-point flux approximation. It is worth to notice that the
recent extension of EDFM called Projection-based EDFM (pEDFM), proposed in Tene
et al. (2016), is also able to handle low permeable fractures. Finally, even if the proposed
benchmark cases are two-dimensional the method can be extended to three dimensions
without any additional constraints.
3.4 Cell-centered, discontinuous-pressure, geometrically-conforming mortar DFM
(mortar-DFM)
The key concept behind the mortar-DFM, as described more thoroughly in Boon et al.
(2016), is the idea that fractures can be considered as interfaces between different subdomains. This has been explored previously by Martin et al. (2005); Frih et al. (2012),
among others. In this context, it is interesting to consider domain decomposition techniques
such as the mortar method to model flow through the fractured porous medium.
The mortar method is generally used to couple equations in different sub-domains by
introducing a so-called mortar variable, defined on the interface. In case of modeling

7

fracture flow, a well-explored choice of the mortar variable is the fracture pressure (Martin
et al., 2005). The method considered here, however, uses as the mortar variable the flux,
which leads to a stronger sense of mass conservation for flows between the matrix and
fractures. One of the main advantages of the close relationship to mortar methods is
the capability to handle non-matching grids. In particular, two sub-domains bordering
a fracture can be meshed independently on both sides, as illustrated in Figure 4. The

Figure 4: The mortar-DFM allows for non-matching grids along fracture interfaces. Fracture and matrix flows are coupled using a mortar variable, defined on a coarser
grid (green dots).
difficulty in mesh generation is then relieved significantly since only the geometry of the
fractures needs to be respected.
By construction, the mortar-DFM is applicable to problems in arbitrary dimensions. The
governing equations in the matrix and the fractures (as well as fracture intersections in 3D)
are identical and thus all fractures, intersections and tips are handled in a unified manner.
Consequently, although only two-dimensional problems are considered in this case study,
the discretization scheme is not at all limited to the presented benchmark problems and
3D cases can easily be considered.
With the use of mixed finite elements, mass is conserved locally in the matrix, fractures,
and fracture intersections. The flux u and pressure p are modeled using the lowest-order
Raviart-Thomas elements and piecewise constants, i.e. RT0 − P0 . Additionally, the mortar variable is given by piecewise constants on a separately generated, lower-dimensional,
mortar grid. This grid matches with the surrounding grids in case of matching grids and
is coarser otherwise (Boon et al., 2016). The resulting mixed finite element formulation
is a saddle-point problem, which may be challenging to solve numerically. To relieve this,
the flux variables may be eliminated through hybridization, which leads to a less expensive
scheme containing solely the cell-center pressures.
Two implementations of the method have been developed, both of which are used in
this benchmark study. The first version, implemented in MATLAB, has the capability
of handling non-matching grids along fractures for two-dimensional problems. It is wellsuited for simpler geometries, containing relatively few fractures, such as those considered
in Benchmarks 1-3. The second version has been implemented for 3D problems and higherorder spaces on matching grids using the open-source finite element library FEniCS (Logg
et al., 2012). This version is more efficient for complex cases such as Benchmark 4.

8

3.5 Discontinuous-pressure, non-conforming primal XFEM (P-XFEM)
The primal XFEM method participating in this benchmarking study is described in detail
in Schwenck (2015), see also Flemisch et al. (2016); Schwenck et al. (2015). The method is
based on the hybrid-dimensional problem formulation investigated in Martin et al. (2005),
where conditions for the coupling between fracture and matrix are derived:
{{um · n}}γ = kf,n /ε Jpm Kγ


ξ0 Jum · nKγ = kf,n /ε {{pm }}γ − pf

(5a)
(5b)

Here, the subscripts “m” and “f” indicate matrix and fracture quantities, while {{·}}γ and
J·Kγ denote the average and the jump of a matrix quantity over the fracture γ, respectively.
The coupling conditions (5) can be used to define a source term for the fracture flow
problem, while they yield an interface problem for the matrix domain. For the discretization
of this interface problem, the methodology presented in Hansbo and Hansbo (2002) is used,
which amounts to applying the eXtended Finite Element Method (XFEM). Together with
an independent standard discretization of the lower-dimensional fracture problem, this
yields a hybrid-dimensional, non-conforming primal XFEM-based method. The XFEM
space is built enriching the standard Lagrangian P1 (or Q1 for quads) finite-element spaces,
whose degrees of freedom are located at the vertices of the full-dimensional grid of the
matrix Ω and the lower-dimensional grid of the fracture γ. A representative example of

Figure 5: Example of meshes, for both fractures and rock matrix, suited for P-XFEM. The
fracture grid vertices are indicated by the blue dots.
matrix and fracture grids is illustrated in Figure 5. Unlike the EDFM method, see Figure 3,
the fracture grid vertices can be placed arbitrarily without taking into account the matrix
grid. On the other hand, the method requires matching fracture branch grids in the form
of vertices placed at the fracture intersections. In particular, special care has to be taken
of intersecting and immersed fractures (Schwenck et al., 2015).
The method is implemented on top of the DUNE framework Bastian et al. (2008) and
the discretization module DUNE-PDELab Bastian et al. (2010). For the enrichment of
the finite-element spaces in the context of XFEM, the modules DUNE-Multidomain and
DUNE-Multidomaingrid are employed Müthing (2015). The simulation code for the XFEM
approach and for the benchmarks studied here is publicly available under https://git.
iws.uni-stuttgart.de/dumux-pub/Flemisch2016a.git.

9

3.6 Discontinuous-pressure, non-conforming dual XFEM (D-XFEM)
The dual XFEM method participating in his benchmark is based on D’Angelo and Scotti
(2012). The method, originally derived for a domain cut by one fracture, was further
developed in Formaggia et al. (2014), Fumagalli and Scotti (2014) to account for intersecting
fractures with different permeabilities. The same equations and coupling conditions as
for the primal XFEM are used, but in a dual formulation where Darcy law and mass
conservation give rise to a saddle-point problem for the fluid mean velocity and pressure,
both in the fracture and in the surrounding medium. Moreover, unlike the previous method,
this method employs triangular/tetrahedral grids. The usual lowest order RT0 −P0 pair for
velocity and pressure is enriched following Hansbo and Hansbo (2002) in the elements of the
porous medium cut by a fracture, or in the elements of a fracture at the intersection with
other fractures. Indeed, triangular/tetrahedral grids are arbitrarily cut by triangulated
lines/surfaces in 2D and 3D respectively. These surfaces can, in turn, intersect each other
in a non-conforming way, as shown in Figure 6.

Figure 6: A portion of the grid cut by two fractures: in the two dimensional case they can
split the elements in two (grey), three (yellow), or four (red) independent parts,
where the restrictions of the basis functions are defined. The fracture grids are
irrespective of the bulk grid and of each other, i.e. the intersection point ip is not
a point of the grid.
In the current implementation of the method no special enrichment is added in the bulk
elements containing the fracture tips. Instead, fractures are artificially extended up to the
boundary of the domain, and in the extension we prescribe the same permeability of the
surrounding porous medium to obtain a “virtual” fracture with no effects on the flow.
The method has been implemented on the basis of the Getfem++ library, http://
download.gna.org/getfem/html/homepage/, which provides support for the computation
of the intersections and the quadrature on sub-elements thanks to an interface with QHull,
http://www.qhull.org/.
3.7 Reference Solutions calculated with mimetic finite differences (MFD)
The reference solutions are computed on very fine grids that discretize both matrix and
fractures by full-dimensional triangular or quadrilateral elements. A mimetic finite differ-

10

ence method, see Brezzi et al. (2005); Flemisch and Helmig (2008), is used to discretize
problem (1). The method is employed as it is implemented in DuMux 2.7 Flemisch et al.
(2011). In particular, a mixed-hybrid approach is used to transform the discrete saddle
point problem in terms of cell pressures and face fluxes into a symmetric positive definite
formulation with face-pressure degrees of freedom.

4 Benchmark Problems
This is the main section which compares the methods described above by means of four
benchmark cases. First, in Section 4.1, we present a well established benchmark for groundwater flow from Swedish Nuclear Power Inspectorate (SKI) (1987) that contains two crossing, highly permeable fractures and a non-straight top surface. The second benchmark case,
considered in Section 4.2, is based on Geiger et al. (2011) and shows a regular fracture network. After that, a small but complex fracture network exhibiting immersed fractures and
intersections at different angles is investigated in Section 4.3. Finally, a case synthesized
from a real application is considered in Section 4.4.
For each benchmark case, a description of the computational domain is provided, including boundary conditions, the geometrical information about the corresponding fracture
network and the associated material parameters such as aperture and permeability. For
some of the cases, the reference solution on the complete domain is visualized. This is
followed by illustrations of the grids used by the participating methods. Since the methods
pose different requirements there, the grid could be chosen arbitrarily for each method,
provided that the number of grid cells or vertices is roughly the same. If a reference solution is available (Benchmarks 1–3), the results of the different methods are compared
by evaluating the errors with respect to the reference in the matrix domain as well as in
the fracture network, indicated by errm and errf , respectively. The errors are calculated
according to the formulas
2
errm
=

errf2 =

1
|Ω|(∆pref )2
1
|γ|(∆pref )2

X

|f | (pm |Km − pref |Kref )2 ,

f =Kref ∩Km

X

|e| (pf |Kf − pref |Kref )2 ,

e=Kref ∩Kf

where |Ω| and |γ| indicate the size of the matrix and fracture domain, respectively, and
∆pref = maxD pref − minD pref . The sum is taken over all intersections of (full-dimensional)
elements Kref of the grid employed for the reference solution with full-dimensional matrix
elements Km in case of errm and lower-dimensional fracture elements Kf in case of errf .
Moreover, the densities and condition numbers of the resulting linear system matrices are
provided. In addition to that, a comparison is performed by means of plots along specific
lines through the domain for some benchmark cases. Each case is concluded by a short
discussion of the results.

11

4.1 Benchmark 1: Hydrocoin
Within the international Hydrocoin project, (Swedish Nuclear Power Inspectorate (SKI)
(1987)), a benchmark for heterogeneous groundwater flow problems was presented. The
domain setup is shown in Figure 7. We point out that we have slightly modified the original
1

50m

5

20 ,30 ,40

9

60 ,70 ,80

-200m
1100m
16
19

17
18

15

10
14 13

12 11

1600m
Figure 7: Geometry of the modeled domain of the Hydrocoin test case 2, Swedish Nuclear
Power Inspectorate (SKI) (1987). Modified node locations are indicated by numbers superscripted with 0 . Boundary conditions are hydraulic head on top and
Neumann no-flow on the other three sides of the domain.
domain such that equi-dimensional and hybrid-dimensional models can be run on exactly
the same domain. This allows for an easier comparison of the solution values over the whole
domain. The exact modifications are described in B.
For this case, we keep the original formulation in terms of the piezometric head and the
hydraulic conductivity instead of pressure and permeability. In particular, the boundary
conditions are Dirichlet piezometric head on the top boundary and Neumann no flow on
the other three boundaries. The hydraulic conductivity is 10−6 m/s in the fracture zones
and 10−8 m/s in the rock matrix respectively.
Table 2 lists the number of degrees of freedom, matrix elements and fracture elements
for all the participating methods. The corresponding grids are visualized in Figure 8.
The original benchmark shows the piezometric head distribution along five horizontal
lines through the modeled domain. Here, we first show in Figure 9 the plot at a depth
of 200 m, as indicated by the dashed line in Figure 7. All participating methods show a
good agreement with the reference solution. Only the EDFM method is a bit off. We
remark that the plots for the methods employing cell-wise constant solution values exhibit
staircase-like patterns corresponding to these values.
Table 3 lists the discretization errors for the different methods, particularly, the error for
the matrix domain and the one along the two fractures. Moreover, it provides the density
of the associated matrix and its condition number for each method. The uniform behavior

12

method
d.o.f. matrix elements fracture elements
Box-DFM
1496
2863 triangles
74
CC-DFM
1459
1416 triangles
43
EDFM
1044
960 quads
84
mortar-DFM 3647
1384 triangles
63
P-XFEM
1667
1320 quads
68
D-XFEM
3514
1132 triangles
160
MFD
889233
424921 mixed
19287
Table 2: Grids for Benchmark 1.

(a) Box-DFM

(b) CC-DFM

(c) EDFM

(d) mortar-DFM

(e) P-XFEM

(f) D-XFEM

Figure 8: Benchmark 1: the grids used by the different methods.
exhibited in Figure 9 is reflected by the error values. Especially the errors in the matrix
domain are within very narrow bounds, while the fracture errors show a larger variation.
The densities of the matrices are also close together. Remarkably high differences can be
observed for the matrix condition numbers. While the ones for Box-DFM, CC-DFM and
EDFM are on the order of 104 , the one for P-XFEM is five orders and the ones for mortarDFM and D-XFEM are even seven orders of magnitude larger, due to their saddle-point
nature.

13

120

reference
Box−DFM
CC−DFM
EDFM
mortar−DFM
P−XFEM
D−XFEM

piezometric head [m]

118
116
114
112
110
108
106
0

0.2

0.4

0.6
0.8
1
arc length [km]

1.2

1.4

1.6

Figure 9: Benchmark 1: head values along a horizontal line at a depth of 200 m.

method
Box-DFM
CC-DFM
EDFM
mortar-DFM
P-XFEM
D-XFEM

errm
9.3e-3
1.1e-2
1.5e-2
1.0e-2
1.2e-2
1.2e-2

errf nnz/size2
3.3e-3
4.5e-3
1.1e-2
2.7e-3
8.3e-3
4.7e-3
7.2e-3
1.5e-3
3.2e-3
6.5e-3
6.9e-3
1.7e-3

k · k2 -cond
5.4e3
3.5e4
3.9e4
9.0e12
2.7e9
6.2e12

Table 3: Discretization errors and matrix characteristics for Benchmark 1.

14

4.2 Benchmark 2: Regular Fracture Network
This test case is based on an article presenting a new dual continuum model, Geiger et al.
(2011), with slightly modified boundary conditions and material properties. The computational domain including the fracture network and boundary conditions is shown in
Figure 10. The matrix permeability is set to K m = I, all fractures have a uniform aperture

y

q̄ N = 0
1
(0.75,0.75)
(0.625,0.625)
p̄D = 1

q̄ N = −1
(0.5,0.5)

0
0

q̄ N = 0

1

x

Figure 10: Benchmark 2: Domain and boundary conditions.
ε = 10−4 . For the fracture permeability we consider two cases: a highly conductive network
with kf,n = kf,t = 104 , as worked out in Subsection 4.2.1, and a case with blocking fractures
by setting kf,n = kf,t = 10−4 , as described in Subsection 4.2.2. The reference solutions are
computed on a grid which resolves every fracture with 10 elements in its normal direction
and becomes coarser away from the fractures. It has a total of 1175056 elements.
The first distinction between the different schemes are given in Table 4, where the number of degrees of freedom, matrix elements and fracture elements for all the participating
methods are listed. The corresponding grids are visualized in Figure 11.
4.2.1 Conductive Fracture Network
First, we consider a highly conductive network by setting kf,n = kf,t = 104 . The pressure
distribution of the corresponding reference solution is shown in Figure 12.
The pressure distributions given by the different methods are first compared along two
lines, one horizontal at y = 0.7 and one vertical at x = 0.5. As shown in Figure 13, all
results are relatively close to the reference solution. Qualitatively, we observe that P-XFEM
produces a more diffuse pressure profile in the vertical fracture.
The results for the conducting fractures are similar to those presented in the first benchmark. In particular, the performance of the methods is comparable as shown by both
the matrix and the fracture errors. In fact, since the degree of sparsity does not differ

15

method
d.o.f. matrix elements fracture elements
Box-DFM
1422
2691 triangles
130
CC-DFM
1481
1386 triangles
95
EDFM
1501
1369 quads
132
mortar-DFM
3366
1280 triangles
75
P-XFEM
1632
961 quads
318
D-XFEM
4474
1250 triangles
126
MFD
2352280
1136456 quads
38600
Table 4: Grids for Benchmark 2.
method
Box-DFM
CC-DFM
EDFM
mortar-DFM
P-XFEM
D-XFEM

errm
6.7e-3
1.1e-2
6.5e-3
1.0e-2
1.7e-2
9.6e-3

errf nnz/size2
1.1e-3
4.7e-3
5.0e-3
2.7e-3
4.0e-3
3.3e-3
7.4e-3
1.8e-3
6.0e-3
7.8e-3
8.9e-3
1.3e-3

k · k2 -cond
7.9e3
5.6e4
5.6e4
2.4e6
6.8e9
1.2e6

Table 5: Discretization errors and matrix characteristics for Benchmark 2 with conductive
fractures.
significantly either, the only notable differences between the methods are the number of
degrees of freedom and the condition numbers, as shown in Table 5. In that context, the
mortar-DFM and D-XFEM are the clear outliers, containing a large number of degrees
of freedom due to the incorporated flux variable and resulting in high condition numbers.
The P-XFEM scheme exhibits the highest condition number, yet we emphasize that it has
significantly fewer degrees of freedom.
4.2.2 Blocking Fracture Network
We now assume a blocking fracture network by setting kf,n = kf,t = 10−4 . The pressure
distribution of the corresponding reference solution is shown in Figure 14. The results
clearly show the pressure discontinuities reminiscent of the low fracture permeability.
Figure 15 compares the results of the different methods along a diagonal line crossing
the whole domain from (0.0, 0.1) to (0.9, 1.0). The discretization errors, sparsity densities,
and condition numbers for the different methods are given in Table 6.
In the case of blocking fractures, the distinction between the different methods is more
apparent. As mentioned above, the Box-DFM and EDFM schemes are unable to capture the
resulting pressure discontinuities. As a result, these methods show large errors in both the
matrix and the fracture domains. The remaining methods, which are capable of handling
discontinuities, do not differ significantly among each other in terms of fracture and matrix
errors. We do note that the condition numbers have improved significantly for the mortarDFM and P-XFEM schemes. Conversely, for CC-DFM and D-XFEM, condition numbers
for the blocking fractures case are similar to those obtained for the permeable fractures
case.

16

(a) Box-DFM

(b) CC-DFM

(c) EDFM

(d) mortar-DFM

(e) P-XFEM

(f) D-XFEM

Figure 11: Benchmark 2: the grids used by the different methods. In the DXFEM grid the
red lines indicate the virtual extension of the fractures up to the boundary.

Figure 12: Benchmark 2 with conductive fractures: pressure reference solution.

17

1.2

1.5

reference
Box−DFM
CC−DFM
EDFM
mortar−DFM
P−XFEM
D−XFEM

1.45
1.4

1.19
1.18

1.3

pressure

pressure

1.35

reference
Box−DFM
CC−DFM
EDFM
mortar−DFM
P−XFEM
D−XFEM

1.25
1.2
1.15

1.17
1.16
1.15

1.1
1.14
1.05
1
0

0.2

0.4
0.6
arc length

0.8

1

0

(a) Horizontal line at y = 0.7.

0.2

0.4
0.6
arc length

0.8

1

(b) Longest vertical fracture at x = 0.5.

Figure 13: Benchmark 2 with conductive fractures: comparison of values along two lines.

Figure 14: Benchmark 2 with blocking fractures: pressure reference solution.

method
Box-DFM
CC-DFM
EDFM
mortar-DFM
P-XFEM
D-XFEM

errm
4.1e-1
5.7e-3
2.9e-1
4.5e-3
2.9e-3
1.0e-2

errf nnz/size2
3.8e-1
4.7e-3
4.6e-3
2.7e-3
3.2e-1
3.3e-3
4.9e-3
1.6e-3
2.2e-2
8.1e-3
1.9e-2
1.3e-3

k · k2 -cond
3.5e3
2.6e4
9.2e3
9.0e2
2.0e4
2.2e6

Table 6: Discretization errors and matrix characteristics for Benchmark 2 with blocking
fractures.

18

3.5

reference
Box−DFM
CC−DFM
EDFM
mortar−DFM
P−XFEM
D−XFEM

pressure

3

2.5

2

1.5

1
0

0.2

0.4

0.6
0.8
arc length

1

1.2

Figure 15: Benchmark 2 with blocking fractures: values along the line (0.0, 0.1) − (0.9, 1.0).

19

4.3 Benchmark 3: Complex Fracture Network
This test case considers a small but complex fracture network that includes permeable
and blocking fractures. The domain and boundary conditions are shown in Figure 16.
The exact coordinates for the fracture positions are provided in C. The fracture network

3
2

p̄D = 4

5

9

8

4

q̄ N = 0

4

10

9

8

q̄ N = 0

10

q̄ N = 0

5
3

7

1

2

6

p̄D = 1

p̄D = 4

7

1
6

p̄D = 1

q̄ N = 0

(a)

(b)

Figure 16: Benchmark 3: Domain and boundary conditions for cases (a) and (b). The red
fractures are conductive, the blue ones are blocking.
contains ten straight immersed fractures, grouped in disconnected networks. The aperture
is ε = 10−4 for all fractures, and permeability is kf,n = kf,t = 104 for all fractures except
for fractures 4 and 5 which are blocking fractures with kf,n = kf,t = 10−4 . Note that we
are considering two subcases a) and b) with a pressure gradient which is predominantly
vertical and horizontal respectively, to better highlight the impact of the blocking fractures.
The corresponding reference solutions are depicted in Figure 17.

(a)

(b)

Figure 17: Benchmark 3: reference solution for cases a) and b)
Table 7 lists the number of degrees of freedom, matrix elements and fracture elements
for all the participating methods. The corresponding grids are visualized in Figure 18.

20

method
d.o.f. matrix elements fracture elements
Box-DFM
1460
2838 triangles
155
CC-DFM
1510
1407 triangles
103
EDFM
1572
1369 quads
203
mortar-DFM
3953
1452 triangles
105
D-XFEM
7180
1922 triangles
199
MFD
1800770
1192504 mixed
7876
Table 7: Grids for Benchmark 3.
The P-XFEM method could not participate in this benchmark example. Its current
implementation requires that each matrix element face is cut by at most one fracture
branch. While it would be possible to construct a matrix grid that satisfies this requirement,
this would contradict the promised advantage of admitting independent fracture and matrix
grids.
4.3.1 Flow from Top to Bottom
Table 8 lists the discretization errors for the first variant, namely, the flow from top to
bottom.
method
Box-DFM
CC-DFM
EDFM
mortar-DFM
D-XFEM

errm
4.4e-2
2.6e-2
3.8e-2
1.0e-2
1.9e-2

errf nnz/size2
3.8e-2
4.6e-3
3.3e-2
2.7e-3
4.5e-2
3.1e-3
1.7e-2
1.4e-3
2.9e-2
8.2e-4

k · k2 -cond
4.5e3
3.8e4
1.2e6
1.1e6
8.1e3

Table 8: Discretization errors and matrix characteristics for the first variant of Benchmark
3.
Even though this is still a synthetic case, we can see that the geometry of the network
starts to be an issue: relatively small intersection angles are present, for instance, between
fractures 1 and 2. Another difficulty consists in the coexistence of permeable and blocking
fractures which intersect each other: on one hand, some of the methods are not well suited to
describe a blocking behavior, on the other hand the coupling conditions at the intersection
become less trivial in these cases. All the participating methods that account explicitly for
the effect of permeability at the fracture intersections have adopted the harmonic average
in the case of a permeable and a blocking fracture crossing each other. The errors reported
in Table 8 show that the methods requiring the continuity of pressure (EDFM and the
Box-DFM) exhibit slightly higher errors in the matrix. However, the difference is not
particularly sharp, since in this sub-case the average pressure gradient is almost parallel to
the blocking fractures.

21

(a) Box-DFM

(b) CC-DFM

(c) EDFM

(d) mortar-DFM

(e) P-XFEM: n/a

(f) D-XFEM

Figure 18: Benchmark 3: the grids used by the different methods. In the DXFEM grid the
red lines indicate the virtual extension of the fractures up to the boundary.
4.3.2 Flow from Left to Right
The discretization errors for the second variant, namely, the flow from left to right, are
summarized in Table 9.
In this second case, since we impose pressure on the sides of the square domain, the solution is more challenging as we can observe from Figure 17 and the gap between continuous
and discontinuous methods increases. However, it should be noted that the errors remain
of the same order of magnitude, indicating that all the methods capture the overall trend of
the solution. The elimination of the fracture intersection cells in the CC-DFM is ill-suited
for cases where fractures of different permeability cross. Therefore, we include a solution
CC-DFM* in which we have not performed the removal for case b). The new results, reported in Table 9, show a far smaller error compared to the CC-DFM with elimination,
but also demonstrate that the elimination significantly reduces the condition number.

22

method
Box-DFM
CC-DFM
CC-DFM*
EDFM
mortar-DFM
D-XFEM

errm
7.5e-2
5.2e-2
1.1e-2
5.8e-2
1.3e-2
2.2e-2

errf nnz/size2
7.0e-2
4.6e-3
7.3e-2
2.7e-3
2.7e-2
2.6e-3
8.9e-2
3.1e-3
2.7e-2
1.4e-3
3.6e-2
8.2e-4

k · k2 -cond
5.6e3
4.5e4
8.1e5
1.2e6
7.3e8
8.1e3

Table 9: Discretization errors and matrix characteristics for the second variant of Benchmark 3.

23

Figure 19: In the left the interpretation of the set of fractures superimposed to the map.
In the right the geometry used in the simulations. The rectified fractures are
depicted in blue.
4.4 Benchmark 4: a Realistic Case
In this last test case we consider a real set of fractures from an interpreted outcrop in
the Sotra island, near Bergen in Norway. The set is composed of 64 fractures grouped
in 13 different connected networks, ranging from isolated fractures up to tens of fractures
each. In the interpretation process two fractures were composed by more than one segment.
However, since the implementation of some methods relay on the fact that one fracture is
represented by a single geometrical object, we substitute them by a single segment. It
is worth to notice that we are changing the connectivity of the system, nevertheless our
goal is to make a comparison of the previous schemes on a complex set of fractures. The
interpreted outcrop and the corresponding set of fractures are represented in Figure 19. The
size of the domain is 700 m × 600 m with uniform scalar permeability equal to 10−14 m2 .
For simplicity all the fractures have the same scalar permeability equal in the tangential
and normal direction to 10−8 m2 , and aperture 10−2 m. We consider no-flow boundary
condition on top and bottom, pressure 1013250 Pa on the left, and pressure 0 Pa on the
right of the boundary of the domain. Due to the high geometrical complexity of the fracture
network not all involved numerical schemes/simulators could be used. Nevertheless, it is
worth to point out that for the others the main difficulty in handling such geometry is an
implementation issue rather than a limitation of the scheme.
Table 10 lists the number of degrees of freedom, the density of the associated matrix, and
its condition number for the different methods. Due to the geometrical difficulties of the
method
d.o.f. matrix elem
Box-DFM
5563 10807 triangles
CC-DFM
8481 7614 triangles
EDFM
3599
2491 quads
mortar-DFM 25258 8319 triangles

frac elem
1386
867
1108
1317

nnz/size2
1.2e-3
4.9e-4
1.4e-3
2.0e-4

k · k2 -cond
9.3e5
5.3e6
4.7e6
2.2e17

Table 10: Discretization and matrix characteristics for Benchmark 4.
network the request of having a similar number of degrees of freedom among the methods

24

is relaxed, as Table 10 indicates. Considering Figure 20, the solutions are reported for the
four methods. We notice that, except for the top right part of the domain in the Box-DFM
method, the solutions are similar and comparable which is an indication of their correctness.
Compared to the previous test cases the mesh generation is the main concern and some of
the methods require a fine tuning to avoid unphysical connections among elements where
the fracture are close. An example can be found in the middle of the domain and reported
in Figure 21. Only EDFM is more robust with respect to this constraint. To present a
more detailed comparison among the methods, Figure 22 represents the pressure solution
along two different lines: for y = 500 m and for x = 625 m. We note that the methods
behave similarly, and the Box-DFM slightly overestimates some peaks. The oscillation of
the methods are related to mesh effects.

25

(a) Box-DFM

(b) CC-DFM

(c) EDFM

(d) mortar-DFM

Figure 20: Representation of the matrix pressures field for the realistic case. The solution
values range between 0 and 101325 Pa.

26

(a) Box-DFM

(b) CC-DFM

(c) EDFM

(d) mortar-DFM

Figure 21: Benchmark 4: Representation of mesh in the middle of the domain. The size of
the picture is approximately 30 m×15 m centered in (360, 350). It is represented
by the small rectangle in the centre of Figure 19 left.

5

12

5

x 10

9

x 10

8
10
7
6
pressure [Pa]

pressure [Pa]

8

6

4

2

0
0

5
4
3
2

Box−DFM
CC−DFM
EDFM
mortar−DFM
100

200

300
400
arc length [m]

500

600

700

(a) y = 500 m

1

Box−DFM
CC−DFM
EDFM
mortar−DFM

0
0

100

200

300
400
arc length [m]

500

600

(b) x = 625 m

Figure 22: Benchmark 4: Pressure solutions of the 4 methods plotted over lines (a) y =
500 m, and (b) x = 625 m.

27

5 Summary and Outlook
Four benchmark cases for single-phase flow in fractured porous media have been proposed
and employed to compare the performances of several state-of-the-art hybrid-dimensional
discrete-fracture-matrix models. If we consider the cases where all the methods are employed within the applicability range for which they were originally developed, the results
are in quite good agreement. In particular, fracture networks exhibiting a larger permeability than the surrounding matrix can be accurately described by all methods. On the
other hand, not all methods are capable of modeling blocking fractures. In this case some
methods fail to predict the correct flow patterns for the corresponding cases. Especially
noteworthy are the large differences in the condition numbers of the associated system matrices. The effect of these differences on the behavior of linear solvers is difficult to quantify
in a comparable manner, since the different methods pose different requirements for such
solvers. In principle, all participating methods should have been able to run all proposed
cases. However, due to implementation restrictions, some methods could not perform the
cases with more complex fracture network geometries. Even if the methodology is general
enough, technical difficulties can become crucial obstacles to tackling realistic scenarios.
All the investigated benchmarks are restricted to simple physics and two-dimensional
computational domains. This should give other researchers developing DFM models the
chance to perform comparison studies for their methods. We encourage the scientific community to contribute their results for the benchmarks to a corresponding Git repository at
https://git.iws.uni-stuttgart.de/benchmarks/fracture-flow.
Further benchmark cases may be developed in the near future. In particular, we are very
interested in enhancing the purely single-phase single-component flow physics by adding
transport, deformation and/or reaction processes. We aim to carry out these efforts in a
broader context by means of international workshops.

A Acknowledgement
The authors warmly thank Luisa F. Zuluaga, from University of Bergen, for constructing
and providing the real fracture network for the example in Section 4.4. The authors wish
to thank also Luca Pasquale and Stefano Zonca.
The second author acknowledges financial support from the GeoStim project from the
Research Council of Norway (project no. 228832) through the ENERGIX program. The
third author was supported by Norwegian Research Council grant 233736. The fourth author acknowledges financial support from the ANIGMA project from the Research Council
of Norway (project no. 244129/E20) through the ENERGIX program.

B Domain modifications for Benchmark 1
Table 11 provides the exact coordinates of the points from Figure 7. In comparison to the
original setup, the plateaus close to the upper left and right corners 1 and 9 have been
omitted. Moreover, the upper ends of the two fractures have been modified according to
Figure 23 which amounts to the changes of nodes 2–4 and 6–8. Finally, the position of
nodes 16–19 has been recalculated with higher precision. The hybrid-dimensional models
do not take into account nodes 2,4,6,8 and 16–19 and combine nodes 11,12 and 13,14, since
the two-dimensional fracture regions have been reduced to two intersecting straight lines.

28

Table 11: Coordinates of the numbered points in the modeled region of the problem depicted in Figure 7.
pt

x (m)

z (m)

pt

x (m)

z (m)

1
0
150
0
2 394.285714286 100.714285714
30
400
100
40 404.444444444 100.555555556
5
800
150
60 1192.66666667 100.916666667
70
1200
100
0
8
1207.6744186 100.959302326
9
1600
150
10
1600
-1000

11
12
13
14
15
16
17
18
19

1505
1495
1007.5
992.5
0
1071.34615385
1084.03846154
1082.5
1069.80769231

-1000
-1000
-1000
-1000
-1000
-566.346153846
-579.038461538
-587.5
-574.807692308

1

2

5
20

8
60

40

3 30 4

9

80

6 70 7

Figure 23: Modifications of the Hydrocoin model domain compared to the original formulation Swedish Nuclear Power Inspectorate (SKI) (1987). The original upper
boundary is drawn with gray thin lines, while thick black lines are used for the
modified boundary. Modified node locations are indicated by numbers superscripted with 0 . The shaded regions show the upper parts of the two slightly
extended equi-dimensional fractures.

C Fracture coordinates for Benchmark 3
The coordinates are listed in Table 12.

References
Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Kornhuber, R., Ohlberger, M.,
Sander, O., 2008. A generic grid interface for parallel and adaptive scientific computing.
Part II: Implementation and tests in DUNE. Computing 82 (2-3), 121–138.
Bastian, P., Heimann, F., Marnach, S., 2010. Generic implementation of finite element
methods in the distributed and unified numerics environment (DUNE). Kybernetika
46 (2), 294–315.
Berkowitz, B., 2002. Characterizing flow and transport in fractured geological media: A
review. Advances in Water Resources 25 (8-12), 861–884.

29

Nf
1
2
3
4
5
6
7
8
9
10

xA
0.0500
0.0500
0.1500
0.1500
0.6500
0.7000
0.6000
0.3500
0.7500
0.1500

yA
0.4160
0.2750
0.6300
0.9167
0.8333
0.2350
0.3800
0.9714
0.9574
0.8363

xB
0.2200
0.2500
0.4500
0.4000
0.8500
0.8500
0.8500
0.8000
0.9500
0.4000

yB
0.0624
0.1350
0.0900
0.5000
0.1667
0.1675
0.2675
0.7143
0.8155
0.9727

Table 12: Benchmark 3: Fracture coordinates
Boon, W. M., Nordbotten, J. M., Yotov, I., 2016. Robust discretization of flow in fractured
porous media, arXiv:1601.06977 [math.NA].
Brezzi, F., Lipnikov, K., Simoncini, V., 2005. A family of mimetic finite difference methods
on polygonal and polyhedral meshes. Math. Mod. Meth. Appl. S. 15 (10), 1533–1551.
Caers, J., Jan. 2013. A special issue on benchmark problems, datasets and methodologies
for the computational geosciences. Computers & Geosciences 50, 1–3.
Class, H., Ebigbo, A., Helmig, R., Dahle, H. K., Nordbotten, J. M., Celia, M. A., Audigane, P., Darcis, M., Ennis-King, J., Fan, Y., Flemisch, B., Gasda, S. E., Jin, M.,
Krug, S., Labregere, D., Naderi Beni, A., Pawar, R. J., Sbai, A., Thomas, S. G., Trenty,
L., Wei, L., 2009. A benchmark study on problems related to CO2 storage in geologic
formations : Summary and discussion of the results (Original paper). Computational geosciences 13 (4), 409–434, http://link.springer.com/article/10.1007%2Fs10596-009-9146x?LI=true.
Dahle, H. K., Eigestad, G. T., Nordbotten, J. M., Pruess, K., May 2010. A model-oriented
benchmark problem for CO2 storage.
D’Angelo, C., Scotti, A., 2012. A mixed finite element method for darcy flow in fractured
porous media with non-matching grids. ESAIM: Mathematical Modelling and Numerical
Analysis 46 (02), 465–489.
de Araujo Cavalcante Filho, J. S., Shakiba, M., Moinfar, A., Sepehrnoori, K., 2015. Implementation of a preprocessor for embedded discrete fracture modeling in an IMPEC compositional reservoir simulator. In: SPE Reservoir Simulation Symposium, 23-25 February,
Houston, Texas, USA. Society of Petroleum Engineers.
De Dreuzy, J.-R., Pichot, G., Poirriez, B., Erhel, J., 2013. Synthetic benchmark for modeling flow in 3d fractured media. Computers & Geosciences 50, 59–71.
Dietrich, P., Helmig, R., Sauter, M., Hötzl, H., Köngeter, J., Teutsch, G., 2005. Flow and
transport in fractured porous media. Springer.

30

Eymard, R., Gallouët, T., Herbin, R., 2000. Finite volume methods. In: Ciarlet, P. G., Lions, J. L. (Eds.), Solution of Equation in Rn (Part 3), Techniques of Scientific Computing
(Part 3). Vol. 7 of Handbook of Numerical Analysis. Elsevier, pp. 713–1018.
Firoozabadi, A., Monteagudo, J. E. P., Jul. 2004. Control-volume method for numerical
simulation of two-phase immiscible flow in two- and three-dimensional discrete-fractured
media. Water Resources Research 40, W07405.
Flemisch, B., Darcis, M., Erbertseder, K., Faigle, B., Lauser, A., Mosthaf, K., Müthing, S.,
Nuske, P., Tatomir, A., Wolff, M., Helmig, R., 2011. DuMux : DUNE for multi-{phase,
component, scale, physics, . . . } flow and transport in porous media. Advances in Water
Resources 34 (9), 1102–1112.
Flemisch, B., Fumagalli, A., Scotti, A., 2016. A Review of the XFEM-Based Approximation
of Flow in Fractured Porous Media. Vol. 12 of SEMA SIMAI Springer Series. Springer
International Publishing, Cham, Ch. Advances in Discretization Methods, pp. 47–76.
Flemisch, B., Helmig, R., 2008. Numerical investigation of a mimetic finite difference
method. In: Eymard, R., Hérard, J. (Eds.), Finite Volumes for Complex Applications V
– Problems and Perspectives. Wiley - VCH, pp. 815–824.
Formaggia, L., Fumagalli, A., Scotti, A., Ruffo, P., 7 2014. A reduced model for Darcy’s
problem in networks of fractures. ESAIM: Mathematical Modelling and Numerical Analysis 48, 1089–1116.
Frih, N., Martin, V., Roberts, J. E., Saâda, A., 2012. Modeling fractures as interfaces with
nonmatching grids. Computational Geosciences 16 (4), 1043–1060.
Fumagalli, A., Pasquale, L., Zonca, S., Micheletti, S., 2016. An upscaling procedure for
fractured reservoirs with embedded grids. Water Resources Research 52 (8), 6506–6525.
Fumagalli, A., Scotti, A., 2014. An efficient XFEM approximation of Darcy flows in fractured porous media. Oil and Gas Sciences and Technologies - Revue d’IFP Energies
Nouvelles 69.4, 555–564.
Geiger, S., Dentz, M., Neuweiler, I., 2011. A novel multi-rate dual-porosity model for
improved simulation of fractured and multi-porosity reservoirs. In: SPE Reservoir Characterisation and Simulation Conference and Exhibition.
Hajibeygi, H., Karvounis, D., Jenny, P., 2011. A hierarchical fracture model for the iterative
multiscale finite volume method. Journal of Computational Physics 230 (24), 8729–8743.
Hansbo, A., Hansbo, P., 2002. An unfitted finite element method, based on Nitsche’s
method, for elliptic interface problems. Computer methods in applied mechanics and
engineering 191 (47-48), 5537–5552.
Helmig, R., Oct. 1997. Multiphase Flow and Transport Processes in the Subsurface: A
Contribution to the Modeling of Hydrosystems, 1st Edition. Springer.
Hoteit, H., Firoozabadi, A., 2008. Numerical modeling of two-phase flow in heterogeneous
permeable media with different capillarity pressures. Advances in Water Resources 31 (1),
56–73.

31

Huang, H., Long, T. A., Wan, J., Brown, W. P., 2011. On the use of enriched finite element
method to model subsurface features in porous media flow problems. Computational
Geosciences 15 (4), 721–736.
Karimi-Fard, M., Durlofsky, L. J., Aziz, K., 2004. An efficient discrete-fracture model
applicable for general-purpose reservoir simulators. SPE Journal 9 (2), 227–236.
Kolditz, O., Shao, H., Wang, W., Bauer, S., 2015. Thermo-hydro-mechanical-chemical
processes in fractured porous media: modelling and benchmarking. Springer.
Li, L., Lee, S. H., 2008. Efficient field-scale simulation of black oil in a naturally fractured
reservoir through discrete fracture networks and homogenized media. SPE Reservoir
Evaluation & Engineering 11, 750–758.
Logg, A., Mardal, K.-A., Wells, G. N., et al., 2012. Automated Solution of Differential
Equations by the Finite Element Method. Springer.
Martin, V., Jaffré, J., Roberts, J. E., 2005. Modeling fractures and barriers as interfaces
for flow in porous media. SIAM Journal on Scientific Computing 26 (5), 1667–1691.
Moinfar, A., Narr, W., Hui, M.-H., Mallison, B. T., Lee, S. H., 2011. Comparison of
discrete-fracture and dual-permeability models for multiphase flow in naturally fractured
reservoirs. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
Moinfar, A., Varavei, A., Sepehrnoori, K., Johns, R. T., 2014a. Development of an efficient
embedded discrete fracture model for 3d compositional reservoir simulation in fractured
reservoirs. SPE Journal 19 (02), 289 – 303.
Moinfar, A., Varavei, A., Sepehrnoori, K., Johns, R. T., April 2014b. Development of an
efficient embedded discrete fracture model for 3d compositional reservoir simulation in
fractured reservoirs. Society of Petroleum Engineers 19 (2), 289–303.
Müthing, S., 2015. A flexible framework for multi physics and multi domain pde simulations.
Ph.D. thesis, Universität Stuttgart.
Neumann, S. P., 2005. Trends, prospects and challenges in quantifying flow and transport
through fractured rocks. Hydrogeol Journal 13, 124–147.
Nordbotten, J., Flemisch, B., Gasda, S., Nilsen, H., Fan, Y., Pickup, G., Wiese, B., Celia,
M., Dahle, H., Eigestad, G., Pruess, K., Jul. 2012. Uncertainties in practical simulation
of CO2 storage. International Journal of Greenhouse Gas Control 9 (0), 234–242.
Panfili, P., Cominelli, A., 2014. Simulation of miscible gas injection in a fractured carbonate reservoir using an embedded discrete fracture model. In: Abu Dhabi International
Petroleum Exhibition and Conference, 10-13 November, Abu Dhabi, UAE. Society of
Petroleum Engineers.
Panfili, P., Cominelli, A., Scotti, A., 2013. Using Embedded Discrete Fracture Models
(EDFMs) to Simulate Realistic Fluid Flow Problems. In: Second EAGE Workshop on
Naturally Fractured Reservoirs, Muscat, Oman.

32

Reichenberger, V., Jakobs, H., Bastian, P., Helmig, R., Jul. 2006. A mixed-dimensional
finite volume method for two-phase flow in fractured porous media. Advances in Water
Resources 29 (7), 1020–1036.
Sahimi, M., 2011. Flow and transport in porous media and fractured rock: from classical
methods to modern approaches. John Wiley & Sons.
Sandve, T. H., Berre, I., Nordbotten, J. M., 2012. An efficient multi-point flux approximation method for discrete fracture–matrix simulations. Journal of Computational Physics.
Schwenck, N., 2015. An XFEM-based model for fluid flow in fractured porous media. Ph.D.
thesis, University of Stuttgart, Department of Hydromechanics and Modelling of Hydrosystems.
Schwenck, N., Flemisch, B., Helmig, R., Wohlmuth, B., 2015. Dimensionally reduced flow
models in fractured porous media: crossings and boundaries. Computational Geosciences
19 (6), 1219–1230.
Singhal, B. B. S., Gupta, R. P., 2010. Applied hydrogeology of fractured rocks. Springer
Science & Business Media.
Swedish Nuclear Power Inspectorate (SKI), 1987. The international hydrocoin project–
background and results. Paris, France: Organization for Economic Co-operation and
Development.
Tatomir, A.-B., 2012. From discrete to continuum concepts of flow in fractured porous
media. Ph.D. thesis, Universität Stuttgart, Holzgartenstr. 16, 70174 Stuttgart.
Tene, M., Al Kobaisi, M. S., Hajibeygi, H., 2016. Multiscale projection-based Embedded
Discrete Fracture Modeling approach (F-AMS-pEDFM). In: ECMOR XIV-15th European Conference on the Mathematics of Oil Recovery. EAGE.

33

