https://doi.org/10.1007/s10596-021-10120-8

ORIGINAL PAPER

Simulation of flow in deformable fractures using a quasi-Newton
based partitioned coupling approach

Patrick Schmidt1 · Alexander Jaust2 ·Holger Steeb1 ·Miriam Schulte2

Received: 24 March 2021 / Accepted: 29 November 2021
© The Author(s) 2022

Abstract
We introduce a partitioned coupling approach for iterative coupling of flow processes in deformable fractures embedded
in a poro-elastic medium that is enhanced by interface quasi-Newton (IQN) methods. In this scope, a unique computational
decomposition into a fracture flow and a poro-elastic domain is developed, where communication and numerical coupling
of the individual solvers are realized by consulting the open-source library preCICE. The underlying physical problem
is introduced by a brief derivation of the governing equations and interface conditions of fracture flow and poro-elastic
domain followed by a detailed discussion of the partitioned coupling scheme. We evaluate the proposed implementation
and undertake a convergence study to compare a classical interface quasi-Newton inverse least-squares (IQN-ILS) with the
more advanced interface quasi-Newton inverse multi-vector Jacobian (IQN-IMVJ) method. These coupling approaches are
verified for an academic test case before the generality of the proposed strategy is demonstrated by simulations of two
complex fracture networks. In contrast to the development of specific solvers, we promote the simplicity and computational
efficiency of the proposed partitioned coupling approach using preCICE and FEniCS for parallel computations of
hydro-mechanical processes in complex, three-dimensional fracture networks.

Keywords Fracture flow · Hydro-mechanics · Hybrid-dimensional modeling · Partitioned coupling · Quasi-Newton
methods · preCICE

Mathematics Subject Classification (2010) 76S05 · 74F10 · 90C53

� Patrick Schmidt
patrick.schmidt@mechbau.uni-stuttgart.de

Alexander Jaust
alexander.jaust@ipvs.uni-stuttgart.de

Holger Steeb
holger.steeb@mechbau.uni-stuttgart.de

Miriam Schulte
miriam.schulte@ipvs.uni-stuttgart.de

1 Institute of Applied Mechanics (CE), University of Stuttgart,
Pfaffenwaldring 7, D-70 569 Stuttgart, Germany

2 Institute for Parallel and Distributed Systems, University
of Stuttgart, Universitätsstraße 38, D-70 569
Stuttgart, Germany

1 Introduction

Modeling of transient flow processes in deformable high-
aspect ratio fractures (length � aperture) embedded in a
poro-elastic medium is a non-trivial task. Complex discreti-
zation of the fracture(-network) geometry and stiff numer-
ical coupling of the surrounding poro-elastic and fracture-
flow domain require specific implementations to guarantee
stability and efficiency of the designed solver. In this work,
we propose a highly efficient, straightforward implementa-
tion of the governing partial differential equations (PDEs) of
each domain using the open-source package FEniCS [1] and
realizing the numerical coupling process, its acceleration via
quasi-Newton methods and inter-solver communication by
the open-source library preCICE [9].

/ Published online: 20 January 2022

Computational Geosciences (2022) 26:381–400

http://crossmark.crossref.org/dialog/?doi=10.1007/s10596-021-10120-8&domain=pdf
http://orcid.org/0000-0001-9834-3084
mailto: patrick.schmidt@mechbau.uni-stuttgart.de
mailto: alexander.jaust@ipvs.uni-stuttgart.de
mailto: holger.steeb@mechbau.uni-stuttgart.de
mailto: miriam.schulte@ipvs.uni-stuttgart.de


Phenomena such as the occurrence of reverse water-
level fluctuations in distant monitoring wells during
hydraulic testing of fractured reservoirs indicate the
importance of hydro-mechanical interaction throughout
fracture-flow processes [17, 48]. Perturbations of the
reservoir’s equilibrium state result in immediate non-local
fracture deformations and time delayed pressure diffusion
with distinct characteristic time scales, which highly
affect the measured pressure transients. Once the fracture
deformation influences the pressure evolution, traditional
[15, 35] and extended [36] diffusion-based models fail and
the necessity of consistent hydro-mechanical simulations
arises [52]. Creeping flow conditions in high-aspect ratio
fractures motivate to simplify the balance equations based
on pressure-driven, Poiseuille-type flow formulations [32,
56] to capture flow processes in deformable fractures.
The resulting fracture-flow equation leads to a reduction
of space dimension of the computational domain by one.
Equivalently, the numerical implementation is introduced
as a hybrid-dimensional model [52, 53], which can be
simplified to the lubrication equation [2] by dimensional
analysis considering the specific case of high-aspect ratio
geometries [52].

The numerical approaches to solve the strongly coupled
hydro-mechanical fracture-flow problem addressed in liter-
ature can be divided in two groups: (i) monolithic and (ii)
partitioned schemes. Monolithic approaches introduce the
fracture-flow domain by zero-thickness interface elements
[22, 45–47] and mostly require direct solution strategies to
guarantee high numerical stability since distinct character-
istic properties of both domains lead to poor conditioning
of the global system [43]. In general, tailored meshing and
integration strategies are required to discretize the fracture
domain by interface elements. The complexity dramati-
cally increases in the presence of intersecting fractures
or regions containing fracture tips. In contrast, partitioned
schemes allow the individual treatment of the subdomain
and, in particular, the use of non-conformal meshes [43].
However, iterative coupling of both domains is required
to ensure global equilibrium for each conducted time step,
which potentially makes this approach computationally
more expensive. Therefore, stabilization and acceleration
methods are used to keep the number of coupling iterations
needed low. In the literature, numerical convergence and
stability of the equilibrium iterations have been achieved
by adapting physics-based preconditioning of unfractured
[26, 27] to the specific case of fractured poro-elastic media
[4, 10, 18, 19, 25]. Nevertheless, parallel communication
within the proposed staggered algorithms is nontrivial, but
required to solve large systems in reasonable time. There-
fore, we propose an easily accessible, efficient and fully
parallelized deformable fracture-flow implementation intro-
ducing a unique computational split of the fracture-flow

and poro-elastic domain. It uses advanced quasi-Newton
schemes for stabilization and acceleration of the coupling.

More specifically, we employ two separate solvers for
(i) flow and mechanical deformation of the porous struc-
ture and (ii) flow in the fracture itself. This allows for the
use of available, highly efficient and parallel solvers for
each of these two subdomains. We avoid assembling the ill-
conditioned monolithic system of equations. The coupling
library preCICE [9] is used to establish and steer the itera-
tive coupling between the two solvers. preCICE provides a
generic coupling solution interface, that is not tailored to a
specific use case. Its library approach makes it simple to use
with a large variety of solvers. It comes with adapter codes
for popular simulation software and a high-level interface
for many popular programming languages. Additionally,
it has been optimized for high computational efficiency
by implementing several acceleration methods, including
a number of interface quasi-Newton methods [41], paral-
lel data communication [29, 31], and sophisticated data
mapping techniques including radial basis function map-
ping [29, 30] for the use with non-matching grids. This sets
preCICE apart from commercial or closed sourced alter-
natives such as MpCCI1 [23], where one cannot access
or change the implementation. Other open-source alter-
natives such as ADVENTURE2 [24], the Data Transfer
Kit3 [49] or OpenPALM4 [8], e.g., are either lacking the
high-level programming interface or combined implementa-
tion of highly-parallel data communication, coupling steer-
ing, sophisticated data mapping, and coupling acceleration.
Combining the efficient iterative solvers in the subdomains
with preCICE to realize the iterative coupling of both
domains yields a fast and efficient overall solution tech-
nique. It eliminates the need to develop a specialized solver
or preconditioner, respectively, for the corresponding mono-
lithic system. Thus, our simulations of hydro-mechanical
processes at high resolution and accuracy can be performed
by simple high-level solvers written in Python based on the
finite-element library FEniCS. We use the Python interface
of preCICE since the recently developed FEniCS-preCICE
adapter [40] had not been available and was lacking some
functionality needed for the applications presented during
the preparation of our work.

After briefly deriving the governing equations, see
Section 2, and introducing the coupling scheme, see
Section 4, we verify and investigate the proposed coupling
approach for an academic test case (single fracture
embedded in a three-dimensional poro-elastic domain) at

1https://www.mpcci.de
2https://adventure.sys.t.u-tokyo.ac.jp/
3https://ornl-cees.github.io/DataTransferKit/
4https://www.cerfacs.fr/globc/PALM WEB/

382 Comput Geosci (2022) 26:381–400

https://www.mpcci.de
https://adventure.sys.t.u-tokyo.ac.jp/
https://ornl-cees.github.io/DataTransferKit/
https://www.cerfacs.fr/globc/PALM_WEB/


high resolution. The results obtained by the staggered
scheme are compared to a monolithic approach to proof
consistency of both strategies in Section 4.1. Afterwards,
we conduct a case study investigating the convergence
behavior of different interface quasi-Newton schemes and
their dependency on coupling specific parameters, see
Section 4.1.1. We demonstrate the ability of this approach
for large problems for three different meshes ranging from
tens of thousands to several million degrees of freedom
(DoF), see Section 4.1.2. Subsequently, the relevance of the
proposed method for more complex applications is shown
by two more challenging test cases. The first case focuses
on flow processes in complex fracture networks on a short
(minutes) time-scale, see Section 4.2, while, in the second
case, flow through fractured porous media on a large time-
scale (days up to years) is investigated, see Section 4.3. In
summary, this work introduces a stable implementation for
fully coupled, parallelized three-dimensional simulations of
flow processes in deformable fractures, using an unique
split of the computational fracture-flow and poro-elastic
domain.

2 Governing equations

We introduce the governing equations for both parts of our
fracture system, the fracture(s) �Fr and the surrounding
poro-elastic BPe domain, resulting in a coupled formulation
capturing solid deformation and fluid flow in fractured
porous media. For the fracture domain, we base our
derivation on the observation that (i) deformation induced
volume changes of high-aspect ratio fractures have a
strong impact on the flow solution and require an implicit
coupling to the surrounding solid domain; (ii) explicit three-
dimensional modeling of flow processes in such fracture
systems is challenging and leads to poor results in case of
insufficient mesh quality.

