Journal of Geophysical Research: Solid Earth
Phase field model of fluid-driven fracture in elastic media:
Immersed-fracture formulation and validation
with analytical solutions
David Santillán1 , Ruben Juanes2,3 , and Luis Cueto-Felgueroso1,2
1Department of Civil Engineering: Hydraulics, Energy and Environment, Technical University of Madrid, Madrid, Spain,
2Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts,
USA, 3Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge,
Massachusetts, USA
Abstract Propagation of fluid-driven fractures plays an important role in natural and engineering
processes, including transport of magma in the lithosphere, geologic sequestration of carbon dioxide, and
oil and gas recovery from low-permeability formations, among many others. The simulation of fracture
propagation poses a computational challenge as a result of the complex physics of fracture and the need
to capture disparate length scales. Phase field models represent fractures as a diffuse interface and enjoy
the advantage that fracture nucleation, propagation, branching, or twisting can be simulated without ad
hoc computational strategies like remeshing or local enrichment of the solution space. Here we propose
a new quasi-static phase field formulation for modeling fluid-driven fracturing in elastic media at small
strains. The approach fully couples the fluid flow in the fracture (described via the Reynolds lubrication
approximation) and the deformation of the surrounding medium. The flow is solved on a lower
dimensionality mesh immersed in the elastic medium. This approach leads to accurate coupling of both
physics. We assessed the performance of the model extensively by comparing results for the evolution of
fracture length, aperture, and fracture fluid pressure against analytical solutions under different fracture
propagation regimes. The excellent performance of the numerical model in all regimes builds confidence
in the applicability of phase field approaches to simulate fluid-driven fracture.
1. Introduction
Fluid-driven fracturing plays an important role in natural and engineering processes. Some examples are the
transport ofmagma in the lithosphere [SpenceandTurcotte, 1985;Rubin, 1995; Parmigiani et al., 2016], the geo-
logic sequestration of carbon dioxide [Cappa and Rutqvist, 2011; Jha and Juanes, 2014], the preconditioning
in rockmasses to control the timing of goaf events [Jeffrey andMills, 2000], the creation of chemically reactive
barriers to inhibit the contaminantmigration in the vadose zone [Murdoch, 2002], the stimulation of geother-
mal reservoirs to increase heat extraction [Legarth et al., 2005; Evans et al., 2005], or the enhanced oil and gas
recovery from unconventional low-permeability reservoirs [Economides and Nolte, 2000; Patzek et al., 2013;
Cueto-Felgueroso and Juanes, 2013], among many others. Since technically recoverable gas resources world-
wide from unconventional reservoirs are estimated at approximately 8000 trillion cubic feet [Kuuskraa et al.,
2013], the last engineering application has a growing interest which stems from the environmental steward-
ship. The created fractures may unintentionally provide access between the reservoir and the surrounding
environment, leading to leakage of fracturing fluid or gas with the potential of groundwater contamina-
tion [Osborn et al., 2011; Vidic et al., 2013]. Both perspectives of the engineering and environmental concerns
require rigorous mathematical models and numerical simulation tools capable of reproducing the complex
physics involved in the process of hydraulic fracturing.
Fluid-driven fracture propagation in elastic media can be simulated through analytical solutions for simple
fracture geometries and injection protocols. The most widely used are the Khristianovic-Geertsma-de Klerk
(KGD)model, which assumes plane strain conditions [GeertsmaandDe Klerk, 1969; Khristianovich and Zheltov,
1955], the Perkins-Kern-Nordgren model, which considers an elliptic-shaped cross-section fracture of con-
stant height underplane straindeformation [PerkinsandKern, 1961;Nordgren, 1972], and thepenny-shapedor
RESEARCH ARTICLE
10.1002/2016JB013572
Key Points:
• We propose a new phase field
approach for fluid-driven fracturing
in elastic media
• The approach fully couples the elastic
deformation of the medium, fracture
apertures, and fluid flow inside the
fractures
• We find excellent agreement between
the proposed model and analytical
solutions of fracture propagation
Correspondence to:
D. Santillán,
david.santillan@upm.es
Citation:
Santillán, D., R. Juanes, and
L. Cueto-Felgueroso (2017), Phase field
model of fluid-driven fracture in elastic
media: Immersed-fracture formulation
and validation with analytical
solutions, J. Geophys. Res.
Solid Earth, 122, 2565–2589,
doi:10.1002/2016JB013572.
Received 20 SEP 2016
Accepted 22 MAR 2017
Accepted article online 27 MAR 2017
Published online 20 APR 2017
Corrected 19 MAY 2017
This article was corrected on 19 MAY
2017. See the end of the full text for
details.
©2017. American Geophysical Union.
All Rights Reserved.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2565
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
radial model, where the fracture propagates symmetrically with respect to the injection well and perpendic-
ular to it. Most analytical formulations are solved in an infinite linear elastic domain, where an incompressible
Newtonian fluid is injected at constant volumetric rate. Detournay [2016] gives an overview and current state
of the art of analytical solutions for the asymptotic behavior of the tip region and the propagation of a
penny-shaped fracture for a constant injection rate.
The propagation of a plane strain fluid-driven fracture in an impermeable elastic medium is controlled by the
competition between two dissipation processes: the viscous dissipation due to flowof the fluid in the fracture
and the toughness dissipation due to fracturing of the medium [Detournay, 2004]. The relative importance
of these dissipation processes leads to two limiting propagation regimes, namely, the viscous-dominated
regime, whose analytical solutions are provided in Adachi andDetournay [2002] and Garagash andDetournay
[2005], and the toughness-dominated regime, whose solution is given in Garagash [2006]. If the medium
is permeable, the storage of fluid must also be considered, and two new scenarios arise: the fluid may be
mostly concentrated in the porous medium—leak-off regime—or in the fracture—storage regime [Hu and
Garagash, 2010]. Analytical solutions with leak off (which assume that the fluid in the medium does not
affect the deformation field) exist for the toughness-dominated regime [Bunger et al., 2005] and for the
viscous-dominated regime [Adachi and Detournay, 2008], among many others available in the literature.
Numerical models of fracture can be classified into two categories: discrete and continuum approaches. The
discrete approach simulates fractures as discontinuities. From a numerical point of view, fractures are prop-
agated by splitting nodes [Fu et al., 2013] or breaking elements [Wangen, 2011] in the case of finite element
models or by splitting nodes and reconnecting springs in the case of spring network models [Hafver et al.,
2014]. Two drawbacks of these approaches are that the discretization must change topology due to fracture
growth and that the fracturepropagation is restricted to followmesh lines. Thesedisadvantages are overcome
by either remeshing techniques [Bouchard et al., 2003] or using advances approaches, such as the cohesive
zone modeling [Carrier and Granet, 2012] or the enriching displacement method [Oliver et al., 2006]. Other
approaches are those based on the boundary element method [Rungamornrat et al., 2005]. They are com-
putationally efficient since the solutions are only solved on the boundaries and the fracture surfaces. The
magnitude and direction of fracture propagation are governed by growth laws related to local stress intensity
factors. The main disadvantages of this method are the difficulties in handling nonlinear or heterogeneous
materials, or topological changes in the fracture trajectories such as joining or branching.
On the other hand, continuous approaches consider the intact and fractured areas as a whole, without the
need of introducing discontinuities. Examples of continuous models include peridynamics, gradient damage
models, or phase field models. Peridynamics is a nonlocal theory in which the solid is assumed to be com-
posed of material points. Each point interacts with all its neighbors within a nonlocal region called horizon
[Ouchi et al., 2015]. Gradient damage models and phase field models exhibit many similarities although dif-
ferences arise in the treatment of the strain localization procedure or the length scale over which damage
spreads [Shojaei et al., 2014]. In phase fieldmodels, fractures are described by incorporating a continuumfield
variable—the phase field—which interpolates between the broken and unbroken regions. The change from
intact to broken regions takes place along a transition band of width controlled by parameters. The main
advantages of this approach are that (i) all the method calculations are conducted on the initial undeformed
topology; (ii) themethod has the ability to simulate complex fracturing processes, such as branching, joining,
propagation or nucleation,without the need for additional criteria; and (iii) heterogeneousmedia are handled
without any additional rule.
Phase fieldmodels of brittle fracture can be classified into two families, although in practical implementations
both families often yield similar numerical schemes [Ambati et al., 2015]. The first family originated within the
physics community on the basis of the work of Aranson et al. [2000], who proposed a continuum field model
formode I propagation inbrittle amorphous solids inspired in theoriginal phase fieldmodels for solidification.
This idea was the seed for further developments. For instance, Karma et al. [2001] proposed an approach for
simulating the dynamic fracture propagation inmode III, andHenry and Levine [2013] developed a phase field
formulation for modeling fracture growth under modes I and III.
The second family emerged in the computational mechanics community, inspired in the work of Ambrosio
and Tortorelli [1990] which proposed a phase field approximation of the Mumford-Shah potential [Mumford
and Shah, 1989] based on Γ convergence. Motivated on these works, Francfort andMarigo [1998] developed
a variational formulation for quasi-static fracture evolution in a brittle material based on the minimization of
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2566
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
the combined elastic energy in the bulkmaterial and the fracture energy. Later, Bourdin et al. [2000] proposed
a numerical implementation of the variational formulation. They defined the phase field as a variable which
distributes the fracture energy over the bulk material. Many other authors have contributed to advancing
this approach.
Analternativequasi-static formulation for brittle fracturebasedoncontinuummechanics and thermodynami-
cal principleswasdevelopedbyMieheetal. [2010a]. The formulationwas later extendedby the introductionof
the so-called localhistoryfield, whichensures that thedamaged regiondoesnotbecomeunbrokenevenwhen
stress disappears (that is the irreversibility of the fracture process) combinedwith an operator split scheme that
successively updates the different fields involved in the formulation and permits simplifying the numerical
implementation of the problem while providing a robust algorithm [Miehe et al., 2010b]. The framework was
also extended to study dynamic problems [Hofacker andMiehe, 2012, 2013]. It also has been coupled to ther-
momechanical problems at large strains [Miehe et al., 2015a] or adapted to simulate ductile fracture coupled
with thermoplasticity at finite strains [Miehe et al., 2015b].
Phase field models have also been extended to model fluid-hydraulic brittle fracturing. An early contribution
was made by Bourdin et al. [2012], who proposed a model for simulating the propagation of a fracture in a
linear elastic impermeable solid. Fracture propagation was limited to the toughness-dominated regime. This
model was employed to study the effect of in situ stresses over the fracture propagation path [Chukwudozie
et al., 2013]. Later,Mikelic et al. [2015a, 2015b] extended the approach to poroelastic media at small strains. In
that work, both the flow equation in the medium and the fracture had the same dimensions and were gov-
erned by Darcy’s law, and a permeability tensor was imposed in the fracture in order to account for the higher
permeability along the crack. Miehe et al. [2015c] and Miehe and Mauthe [2016] proposed a new approach
for fracture propagation in a poroelastic medium at finite strains. The phase field evolution was driven by
the effective stress in the solid skeleton. Moreover, the flow in the fractures was governed by Poiseuille flow
modeled through a transition rule for Darcy flow combined with an anisotropic permeability tensor.
Although thepreviousmethods lead to straightforward computational frameworks, the storageof fluidwithin
the fracture is approximated by the storage within the solid. The smearing effect introduced by the phase
field leads to an artificial larger strain in the transition region and corresponds an artificial large storage fluid
volume. Whereas this approximation may be acceptable in poroelastic solids in which the leak-off rate from
the fracture to the solid may be of the same order of magnitude as the injection rate, it can result in large
errors when applied to impermeable elastic solids.
Here we propose a new quasi-static formulation for modeling fluid-driven fracturing in elastic media at small
strains, which describes fluid flow in the fracture and its coupling with the elastic deformation of the sur-
rounding medium. We simulate the flow in the fracture through the Reynolds lubrication equation. The fluid
pressure in the fracture, the aperture, and the deformation of themediumare solved in a fully coupled fashion
employing two separate meshes which are coincident in space. This strategy allows for a careful description
of the storage within the fracture. To assess the performance of phase field models for hydraulic fracturing
simulation, we compare the results of the evolution of the fracture length, aperture, and fluid pressure against
analytical solutions for several well-understood fracture propagation regimes.
The paper is organized as follows. Section 2 provides the governing equations for the fluid flow in the fracture
and the rockmechanicswith damage. A fewnotes on the numerical implementation of themodel are given at
the end of the section. A detailed description of the analytical benchmarkmodels is included in section 3 and
Appendix A. The case studywhere the numerical and analytical approaches are applied is described in section 4,
together the results and analysis of the simulations. Finally, some conclusions are provided in section 5.
2. Governing Equations
In this section, we begin by describing the geometrical representation of the fracture/matrix domains,
together with the fracture flowmodel, mechanical equilibrium equations and stress-strain relations. We then
derive our damage propagation framework from variational principles and exploit this approach to introduce
the phase field methodology for the simulation of quasi-static, brittle, fluid-driven fracture propagation in
elastic materials. Finally, we summarize the analytical solutions to the problem that will be used to evaluate
the numerical model results.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2567
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 1. Schematic of the fluid-pressurized fracture domain, ΩF , and the elastic medium domain, ΩR.
2.1. Geometry
Consider a domainΩ ⊂ ℜ𝛿 with spatial dimension 𝛿 ∈ {2, 3}, as depicted in Figure 1. Themedium comprises
an impermeable elastic subdomain,ΩR, and a pressurized fracture subdomain,ΩF ; i.e.,Ω = ΩR∪ΩF . A slightly
compressibleNewtonianfluid is injected inside the fracture at constant rateQ.While nomass exchangeoccurs
between subdomains, the fluid pressure exerts a force over the elastic subdomain that can propagate the
fracture. The boundary ofΩi is denoted by 𝜕Ωi, where i = R, F. The boundary ofΩR is divided into two subsets,
𝜕DΩR and 𝜕NΩR, where Dirichlet and Neumann boundary conditions are imposed, respectively.
2.2. Fluid Flow in the Fracture
The fracture subdomainΩF is represented as a lower dimensional manifold, and fluid flowwithin the fracture
is described through lubrication theory. The fluid motion is therefore governed by the Reynolds equation:
𝜕w
𝜕t
+ Cfw
𝜕pf
𝜕t
= ∇s ⋅
(
K(w)
(
∇spf − 𝜌f g∇sz
))
+ q, (1)
wherew is the fracture aperture, t is time, Cf is the compressibility of the fluid, s is the longitudinal coordinates
along the fracture, pf is the fluid pressure, 𝜌f is the fluid density, g is the gravity acceleration, z is depth, q
is the source or sink flow rate, and K(w) is the fracture mobility, given by the local cubic law [Witherspoon
et al., 1980]:
K(w) = w
3
12𝜇f
, (2)
where 𝜇f is the fluid dynamic viscosity.
2.3. Variational Approach to Rock Damage
We adopt the Griffith interpretation to model quasi-static brittle fracture propagation, which states that the
elastic energy release due to the propagation of the fracture is balanced by the newly created surface energy
[Griffith, 1920; Francfort and Marigo, 1998]. The total potential energy Ψ involved in the process of fracture
propagation has three contributions: the energy dissipated in the fracture process, Ψd ; the energy stored in
the bulk of the solid,Ψe; and the external sources of energy,Ψs,
Ψ = Ψd + Ψe − Ψs. (3)
The energy dissipated in the fracture process,Ψd , is the work required to create a unit fracture area, i.e.,
Ψd = ∫
𝜕ΩF
gc d𝜕, (4)
where gc is the Griffith critical energy release rate for mode I failure. This formulation requires knowing the
fracture surface, which is unknown a priori. In the phase field approach, this shortcoming is remedied by the
regularization of the fracture surface, i.e., by using adiffuse fracture surface definedwith a phase field variable,
d. This variable interpolates between the broken state (d = 1) and the intact one (d = 0). The fracture surface
is then defined through and integral extended to the whole domain [Miehe et al., 2010a]:
Γ𝓁 = ∫Ω 𝛾𝓁 dΩ, (5)
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2568
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
where Γ𝓁 is the regularized fracture functional and 𝛾𝓁 is the fracture density given by
𝛾𝓁 =
1
2𝓁
d2 + 𝓁
2
|∇d|2 , (6)
and 𝓁 is a length scale parameter. The dissipated energy,Ψd , reads
Ψd = ∫Ω gc𝛾𝓁 dΩ, (7)
which allows for the energy dissipated during the fracture process to be computed without prior knowledge
of the discrete fracture surface.
The energy stored in the bulk of the solid,Ψe, evolves due to stress redistribution during deformation and due
to the fracture process. The loss of stiffness of the solid due to fracture growth implies that the solid cannot
store the same amount of elastic energy, and it is therefore dissipated. From amathematical point of view, an
energy transfer fromΨe toΨd occurs. The transfer, or degradation of the elastic energy, is computed through
the introduction of a degradation function, g, that depends on the phase field variable, d. While other options
have been proposed, in this study we employ a quadratic degradation function [Vignollet et al., 2014], since
the impact of the choice of the function diminishes after the fracture has formed [Kuhn et al., 2015], and the
Γ convergence has so far only been proved for the quadratic polynomial [Chambolle, 2004]. The function has
the following form:
g(d) = (1 − d)2, (8)
which assumes that the phase field variable is directly related to fracture growth.
The energy stored in the undamaged bulk of the solidΨe0 is given by:
Ψe0 = ∫Ω 𝜓
e
0 (𝜀)dΩ, (9)
where 𝜓e0 is the undamaged elastic energy density,
𝜓e0 (𝜀) =
𝜆
2
tr(𝜀)2 + 𝜇tr(𝜀2), (10)
where 𝜆 and 𝜇 are the Lamé constants and tr(𝜀) is the trace of the strain tensor. The degradation of the elas-
tic energy density can be modeled mathematically in two ways. The isotropic damage formulation assumes
degradation of the whole elastic energy as the fracture propagates [Miehe et al., 2010b]:
Ψe = ∫Ω g(d)𝜓
e
0 (𝜀)dΩ = ∫Ω 𝜓
e(𝜀, d)dΩ. (11)
The corresponding Cauchy stress tensor incorporating damage is therefore given by:
𝜎(u, d) ∶= 𝜕𝜓
e
𝜕𝜀
= g(d)
𝜕𝜓e0
𝜕𝜀
. (12)
The above isotropic approach leads to a symmetric response of the model for fracturing under tension and
compression, which may result in unrealistic fracture evolution patterns for rocks and other materials. An
alternativemethodology is the so-calledanisotropicdamage formulation [Mieheetal., 2010b],whichallows for
fracturing in tension only. The formulation is suitable for hydraulic fracturing under mode I (tensile opening),
and itwould need tobe revisited for fracture propagationundermode II (in-plane shear) ormode III (antiplane
shear). The undamaged elastic energy density function, 𝜓e0 , is split into a positive part, 𝜓
e+
0 , due to tension,
and a negative part, 𝜓e−0 , due to compression, as follows:
𝜓e0 (𝜀) = 𝜓
e+
0 (𝜀) + 𝜓
e−
0 (𝜀), (13)
where tension/compression components are defined as follows:
𝜓e±0 (𝜀) =
𝜆
2
(⟨
𝛿∑
i=1
𝜀i
⟩
±
)2
+ 𝜇
𝛿∑
i=1
(⟨𝜀i⟩±)2 , (14)
where ± refers to the tension/compression parts of the undamaged elastic energy, 𝛿 is the number of space
dimensions, and ⟨x⟩± = (x ± |x|) ∕2.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2569
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
This approach is based on the spectral decomposition of the strain tensor:
𝜀±(u) =
𝛿∑
a=1
⟨𝜀a⟩± na ⊗ na, (15)
where na is the principal strain direction associated to the principal strain 𝜀a. The eigenvalue bases na ⊗ na
are computed in terms of the principal strains 𝜀i and the strain tensor 𝜀 as follows [Miehe, 1998]:
na ⊗ na =
1∏𝛿
b≠a
(
𝜀a − 𝜀b
) 𝛿∏
b≠a
(
𝜀 − 𝜀b1
)
. (16)
To consider degradation of the tensile energy alone, the damaged elastic energyΨe is given by
Ψe = ∫Ω
[
g(d)𝜓e+0 (𝜀) + 𝜓
e−
0 (𝜀)
]
dΩ = ∫Ω 𝜓
e(𝜀,d)dΩ, (17)
and the resulting stress-strain constitutive relation yields:
𝜎(u, d) ∶= g(d)
𝜕𝜓e+0
𝜕𝜀
+
𝜕𝜓e−0
𝜕𝜀
. (18)
Alternatively,
𝜎+(u, d) = g(d)
𝛿∑
a=1
[
𝜆
⟨
𝛿∑
i=1
𝜀i
⟩
+
+ 2𝜇 ⟨𝜀a⟩+
]
na ⊗ na (19)
and
𝜎−0 (u) =
𝛿∑
a=1
[
𝜆
⟨
𝛿∑
i=1
𝜀i
⟩
−
+ 2𝜇 ⟨𝜀a⟩−
]
na ⊗ na. (20)
The external energy functional, Ψs, accounts for the exchanged energy between the elastic domain ΩR and
the surrounding environment. It is given by
Ψs = ∫Ω f̄ ⋅ udΩ + ∫𝜕NΩ t̄ ⋅ ud𝜕, (21)
where f̄ is the body force per unit volume, t̄ is the vector of applied forces, and u is the displacement field. The
second term on the right-hand side of equation (21) includes the force introduced by the pressure inside the
fracture, which is applied on a surface that is a priori unknown:
∫
𝜕NΩF
t̄ ⋅ ud𝜕 = −∫
𝜕NΩF
pf n̄ ⋅ ud𝜕, (22)
where n̄ is the normal unit vector to the fracture surface. Equation (22) can be interpreted as a flux of energy
through the surface of the fracture. Using the divergence theorem, we rewrite equation (22) as
∫
𝜕NΩF
pf n̄ ⋅ ud𝜕 = ∫ΩR ∇ ⋅
(
pfu
)
dΩ − ∫
𝜕NΩR
pf n̄ ⋅ ud𝜕, (23)
where the pressure exerted on the fracture is extended as a function defined over the entire domain. The
practical implementation of this fictitious pressure approach is discussed in section 2.4 below. As formulated,
the integral (23) is defined over the elastic medium, whose domain changes with time due to the fracture
propagation. To avoid this difficulty, we introduce the degradation function, g(d), and redefine the integral
over the entire domain [Mikelic et al., 2015b]:
∫ΩR ∇ ⋅
(
pfu
)
dΩ = ∫Ω g(d)∇ ⋅
(
pfu
)
dΩ. (24)
With the above definitions, the external energy functional finally reads:
Ψs = ∫Ω f̄ ⋅ udΩ − ∫Ω g(d)∇ ⋅
(
pfu
)
dΩ + ∫
𝜕NΩR
pf n̄ ⋅ ud𝜕, (25)
where we have assumed that there are no tractions applied at the domain outer boundary.
We derive the strong form of the problem via minimization of the total potential energy Ψ with respect to
the displacement and phase fields, u and d. Herein, we have adopted the anisotropic formulation for the
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2570
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
degradation of the elastic energy, equation (17). The Fréchet derivative of Ψ with respect to u provides the
elasticity, or equilibrium, equations:
𝜕Ψ
𝜕u
− ∇ ⋅
(
𝜕Ψ
𝜕∇u
)
= 𝜕Ψ
𝜕u
− ∇ ⋅
(
𝜕Ψ
𝜕𝜀
𝜕𝜀
𝜕∇u
)
= 0, in Ω, (26)
𝜕Ψ
𝜕∇u
n̄ = t̄ + pf n̄, in 𝜕NΩR. (27)
In terms of the positive and negative stress contributions, the equilibrium equations for the damage model
are given by:
− ∇ ⋅
{
g(d)𝜎+0 + 𝜎
−
0
}
− pf∇g(d) − f̄ = 0, in Ω, (28)
𝜎n̄ = t̄, in 𝜕NΩ. (29)
The Fréchet derivative ofΨwith respect to d,
𝜕Ψ
𝜕d
− ∇ ⋅
(
𝜕Ψ
𝜕∇d
)
= 0, in Ω, (30)
𝜕Ψ
𝜕∇d
⋅ n̄ = 0, in 𝜕NΩ (31)
provides the Euler equations of the phase field problem:
gc
𝓁
(
d − 𝓁2∇2d
)
= 2(1 − d)
(
𝜓e+0 (𝜀) + pf∇ ⋅ u + u ⋅ ∇pf
)
, in Ω, (32)
∇d ⋅ n̄ = 0, in 𝜕NΩ. (33)
The total energy density in equation (32),𝜓e+0 (𝜀)+pf∇ ⋅u+u ⋅∇pf , determines the amount of phase field vari-
able d at the current time. The above equation does not account for the deformation history; thus, d becomes
zero when stress disappears. Introducing the maximum local history field H+ permits accounting for the irre-
versibility of the process in a straightforward way [Miehe et al., 2010b]. H+ may be considered as a measure of
the maximum historic tensile strain, defined as:
H+(u, pf , t) = maxs∈[0,t]
(
𝜓e+0 (𝜀) + pf∇ ⋅ u + u ⋅ ∇pf
)
. (34)
The Euler equation of the phase field is then written as follows:
gc
𝓁
(
d − 𝓁2∇2d
)
= 2(1 − d)H+. (35)
2.4. Numerical Implementation of the Model
The system of governing equations comprises three partial differential equations (PDEs), modeling fluid flow
inside the fracture, mechanical equilibrium, and evolution of the phase field variable. Neglecting gravity, the
Reynolds lubrication equation couples the evolution of fluid pressures and fracture apertures:
𝜕w
𝜕t
+ Cfw
𝜕pf
𝜕t
= ∇s ⋅
(
K(w)
(
∇spf
))
+ q. (36)
The second model equation is the mechanical equilibrium equation,
− ∇ ⋅ 𝜎 − pf∇g(d) = 0, (37)
where 𝜎 is the damaged stress tensor, computed by equation (18). Finally, the third equation is the phase field
equation,
gc
𝓁
(
d − 𝓁2∇2d
)
= 2(1 − d)H+, (38)
where H+ is given by equation (34).
Equations (36)–(38) constitute a coupled nonlinear systemof equations. We adopt a finite volume scheme on
a fixed regular quadrilateral mesh for the spatial discretization of the system of governing equations. As the
fracture flow equation, equation (36), has a lower dimensionality in space, it is solved on a compatible regular
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2571
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 2. Propagation of the fracture. Example of contour plot of
phase field variable at time tn
one-dimensional mesh embedded in
the 2-D domain. The fracture domain is
defined as a lower dimensionality entity
where the gradient of the phase field
is zero and the phase field is above a
threshold df : we adopt df = 0.90.
We explain the strategy through Figure 2,
which draws a surface plot of d at time tn.
Themesh of the problem and the integra-
tion points are depicted in black. The frac-
ture domain at time tn−1 is represented
by the solid white line. To determine the
fracture domain at time tn, we first check if
d is larger than df in the neighboring cells
of the tip. If this is the case, the fracture
propagates following the line that starts
at the previous tip, follows the direction
∇d = 0, and passes through the integration points where d> df . In Figure 2 the fracture propagates one cell,
up to the point in the upper right diagonal of the tip.
We discretize equation (36) in time using a backward Euler scheme and adopt a simple two-point flux
approximation in space. The aperture at point x along the fracture is evaluated through the line integral
[Mikelic et al., 2015b]
w(x) = ∫
b
a
u(r(𝜂)) ⋅ ∇d(r(𝜂)) ‖‖r′(𝜂)‖‖d𝜂, (39)
where r(𝜂) is the line normal to the fracture that passes through x. The line integral (39) has been truncated
to the interval [r(a), r(b)], whose limits denote points that are far enough from the fracture so that ∇d(r(a))
and∇d(r(b)) are negligible.
The equilibrium equation, equation (37), includes a nonlinear stress-strain relation due to the spectral decom-
position of the strain tensor, equation (15), and the presence of the phase field variable through the degra-
dation function, equation (18). This source of nonlinearity can be handled effectively by using the so-called
hybrid formulation for the stress-strain relationship [Ambati et al., 2015], combined with a staggered scheme
for updating the displacement field, u, and the phase field, d [Miehe et al., 2010b].
The hybrid formulation for the spectral decomposition of the strain tensor contains features of the isotropic
and the anisotropic formulations. The idea behind this formulation is to retain a linear equilibrium equation
within a staggered scheme. The formulation requires the stress-strain relationship (12), combined with the
usual phase field equation (38), where the evolution of d is driven by the tensile energy evolution Ψe+0 only.
To avoid the interpenetration of rock masses on each side of the fracture, equation (38) is subjected to
the constraint:
∀ x ∶ Ψe+0 < Ψ
e−
0 ⇒ d ∶= 0. (40)
The above procedure switches the value of the phase field variable to zero (undamaged solid) in those areas
subject to compression. When the tensile part of the energy is again larger than the compressive one, the
phase field variable recovers its previous value (damaged solid) and is again driven by H+.
Ambati et al. [2015] provide numerical evidence that the hybrid formulation leads to physical results for frac-
ture evolution that are qualitatively and quantitatively similar to those of the anisotropic formulation, while
requiring a much lower computational effort, in fact, comparable to that of the isotropic model.
Fluid pressures along the fracture contribute to themechanical equilibrium equation, equation (37), through
the body force pf∇g(d). Since the degradation function g(d) decays quickly away from the fracture zone, this
term localizes theeffect of thepressure at the fracture andmodels the force exertedby thefluidpressureof the
fracture walls. The conversion from surface to volume integrals requires, however, the ad hoc reconstruction
of a pressure field extended to the full domain, not just at the fracture. In the present case of an impervious
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2572
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
medium, such pressure is a “fictitious” one, definedwith the sole purpose of converting the impact of surface
forces into body forces. In the present study, points inside the elastic domain are assigned the corresponding
pressure of the closest point at the fracture. The degradation functiong(d)has the effect of localizing the pres-
sure contribution, and therefore, the impact of different choices on the definition of the extended pressure is
small as the fracture width decreases. We present a validation of the present approach in section 3.1.
We adopt the staggered scheme proposed byMiehe et al. [2010b] for computing (u, d), which is a simple and
robust alternative to themonolithic scheme. When combined with the hybrid formulation, it provides results
of comparable quality to thoseof the anisotropic formulationwith amonolithic solution scheme [Ambati etal.,
2015]. However, the scheme requires sufficiently small loading increments to properly capture the nonlinear-
ity and strong coupling of the problem. A single step of the staggered scheme in the time interval [tn, tn+1]
comprises three substeps:
1. Assuming that the displacements, u, fracture pressures, pf , phase field, d, andmaximumhistorical energies,
H+, are known at time tn, update H
+ with the values un.
2. With the udpdated energy field, H+, update the phase field variable, d.
3. Finally, advance the displacement and pressure fields u and pf , using a fully coupled approach with frozen
phase field d.
3. Benchmark Analytical Models
We assess the performance of our model by comparing numerical results with known analytical solu-
tions describing asymptotic propagation regimes. To validate the fictitious pressure approach adopted here
(equations (23) and (24)), we first consider the static case of a fracture subjected to a constant pressure force
that is small enough to avoid growth. In the absence of fracture growth, the comparison focuses on the effect
of the fracture pressure on the elastic deformations through the conversion of surface to volume forces.
We subsequently consider three fluid-driven fracture propagation cases. The analytical solutions are obtained
for an infinite elastic medium under plane strain conditions: the KGD fracture (Figure 3). A far-field compres-
sion 𝜎0 acts in the direction perpendicular to the fracture, and an incompressible Newtonian fluid fills the
fracture at constant pressure in the static case or it is injected at constant volumetric rate Q at the center of
the fracture in the dynamic cases.
Further assumptions underlying the fluid-driven fracture propagation models include the following:
1. The fracture is fully filled with fluid at all times.
2. The fracture is inmobile equilibrium, and its quasi-static propagation is described in the framework of linear
elastic fracture mechanics.
3. The flow of the fluid is unidirectional and laminar, and it is governed by the Reynolds lubrication equation.
Sections3.1 and3.2 summarize theanalytical results applicable to thedifferent tests andpropagation regimes.
These analytical solutions are compared with numerical ones in section 4.
3.1. Static Case
The first benchmark is the analytical solution of Sneddon [1946], who considered an infinite 2-D elastic domain
with a fracture of constant length 2l filled with a fluid at pressure pf . The fracture aperture profile, w(x),
is given by:
w(x) = 4pl
E′
√
1 −
(x
l
)2
, (41)
where p = pf −𝜎0 is the net pressure, E′ = E∕(1−𝜈2), and x is the position along the fracture. This fairly simple
model allows us to validate equation (23), where pressure forces are accounted for as body forces using the
phase field. Moreover, we test the convergence and mesh independence of our model results by conducting
amesh refinement study. In particular, we study the effect of mesh size on the fracture aperture. Wemeasure
the discrepancy between the analytical and numerical apertures through the L2 error norm, given by:
L2 = ‖Vex − Vnum‖2 = { Ni∑
i=1
∫Ω
(
Vexi − V
num
i
)2
ds
}1∕2
, (42)
where Vexi is the exact value of the variable at location i and V
num
i is the numerical value.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2573
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 3. Schematic of a two-dimensional, fluid-driven fracture in an infinite elastic medium under plane strain
conditions. The fracture is described as a lower dimensional manifold, and fluid flow through the fracture is modeled
through the lubrication approximation, equation (1).
3.2. Regimes of Fracture Propagation
The propagation of a fluid-filled fracture in an elastic permeable medium is governed by two competing dis-
sipative processes related to fluid viscosity andmedium toughness and two competing storage mechanisms
related to storage in the fracture and in the medium [Hu and Garagash, 2010]. In the case of impermeable
media, the fracture is the only storagemechanism. We have chosen three analytical solutions, corresponding
to the following scenarios:
1. Propagation of a fracture in an impermeable medium where the energy expended in fracturing the rock is
much larger than the viscous dissipation. This regime is called toughness dominated.
2. Propagation of a fracture in an impermeable medium where the dissipation of energy in propagating the
fracture is negligible compared to the viscous dissipation. This regime is known as viscous dominated.
3. Propagation of a fracture in a permeablemediumwhere the storage of fluid in the porousmedium ismuch
larger than the storage in the fracture, and thedissipationof energy in propagating the fracture is alsomuch
larger than the viscous dissipation. This regime is referred to as leak-off toughness dominated.
The analytical solutions are included in Appendix A.
4. Numerical Results and Discussion
The geometry of the numerical model is depicted in Figure 4. The injection point is located at the center of
the domain, and the fracture propagates following theminimum stress path, i.e., the horizontal direction. The
dimensions of the domain are 90 m in the horizontal direction (x axis) and 60 m in the vertical direction
(y axis). The mechanical properties of the solid are listed in Table 1. We assume plane strain conditions, and
the injected fluid is incompressible and Newtonian.
Figure 4. Model setup and boundary conditions. The injection point is located at the center of the domain, and the
fracture propagates following the minimum stress path, i.e., the horizontal direction. Given that the problem is
symmetric, we model one quarter of the domain only. For the x axis, the boundary conditions are zero vertical
displacement, v = 0, and zero derivative of the horizontal displacement with respect to the vertical coordinate,
𝜕u∕𝜕y = 0. For the y axis, the boundary conditions are zero horizontal displacement, u = 0, and zero derivative of the
vertical displacement with respect to the horizontal coordinate, 𝜕v∕𝜕x = 0. The domain is discretized using a regular
square fixed mesh composed of 1.5 million cells of size 𝛿x = 𝛿y = 0.03 m.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2574
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Table 1. Mechanical Properties of the Elastic Medium
Young’s modulus E 17 GPa
Poisson’s ratio 𝜈 0.2
Griffith energy gc 120 Pa m
The symmetry of the problem allows us
to model one quarter of the domain
only (Figure 4). Hence, the mechanical
boundary conditions are zero displace-
ments except in the axes of symmetry.
For the x axis, the boundary conditions
are zero vertical displacement, v=0, and
zero derivative of the horizontal displacementwith respect to the vertical coordinate, 𝜕u∕𝜕y=0. For the y axis,
the boundary conditions are zero horizontal displacement, u=0, and zero derivative of the vertical displace-
ment with respect to the horizontal coordinate, 𝜕v∕𝜕x = 0. The domain is discretized using a regular square
fixed mesh composed of 1.5 million cells of size 𝛿x= 𝛿y=0.03 m. As the analytical solutions are obtained for
an infinite domain, we confirmed that the computational domain is sufficiently large, so that boundary effects
do not pollute the results.
4.1. Static Case
We compute the aperture of a fracture of half-length l = 2mfilledwith fluid at net pressure pf = 2.25×105 Pa
with several mesh sizes ranging from 0.100 to 0.017 m. The fracture is introduced by setting a large-enough
Figure 5. Static case: the Sneddon’s crack. We compute the aperture of one fracture of half-length l = 2 m filled with
fluid at pressure pf = 2.25 × 105 Pa and compare the results with the analytical solution of Sneddon [1946]. (a) We plot in
log-log scale the L2 error norm for the fracture aperture against the mesh size 𝛿x . (b) We plot the fracture profiles for the
analytical solution and the numerical model with mesh size 𝛿x = 𝛿y = 0.03 m, and we zoom the fracture tip. The phase
field model results agree well with the analytical solution except near the tips (1 − 𝜉 → 0), where the phase field model
introduces a smearing effect due to the regularization of the fracture surface. (c) Here we depict the deviation at the tip
Δwtip for several mesh sizes. We compute Δwtip as the difference between the fracture aperture at the tip computed
with both the analytical solution and the numerical model. (d) Setting the mesh to a size 𝛿x = 𝛿y = 0.03 m, as the value
of 𝓁 increases the L2 error norm for the fracture aperture decreases and Δwtip slightly increases.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2575
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
value of H+ in the damaged region that the phase field is close to 1 in the fracture. We set the length
scale parameter 𝓁 = 0.35 m. We also compute the aperture profiles with the Sneddon’s analytical solution
[Sneddon, 1946].
We plot the L2 error norm of the fracture aperture against the mesh size in log-log scale in Figure 5a. The
numerical results converge to the exact solution as the element size decreases. The aperture converges quar-
tically for the largest mesh sizes up to 𝛿x = 𝛿y = 0.043 m. For smaller sizes, the convergence order is lower
than one. The convergence order of the aperture differs from the order of the numerical scheme—equal to
2—since the aperture is a variable computed from the solution of the numericalmodel. These results suggest
that a good choice of the mesh size is 𝛿x = 𝛿y = 0.03 m. The accuracy of the numerical solution, in terms of
the root-mean-square error, is 1.3 × 10−5 for this size.
To illustrate typical profiles of fracture aperture,weplot the aperture computedwith the numericalmodel and
the analytical solution (Figure 5b). The mesh size is 𝛿x = 𝛿y = 0.03 m. The numerical model agrees well with
the analytical solution except near the fracture tips, where the phase field model does not match the exact
solution due to the smearing effect introduced by the regularization of the fracture, as can be appreciated in
the zoom of the tip.
To perform a quantitative assessment of the deviation near the tip owing to the diffuse modeling of the frac-
ture, we compute the difference between the exact and the numerical aperture at the tip for several mesh
sizes (Figure 5c). The deviation has the same order of magnitude as the L2 norm, and it converges to the exact
solution as the size decreases.
The distance occupied by the deviation—distance between the position of the exact tip and the point where
the numerical aperture is zero—is controlled by the length scale parameter 𝓁. The deviation is proportional
to the transition between the broken and undamaged regions, which increases with the value of 𝓁. For a
given mesh size, the accuracy of the numerical model also depends on 𝓁. The evolution of the L2 error for
the fracture aperture with the value of 𝓁 is depicted in Figure 5d. Larger values of 𝓁 provides more accurate
apertures, but the deviation slightly increases due to the smearing effect. In the following numerical models
we set 𝓁 = 0.35 m.
4.2. Fracture Propagation Cases
4.2.1. Toughness-Dominated Regime
To perform a simulation in the toughness-dominated regime, we inject fluid with dynamic viscosity
𝜇′
f
= 10−6 Pa s at a rate Q = 10−3 m2/s. The dimensionless viscosity is equal to 1.4 × 10−4, which is sig-
nificantly less than0 = 3.4 × 10−3, indicating that the fracture propagates in the desired regime. We start
the simulation with an initial fracture of length 2.2 m. We do not provide initial fracture aperture or fluid pres-
sure consistent with the analytical solutions since the net pressure depends on the far-field stress induced
by the boundaries, which are initially unknown. Some discrepancies at early times arise that after some time
steps disappear.
To illustrate typicalmodel results, we plot the phase field, vertical displacements, and principal tensile stresses
at time t = 10 s, at which the fracture has grown to a length of approximately 10 m (Figure 6). Since phase
fieldmodels regularize damage and define the fracture as a diffuse interface, a small transition exists between
the fractured (d = 1) and intact (d = 0) regions. This transition, whose width is controlled by the length scale
parameter 𝓁, is fully resolved by the computational grid (Figure 6a).
The vertical displacement field is shown in Figure 6b. The displacement is largest at the center of the frac-
ture, and it decreases toward the fracture tips. However, the computed displacements in the transition
region are artificially large as a result of the mathematical regularization of the fracture. Therefore, frac-
ture apertures are not directly provided by the displacement fields but, rather, require the evaluation of the
integral equation (39); that is, in the present diffuse interface model, fracture apertures are recovered from
the displacements as a nonlocal quantity. If we were to estimate apertures directly from the displacements
(Figure 6b), the aperture of the central point of the fracture is approximately 0.8mm. In contrast, the aperture
computed using equation (39) is smaller, approximately 0.6 mm.
Stress concentration occurs at both fracture tips, while inside the crack the stress is close to zero (Figure 6c).
This pattern is consistent with predictions from the classical theory of linear elastic fracture mechanics.
Stress is calculated through a stress-strain constitutive relation that involves the phase field variable. In turn,
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2576
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 6. Toughness-dominated regime: Contour plots of (a) phase field, (b) vertical displacement, and (c) principal
tensile stress, at time t = 10 s.
the evolution of the phase field is governed by the stress-strain constitutive relation over the undamaged
solid, which provides larger stresses than the ones shown in Figure 6c.
To perform a quantitative assessment of the computational model, we plot the evolution in time of the frac-
ture half-length l(t), fracture midpoint aperture w(x = 0, t), and net fluid pressure at the fracture midpoint,
p(x= 0, t) (Figures 7a–7c). The results from the numerical model agree well with the analytical solution pre-
sented in section A1. We limite the simulation to the point where the half-length was approximately 14 m
or about 30% of the half domain, to avoid boundary effects. The agreement between simulated fracture
apertures and analytical solution is excellent except at very early times, when the aperture computed by the
numerical model is smaller than the analytical one (Figure 7b). This effect is due to the time needed by the
numerical model to increase the fracture apertures to the point where fracture propagation starts. In that
sense, the deviation from the analytical results can be attributed to an onset time.We providemore details on
this issue in the next section, where we discuss the viscous-dominated case. The far-field stress 𝜎0 is constant
in time and space in the analytical model. However, it is not in the numerical approach due to the boundary
condition of zero normal displacements. We compute the net pressure as the difference between the fluid
pressure and the average vertical stresses along the horizontal boundaries. The resulting pressure is slightly
higher than the analytical solution, but both evolutions exhibit the same power law scaling except at the
beginning of the simulation (Figure 7c).
An important characteristic of the toughness-dominated regime is the pressure distribution along the frac-
ture. For a given time, the pressure is almost uniform in the fracture except at the fracture tip, where a
significant pressure drop occurs (Figure 8a). The first-order analytical solution provides a uniform value of the
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2577
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 7. Toughness-dominated regime: Evolution of (a) half-length of the fracture, (b) aperture at the center of the
fracture, and (c) net pressure at the center of the fracture.
pressure in the fracture (Figure 8a). The numerical model reproduces the pattern except at the fracture tip
where a significant pressure drop occurs. The near-tip drop in pressure is owed to a dominance of viscous
effects near the tip, where the simulated aperture is lower than the analytical one (Figure 8b) and viscous
dissipation dominates.
The simulated pressure is slightly higher than the analytical solution, but the difference decreases as the frac-
ture grows. The error is presumably owing to the far-field stress induced by the boundaries whosemagnitude
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2578
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 8. Toughness-dominated regime. (a) We plot the net pressure along the fracture computed with the analytical and numerical models and the vertical
stress at the boundary at two different times. (b) Here we depict the dimensionless aperture against the dimensionless coordinate 𝜉 = x∕l in log-log scale at the
same two time steps.
becomes lower as the fracture grows and the spatial distribution is more uniform with time (Figure 8a). The
numerical solution underpredicts the near-tip fracture aperture, but the error decreases as the fracture grows
(Figure 8b).
4.2.2. Viscous-Dominated Regime
To perform a simulation in the viscous-dominated regime, the injection flow rate is Q = 2 × 10−4 m2/s and
the dynamic viscosity of the fluid is 𝜇′
f
= 0.1 Pa s. The initial fracture has a length of 2.2 m. Moreover, as the
net pressure depends on the far-field stress induced by the boundaries, we impose an initial fracture aperture
and fluid pressure lower than those from the analytical solution. After a few initial time steps, the numerical
solution converges to the analytical one.
As in the previous section, we plot the evolution of the fracture half-length, aperture at the central point
of the fracture, and fracture pressure at the center of the fracture (Figures 9a–9c). The numerical approach
exhibits a good agreement with the analytical solution, except for minor differences at the early stages of the
simulation. Initially, the pressure computed with the numerical model is lower than the one provided by the
analytical formulation. Once the pressure increases and the fracture starts to grow, the pressure predicted
by the numerical model agrees well with the analytical solution. Although the numerical solution slightly
underpredicts the fracture pressure, both solutions show the same power law scaling after the early stages of
the simulation (Figure 9c).
The reason for the discrepancy during the early stages arises from the initial value of the pressure, which is
not sufficiently high in the numerical simulation. In Figure 10 we plot the fracture profiles at early times and
show that during this period the fracture length does not increase. During that early period, the numerically
simulated fracture “inflates” and the fluid pressure increases up to the point where the fracture starts to prop-
agate. This onset timeduringwhich the fracture inflates prior to propagation explains, in part, the discrepancy
between numerical and analytical results at early times.
The pressure distribution inside the fracture is shown in Figure 11a at two different times. In contrast with
the toughness-dominated regime, the viscous-dominated regime is characterized by appreciable fluid pres-
sure loss within the fracture which even leads to negative values of pressure. The agreement between
our numerical model and the analytical solution is excellent except at the fracture tip, where the numer-
ical solution provides lower values of the pressure because the simulated near-tip fracture aperture is
smaller than the analytical one (Figure 11b). The far-field stress induced by the boundaries is smaller
than in the toughness-dominated regime, which leads to better simulations of fluid pressure and fracture
aperture.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2579
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 9. Viscous-dominated regime: Evolution of (a) half-length of the fracture, (b) aperture at the center of the
fracture, and (c) net pressure at the center of the fracture.
4.2.3. Leak-Off Toughness-Dominated Regime
To perform a simulation in the leak-off toughness-dominated regime, we inject a fluid with dynamic viscosity
𝜇′
f
= 10−6 Pa s at a rate Q = 10−2 m2/s. The leak-off coefficient is 5 × 10−4 m/s1∕2. The initial fracture is 2.2 m
length, and it is reached in the analytical solution at time 0.65 s. We impose consistent initial conditions as
the Carter’s model for leak-off flow rate depends on time. We calculate the initial values through an iterative
process. After 0.5 s, the dimensionless time is 𝜏 = 1.05, which means that the transition between the storage
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2580
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 10. Viscous-dominated regime: Fracture profiles at early times. Fluid injection induces the fracture opening, but
the fracture length remains constant.
and leak-off toughness-dominated regimes occurs at that moment. Since the analytical formulation indi-
cates that the initial fracture is approximately reached at time 0.65 s, the numerical simulation starts in the
leak-off regime.
The evolution of the fracture half-length, aperture at the central point of the fracture, and fracture pressure
at the center of the fracture in this regime is shown inFigures 12a–12c. Aswas the case in theprevious regimes,
the phase field model prediction of the evolution of these three variables agrees well with the analytical
solution. Small differences exist only during the very early stage of the simulation. It is interesting to compare
the behavior of fracture propagationbetween this regime and the toughness-dominated regimewithout leak
off. While the half-length at the end of the simulation is similar for both simulations, this is the result of simu-
lations which differ by an order of magnitude in the injection flow rate: 10−3 m2/s for the impermeable solid
and 10−2 m2/s. Despite the tenfold larger injection rate, the simulation with leak off results in both smaller
fracture aperture and fracture fluid pressure. It is important to note that the analytical model is developed for
an elastic solid inwhich the presence of a fluidwithin the pores does not alter the stress field, i.e., the Biot coef-
ficient is zero. This feature has also been considered in the numerical model, permitting a direct comparison
of the numerical results with the analytical solution.
4.3. Practical Remarks
Wehave simulated the propagation of fractures in homogeneous elastic solidswhere cracks followplanar tra-
jectories. Phasefieldmodels enjoy theadvantage that fracturenucleation, propagation, branching, or twisting
Figure 11. Viscous-dominated regime. (a) We plot the net pressure along the fracture computed with the analytical and
numerical models and the vertical stress at the boundary at two different times. (b) Here we depict the dimensionless
aperture against the dimensionless coordinate 𝜉 = x∕l in log-log scale at the same two time steps.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2581
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 12. Leak-off toughness-dominated regime: Evolution of (a) half-length of the fracture, (b) aperture at the center
of the fracture, and (c) net pressure at the center of the fracture.
can be simulated without ad hoc computational strategies. These capabilities have not been shown in the
simulations of the previous section.
The method is particularly well suited to simulate fracture propagation under the toughness-dominated
regime and when the injection pressure is prescribed. In this case, the pressure along the fracture is uniform
and equal to the injection pressure and the injection flow rate can be computed a posteriori. This scenario
is particularly suitable for determining fracture propagation patterns and constitutes a potential direction in
which our method could be used.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2582
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Figure 13. Toughness-dominated regime. (a) We plot one realization of 2-D lognormal Young’s modulus with a Gaussian
correlation function. 𝜇E is 17 GPa, 𝜎E is 6.8 GPa, and the correlation length is 0.05 m. (b) Here we depict the contour plot
of the phase field. The fracture propagates in an elastic solid whose Young’s modulus is depicted in the left panel.
We illustrate this with a simulation in which a fracture propagates under the toughness-dominated regime
in an elastic solid where both the Young’s modulus and the Griffith critical energy release rate vary spatially
following isotropic lognormal fields. Both fields are assumed to have the same Gaussian correlation function
𝜌(r2) where r is spatial distance. We assume that at any given location x, both the log-Griffith energy ln gc(x)
and the log-modulus field ln E(x) come from the same local value F(x), which is an isotropic normal field with
Gaussian correlation function.
We showone realization of the Young’smodulus in Figure 13a. Themean value of E is 17 GPa and the standard
deviation 6.8 GPa; the mean value of gc is 120 Pa m and the standard deviation 48 Pa m. The correlation
length is 0.05 m. An initial horizontal fracture of length 0.2 m is located at the center of the left side. We
minimize boundary effects by embedding this domain inside a bigger square domain with sides of 60 m. We
inject fluid in the initial fracture at prescribed pressure. The fracture pattern (Figure 13b) tends to follow the
weakest points and results in a complex nonlinear trajectory that eventually branches into two fractures. This
simulation provides an example of the capabilities of the approach to capture complex fracture patterns.
5. Conclusions
The simulation of fluid-driven fractures in geologic porous media requires a modeling framework that is
able to incorporate the controlling mechanisms of fluid flow and rock fracture mechanics. Given the extreme
range of spatial and temporal scales involved in these processes for realistic applications, striking a balance
between physical fidelity and numerical tractability is a significant modeling challenge. Phase field method-
ologies, based on the variational approach to fracture, have recently emerged as promising tools tomeet that
challenge, but their ability to reproduce the complex processes involved in hydraulic fracturing in realistic
environments is still being debated.
The objective of the present study was to assess the performance of phase field modeling in simple injec-
tion scenarios where analytical solutions are available. Here we presented a phase field strategy tomodel the
propagation of fluid-filled fractures in brittle elastic media and validated our approach through comparison
with analytical solutions. The model considers the Reynolds lubrication equation to describe fluid flow in the
fracture and linear elastic fracturemechanics.We implemented a fully coupled approach to simulate the inter-
action between the pressure field in the fracture and the displacement field in the solid. We adopted a robust
sequential algorithm to update damage in the solid, along with a local history field that guarantees the irre-
versibility of the damage (fracture) process and an anisotropic degradation of the elastic energywhich honors
that fractures propagate only under tension.
We considered injection scenarios for which different propagation regimes have been identified: the
toughness- and viscous-dominated regimes and the leak-off toughness regime. We compared simulated
results of the evolution of fracture lengths, fluid pressures, and fracture apertures against their analytical
counterparts and found good agreement in all three propagation regimes. The present validation exercise
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2583
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
implies that the phase field approach captures the key controllingmechanisms of coupled flow, rockmechan-
ical deformation, and brittle fracture mechanics and therefore is a suitable tool to understand more complex
processes of fluid-driven fracturing of geologic porous media.
Appendix A: Fracture Propagation. Analytical Solutions
The solution of a plane strain hydraulic fracture in an infinite elastic solid depends on the injection rateQ, the
fluid dynamic viscosity 𝜇f , and four material parameters: Young’s modulus E, Poisson’s ratio 𝜈, toughness KIc,
and leak-off coefficient Cl . The formulations use the following effective material parameters:
E′ = E
1 − 𝜈2
, 𝜇′f = 12𝜇f , K
′ = 4
√
2
𝜋
KIc, C
′ = 2Cl, (A1)
where the toughness KIc is related to the Griffith critical energy release rate gc through
KIc =
√
1 − 𝜈2
gcE
. (A2)
The pressure inside the fracture is denoted as pf , the net pressure p is defined as p = pf (x, t) − 𝜎0, the fracture
opening isw(x, t), and the fracture half-length is l(t), where x is the position along the fracture. These variables
are expressed in the following form:
w(x, t) = 𝜖(t)L(t)Ω(𝜉, t), (A3)
p(x, t) = 𝜖(t)E′Π(𝜉, t), (A4)
l(t) = L(t)𝛾(t), (A5)
whereΩ is the dimensionless opening,Π is the dimensionless net pressure, 𝛾 is the dimensionless half-length,
𝜉 = x∕l(t) ∈ [0, 1] is the scaled coordinate, and L(t) is the fracture length scale; for an impermeable solid
𝜖(t) = V(t)
L2(t)
, (A6)
and for a permeable one,
𝜖(t) = C
′2
Q
, (A7)
where 𝜖(t) is a characteristic fracture aspect ratio—a small dimensionless parameter which depends on the
volume of injected fluid V(t), among other variables.
The fluid leak off per unit fracture length is studied within Carter’s model [Carter, 1957], which is described by
an inverse square root of time law of the form
q(x, t) = C
′√
t − t0(x)
, t − t0(x)> 0, (A8)
where t0 is the time at which the fracture first propagates to a given point of coordinate x.
A1. Toughness-Dominated Regime
Garagash [2006] provides an explicit solution for a fracture propagating in an impermeable elastic media in
the toughness-dominated regime. The fracture length scale Lk(t) is given by
Lk =
(
E′Q0t
K ′
)2∕3
, (A1.1)
where the subscript k denotes the regime. The dimensionless viscosity is defined as
 = 𝜇′Q