We overcome these challenges by lower dimensional
modeling of the fracture flow domain where the fracture
aperture is implicitly treated as a function on a two-
dimensional manifold. The value of this function is given by
an initial opening value and the poro-elastic deformation of
the surrounding porous medium.

2.1 Flow in a deformable fracture�Fr

In the following, we derive a hybrid-dimensional model
[52, 53] governing the flow processes in deformable high-
aspect ratio fractures by evaluation of the balance of mass
and momentum within the fracture domain �Fr. Equilibrium
conditions of the fluid and the biphasic porous medium
are enforced in terms of fracture aperture δ and fluid
pressure p.

2.1.1 Balance of Momentum

In our fracture model, the balance of momentum equation
can be simplified drastically compared to the full Navier-
Stokes equations for viscous and compressible fluids.
Studies on flow processes in hydraulically transmissive
fractures with low contact areas motivate the simplification
to a pressure-driven Poiseuille-type flow between two
parallel plates [52, 56] under creeping flow conditions.
Assuming a quasi-static, deformation dependent fracture
aperture

δ = u+nFr+ + u−nFr− , (1)

the geometrical characteristics of high-aspect-ratio fractures
lead to predominant flow within the fracture plane. In Eq. 1,
we construct the effective fracture aperture δ based on the
jump u± in normal deformation over the fracture surface,
where we introduced the fracture normal vector nFr± (see
Fig. 2). Integration of the velocity profile of the respective
Poiseuille-type (assuming no-slip boundary conditions at
the fracture surfaces) yields the relative fluid velocity

ŵf = −δ2(x, t)
12 ηfR

ˆgrad p̂ = −ksFr(x, t)

ηfR
ˆgrad p̂, (2)

where ηfR denotes the fluid’s effective dynamic viscosity
and ksFr(x, t) := δ2(x, t)/12 the space and time resolved
effective fracture permeability5. Variables and mathematical
operations defined by means of the fracture domain �Fr

are denoted with �̂ to avoid confusion with quantities and
higher dimensional operations defined in the poro-elastic
domain BPe (Fig. 1).

2.1.2 Balance of Mass

Mass conservation has to take into account volumetric
changes and fluid compressibility βf. For a given fluid
compressibility βf = 1/Kf with respect to the fluid’s bulk
modulus Kf, the local mass balance reads

∂(ρfR δ)

∂t
+ ˆdiv

(
ŵf ρfR δ

)
= wN

f , (3)

where wN
f

is the leak-off triggered seepage velocity of the
porous domain normal to the fracture domain.

2.1.3 Governing Equations

The governing equation for the flow in deformable fractures
is obtained by combining (2) with the balance of mass in

5Geometrical characteristics such as the fracture surface roughness can
be governed by adaption of the proportionality factor 1/12 [39].

383Comput Geosci (2022) 26:381–400



Fig. 1 Three dimensional representation of an embedded fracture
characterized by its surface �Fr in a deformable poro-elastic body BPe

with its boundary �Pe. A local coordinate system êi at the fracture level
is introduced, where ê3 is pointing in the direction of the fracture sur-
face normal nFr. The volume V Fr and aperture δ changes of the fracture

between the reference configuration X at time t0 and the current con-
figuration x at time t of a material point P are implicitly coupled to
the poro-elastic deformations expressed by the unique motion function
χ(t)

Eq. 3. To close the set of governing equations, we introduce
a linear equation of state for the fluid pressure p ∝ ρfR

p = Kf

[
ρfR

ρ
fR
0

− 1

]

where ρ
fR
0 is the effective density at t = 0. This results in

the governing equation

δ
∂p̂

∂t︸︷︷︸
I)

− δ3

12 ηfR
ˆgrad p̂ · ˆgrad p̂

︸ ︷︷ ︸
II)

− δ2

12 ηfR βf
ˆgrad δ · ˆgrad p̂

︸ ︷︷ ︸
III)

− 1

12 ηfRβf
ˆdiv

(
δ3 ˆgrad p̂

)

︸ ︷︷ ︸
IV)

+ 1

βf

∂δ

∂t︸ ︷︷ ︸
V)

= wN
f /βf

︸ ︷︷ ︸
VI)

.

(4)

Note that Eq. 4 establishes a lower-dimensional non-
linear partial differential equation for the pressure p̂ in the
fracture, i.e., for a scalar unknown function. It consists of
a transient term I), a quadratic term II), a convection term
III), a diffusion term IV), a volumetric coupling term V)
and a leak-off term VI). The leak-off term q̂lk = wN

f
/βf is

related to the leak-off triggered seepage velocity wN
f
.

The authors of [52] have shown via a dimensional
analysis, that the terms II) and III) can be neglected in
the regime of high aspect-ratio fractures due to their minor

contribution to the overall solution. Thus, the final form of
the fracture flow equation reads

δ
∂p̂

∂t
− 1

12 ηfRβf
ˆdiv

(
δ3 ˆgrad p̂

)

+ 1

βf

∂δ

∂t
= wN

f /βf in �Fr.
(5)

The reduced form introduced by this governing equation
still consists of the pressure diffusion term II) and the
volumetric coupling term V) which require an implicit, non-
trivial coupling to the deformation state of the biphasic
material. In case of fracture intersections, the lower
dimensional formulation induces moderate inaccuracies of
the flow field in limited regions close to the intersection
[37]. In case of high aspect ratio fractures with apertures
in the micrometer range, the deviating flow field in the
intersection area on the overall solution is expected to be
negligible and is, therefore, not considered in this work.
Still, we compare the two resulting apertures at fracture
intersections and choose the larger value considering that
the fluid follows the path of highest permeability .

2.2 Partial Differential Equations in the Poro-Elastic
DomainBPe

Our fracture equation derived above calculates the pressure
evolution for a given deformation and seepage velocity
of the porous matrix. In this section, we present the
respective biphasic poro-elastic formulation [6, 38, 54] for
the porous medium. As the main emphasis of this work is
on the hybrid-dimensional fracture flow formulation and its

384 Comput Geosci (2022) 26:381–400



coupling to the porous medium, we introduce and explain
the governing equations without any explicit derivation of
the balance equations. A more detailed derivation of the
poro-elastic formulation can be found in [43].

2.2.1 Governing Equations

The governing equations of the poro-elastic mixture consist
of a solid and a pore-fluid constituent that are studied
under quasi-static conditions. The unknown functions are
the pore-fluid pressure p and the solid displacement us

given by the following equations:

−div(σ s
E − α p I) = ρ b,

1

M

∂p

∂t
− kf

γ
fR
0

div gradp = −α div
∂us

∂t
,

⎫⎪⎪⎪⎬
⎪⎪⎪⎭

in BPe

with

σ s
E = 3K vol(εs) + 2G dev(εs) + (1 − α) p I.

(6)

The parameters are the Biot parameter α = 1 − K/Ks,
where K is the dry bulk modulus of the solid skeleton, Ks

the bulk modulus of the compressible grains forming the
skeleton and b are the body forces. The (local) storativity
or inverse storage capacity is defined by 1/M = φ0/K

f +
(α − φ0)/K

s with the bulk modulus Kf of the pore fluid
and the porosity in the initial configuration φ0 = φ(t0) of
the mixture relating the partial and effective densities of
the constituents. kf is the Darcy permeability or hydraulic
conductivity, γ

fR
0 the effective weight at t = t0, εs the

elastic strain (split into a volumetric vol(εs) and a deviatoric
dev(εs) part) and G the dry shear modulus of the skeleton
[50, 54].

Compressibility of the porous material is modeled
following the assumption of linear poroelasticity [6]. Fluid-
flow processes and solid deformations within the fracture
and porous domain are now governed by Eqs. 5 and 6.

2.3 Coupling conditions

For a consistent coupling of the fracture and the poro-
elastic domain, the system needs to be closed by coupling
conditions along the fracture surface. The respective
equilibrium conditions along the fracture surface �Fr as
displayed in Fig. 2 extend (5) and (6) to a well-defined
system of equations for the whole simulation domain
comprising fractures and surrounding porous matrix. The
equilibrium state of the fluid flow is reached, once
equivalent exchange between the fracture zone and the
poro-elastic matrix becomes apparent and no-flow boundary
conditions at the fracture tips x̂Tip are met. Flow exchange
between both domains is governed by leak-off flow q̂lk

within the fracture domain and the normal seepage velocity
wN
f
within the poro-elastic matrix (Fig. 2).

The equilibrium conditions of the mechanical state are
met once surface traction tFr = −σ · nFr± , defined by the
total stresses σ , and the fluid pressure p̂ acting normal to
the fracture surface are balanced. Despite this definition of
the coupled system’s equilibrium state via surface terms,
the physical nature of the coupling is volumetric. This is
due to the fluid pressure p̂ being relevant for both, the fluid
and mechanical boundary conditions. It is defined by the
evolution (5) consisting of the volumetric coupling term V).
In conclusion, the necessary boundary conditions along the
coupling interface are

q̂lk = wN
f /βf on �Fr and ŵs = ŵf = 0 at x̂Tip, (7a)

tFr = −p̂ nFr± on �Fr and p̂ = p on �Fr. (7b)

Fig. 2 Representation of flow
and mechanical coupling
conditions along the fracture
surface �Fr. Fracture normal
vectors nFr± are introduced and
can vary along the fracture
surface dependent on the
evaluation position x̂. Flow
coupling conditions are
characterized by the outflow
q̂lk = wN

f
/βf and normal

seepage velocity wN
f
in

combination with a no-flow
requirement at the fracture tips
x̂Tip. Mechanical coupling
conditions are described with
respect to the fracture surface
traction tFr and the acting fluid
pressure p̂ [42]

385Comput Geosci (2022) 26:381–400



2.4Weak formulation of governing equations

We employ a Bubnov-Galerkin finite element scheme to
introduce a spatial discretization to governing equations
(4) and (6). Therefore, the weak form for each domain is
introduced based on the principal of virtual work [e.g., 3].

2.4.1 Deformable fracture�Fr

The weak form governing flow processes within the
fracture domain �Fr is introduced considering a trial
function p̂t and a test function wp̂ for the fluid pressure
p̂, both required to be smooth enough for the applied
mathematical operations. Further, the trial function p̂t must
be suitable to reproduce the applied boundary conditions
and the test function wp̂ needs to vanish in regions where
Dirichlet boundary conditions are defined. Reduction of the
smoothness demand on the trial function p̂t in the second
term of Eq. 4 is achieved by Green’s first identity resulting
in the implemented weak form

∫

�Fr

[
δ wp̂

∂p̂t

∂t
+ δ3

12 ηfRβf
ˆgradwp̂ · ˆgradp̂t

+ wp̂

1

βf

∂δ

∂t

]
da+

∫

LFr
Ne

δ

βf
wp̂ ŵf ·nl

Ne dl=
∫

�Fr
wp̂ q̂lk da.

(8)

In the presence of a lower-dimensional formulation
evaluation of Green’s first identity introduces a Neumann
boundary term along a line segment LFr

Ne in Eq. 8 which
requires a brief discussion. Contrary to the leak-off term,
which considers the fluid exchange with the surrounding
poro-elastic domain regarding the normal seepage velocity,
the introduced Neumann boundary condition defined on
LFr

Ne considers the fluid exchange in longitudinal flow
direction, i.e., relevant once the fracture domain is
intersected by a borehole. The line segment LFr

Ne then
results from the intersection of the fracture domain with the
intersecting geometry, where nl

Ne introduces the outward
pointing normal corresponding to the intersection area. On
a side note, in many applications the explicit dependency
of the boundary term on the fracture aperture δ vanishes,
since it is implicitly considered by the Neumann boundary
condition, i.e., when volumetric flow rate protocols are
applied.

2.4.2 Poro-Elastic DomainBPe

The weak form of the governing equations describing the
hydro-mechanical response within the poro-elastic domain
BPe is introduced considering the testwus and trial functions
us,t for the solid deformation us and the test wp and trial pt

functions for the pore fluid pressure p. Their requirements
regarding mathematical smoothness and fulfillment of
boundary conditions are similar to those required from the
test and trial functions introduced for the fracture domain.
Evaluation of Green’s first identity leads to the implemented
weak form

∫

V Pe

(
σ

E,s
t − α pt I

)
: gradwus dv =

∫

�s
Ne

t̄ · wusda,

∫

V Pe

[
wp

1

M

∂pt

∂t
− kf

γ
fR
0

gradwp · gradpt

+ wp α div
∂us,t
∂t

]
dv =

∫

�
f

Ne

wp w̄N
f da.

(9)

In Eq. 9 the seepage velocity w̄N
f

represents normal
flow conditions on the total Neumann boundary region of
the poro-elastic domain �

f

Ne. The total Neumann boundary
region includes the fracture surface specific region where
the normal seepage velocity wN

f
is evaluated to describe the

fluid exchange between fracture and poro-elastic domain.
Further, we employ a mixed element formulation for the
poro-elastic domain, namely Taylor-Hood elements, by
introducing linear, first order elements for the pore fluid
pressure and quadratic, second order elements for the solid
deformation to guarantee numerical stability [7].

2.4.3 Fluid exchange boundary condition

Considering the fluid exchange of the biphasic poro-elastic
domain with the fracture domain, the normal seepage
velocity is constructed based on the residual state of the
poro-elastic domain obtained for a prescribed fluid pressure
along the fracture surface matching the fluid pressure
state of the fracture domain. Under the consideration of
Green’s theorem, the reconstruction is given by inserting the
obtained solution fields into the right hand side of

∫

�Fr
wp wN

f da =
∫

V Pe

[
wp

1

M

∂pt

∂t

− kf

γ
fR
0

gradwp · gradpt + wp α div
∂us,t
∂t

]
dv

(10)

When considering a finite element discretization we refer to
[e.g., 55, sec. 1.2.] for the node wise evaluation of Eq. 10.
The obtained normal seepage velocity can then be applied
within the fracture domain in terms of the right hand side
of Eq. 8 where we account for the leak-off term by q̂lk =
wN
f

/βf like introduced in Eq. 7a.

386 Comput Geosci (2022) 26:381–400



2.4.4 Treatment of the Poro-Elastic Response on Different
Time Scales

Investigations of systems involving (crystalline) rock
characterized by a low permeability, perturbations in the low
frequency range (� 100 Hz), under undrained conditions,
and with a low viscous pore fluid motivate the application
of Gassmann’s effective low-frequency result [16, 33].
Considering pore pressure effects by effective parameters
in a single-phase, solid formulation introduces two major
numerical advantages since a) the number of degrees of
freedom considered in the rock domain is reduced and b)
numerical instabilities due to distinct diffusion times of
the rock and the fracture domain can be avoided. Still,
investigation periods must be considerably small compared
to the characteristic diffusion time of the rock matrix to
ensure the negligible contribution of leak-off effects. The
effective parameters based on Gassmann’s effective low-
frequency result are determined by

Keff =
φ0

(
1

Ks − 1
Kf

)
+ 1

Ks − 1
K

φ0
K

(
1

Ks − 1
Kf

)
+ 1

Ks

(
1

Ks − 1
K

) ,

Geff = G.

(11)

Here, Keff is the effective bulk and Geff the effective
shear modulus. The poro-elastic effects are governed by the
effective Gassman modulus Keff. In this work, we assume
isotropic linear elasticity of the surrounding bulk material
(with effective elastic parameters Keff and Geff) to avoid
numerical instabilities in the pore pressure solution once the
introduced coupling conditions (7) are fulfilled.

3 Partitioned coupling

The coupling of the poro-elastic medium and the fracture
flow is realized using partitioned coupling methods instead
of a monolithic solver, i.e., we solve the equations for
the fracture domain �Fr (5) and poro-elastic domain BPe

(6) separately. The coupling conditions eq. (7a) and eq.
(7b) are then enforced via suitable boundary conditions for
the subdomains, see Fig. 2, and iterative exchange of the
respective boundary values within a time step.

We use a fracture solver F that maps a seepage velocity
wN
f
and an aperture δ to a pressure p by executing a discrete

time step for eq. (5). The corresponding porous medium
solver S maps the pressure p at the fractures back to a
seepage velocity wN

f
and an aperture δ by executing a

time step for eq. (6) while deriving wN
f

and δ via eq. (10)

and the difference between displacements u± at both sides
of the fracture. This can be considered equivalent to the
typical approach in classical fluid-structure interactions
(FSI) with elastic solid structures as Dirichlet-Neumann
domain decomposition approach [13, 14, 34]. Note that
there are two major differences compared to classical FSI:
(i) the fluid domain is modeled with a lower-dimensional
simplified equation and, thus, the transferred aperture and
seepage velocities are not boundary conditions in the strict
sense; (ii) the fluid exchange between porous matrix and
fractures and the compressibility of the porous structure and
the fluid might be a factor that makes our problem less prone
to instabilities than FSI, in particular with incompressible
fluids. A potential third point are the relatively large fracture
aperture changes that appear in the first time step for
most applied boundary conditions. This makes the coupled
problem hard to solve in the very beginning which is also
the case for the studied test cases in this work.

By iteratively solving one of the fixed-point equations

S ◦ F(δ,wN
s ) = (δ,wN

s ) or (12a)

(
F(δ,wN

s )

S(p)

)
=

(
p

(δ,wN
s )

)
, (12b)

we, thus, ensure fulfilment of all equations and coupling
conditions at the new time step. In this formulation, input
and output of F and S refer to the end point of the respective
time step.

One of the main advantages of the approach presented
here is that it allows for “black box” coupling, i.e., we
can reuse existing solvers that have been tailored for
the problems in the subdomains. If we can provide an
efficient iterative scheme to solve eq. (12), we do not
have to develop a dedicated solver for the overall highly
ill-conditioned system of equations that arises from the
monolithic approach, either. Note that, despite of the
equilibrium state definition for the coupled system via the
fracture surface, the nature of the coupling from the point
of view of the fracture domain �Fr is volumetric. This is
due to the involvement of the fluid pressure p̂ in both, the
fluid and mechanical boundary conditions, which is defined
by evolution equation (4) incorporating the volumetric
coupling term V).

3.1 Iterative Coupling

We present different options to solve Eq. 12 within each
time step. These options have been presented before in [9,
11, 12, 41] and evaluated for classical FSI problems. In this

387Comput Geosci (2022) 26:381–400



Fig. 3 Sketch of the used coupling schemes. For the sake of readability, we display the case without the normal seepage velocity wN
s .

work, we analyze their potential for the lower-dimensional
fracture flow problem. The cheapest version in terms of cost
per time step, but in general known to generate unstable time
stepping is the explicit couplingwere each solver is executed
only once per time step, i.e.,

(δ(n+1),wN,(n+1)
s ) = S ◦ F(δ(n),wN,(n)

s )︸ ︷︷ ︸
=: p(n+1)

or (13a)

(
p(n+1)

(δ(n+1),wN,(n+1)
s )

)
=

(
F(δ(n),wN,(n)

s )

S(p(n))

)
, (13b)

where the superscript (n) denotes the discrete solution at
time tn. Note that the first option in Eq. 13b implies a
sequential one-after-the-other execution of the two solvers
F and S, where S already uses the “new” pressure p(n+1)

as an input, whereas, in the second option, both solvers
can be executed simultaneously. Explicit coupling cannot
capture the strong physical interaction between the fracture
and porous matrix and is not considered further.

Implicit coupling can be achieved via fixed point
iterations

(δ(n+1),i+1,wN,(n+1),i+1
s )

=S ◦ F(δ(n+1),i ,wN,(n+1),i
s )︸ ︷︷ ︸