E′
(
E′
K ′
)4
. (A1.2)
A threshold value of the viscosity for this propagation regime is0 = 3.4×10−3, and this regime implies that < 0.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2584
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
The solution of the problem is given by a series expansion in the small parameter,
 (, 𝜉) =
∞∑
j=0
jj(𝜉) (A1.3)
with j(𝜉) = {Ω̄kj(𝜉),Πkj(𝜉), 𝛾kj}, and Ω̄kj(𝜉) = Ωkj(𝜉)∕𝛾k(t), combined with equations (A3)–(A5). The
zeroth-order solution, i.e., the first term in the expansion or j = 0, is the zero-viscosity solution. The first-order
solution, j = 0, 1, provides the correction of order to the zero-viscosity solution and is denoted as the small
viscosity solution or large toughness solution. Garagash [2006] points out that the first-order solution provides
an excellent approximation for a wide range of the viscosity parameter. Moreover, he improved this solution
and proposed a new formulation in terms of the small parameter 𝛿()
𝛿() = √
1 +∕Si
, (A1.4)
whereSi = 1∕30. The improved first-order solution has the form
 (, 𝜉) = 0(𝜉) + 𝛿()1(𝜉). (A1.5)
The zeroth-order dimensionless variables are given by
Ω̄k0(𝜉) =
𝜋1∕3
2
√
1 − 𝜉2, (A1.6)
Πk0 =
𝜋1∕3
8
, (A1.7)
𝛾k0 =
2
𝜋2∕3
(A1.8)
and the first-order variables by
Ω̄k1(𝜉) =
8
3𝜋2∕3
⎛⎜⎜⎜⎝2𝜋 − 4𝜉 arcsin 𝜉 −
(5
6
− ln 2
)√
1 − 𝜉2 − 3
2
ln
|||||||||
(
1 +
√
1 − 𝜉2
)1+√1−𝜉2
(
1 −
√
1 − 𝜉2
)1−√1−𝜉2
|||||||||
⎞⎟⎟⎟⎠ , (A1.9)
Πk1 =
8
3𝜋2∕3
(
1
24
+ ln
(
4
√
1 − 𝜉2
)
− 3𝜉 arccos 𝜉
4
√
1 − 𝜉2
)
, (A1.10)
𝛾k1 = −
32(1 + 6 ln 2)
9𝜋5∕3
. (A1.11)
A2. Viscous-Dominated Regime
Anexplicit analytical solution for theproblemof fluid-driven fracturing in an impermeable elasticmedia in the
viscous-dominated regime was found by Adachi and Detournay [2002] and Garagash and Detournay [2005].
The fracture length scale Lm(t) is given by
Lm =
(
E′Q3t4
𝜇′
)1∕6
, (A2.1)
where the subscriptm denotes the regime. The dimensionless toughness is given by
 = K ′