=: p(n+1),i+1

or (14a)

(
p(n+1),i+1

(δ(n+1),i+1,wN,(n+1),i+1
s )

)

=
(

F(δ(n+1),i ,wN,(n+1),i
s )

S(p(n+1),i )

)
.

(14b)

The first equation refers to a serial-implicit, cf. Fig. 3a,
and the second equation to a parallel-implicit coupling, cf.
Fig. 3b. In the figures, an additional element, the so-called
acceleration scheme, is shown. This numerical component
generates an improved next iterate based on the output of the
respective fixed-point iteration either by under-relaxation

or interface quasi-Newton (IQN) methods. To simplify the
notation in the following, we use the general formulation

H(a) = a (15)

with a vector of unknowns a ∈ R
m at the fractures and the

fixed-point operator H : R
m �→ R

m for both variants of
fixed-point equations in Eq. 12.

We write the general accelerated version of the fixed-
point iterations in eq. (14b) as

ai+1 = A(H(ai)) ,

with the acceleration operator A. A straight-forward option
to realize the acceleration is to use underrelaxation,

ai+1 = ωH(ai) + (1 − ω)ai , ω ∈]0; 1] , (16)

e.g., based on a dynamic Aitken scheme [28]. A more
advanced approach is to reuse past iterates to establish
quasi-Newton iterations [11, 12, 20]6:

ai+1 = ãi + 
ãi with 
ãi = J−1
R̃

(
ãi)R̃(
ãi),

with ãi := H(ãi), the modified residual R̃(ã) := ã −
H−1(ã), and J−1

R̃
(ãi ) as an approximation of the inverse of

the Jacobian of R̃.
To compute J−1

R̃
(ãi ), we collect input-output information

from past iterates of H in tall and skinny matrices Vi, Wi ∈
R

m×i , m � i,

Vi =
[
R̃(ã1) − R̃(ã0), R̃(ã2) − R̃(ã1), . . . , R̃(ãi ) − R̃(ai−1)

]
,

Wi =
[
ã1 − ã0, ã2 − ã1, . . . , ãi − ãi−1

]
.

The matrices Vi andWi define the multi-secant equations
for the inverse Jacobian

J−1
R̃

(ãi ) Vi = Wi . (17)

6To avoid linear dependencies between information from previous
iterations, modified Newton iterations starting from the result of the
pure fixed-point iteration are used (for details, see [51]).

388 Comput Geosci (2022) 26:381–400



To get the classical interface quasi-Newton inverse least-
squares (IQN-ILS) method [11, 12], we close (17) by

‖J−1
R̃

(ãi )‖F → min .

For time-dependent problem, we have to solve the
coupling problem for every time step n. We can extent the
notation of Vi and Wi to V n

i and Wn
i to emphasize the

collection of differences in the current time step. In this
case, the convergence of the quasi-Newton method can be
improved by additionally consider information of previous
time steps, i.e. to use V n−1

i , V n−2
i . . . and Wn−1

i , Wn−2
i . . . .

This is referred to as reuse of time steps method. However,
we cannot store information from an infinite amount of
time due to memory restrictions and since information can
be outdated. Thus, we define a reuse parameter m that
defines for how long information is retained, i.e., we keep
V n−1

i , V n−2
i , . . . , V n−m

i and Wn−1
i , Wn−2

i , . . . , Wn−m
i . If

there are many coupling iterations per time step and, thus,
the matrices V n

i and Wn
i are very large, this can lead to

excessive memory requirements as well. To avoid storing
the information of too many time steps, we also define the
iteration reuse parameter M that limits how many data pairs
over previous coupling iterations and time steps may be kept
in total.

An alternative to IQN-ILS is based on the multi-vector
(MV) approach [5].

The multi-secant equation (17) is closed by

‖J−1
R̃

(ãi ) − J−1
R̃

(ãprev)‖F → min

where J−1
R̃

(ãprev is the approximation of the inverse
Jacobian from the previous time step. The method is also
referred to as interface quasi-Newton inverse multi-vector
Jacobian (IQN-IMVJ) method [41].

In contrast to the IQN-ILS method, we do not need
to store information from previous time steps explicitly,
since the information is kept by the explicit incorporation
of J−1

R̃
(ãprev) in the minimization condition. However, this

also leads to increased memory and runtime requirements
for an increasing number of time steps, as the amount of
information stored in J−1

R̃
(ãi ) and, thus, its size increases.

Therefore, we use a IMVJ flavor with periodic restart to
reduce runtime complexity and storage requirements.

After number η time steps, the method computes a
singular value decomposition (SVD) of the approximated
Jacobian and drops all information connected to singular
values smaller than a user-specified threshold εSVD. We
refer to this restart method as RS-SVD (restarted via
singular value decomposition). For a detailed derivation and
analysis of this restarted method, we refer to [41]. Main
advantages over the ILS method are the low number of
parameters and that it showed less dependency on the choice
of coupling parameters.

Both quasi-Newton methods may suffer from a lack
of stability, if the columns in Vi become (nearly)
linearly dependent. Therefore, additional stabilization of
the numerical method is achieved by filtering to remove
such nearly linearly dependent columns. In this work, we
use a QR filter, also called QR2 filter, which constructs
the QR decomposition QR = Vi and drops columns
that do not add sufficiently much new information to the
problem based on a user-defined filter limit εF. A detailed
description of this filtering technique and a comparison with
other filter techniques can found in literature such as [41,
Scheufele2017].

All grids used in this work are generated such that
the location of the degrees of freedom in both domains
matches on the coupling interface, i.e., the grids are
matching. In general, preCICE also allows for non-matching
meshes by mapping data between the meshes. A variety
of different mapping methods such as low order nearest-
neighbor projection or higher-order approaches using radial
basis functions are available for the mapping. However, the
introduction of data mapping increases the parameter space
for the evaluation of the numerical stability and efficiency
of the investigated coupling approach. Therefore, we limit
this work to matching grids and leave the study of different
data mapping methods and their influence on the coupling
for future work. Implementation details of the partitioned
coupling methods, including the efficient parallelization,
data-mapping techniques etc., are out of the scope of this
work. Instead, we refer to the reference paper of preCICE
[9] and the references therein.

4 Numerical results

Computational efficiency and flexibility of the proposed
partitioned coupling algorithm allow numerical investiga-
tions of non-trivial flow processes in deformable, arbitrarily
fractured porous media in three dimensions for realistic ini-
tial aperture openings in the micrometer range. The capacity
of the proposed method is shown throughout a number of
numerical studies. First, the implementation of the parti-
tioned scheme is verified by a comparison to a reference
solution (computed by a monolithic approach) for a simple
boundary value problem consisting of a single embedded
fracture. Afterwards, the efficiency and robustness of the
method is explored carrying out a study on the convergence
behavior of the interface quasi-Newton schemes and their
dependence on coupling parameters such as reuse m, restart
η, and filter limit εF. The mesh dependency of the solution
is investigated via parallel computations on meshes ranging
from tens of thousand to several million degrees of free-
dom. The work is closed by demonstrating the potential of
the approach to answer relevant questions in a broad range

389Comput Geosci (2022) 26:381–400



of fields including the inverse analysis of pumping opera-
tions and investigations related to nuclear waste disposal;
two fields that clearly require computations with distinct
boundary conditions and time scales.

The implementation of the proposed strategy uses the
open-source libraries preCICE and FEniCS to couple flow
processes in the fracture domain �Fr

s governed by eq. (5)
with responses of the poro-elastic domain BPe

s defined by
eqs. (6), resulting in a non-linear system. In both domains,
the governing equations are solved by standard continuous
Galerkin methods [3, e.g.,] and the grids at the coupling
interface match.

4.1 Verification of the partitioned implementation

Closed form solutions for non-linear flow processes in
deformable fractures are not known to exist. Thus, the
partitioned approach is verified by comparison to a
monolithic scheme.

Focusing on the hydro-mechanical interaction within the
fracture domain, the surrounding bulk material is assumed
to be linear-elastic where poro-elastic effects are transferred
to the material parameters using Gassmann’s effective low-
frequency result, see Eq. (11). The numerical stiffness
of the coupled system of governing equations mainly
scales with (i) the compressibility of the fluid βf and (ii)
the initial fracture aperture δ0. For the limiting case of
βf −→ 0 and δ0 −→ 0, the fracture flow formulation
(5) reduces to the volumetric coupling term V) and the

leak-off term VI) leading to a non-unique solution in case
of partitioned coupling approaches due to the reduction to
Neumann boundary condition terms. Nevertheless, such a
combination of compressibility and initial fracture aperture
is non-physical and has no relevance for experimental
investigations. As we use implicit coupling based on
sophisticated quasi-Newton iterations, we did not observe
numerical instabilities, but an increase in equilibrium
iterations for decreasing values of initial fracture aperture
δ0 and fluid compressibility βf, where the number of
coupling iterations increased stronger with decreasing initial
fracture aperture than with decreasing fluid compressibility.
For the evaluation of the proposed coupling strategy,
we therefore consider a numerically challenging set of
parameters in terms of a low compressible fluid with
properties comparable to those of water under 20◦C and
an initial fracture aperture in the micrometer range. The
chosen parameters have a high relevance for hydraulic
testing operations under in-situ conditions [42; 44, e.g.,].

The parameters used in the investigated boundary value
problem are presented in Tab. 1 and the geometrical set
up is given in Fig. 4. Throughout all numerical studies,
Dirichlet deformation boundary conditions are applied to
the poro-elastic domain BPe by setting deformations on the
outer surfaces in normal direction equal to zero. Within the
fracture domain �Fr

s , pressure Dirichlet conditions p̂s
i are

applied at the intersection of fracture and well. In Fig. 4,
we compare the monolithic and partitioned solutions. We
exploit the radial-symmetric characteristic of the boundary

Table 1 Parameters defining the validation boundary value problem

Quantity Value Unit Quantity Value Unit

Numerical parameters

solid depth dBPe
s

10.0 [m] solid width wBPe
s