E′
(
E′
𝜇′Q
)1∕4
. (A2.2)
A threshold value of the toughness for this propagation regime is0 = 1.42, with a range of validity < 0.
As in the toughness-dominated regime, the solution of the problem is given by a series expansion in the
parameter (),
 (, 𝜉) =
∞∑
j=0
()jj(𝜉), (A2.3)
with j(𝜉) = {Ω̄mj(𝜉),Πmj(𝜉), 𝛾mj}, and Ω̄mj(𝜉) = Ωmj(𝜉)∕𝛾m(t), combined with equations (A3)–(A5). The
parameter () will be defined later. The zeroth-order solution, j = 0, is the zero-toughness solution.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2585
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
The first-order solution, j = 0, 1, provides the correction of order () to the zero-toughness solution and is
denoted as the small toughness solutionor large viscosity solution. As shown inGaragashandDetournay [2005],
the first-order solution provides an excellent approximation.
The zeroth-order dimensionless variables are given by Adachi and Detournay [2002]:
Ω̄m0(𝜉) = A00(1 − 𝜉2)2∕3 + A01(1 − 𝜉2)5∕3 + B01
[
4
√
1 − 𝜉2 + 2𝜉2 ln
||||||
1 −
√
1 − 𝜉2
1 +
√
1 − 𝜉2
||||||
]
, (A2.4)
Πm0(𝜉) = A00
1
3𝜋
B
(1
2
,
2
3
)
2F1
(
−1
6
, 1; 1
2
, 𝜉2
)
+ A01
10
7 2
F1
(
−7
6
, 1; 1
2
, 𝜉2
)
+ B01(2 − 𝜋 |𝜉|), (A2.5)
𝛾m0 = 0.61524, (A2.6)
where A00 =
√
3, A01 ≈ −0.15601, B01 ≈ 6.632 × 10−2, B(⋅, ⋅) is the Euler’s beta function, and 2F1(⋅, ⋅; ⋅, ⋅) is the
Gauss hypergeometric function.
The first-order dimensionless variables read as follows [Garagash and Detournay, 2005]:
Ω̄m1(𝜉) = A10(1 − 𝜉2)h + A11(1 − 𝜉2)h+1C
h+1∕2
2 (𝜉) + B11
[
4
√
1 − 𝜉2 + 2𝜉2 ln
||||||
1 −
√
1 − 𝜉2
1 +
√
1 − 𝜉2
||||||
]
, (A2.7)
Πm1(𝜉) = A10
Γ(1 + h)
2𝜋1∕2Γ
(
1
2
+ h
) 2F1 (12 − h, 1; 12 , 𝜉2)
+ A11
−1
B
(
1
2
+ h, 3
2
) [1
2 2
F1
(
−3
2
− h, 1; 1
2
, 𝜉2
)
− (1 + h)𝜉22F1
(
−1
2
− h, 2; 3
2
, 𝜉2
)]
+ B11(2 − 𝜋 |𝜉|),
(A2.8)
𝛾m1 = −
1
2
𝛾3m0
[
A10B
(1
2
, 1 + h
)
− A11
1∕2 + h
5∕2 + h
B
(1
2
, 2 + h
)
+ B11
4𝜋
3
]
, (A2.9)
whereh≈ 0.138673,C(𝜅)2 (⋅) is theGegenbauer polynomial of degree 2 and index𝜅,A10 = 2
−h,A11≈ −4.0042 ×
10−2, B11≈ −4.5476 × 10−2, and Γ(⋅) is the Euler gamma function.
The Gegenbauer polynomial of degree n and index 𝜅, C(𝜅)n (⋅), has the form
C(𝜅)n (x) =
n∕2∑
k=0
(−1)k Γ(n − k + 𝜅)
Γ(𝜅)k!(n − 2k)!
(2x)n−2k. (A2.10)
Finally, the parameter () is defined as
() = B1b, (A2.11)
where b and B1 are computed from the constants 𝛾m0, h, and 𝛽1 ≈ 0.03719 according to
b = 4 − 6h ≈ 3.16796 (A2.12)
and
B1 =
(2
3
)2h−1
𝛾3h−20 𝛽1 ≈ 0.1076. (A2.13)
A3. Leak-Off Toughness-Dominated Regime
Bunger et al. [2005] consider the problem of the propagation of a fracture by a zero fluid viscosity in a perme-
able medium. The authors studied the so-called small and large time asymptotic behaviors, which correspond
to the storage toughness-dominated regime and leak-off toughness-dominated regime, respectively. During
fracture growth, the propagation evolves in time and transitions from the storage toughness-dominated
regime to the leak-off toughness-dominated regime.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2586
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
The fracture length scale Lk,k̃(t) is given by
Lk,k̃ =
(
K ′Q
E′C′2
)2
, (A3.1)
where the subscript k or k̃ denotes the regime. The dimensionless viscosityk,k̃ is defined as
k,k̃ = 𝜇
′QE
′3
K ′4
, (A3.2)
and the dimensionless time 𝜏 is defined as
𝜏 = t E
′4C
′2
K ′4Q2
. (A3.3)
The solution obtained by the authors is valid whenk,k̃ ≪ 1.
For 𝜏 ≪ 1, the fracture propagates under the storage toughness-dominated regime. The dimensionless
variable 𝛾k is then given by
𝛾k(𝜏) = 𝜏2∕3
n∑
i=0
𝛾ki𝜏
−i∕6, (A3.4)
where
𝛾k0 =
2
𝜋2∕3
, (A3.5)
𝛾k1 = −
16I2∕3
3𝜋4∕3
, (A3.6)
𝛾k2 =
32I2∕3
9𝜋2
(
I2∕3 − 4I5∕6
)
, (A3.7)
𝛾k3 = −
512I2∕3
81𝜋8∕3
[
−8I5∕6 + I2∕3
(
2 − I2∕3 + 3I5∕6
)]
, (A3.8)
and I𝜑 is an integral solved in terms of the Gamma function Γ(⋅) as
I𝜑 = ∫
1
0
𝜂𝜑√
1 − 𝜂
d𝜂 =
Γ(𝜑 + 1)Γ(1∕2)
Γ(𝜑 + 3∕2)
. (A3.9)
For 𝜏 ≫ 1, the propagating regime has evolved to the leak-off toughness-dominated regime. The dimension-
less variable 𝛾k̃ is now given by
𝛾k̃(𝜏) =
√
𝜏
n∑
i=0
𝛾k̃i𝜏
−i∕4, (A3.10)
where
𝛾k̃0 =
1
𝜋
, (A3.11)
𝛾k̃1 = −
1
4
√
2𝜋I1∕4
, (A3.12)
𝛾k̃2 =
3
128I1∕4
, (A3.13)
𝛾k̃3 = −
3
√
𝜋(1 + 3I1∕4)
1024
√
2I−1∕4I
2
1∕4
. (A3.14)
The dimensionless opening and pressure are given as
Ωk,k̃ = 2−1∕2𝛾 (1∕2)
√
1 − 𝜉2 (A3.15)
and
Πk,k̃ = 2−5∕2𝛾−1∕2. (A3.16)
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2587
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
References
Adachi, J. I., and E. Detournay (2002), Self-similar solution of a plane-strain fracture driven by a power-law fluid, Int. J. Numer. Anal. Methods
Geomech., 26(6), 579–604.
Adachi, J. I., and E. Detournay (2008), Plane strain propagation of a hydraulic fracture in a permeable rock, Eng. Fract. Mech., 75(16),
4666–4694.
Ambati, M., T. Gerasimov, and L. De Lorenzis (2015), A review on phase-field models of brittle fracture and a new fast hybrid formulation,
Comput. Mech., 55, 383–405.
Ambrosio, L., and V. M. Tortorelli (1990), Approximation of functional depending on jumps by elliptic functional via Γ-convergence,
Commun. Pure Appl. Math., 43, 999–1036.
Aranson, I. S., V. A. Kalatsky, and V. M. Vinokur (2000), Continuum field description of crack propagation, Phys. Rev. Lett., 85, 118–121.
Bunger, A. P., E. Detournay, and D. I. Garagash (2005), Toughness-dominated hydraulic fracture with leak-off, Int. J. Fract., 134, 175–190.
Bouchard, P. O., F. Bay, and Y. Chastel (2003), Numerical modelling of crack propagation: Automatic remeshing and comparison of different
criteria, Comput. Methods Appl. Mech. Eng., 192, 3887–3908.
Bourdin, B., G. A. Francfort, and J. J. Marigo (2000), Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids, 48, 797–826.
Bourdin, B., C. P. Chukwudozie, and K. Yoshioka (2012), A variational approach to the numerical simulation of hydraulic fracturing, in SPE
Annual Technical Conference and Exhibition, Soc. of Pet. Eng., SPE-159154-MS, San Antonio, Tex.
Cappa, F., and J. Rutqvist (2011), Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep
underground injection of CO2, Int. J. Greenhouse Gas Control, 5(2), 336–346.
Carrier, B., and S. Granet (2012), Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model,
Eng. Fract. Mech., 79, 312–328.
Carter, E. D. (1957), Optimum fluid characteristics for fracture extension, in Drilling and Production Practice, edited by G. C. Howard and
C. R. Fast, pp. 261–270, Am. Pet. Inst., New York.
Chambolle, A. (2004), An approximation result for special functions with bounded deformation, J. Math. Pures Appl., 83, 929–954.
Chukwudozie, C., B. Bourdin, and K. Yoshioka (2013), A variational approach to the modeling and numerical simulation of hydraulic
fracturing under in-situ stresses, in Proceedings of the 38th Workshop on Geothermal Reservoir Engineering, Stanford Geothermal
Program, Stanford, Calif.
Cueto-Felgueroso, L., and R. Juanes (2013), Forecasting long-term gas production from shale, Proc. Natl. Acad. Sci. U.S.A., 110(49),
19660–19661.
Detournay, E. (2004), Propagation regimes of fluid-driven fractures in impermeable rocks, Int. J. Geomech., 4(1), 35–45.
Detournay, E. (2016), Mechanics of hydraulic fractures, Annu. Rev. Fluid Mech., 48, 311–339.
Economides, M. J., and K. G. Nolte (2000), Reservoir Stimulation, 3rd ed., Wiley, Chichester, West Sussex, U. K.
Evans, K. F., H. Moriya, H. Niitsuma, R. H. Jones, W. S. Phillips, A. Genter, J. Sausse, R. Jung, and R. Baria (2005), Microseismicity and
permeability enhancement of hydrogeologic structures during massive fluid injections into granite at 3 km depth at the Soultz HDR site,
Geophys. J. Int., 160(1), 388–412.
Francfort, G. A., and J. J. Marigo (1998), Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids, 46, 1319–1342.
Fu, P., S. M. Johnson, and C. R. Carrigan (2013), An explicitly coupled hydro-geomechanical model for simulating hydraulic fracturing in
arbitrary discrete fracture networks, Int. J. Numer. Anal. Methods Geomech., 37, 2278–2300.
Garagash, D. I. (2006), Plane-strain propagation of a fluid-driven fracture during injection and shut-in: Asymptotics of large toughness,
Eng. Fract. Mech., 73(4), 456–481.
Garagash, D. I., and E. Detournay (2005), Plane-strain propagation of a fluid-driven fracture: Small toughness solution, J. Appl. Mech., 72(6),
916–928.
Geertsma, J., and P. De Klerk (1969), A rapid method of predicting width and extent of hydraulically induced fractures, J. Pet. Technol., 246,
1571–1581.
Griffith, A. A. (1920), The phenomena of rupture and flow in solids, Philos. Trans. R. Soc., A221, 163–198.
Hafver, A., E. Jettestuen, J. Feder, P. Meakin, and A. Malthe-Sørenssen (2014), A node-splitting discrete element model for fluid-structure
interaction, Physica A, 416, 61–79.
Henry, H., and H. Levine (2013), Dynamic instabilities of fracture under biaxial strain using a phase field model, Phys. Rev. Lett., 93, 105504.
Hofacker, M., and C. Miehe (2012), Continuum phase field modeling of dynamic fracture: Variational principles and staggered
FE implementation, Int. J. Fract., 178, 113–129.
Hofacker, M., and C. Miehe (2013), A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns,
Int. J. Numer. Methods Eng., 93, 276–301.
Hu, J., and D. I. Garagash (2010), Plane-strain propagation of a fluid-driven crack in a permeable rock with fracture toughness, J. Eng. Mech,
136(9), 1152–1166.
Jeffrey, R. G., and K. W. Mills (2000), Hydraulic fracturing applied to inducing longwall coal mine goaf falls, in Pacific Rocks
2000—Proceedings of the 4th North American Rock Mechanics Symposium, edited by J. Girard et al., pp. 423–430, Balkema, Rotterdam,
Netherlands.
Jha, B., and R. Juanes (2014), Coupled multiphase flow and poromechanics: A computational model of pore-pressure effects on fault slip
and earthquake triggering,Water Resour. Res., 50, 3776–3808, doi:10.1002/2013WR015175.
Karma, A., D. A. Kessler, and H. Levine (2001), Phase-field model of mode III dynamic fracture, Phys. Rev. Lett., 87, 45501.
Khristianovich, S. A., and Y. P. Zheltov (1955), Formation of vertical fractures by means of highly viscous liquid, in Proceeding 4th World
Petroleum Congress, vol. 2, pp. 556–579, World Petroleum Congress, Rome.
Kuhn, C., A. Schlüte, and R. Müller (2015), On degradation functions in phase field fracture models, Comput. Mater. Sci., 108, 374–384.
Kuuskraa, V., S. H. Stevens, and K. D. Moodhe (2013), Technically Recoverable Shale Oil and Shale Gas Resources: An Assessment of 137 Shale
Formations in 41 Countries Outside the United States, U.S Energy Inform. Administr., Washington, D. C.
Legarth, B., E. Huenges, and G. Zimmermann (2005), Hydraulic fracturing in a sedimentary geothermal reservoir: Results and implications,
Int. J. Rock Mech. Min. Sci., 42(7), 1028–1041.
Miehe, C. (1998), Comparison of two algorithms for the computation of fourth-order isotropic tensor functions, Comput. Struct., 66, 37–43.
Miehe, C., F. Welschinger, and M. Hofacker (2010a), Thermodynamically consistent phase-field models of fracture: Variational principles and
multi-field FE implementations, Int. J. Numer. Methods Eng., 83, 1273–1311.
Miehe, C., M. Hofacker, and F. Welschinger (2010b), A phase field model for rate-independent crack propagation: Robust algorithmic
implementation based on operator splits, Comput. Methods Appl. Mech. Eng., 199, 2765–2778.
Miehe, C., L. M. Schänzel, and H. Ulmer (2015a), Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface
and failure criteria for brittle crack propagation in thermo-elastic solids, Comput. Methods Appl. Mech. Eng., 294, 449–485.
Acknowledgments
D.S. acknowledges funding from
the Spanish Ministry of Economy
and Competitiveness (grant
CTM2014-54312-P). R.J. was supported
by the U.S. Department of Energy
through a DOE Mathematical
Multifaceted Integrated Capability
Center (grant DE-SC0009286). L.C.F.
gratefully acknowledges funding
from the Spanish Ministry of
Economy and Competitiveness
(grants RyC-2012-11704 and
CTM2014-54312-P). D.S. specially
thanks Department of Civil & Environ-
mental Engineering, Massachusetts
Institute of Technology for hospital-
ity during his sabbatical stay from
March to June, 2015, and “Fun-
dación José Entrecanales Ibarra” for
its support. The authors would like
to thank the anonymous reviewers
for their insightful comments and
suggestions that have contributed to
improve this paper. No data were used
in producing this manuscript.
SANTILLÁN ET AL. PHASE FIELDMODEL FLUID-DRIVEN FRACTURE 2588
Journal of Geophysical Research: Solid Earth 10.1002/2016JB013572
Miehe, C., M. Hofacker, L. M. Schänzel, and F. Aldakheel (2015b), Phase field modeling of fracture in multi-physics problems. Part II. Coupled
brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids, Comput. Methods Appl. Mech. Eng., 294, 486–522.
Miehe, C., S. Mauthe, and S. Teichtmeister (2015c), Minimization principles for the coupled problem of Darcy-Biot-type fluid transport in
porous media linked to phase field modeling of fracture, J. Mech. Phys. Solids, 82, 186–217.
Miehe, C., and S. Mauthe (2016), Phase field modeling of fracture in multi-physics problems. Part III. Crack driving forces in
hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media, Comput. Methods Appl. Mech. Eng., 304, 619–655.
Mikelic, A., M. F. Wheeler, and T. Wick (2015a), A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous
medium,Multiscale Model. Simul., 13(1), 367–398.
Mikelic, A., M. F. Wheeler, and T. Wick (2015b), A quasi-static phase-field approach to pressurized fractures, Nonlinearity, 28(5), 1371–1399.
Mumford, D., and J. Shah (1989), Optimal approximations by piecewise smooth functions and associated variational problems, Commun.
Pure Appl. Math., 42, 577–685.
Murdoch, L. C. (2002), Mechanical analysis of idealized shallow hydraulic fracture, J. Geotech. Geoenviron. Eng., 128(6), 488–495.
Nordgren, R. P. (1972), Propagation of a vertical hydraulic fracture, Soc. Pet. Eng. J., 12, 306–314.
Oliver, J., A. E. Huespe, and P. J. Sánchez (2006), A comparative study on finite elements for capturing strong discontinuities: E-FEM vs X-FEM,
Comput. Methods Appl. Mech. Eng., 195, 4732–4752.
Osborn, S. G., A. Vengosh, N. R. Warner, and R. B. Jackson (2011), Methane contamination of drinking water accompanying gas-well drilling
and hydraulic fracturing, Proc. Natl. Acad. Sci. U.S.A., 108(20), 8172–8176.
Ouchi, H., A. Katiyar, J. York, J. T. Foster, and M. M. Sharma (2015), A fully coupled porous flow and geomechanics model for fluid driven
cracks: A peridynamics approach, Comput. Mech., 55, 561–576.
Parmigiani, A., S. Faroughi, C. Huber, O. Bachmann, and Y. Su (2016), Bubble accumulation and its role in the evolution of magma reservoirs
in the upper crust, Nature, 532(7600), 492–495.
Patzek, T. W., F. Male, and M. Marder (2013), Gas production in the Barnett Shale obeys a simple scaling theory, Proc. Natl. Acad. Sci. U.S.A.,
110(49), 19,731–19,736.
Perkins, T. K., and L. R. Kern (1961), Widths of hydraulic fractures, J. Petrol. Technol., 13, 937–949.
Rubin, A. M. (1995), Propagation of magma-filled cracks, Annu. Rev. Earth Planet. Sci., 23, 287–336.
Rungamornrat, J., M. F. Wheeler, and M. E. Mear (2005), A numerical technique for simulating nonplanar evolution of hydraulic fractures, in
Proceedings—SPE Annual Technical Conference and Exhibition, pp. 3813–3821, SPE, Dallas, Tex.
Shojaei, A., A. D. Taleghani, and G. Li (2014), A continuum damage failure model for hydraulic fracturing of porous rocks, Int. J. Plast., 59,
199–212.
Sneddon, I. (1946), The distribution of stress in the neighbourhood of a crack in an elastic solid, Proc. R. Soc. A, 187, 229–260.
Spence, D. A., and D. L. Turcotte (1985), Magma-driven propagation of cracks, J. Geophys. Res., 90(B1), 575–580.
Vidic, R. D., S. L. Brantley, J. M. Vandenbossche, D. Yoxtheimer, and J. D. Abad (2013), Impact of shale gas development on regional water
quality, Science, 340(6134), 1235009.
Vignollet, J., S. May, R. De Borst, and C. V. Verhoosel (2014), Phase-field models for brittle and cohesive fracture,Meccanica, 49, 2587–2601.
Wangen, M. (2011), Finite element modeling of hydraulic fracturing on a reservoir scale in 2D, J. Petrol. Sci. Eng., 77, 274–285.
Witherspoon, P. A, J. Wang, K. Iwai, and J. Gale (1980), Validity of cubic law for fluid flow in a deformable rock fracture,Water Resour. Res.,
16(6), 1016–1024.
Erratum
In the originally published version of this article, the acknowledgments were incomplete. The text has since
been updated, and this version may be considered the authoritative version of record.
SANTILLÁN ET AL. PHASE FIELD MODEL FLUID-DRIVEN FRACTURE 2589