10.0 [m]

solid height hBPe
s

10.0 [m] fracture radius r�Fr
s

1.25 [m]

well diameter dw 0.28 [m] solid vertices part. 8.3 · 104 [-]

fracture vertices part. 3.8 · 103 [-] solid vertices mono. 4.4 · 104 [-]

fracture vertices mono. 1.0 · 103 [-] time step size 
t 0.1 [sec]

Rock parameters

dry frame bulk modulus K 8.0 [GPa] grain bulk modulus Ks 33.0 [GPa]

shear modulus G 15.0 [GPa] initial porosity φ0 0.01 [-]

intrinsic permeability ks 1.1 · 10−19 [m2] fluid compressibility βf 0.45 [1/GPa]

effective bulk modulus Keff 29.1 [GPa] effective shear modulus Geff 15.4 [GPa]

Fracture parameters

initial aperture δ0 50.0 [µm] effective fluid viscosity ηfR 0.001 [Pa·s]
fluid compressibility βf 0.45 [1/GPa] injection pressure p̂s

i 200.0 [kPa]

Coupling Parameters

coupling serial-implicit quasi-Newton ILS

filter QR2 initial relaxation ω 0.001 [-]

rel. convergence 10−6 [-]

390 Comput Geosci (2022) 26:381–400



Fig. 4 Left: Sketch of a single
fracture embedded in a solid
matrix used for comparisons of
the monolithic and the
partitioned scheme. Top Right:
Pressure results obtained from
both methods plotted over the
fracture radius for time steps ti
with i = 1, 5, 10, 15, 20, 25.
Bottom Right: Relative error
plots based on the deviation
between solutions obtained from
the monolithic and the
partitioned scheme.

value problem and plot the pressure over the fracture radius
at six different times. We observe very good agreement of
the pressure for the monolithic and partitioned approach for
all time steps. Additionally, we give the deviations between
both strategies evaluated by a relative error

error =
∑Nc

i=1|p̂m
i − p̂

p
i |∑Nc

i=1|p̂m
i | · 100 (18)

The pressure solution of both approaches has been
interpolated to Nc = 50 control points equally distributed
over the fracture radius. The monolithic approach is based
on quad meshes, while the partitioned approach is based on
tetrahedral (porous matrix) and triangular (fracture) meshes.
3 808 degrees of freedom (DoF) were used in the fracture
and 83 109 DoF in the solid domain for the partitioned
approach and 44 418 DoF for the monolithic approach,
where the fracture flow domain consists of 952 DoF. Due
to this, the elements and vertices of the monolithic and the
partitioned coupling strategy do not match and a comparison
of the obtained solutions can not be expected to result in
perfect agreement. Nevertheless, the results given by Fig. 4
show reasonably small errors below 2.5% and decrease over
time.

From a physical perspective, the obtained results are
sound and the expected response in form of reverse water

level fluctuations can be observed at an early stage,
where negative pressure values are induced by non-local
volumetric deformations of the fracture. Summarizing, the
verification indicates the convergence of both coupling
approaches towards the same solution.

4.1.1 Investigation of interface quasi-Newton schemes

To verify the suitability of the partitioned coupling
approach, we run a parameter study with the same settings
as in Tab. 1 for the physical setup and coupling parameters
as given in Tab. 2. All simulations are run in serial-implicit
mode, i.e., we iterate in every time step until we fulfill the
first fixed-point equation in eqs. 12. The initial relaxation
is ω = 0.1 and we use the QR2 filter. In order to keep
the simulation feasible, we use a somewhat coarse mesh
with 3 292 degrees of freedom in the fracture domain and
68 853 degrees of freedom in the porous-medium domain.
If the reuse parameter is set to m = ∞, we allow the ILS
method to keep information from all time steps. In this case,
information from the current or previous time step is only
dropped due to the QR2 filter.

In (12), we have formulated the coupling in terms of
the aperture δ. We refer to this as “pre-accumulated”
case here as it uses δ = u+nFr+ + u−nFr− in the
coupling. An alternative way, that we have studied is the

Table 2 Coupling parameters used in the study of quasi-Newton methods

Quantity Value Unit Quantity Value Unit

Coupling Parameters

coupling {serial-implicit, parallel-implicit} quasi-Newton {ILS, IMVJ}
filter QR2 initial relaxation ω 0.1 [-]

filter limit εF {10−4, 10−3, 10−2, 10−1} relative convergence 10−3 [-]

ILS parameters

Reuse parameter m {0, 8, 16,∞} [-]

IMVJ parameters

SVD threshold εSVD {10−4, 10−3, 10−2, 10−1} [-] Restart parameter η 8 [-]

391Comput Geosci (2022) 26:381–400



Fig. 5 Average number of
coupling iterations per time step.
The serial-implicit coupling and
the IQN-ILS quasi-Newton
method are used. The initial
relaxation parameter is ω = 0.1,
the filter limit and the reuse
parameter are varied.

explicit use of the displacements u± instead. We refer to
this as “post-accumulated” approach. In the latter case,
summation of the displacements to obtain the aperture is
done in the fluid solver and the amount of data has to
be exchanged between the solvers is increased. The latter
case was straightforward to couple in preCICE as it could
rely on already implemented features. We do not expect,
that the different approaches affect the obtained solution,
but rather influence the coupling iteration convergence.
Minimal deviations can be expected as the definition of the
coupling residual differs slightly as it depends on δ in the
pre-accumulated case and on u± otherwise.

In Figs. 5 and 6, the average number of coupling
iterations per time step using the serial-implicit approach
is given for the IQN-ILS and the IQN-IMVJ quasi-Newton
method, while the filter limit and the quasi-Newton method-
specific parameters are varied. The coupling works for all
parameter settings investigated.

The IQN-ILS method, see Fig. 5, requires the largest
number of coupling iterations, when no information from
previous time steps is kept, i.e., the reuse parameter is 0,
and the filter limit is too large, i.e., 10−1. Consequently,
the largest time-averaged number of coupling iterations is
observed, when the reuse parameter is m = 0 and the filter

limit is εF = 10−1. All other cases need an average of
approximately 4.5 coupling iterations per time step, which
is reasonably small.

For the IQN-IMVJ method, see Fig. 6, we observe only
weak dependence on the SVD truncation threshold εSVD,
but stronger dependence on the filter limit. In contrast to the
IQN-ILS method, the average number of iterations increases
for larger (εF = 10−1) and smaller filter limits (εF = 10−4).
Bad choices of coupling parameters lead to approximately
6 coupling iterations per time step, while, in the ideal cases,
only 4 iterations per time step are needed. This is even less
than for the best cases of the ILS method.

The effect of the coupling parameters is very similar
to what is reported in [41], where a study of different
quasi-Newton methods was carried out for fluid-structure
interaction test cases. The authors also observed, that the
IQN-ILS method with m = 0 performs worst and it
depends stronger on the actual choice of parameters than the
IQN-IMVJ method.

Using the pre- or post-accumulation approach has barely
any effect in the current setting. In the observed time
frame, the time-averaged number of coupling iterations
are nearly identical. However, using the aperture δ in
the coupling instead of the displacements u± reduces the

Fig. 6 Average number of
coupling iterations per time step.
The serial-implicit coupling and
the IQN-IMVJ quasi-Newton
method are used. The initial
relaxation parameter is ω = 0.1,
the filter limit and the reuse
parameter are varied.

392 Comput Geosci (2022) 26:381–400



Fig. 7 Average number of
coupling iterations per time step.
The parallel-implicit coupling
and both quasi-Newton methods
are used. The initial relaxation
parameter is ω = 0.1 and the
filter limit and quasi-Newton
specific parameters are varied.

amount of data, that needs to be exchanged. Additionally,
we observed strongly improved coupling stability for
the fracture networks using the pre-accumulated coupling
approach and, thus, it has been used in all following test
cases.

In Fig. 7, we show the average number of coupling
iterations per time step when a parallel-implicit coupling
method is used. The coupling behavior is really similar
to the serial-implicit approach, see Figs. 5 and 6. The
IQN-IMVJ method needs less iterations than the IQN-
ILS method and shows less dependency on the coupling
parameters. For the IQN-ILS, one needs most coupling
iterations, again, for very small filter limits and especially
when m = 0. Surprisingly, the parallel-implicit even beats
the serial-implicit coupling in terms of coupling iterations
needed. This was not expected as the parallel-implicit
approach is weaker and thus normally needs more iterations
to recover the strong coupling behavior of the underlying
physical problem. It is not clear what causes this and should
be investigated further.

4.1.2 Mesh convergence study

For the numerical investigation of mesh convergence,
we use three meshes for the fluid and the poro-elastic
subdomain with different resolution: (i) (68 853, 3 292)

degrees of freedom (coarse, the same mesh as in the
parameter study), (ii) (447 816, 24 188) degrees of freedom
(medium), and (iii) (6 586 770, 323 976) degrees of freedom
(fine) in the porous matrix and fluid domain.

In Fig. 8, we present the pressure over the fracture radius
at the final time t = 2.5 for the different meshes. In
all cases, a reasonable pressure curve is obtained. On the
coarsest mesh, the pressure is highest. On the mediummesh,
the pressure is clearly lower than on the coarse mesh. Thus,
the predicted pressure tends to get lower for higher mesh
resolution. This is confirmed by the solution on the finest
mesh where the pressure is again lower than on the medium
mesh. At the same time, the pressure difference between
the fine and the medium mesh is smaller than between the
medium and the coarse mesh although the difference in grid
points increased severely. This indicates that the simulations
converges toward a grid-converged solution.

The simulation setup was tested on shared-memory and
distributed-memory systems with up to 256 cores. No
adjustments had to be made to the code as the parallelization
is handled internally via FEniCSand preCICE. However,
parallel efficiency is not the focus of this work, but the
simplicity of the current approach to realize solvers and
couplings that can be executed on parallel computers.
Therefore, we do not report scaling results and instead leave
it for future work.

Fig. 8 Mesh convergence study
for different mesh sizes. We plot
the pressure over the fracture
radius at time t = 2.5 similar to
Fig. 4.

393Comput Geosci (2022) 26:381–400



4.2 Injection and production in an arbitrarily
fractured reservoir

Transient flow and pressure data obtained by experimental
field operations on fractured reservoirs provide information
of their storage capacity, when evaluated by best numeri-
cal fits. Computational effort might be reduced for specific
investigations on circular fractures by using two dimen-
sional radial-symmetric models, but most field settings
require consideration of several interacting fractures and
three-dimensional modeling.

Due to the low permeability of the surrounding
bulk matrix, experimental pumping operations are often
performed on fractured granite reservoirs. Such a problem
setting might lead to instabilities throughout the numerical
analysis, since characteristic pressure diffusion times of
fracture network and granite bulk material greatly differ.
Based on the short experimental execution time, outflow
into the surrounding bulk material can be neglected and
effective material parameters can be introduced based on
Gassmann’s solution defined by Eq. (11). This reduces
the matrix response to linear-elastic behavior. Hence, the
hydro-mechanically coupled system consists of the hybrid-
dimensional flow Eq. (5) and the linear-elastic response of
Eq. (6), similar to the problem solved in Section 4.

Here, we demonstrate the capability of the proposed
method to solve flow problems in fracture networks
embedded in a low permeable porous bulk material.
Therefore, we test our approach on an arbitrarily generated
fracture network containing 17 fractures by inducing
injection and production Dirichlet pressure boundary
conditions in the fracture domain �Fr

PN, see Fig. 9. The
regions of injection and production are highlighted. The

parameters describing the boundary value problem and the
coupling parameters are given in Tab. 3.

The applied coupling strategy shows a convergence
behavior, which is characteristic for quasi-Newton schemes.
Convergence in the first, critical time step is reached within
27 iterations, before the required number of iterations
reduces to 14 in the second time step and reaches its
minimum of 4–5 iterations for later time steps. The
characterization of the tested network is carried out by
investigations on the preferential flow path through the
network connecting the injection with the production well.
The chosen time step size resolves the fracture aperture
evolution and allows to study transient hydro-mechanical
effects such as the inverse pressure response at an early
stage of the simulation. Nevertheless, the results are
evaluated at time t = 900 sec to focus on a solution close to
the quasi-static equilibrium.

In Fig. 9, the aperture changes of the fracture network
indicate strong hydro-mechanical interaction showing, that
fractures with dominant opening behavior are closing
fractures with similar orientations by reallocation of the
surrounding bulk material. This phenomenon is evident
when looking at the aperture change distribution of fractures
9 and 10 and to some extend for fractures 1 and 2, or 15 and
17, respectively.

The pressure distribution in Fig. 10 shows a smooth
pressure field, where pressure drops between fractures
are highest, when large fractures are connected by small
fractures which is due to their lower cross-section area
and higher geometrical stiffness. The phenomenon is
demonstrated best by the pressure drop between large
fractures 3 and 8 interconnected by the small fractures 5
and 7.

Fig. 9 Left: Sketch of the connected fracture network �Fr
Nw embedded in a solid matrix BPe

Nw. Right: Change of fracture aperture after t = 900
secs including numbering of the embedded fractures.

394 Comput Geosci (2022) 26:381–400



Table 3 Parameters defining the injection-production boundary value problem

Quantity Value Unit Quantity Value Unit

Numerical parameters

solid depth dBPe
Nw

50.0 [m] solid width wBPe
Nw

50.0 [m]

solid height hBPe
Nw

50.0 [m] small fracture radius rsm
�Fr
Nw

2.0 [m]

large fracture radius rla
�Fr
Nw

4.5 [m] solid DoFBPe
Nw

9.5 · 105 [-]

fracture DoF�Fr
Nw

5.0 · 104 [-] time step size 
t 5.0 [sec]

Rock parameters

dry frame bulk modulus K 8.0 [GPa] grain bulk modulus Ks 33.0 [GPa]

shear modulus G 15.0 [GPa] initial porosity φ0 0.01 [-]

intrinsic permeability ks 1.1 · 10−19 [m2] fluid compressibility βf 0.45 [1/GPa]

effective bulk modulus Keff 29.1 [GPa] effective shear modulus Geff 15.4 [GPa]

Fracture parameters

initial aperture δ0 75.0 [µm] effective fluid viscosity ηfR 0.001 [Pa·s]
fluid compressibility βf 0.45 [1/GPa] injection pressure p̂Nw

i 200.0 [kPa]

production pressure p̂Nw
p −200.0 [kPa]

Coupling Parameters

coupling serial-implicit quasi-Newton ILS

filter QR2 initial relaxation ω 0.001 [-]

rel. convergence 10−5 [-]

The calculated flow solution shown in Fig. 10 visualizes
the preferential flow path of the system through fractures 2,
4, 6, 9, 15 and 17 and confirms that the dominant opening
of fracture 9 leads to the reduction of the flow through
fracture 10. Regions of high flow rates are small connecting
fractures and regions close to the injection or production
area, which is consistent with the investigated pressure and
aperture distributions. The study is representative for hydro-
mechanical investigations on tested networks embedded in
low permeable rock at a short time scale. It emphasizes
the numerical capacity of the proposed approach and its
relevance for detailed investigations of research questions in
the field of experimental pumping operations.

4.3 Flow through fractured porous media

Approximations of preferential flow patterns through
fractured poro-elastic media are of high interest in the
field of nuclear waste disposal to reduce the risk of
potential pollution by leak-off of contaminated matter. In
contrast to the previous test case, an entirely different time
scale is required, since investigation periods might last up
to a million years. Discrete fracture networks influence
the effective transport characteristics of a reservoir and
even slight hydro-mechanically induced changes of the
fractures’ permeability have an immediate impact on its
final characteristic diffusion time.

Fig. 10 Left: Pressure state at time t = 900 secs highlighting positions of pressure injection p̂NW
i and production p̂NW

p . Right: Post-processed
fluid flow field obtained by inserting pressure and aperture solutions into the balance of momentum (2) at time t = 900 secs.

395Comput Geosci (2022) 26:381–400



Fig. 11 Left: Sketch of an arbitrarily chosen fracture network �Fr
PN embedded within a poro-elastic matrix highlighting the fluid pressure boundary

conditions p0 and p1 applied to the poro-elastic domain BPe
PN. Right: Change of fracture aperture after t = 1.65 days including numbering of the

embedded fractures.

The following boundary value problem investigates
fractures embedded in a poro-elastic matrix with a low
permeability. Here, we solve the fully coupled system
consisting of the hybrid-dimensional flow eq. (5) and the
set of poro-elastic governing eqs. (6), where fluid exchange
between fracture and poro-elastic domain is considered by
evaluating eq. (10). We induce a pressure gradient on the
poro-elastic domain by prescribing the fluid pressure p0

and p1 at the top and bottom, see Fig. 11. The parameters
defining the boundary value problem and the numerical

coupling are given in Tab. 4, the geometrical set up is
shown in Fig. 11. The embedded fractures can be grouped
into three single fractures (fractures 3, 15 and 17), a small
fracture network consisting of fractures 1 and 2 and a large
fracture network formed by the remaining fractures (4–14
and 16). The position, shape and orientation of the fractures
is chosen arbitrarily to demonstrate the flexibility of the
method.

The applied coupling parameters result in a stable
convergence behaviour. Convergence is reached within 31

Table 4 Parameters defining the fractured porous media boundary value problem

Quantity Value Unit Quantity Value Unit

Numerical parameters

solid depth dBPe
PN

30.0 [m] solid size wBPe
PN

30.0 [m]

solid height hBPe
PN

20.0 [m] smallest fracture are Asm
�Fr
PN

240.5 [m2]

largest fracture area Ala
�Fr
PN

4.9 [m2] solid DoFBPe
PN

6.0 · 105 [-]

fracture DoF�Fr
PN

3.2 · 104 [-] time step width 
t 200.0 [sec]

Rock parameters

dry frame bulk modulus K 8.0 [GPa] grain bulk modulus Ks 33.0 [GPa]

shear modulus G 15.0 [GPa] initial porosity φ0 0.01 [-]

intrinsic permeability ks 1.0 · 10−17 [m2] fluid compressibility βf 0.45 [1/GPa]

effective bulk modulus Keff 29.1 [GPa] effective shear modulus Geff 15.4 [GPa]

applied pressure pPN
0 200.0 [kPa] applied pressure pPN

1 0.0 [kPa]

Fracture parameters

initial aperture δ0 75.0 [µm] effective fluid viscosity ηfR 0.001 [Pa·s]
fluid compressibility βf 0.45 [1/GPa]

Coupling Parameters

coupling serial-implicit quasi-Newton ILS

filter QR2 initial relaxation ω 0.001 [-]

rel. convergence 10−5 [-]

396 Comput Geosci (2022) 26:381–400



Fig. 12 Left: Post-processed fluid-flow field in the poro-elastic
domain BPe

PN after t = 1.65 days obtained by inserting the pressure
solution into the governing equation (6) and neglecting of time depen-
dent terms by assuming quasi static conditions. Right: Post-processed

fluid-flow field within the fracture domain �Fr
PN after t = 1.65 days

obtained by inserting the pressure and aperture solutions into the
balance of momentum (2).

iterations in the first, critical time step, 15 iterations in the
following step and 4–5 iterations for later time steps.

Investigations on the impact of hydraulically highly con-
ductive deformable fractures on the transport characteristic
of the tested reservoir are evaluated by means of the pref-
erential flow pattern in the poro-elastic and the fracture
domain. Results are displayed for a stage close to the con-
verged quasi-static equilibrium after t = 1.65 days, see
Figs. 11 and 12. Similar to the findings of the previous study,
hydro-mechanical interactions are evident in the fracture
aperture change distribution displayed, where the fracture
pairs 7 and 8, 9 and 14 and 4 and 8 have the strongest inter-
actions, see Fig. 11. The interaction results in opening or
closing of the involved fractures, which has an immediate
impact on the local conductivity of the network.

The post-processed flow solutions in the poro-elastic
BPe
PN and in the fracture domain �Fr

PN are presented in Fig.
12. Depending on the fracture orientation and the fracture
connectivity, embedded fractures have a distinct impact on
the flow through the poro-elastic medium. Related to the
orientation orthogonal to the pressure gradient and the lack
of connectivity, the embedded single fractures 3, 15 and 17
have a minor contribution. At the same time, the fluid is
strongly attracted by the small and large fracture network.
The fluid mainly enters the networks through fractures 1,
4, 5 and 7 which are closest to the pressure boundary p0

and is released back into the poro-elastic matrix through
fractures 13 and 16 which are closest to the applied pressure
boundary p1. The impact of the fracture network on the
transport characteristic of the studied poro-elastic domain
is evident in terms of the distinct difference of flow rates
in both domains. The hydro-mechanical interaction between
the fracture and the poro-elastic domain has shown the

potential of the method to consider flow through deformable
fractures embedded in a hydro-mechanically interacting
poro-elastic medium throughout long term investigations
with a relevance for fields such as nuclear waste disposal.

5 Conclusion and outlook

We proposed a new partitioned coupling approach for
hydro-mechanical flow processes in deformable fractures
embedded in a poro-elastic medium. Implicit coupling of
the decomposed fracture flow and poro-elastic domain
under consideration of the introduced interface conditions
was realized by an iterative approach. The latter solves
the underlying fixed-point problem using interface quasi-
Newton methods. It was implemented using the open-source
computing platform FEniCSto solve the individual systems
of PDEs and the open-source coupling library preCICEto
realize the implicit coupling.

We showed, that the proposed coupling strategy enables
straight-forward usage of parallel computations throughout
solution and coupling steps. This allows to study complex
fracture systems in three dimensions with a high computa-
tional resolution. Evaluation of the proposed implementa-
tion against solutions obtained from a monolithic approach
showed good agreement in terms of the transient pressure
evolution in a single deformable fracture. Throughout a
coupling convergence study, we showed the slightly bet-
ter performance of the advanced IQN-IMVJ quasi-Newton
scheme in comparison to the performance of the classical
IQN-ILS quasi-Newton scheme, which is in good agree-
ment with results on classical fluid-structure interaction
problems in [41].

397Comput Geosci (2022) 26:381–400



The generality of the proposed strategy and its relevance
for research topics such as modeling of injection and
production in fractured reservoirs and the flow through a
fractured poro-elastic domain was demonstrated throughout
two numerical studies of complex fracture networks in three
dimensions. We emphasized the advantage of the fracture
and poro-elastic domain decomposition in terms of the
creation of complex networks and straight forward post-
processing of the numerical solutions in each computational
region to identify preferential flow paths through the
fracture network and the poro-elastic medium, respectively.

Future work can focus on the extension of the physical
models to include temperature, e.g., extension of the
partitioned coupling schemes, and the investigation of other
discretization methods in for the subproblems. Considering
additional physics in the model opens new application such
as heat related energy production relevant for geothermal
applications. The partitioned coupling schemes, especially
in the black-box setting of preCICE, can be improved by
developing more sophisticated start-up strategies to reduce
the number of coupling iterations in the first time steps and
by new data mapping and communication concepts. The
mixed-dimensional modeling leads to several challenges on
how to communicate data between the different models,
especially at fracture intersections, where the dominant
deformation has to be identified.

Acknowledgments Holger Steeb and Patrick Schmidt gratefully
acknowledge the funding provided by the German Federal Ministry
of Education and Research (BMBF) for the GeomInt (I &
II) project (Grant Numbers 03A0004E and 03G0899E) in the
BMBF Geoscientific Research Program “Geo:N Geosciences for
Sustainability”. Alexander Jaust, Miriam Schulte, and Holger Steeb
thank the German Research Foundation (DFG) for supporting this
work under Grant No. SFB 1313 (Project No. 327154368). We thank
the preCICE developers for their support, especially B. Uekermann.

Funding Open Access funding enabled and organized by Projekt
DEAL. Holger Steeb and Patrick Schmidt gratefully acknowledge the
funding provided by the German Federal Ministry of Education and
Research (BMBF) for the GeomInt (I & II) project (Grant Numbers
03A0004E and 03G0899E) in the BMBF Geoscientific Research
Program “Geo:N Geosciences for Sustainability”. Alexander Jaust,
Miriam Schulte, and Holger Steeb thank the DFG for supporting this
work under Grant No. SFB 1313 (Project No. 327154368).

Availability of data andmaterial The data sets and codes are available
via the data repository (DaRUS) of the University of Stuttgart [21]
which has the identifier https://doi.org/10.18419/darus-1778 assigned.

Code availability See Availability of data and material.

Declarations

Ethics approval Not applicable

Consent to participate Not applicable

Consent for publication Not applicable

Conflicts of interest/Competing interests Not applicable

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as
long as you give appropriate credit to the original author(s) and the
source, provide a link to the Creative Commons licence, and indicate
if changes were made. The images or other third party material in this
article are included in the article’s Creative Commons licence, unless
indicated otherwise in a credit line to the material. If material is not
included in the article’s Creative Commons licence and your intended
use is not permitted by statutory regulation or exceeds the permitted
use, you will need to obtain permission directly from the copyright
holder. To view a copy of this licence, visit http://creativecommons.
org/licenses/by/4.0/.

References

1. Alnæs, M.S., Blechta, J., Hake, J., Johansson, A., Kehlet, B.,
Logg, A., Richardson, B.C., Ring, J., Rognes, M.E., Wells, G.N.:
The FEniCS Project Version 1.5. Arch. Numer. Softw. 3.100
(2015)

2. Batchelor, G.K.: An introduction to fluid dynamics. Cambridge
mathematical library. Cambridge University Press (2000)

3. Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K.: Nonlinear
finite elements for continua and structures. Wiley (2013)

4. Berge, R.L., Berre, I., Keilegavlen, E., Nordbotten, J.M.,
Wohlmuth, B.: Finite Volume Discretization for Poroelastic Media
with Fractures Modeled by Contact Mechanics. Int. J. Numer.
Methods Eng. 121.4, 644–663 (2020)

5. Bogaers, A. E. J.., Kok, S., Reddy, B. D.., Franz, T.: Quasi-Newton
methods for implicit black-box FSI coupling. Computer Methods
in Applied Mechanics and Engineering 279, 113–132 (2014)

6. Biot, M.A.: General Theory of Three-Dimensional Consolidation.
J. Appl. Phys. 12.2, 155–164 (1941)

7. Brezzi, F., Fortin, M.: Mixed and hybrid finite element methods.
vol. 15. Springer Science & Business Media (2012)

8. Buis, S., Piacentini, A., Déclat, D.: PALM: A Computational
Framework for Assembling Highperformance Computing Appli-
cations. Concurr. Comput. Practice Exper. 18.2, 231–245 (2006)

9. Bungartz, H.-J., Lindner, F., Gatzhammer, B., Mehl, M.,
Scheufele, K., Shukaev, A., Uekermann, B.: PreCICE - A Fully
Parallel Library for Multi-Physics Surface Coupling. Comput.
Fluids, vol. 141. Advances in Fluid- Structure Interaction,
pp. 250–258 (2016)

10. Castelletto, N., White, J.A., Tchelepi, H.A.: Accuracy and
Convergence Properties of the Fixed-Stress Iterative Solution of
Two-Way Coupled Poromechanics. Int. J Numer. Anal. Methods
Geomechan. 39.14, 1593–1618 (2015)

11. Degroote, J., Bruggeman, P., Haelterman, R., Vierendeels, J.:
Stability of a coupling technique for partitioned solvers in FSI
applications. Computers & Structures 86(23), 2224–2234 (2008)

12. Degroote, J., Haelterman, R., Annerel, S., Bruggeman, P.,
Vierendeels, J.: Performance of partitioned procedures in fluid-
structure interaction. Computers & Structures 88(7), 446 - 457
(2010)

13. Deparis, S., Discacciati, M., Fourestey, G., Quarteroni, A.:
Fluid-structure algorithms based on Steklov-Poincaré operators.
Computer Methods in Applied Mechanics and Engineering
195(41), 5797-5812 (2006). John H. Argyris Memorial Issue. Part
II

14. Farhat, C.: CFD-Based Nonlinear Computational Aeroelastic-
ity. In: Encyclopedia of Computational Mechanics. John Wiley

398 Comput Geosci (2022) 26:381–400

https://doi.org/10.18419/darus-1778
http://creativecommons.org/licenses/by/4.0/
http://creativecommons.org/licenses/by/4.0/


& Sons (2004). https://doi.org/10.1002/0470091355.ecm063.
https://onlinelibrary.wiley.com/doi/abs/10.1002/0470091355.
ecm063

15. Fetter, C.: Applied Hydrogeology. 4th edn. Prentice Hall (2001)
16. Gassmann, F.: Über Die Elastizität Poröser Medien. Vierteljahrss-

chrift Naturforsch. Gesellschaft Zürich 96, 1–23 (1951)
17. Gellasch, C.A., Wang, H.F., Bradbury, K.R., Bahr, J.M., Lande,

L.L.: Reverse Water- Level Fluctuations Associated with Fracture
Connectivity. Groundwater 52.1, 105–117 (2014)

18. Girault, V., Wheeler, M.F., Ganis, B., Mear, M.E.: A Lubrication
Fracture Model in a Poro-Elastic Medium. Math. Models Methods
Appl. Sci. 25.04, 587–645 (2015)

19. Girault, V., Kumar, K., Wheeler, M.F.: Convergence of Iterative
Coupling of Geomechanics with Ow in a Fractured Poroelastic
Medium. Comput. Geosci. 20.5, 997–1011 (2016)

20. Haelterman, R., Degroote, J., Van Heule, D., Vierendeels, J.: The
Quasi-Newton Least Squares Method: A new and fast secant
method analyzed for linear systems. SIAM Journal on Numerical
Analysis 47(3), 2347-2368 (2009)

21. Jaust, A., Schmidt, P.: Replication Data for Simulation of ow
in deformable fractures using a quasi-Newton based partitioned
coupling approach. Version V1 (2021)

22. Jin, L., Zoback, M.D.: Fully Coupled Nonlinear Fluid Ow and
Poroelasticity in Arbitrarily Fractured Porous Media: a Hybrid-
Dimensional Computational Model. J. Geophys. Res. Solid Earth
122.10, 7626–7658 (2017)

23. Joppich, W., Kürschner, M.: MpCCI — a tool for the simulation
of coupled applications. Concurr. Comput. Practice Exper. 18.2,
183–192 (2006)

24. Kataoka, S., Minami, S., Kawai, H., Yamada, T., Yoshimura, S.: A
Parallel Iterative Partitioned Coupling Analysis System for Large-
Scale Acoustic Fluid-Structure Interactions. Comput. Mech. 53.6,
1299–1310 (2014)

25. Keilegavlen, E., Berge, R., Fumagalli, A., Starnoni, M., Stefans-
son, I., Varela, J., Berre, I.: PorePy: An Open-Source Simulation
Tool for Flow and Transport in Deformable Fractured Rocks.
Comput. Geosci. 25.1, 243–265 (2021)

26. Kim, J., Tchelepi, H., Juanes, R.: Stability and Convergence
of Sequential Methods for Coupled Ow and Geomechanics:
Drained and Undrained Splits. Comput. Methods Appl. Mech.
Eng. 200.23, 2094–2116 (2011)

27. Kim, J., Tchelepi, H., Juanes, R.: Stability and Convergence of
Sequential Methods for Coupled Ow and Geomechanics: Fixed-
Stress and Fixed-Strain Splits. Comput. Methods Appl. Mech.
Eng. 200.13, 1591–1606 (2011)

28. Küttler, U., Wall, W.: Fixed-point fluid-structure interaction
solvers with dynamic relaxation. Computational Mechanics 43,
61-72 (2008). 10.1007/s00466-008-0255-5

29. Lindner, F.: Data transfer in partitioned multiphysics simulations :
interpolation & communication. PhD thesis University of Stuttgart
(2019)

30. Lindner, F., Mehl, M., Uekermann, B.: Radial basis function inter-
polation for Black-Box Multi-Physics simulations. In: Conference
Proceedings at the ECCOMAS Coupled Problems (2017)

31. Lindner, F., Totounferoush, A., Mehl, M., Uekermann, B.,
Pour, N.E., Krupp, V., Roller, S., Reimann, T., Sternel, D.C.,
Egawa, R., Takizawa, H., Simonis, F.: ExaFSA Parallel Fluid-
Structure-Acoustic Simulation. Software for Exascale Computing
- SPPEXA 2016-2019. In: Bungartz, H.-J., Reiz, S., Uekermann,
B., Neumann, P., Nagel, W.E. (eds.), pp. 271–300. Springer
International Publishing, Cham (2020)

32. Louis, C.: A study of groundwater ow in jointed rock and its in
uence on the stability of rock masses. Imperial College. Rock
Mech. Res. Rep. 10, 1–90 (1969)

33. Mavko, G., Mukerji, T., Dvorkin, J.: The rock physics handbook:
Tools for seismic analysis of porous media. Cambridge University
Press (2009)

34. Monge, A., Birken, P.: On the convergence rate of the
Dirichlet-Neumann iteration for unsteady thermal fluid-structure
interaction. Computational Mechanics 62(3), 525-541 (2018)

35. Muskat, M.: The Ow of Homogeneous Fluids through Porous
Media. Soil Sci. 46.2, 169 (1938)

36. Ortiz, A.E.R., Renner, J., Jung, R.: Hydromechanical analyses of
the hydraulic stimulation of borehole Basel 1. Geophys. J. Int.
185.3, 1266–1287 (2011)

37. Quintal, B., Caspari, E., Holliger, K., Steeb, H.: Numerically
quantifying energy loss caused by squirt ow. Geophys. Prospect.
67.8, pp. 2196–2212 (2019)

38. Renner, J., Steeb, H.: Modeling of Fluid Transport in Geothermal
Research. Handbook of Geomathematics. Berlin, Heidelberg:
Springer Berlin Heidelberg, pp. 1443–1500 (2015)

39. Renshaw, C.E.: On the Relationship between Mechanical and
Hydraulic Apertures in Rough-Walled Fractures. J. Geophys Res.
Solid Earth 100.B12, 24629–24636 (1995)

40. Rodenberg, B., Desai, I., Hertrich, R., Jaust, A., Uekermann, B.:
FEniCS-preCICE: Coupling FEniCS to other Simulation Software
(2021)

41. Scheufele, K., Mehl, M.: Robust Multisecant Quasi-Newton Vari-
ants for Parallel Fluid-Structure Simulations—and Other Mul-
tiphysics Applications. SIAM Journal on Scientific Computing
39(5), S404-S433 (2017)

42. Schmidt, P., Dutler, N., Steeb, H.: Importance of Fracture Defor-
mation throughout Hydraulic Testing under In-Situ Conditions.
Geophys. J. Int. (2021)

43. Schmidt, P., Steeb, H.: Numerical Aspects of Hydro-Mechanical
Coupling of Fluid-Filled Fractures Using Hybrid-Dimensional
Element Formulations and Non-Conformal Meshes. GEM - Int. J.
Geomathematics 10.1, 14 (2019)

44. Schmidt, P., Steeb, H., Renner, J.: Investigations into the
Opening of Fractures during Hydraulic Testing Using a Hybrid-
Dimensional Ow Formulation. Environ. Earth Sci. 80.497 (2021)

45. Segura, J.M., double, I.Carol.: Coupled HM Analysis Using Zero-
Thickness Interface Elements with Nodes—Part II: Verification
and Application. Int. J. Numer. Anal. Methods Geomech. 32.18,
2103–2123 (2008)

46. Segura, J.M.: Coupled HM analysis using zero-thickness interface
elements with double nodes—Part II: Verification and application.
Int. J. Numer. Anal. Methods Geomech. 32.18, 2083–2101 (2008)

47. Settgast, R.R., Fu, P., Walsh, S.D., White, J.A., Annavarapu, C.,
Ryerson, F.J.: A Fully Coupled Method for Massively Parallel
Simulation of Hydraulically Driven Fractures in 3-Dimensions.
Int. J. Numer. Anal. Methods Geomech. 41.5, 627–653 (2017)

48. Slack, T.Z., Murdoch, L.C., Germanovich, L.N., Hisz, D.B.:
Reverse Water-Level Change during Interference Slug Tests
in Fractured Rock. Water Resour. Res. 49.3, 1552–1567
(2013)

49. Slattery, S.R., Wilson, P.P.H., Pawlowski, R.P.: The Data Transfer
Kit: A geometric rendezvous-based tool for multiphysics data
transfer. In: Proceedings of the 2013 International Conference
on Mathematics and Computational Methods Applied to Nuclear
Science and Engineering - M and C 2013 (2013)

50. Steeb, H., Renner, J.: Mechanics of Poroelastic Media a Review
with Emphasis on Foundational State Variables. Transport Porous
Media 130, 437–461 (2019)

51. Uekermann, B.: Partitioned Fluid-Structure Interaction on Mas-
sively Parallel Systems. Institut für Informatik, Technische
Universität München. https://mediatum.ub.tum.de/doc/1320661/
document.pdf. https://doi.org/10.14459/2016md1320661 (2016)

399Comput Geosci (2022) 26:381–400

https://onlinelibrary.wiley.com/doi/abs/10.1002/0470091355.ecm063
https://onlinelibrary.wiley.com/doi/abs/10.1002/0470091355.ecm063
https://mediatum.ub.tum.de/doc/1320661/document.pdf
https://mediatum.ub.tum.de/doc/1320661/document.pdf
https://doi.org/10.14459/2016md1320661


52. Vinci, C., Renner, J., Steeb, H.: A Hybriddimensional
Approach for an Efficient Numerical Modeling of the Hydro-
Mechanics of Fractures. Water Resour. Res. 50.2, 1616–1635
(2014)

53. Vinci, C., Steeb, H., Renner, J.: The Imprint of Hydro-Mechanics
of Fractures in Periodic Pumping Tests. Geophys. J. Int. 202.3,
1613–1626 (2015)

54. Wang, H.F.: Theory of linear poroelasticity. Princeton University
Press Princeton & Oxford (2000)

55. Widlund, O., Toselli, A.: Domain Decomposition Methods - Algo-
rithms and Theory. English (US). Computational Mathematics,
vol. 34. Springer (2004)

56. Witherspoon, P.A., Wang, J.S.Y., Iwai, K., Gale, J.E.: Validity of
Cubic Law for Fluid Ow in a Deformable Rock Fracture. Water
Resour. Res. 16.6, 1016–1024 (1980)

Publisher’s note Springer Nature remains neutral with regard to
jurisdictional claims in published maps and institutional affiliations.

400 Comput Geosci (2022) 26:381–400


	Simulation of flow in deformable fractures using a quasi-Newton based partitioned coupling approach
	Abstract
	Introduction
	Governing equations
	Flow in a deformable fracture Fr
	Balance of Momentum
	Balance of Mass
	Governing Equations

	Partial Differential Equations in the Poro-Elastic Domain BPe
	Governing Equations

	Coupling conditions
	Weak formulation of governing equations
	Deformable fracture Fr
	Poro-Elastic Domain BPe
	Fluid exchange boundary condition
	Treatment of the Poro-Elastic Response on Different Time Scales


	Partitioned coupling
	Iterative Coupling

	Numerical results
	Verification of the partitioned implementation
	Investigation of interface quasi-Newton schemes
	Mesh convergence study

	Injection and production in an arbitrarily fractured reservoir
	Flow through fractured porous media

	Conclusion and outlook
	Declarations
	References