# A Multi-Resolution Approach
to Hydraulic Fracture Simulation

Andre Costa  Department of Mechanical Engineering and Material Science, Duke
University, Durham, 27708, NC, USA  Matteo Cusini  Atmospheric, Earth and
Energy Division, Lawrence Livermore National Laboratories, Livermore, 94550,
CA, USA  Tao Jin  Department of Mechanical Engineering, University of Ottawa,
Ottawa, K1N6N5, ON, Canada  Randolph Settgast  Atmospheric, Earth and Energy
Division, Lawrence Livermore National Laboratories, Livermore, 94550, CA, USA
John Dolbow  Department of Mechanical Engineering and Material Science, Duke
University, Durham, 27708, NC, USA  Corresponding author: John Dolbow,
<EMAIL_ADDRESS>
###### Abstract

We present a multi-resolution approach for constructing model-based
simulations of hydraulic fracturing, wherein flow through porous media is
coupled with fluid-driven fracture. The approach consists of a hybrid scheme
that couples a discrete crack representation in a global domain to a phase-
field representation in a local subdomain near the crack tip. The multi-
resolution approach addresses issues such as the computational expense of
accurate hydraulic fracture simulations and the difficulties associated with
reconstructing crack apertures from diffuse fracture representations. In the
global domain, a coupled system of equations for displacements and pressures
is considered. The crack geometry is assumed to be fixed and the displacement
field is enriched with discontinuous functions. Around the crack tips in the
local subdomains, phase-field sub-problems are instantiated on the fly to
propagate fractures in arbitrary, mesh independent directions. The governing
equations and fields in the global and local domains are approximated using a
combination of finite-volume and finite element discretizations. The efficacy
of the method is illustrated through various benchmark problems in hydraulic
fracturing, as well as a new study of fluid-driven crack growth around a stiff
inclusion.

##  1 Introduction

A wide range of approaches for model-based simulations of hydraulic fracturing
have been developed over the past several decades [1, 2]. These range from
production-level reservoir modeling tools such as ResFrac[3, 4], to GEOS[5, 6,
7], PyFrac[8], and others. Many of the models and associated codes assume
fracture networks that remain planar, but in recent years strides have been
made towards modeling cracks that evolve in arbitrary ways in response to
fluid-driven loads. High-resolution models for complex fracture evolution
generally fall into two categories: sharp interface models that explicitly
model the fracture surface, and diffuse crack approaches that effectively
smear the geometry over the underlying grid or mesh. Techniques that represent
the crack as a sharp interface can be advantageous when the fracture
configuration is relatively simple, but representing complex geometric
evolution can be challenging [9, 10, 11]. By contrast, diffuse crack models
offer more flexibility for representing complex fracture evolution, but
introduce other challenges such as the lack of a well-defined fracture surface
and increased computational expense[12]. In this work, we introduce a multi-
resolution scheme for hydraulic fracturing simulation that makes use of both
sharp and diffuse crack representations within a single framework. The
objective is to establish a methodology that makes use of the advantageous
aspects of both sharp and diffuse crack models while circumventing some of the
drawbacks.

Over the past several decades, the phase-field model for fracture [13, 14, 15]
has emerged as a promising approach for constructing robust simulations of
complex crack evolution. The method has shown considerable success for
simulating fracture evolution in quasi-brittle materials, and there have been
several recent efforts to extend the approach to hydraulic fracturing. In what
follows, we review some prior works of particular relevance to the current
manuscript. For additional references in this topic, we refer the reader to
the recent review by Heider [12].

The first attempts towards a phase-field model for hydraulic fracture began
with extensions of the traditional phase-field model to pressurized cracks, as
in Bourdin et al. [16] and Wheeler et al. [17]. Subsequently, fluid flow in
the fractures, and also poromechanics were considered. Miehe et al. [18, 19]
developed a thermodynamically consistent framework, from minimization
principles, to couple poromechanics, fluid-flow and phase-field fracture. The
flow problem was modeled via the Darcy’s equation, containing a permeability
coefficient that used the phase-field variable and the crack opening to mimic
the cubic relationship from the lubrication theory in the crack region.
Mikelic et al. [20, 21] developed a model that separated the domain into
fracture and reservoir, by using the phase-field variable as an indicator
function. They also considered the flow inside the fracture as a Darcy flow,
but their model treated the fracture as a three-dimensional entity, which led
to a different permeability tensor compared to [18, 19]. Yet another approach
concerns the work of Wilson and Landis [22], who proposed a model that
included fluid velocities as primary variables. This allowed for a more
detailed description of the flow within the fracture, which was modeled by a
Brinkman-type equation [23]. The phase-field parameter acted as an indicator
of the flow regime, between Darcy flow (away from cracks) and Stokes flow
(inside cracks). Finally, the recent work of Chukwudozie et al. [24] presented
a different model, wherein the lubrication theory equations were included in
the weak form by means of a $\Gamma$-convergent regularization.

The use of a phase-field to represent a fracture network in a diffuse manner
certainly facilitates the representation of complex geometric evolution,
including crack branching and merging. However, it also requires the use of
meshes or grids that are capable of resolving the regularization length,
making these approaches computationally expensive. In the specific case of
hydraulic fracturing, another challenge concerns the crack opening or
aperture, a field that is tightly coupled with the fluid pressure within
fractures. In a phase-field setting, due to the lack of an explicit crack
surface, extracting the aperture or accounting for its effects requires
additional considerations. All of the aforementioned works present some way to
account for the aperture within a diffuse setting, but the robustness of these
approaches remains unclear [2]. For a review of the most frequently used
methods to calculate the crack aperture from phase-field simulations, see the
recent work of Yoshioka et al. [25].

Outside of the context of hydraulic fracturing, some researchers in the phase-
field community have developed “hybrid” approaches, wherein the phase-field
formulation was combined with a sharp crack representation. The motivation for
these approaches varies, from “cutting” the mesh to remove artificial traction
transmission and circumvent element distortion [26] to reducing the overall
computational cost [27, 28]. In the work of Giovanardi et al. [27], phase-
field subproblems in the vicinity of crack tips were used to propagate a
global, discrete crack. The eXtended Finite Element Method (XFEM)[29] was used
to place fracture discontinuities in the displacement field within the
background global mesh. More recently, Muixi et al. [28] created an approach
that uses the phase-field method only at the crack tips, and XFEM in the rest
of the domain. In contrast to [27], there is no overlap of the representations
in crack tip areas.

The success of these hybrid approaches for purely mechanical cases opens the
door for their extension to hydraulic fracturing. Such approaches are
appealing because in principle they can circumvent the need for a complicated
reconstruction of the crack opening from the phase-field. This area is
relatively unexplored, although there have been some recent efforts that are
similar in spirit, such as the recent work of Sun et al. [30]. They developed
a Finite Element-Meshfree method to represent the crack surfaces in a discrete
fashion. The computed displacement field was used to obtain a driving force
which was employed within a phase-field evolution equation near the crack
tips. This approach eliminated the need for the reconstruction of crack
openings from the diffuse crack representation, but it also largely decoupled
the phase field from the equations governing the force balance near the crack
tips.

The approach presented in this work, which we refer to as a multi-resolution
method, extends the concept of a hybrid phase-field method to hydraulic
fracturing. It discretizes the problem at a global level using a sharp
interface approach based on the Embedded Finite Element Method of Cusini et
al. [31]. It then approaches the simulation of fracture evolution by coupling
the global fields with a phase-field fracture problem posed over a subdomain
in the vicinity of the crack front. This approach has several advantages. The
crack aperture and flow inside the fracture are handled at the global scale
using techniques that work well and are efficient when the crack geometry is
known. Then the phase-field model is employed in the subdomain to effectively
update the crack geometry. This framework lends itself to incorporation with a
wide range of existing hydraulic fracturing solvers, as the phase-field sub-
problem is agnostic about the type of numerical treatment used in the global
domain.

The paper is structured as follows. In Section 2, we present the governing
equations and constitutive assumptions for both hydraulic fracturing and
phase-field for fracture. In Section 3, we propose our multi-resolution
framework and present numerical schemes to discretize both the hydrofracture
and the phase-field subproblems. In Section 4, we apply our method to study
hydraulic fracturing in some simple scenarios, in order to verify its
accuracy. These include the well-known KGD[32, 33] problems and a case of non-
planar fracture propagation around a stiff inclusion. Finally, in Section 5,
we provide a summary and some concluding remarks.

##  2 Model formulation

In this section, we describe the models that are used to develop our multi-
resolution scheme. We start by presenting the governing equations to model
hydraulic fracture in poroelastic rock that forms the basis for the solver at
the global scale. We then describe a phase-field model for fracture that forms
the basis for the solver used in local subdomains near the crack tips.

###  2.1 Governing equations for hydraulic fracture

Figure 1: Schematic of a poroelastic rock with a fracture, inspired by [34].

We consider a model that couples flow and elastic deformation in a porous
media with evolving fracture surfaces. Consider the domain $\Omega$ consisting
of a porous rock, that is fully saturated with a single-phase Newtonian fluid.
The external boundary is composed of both traction $\partial\Omega_{t}$ and
displacement $\partial\Omega_{u}$ surfaces, viz.
$\partial\Omega=\partial\Omega_{t}\cup\partial\Omega_{u}$.

The domain contains fractures $\Gamma$, as shown in Figure 1. For simplicity,
we assume quasi-static conditions and small-strain kinematics. Neglecting
inertial effects, the balance of linear momentum reads

$\nabla\cdot\boldsymbol{\sigma}+\textbf{b}=\boldsymbol{0}\ \text{on
}\Omega\setminus\Gamma,$ (1)

where $\boldsymbol{\sigma}$ denotes the total Cauchy stress and b the body
force. The mechanical boundary conditions are given by

$\boldsymbol{\sigma}\cdot\textbf{n}=\textbf{t}\text{ on }\partial\Omega_{t},$
(2) $\textbf{u}=\overline{\textbf{u}}\text{ on }\partial\Omega_{u},$ (3)

where n denotes the normal direction at any point on $\partial\Omega$, t are
the applied tractions, and $\overline{\textbf{u}}$ denotes the prescribed
displacements.

Fluid flow within the cracks gives rise to pressure loads on the crack
surfaces, translating into the boundary condition

$\boldsymbol{\sigma}^{+}\cdot\textbf{n}_{\Gamma}=-\boldsymbol{\sigma}^{-}\cdot\textbf{n}_{\Gamma}=-p_{f}\textbf{n}_{\Gamma}\text{
on }\Gamma,$ (4)

where $p_{f}$ is the pressure inside the crack and $\textbf{n}_{\Gamma}$ is
the normal vector to $\Gamma$ at any point. In this work, we assume that
fractures are open and neglect contact conditions on crack faces. For
additional considerations to account for closed fractures, see the model
described in Cusini et al.[31].

In terms of the fluid flow, the fluid velocity $\textbf{v}_{m}$ in the matrix
is governed by the mass balance

$\dfrac{\partial(\rho\phi)}{\partial
t}+\nabla\cdot(\rho\textbf{v}_{m})=Q_{m}+Q_{mf}\ \text{on
}\Omega\setminus\Gamma.$ (5)

where $\rho$ denotes the fluid density and $\phi$ is the porosity, and $Q_{m}$
is a prescribed source term in the matrix. The source term $Q_{mf}$ accounts
for the exchange of fluid between the matrix and the fractures. The boundary
is partitioned into pressure $\partial\Omega_{p}$ and flux
$\partial\Omega_{q}$ portions, such that
$\partial\Omega=\partial\Omega_{p}\cup\partial\Omega_{q}$, and the following
boundary conditions are applied,

$\rho\textbf{v}_{m}\cdot\textbf{n}=\textbf{q}\text{ on }\partial\Omega_{q},$
(6)

$p_{m}=\overline{p}_{m}\text{ on }\partial\Omega_{p},$ (7)

where $p_{m}$ denotes the pore-pressure and $\overline{p}_{m}$ is a prescribed
pressure.

In a comparable manner, the fluid velocity $\textbf{v}_{f}$ within fractures
is governed by the mass balance

$\dfrac{\partial(\rho w)}{\partial t}+\nabla_{\Gamma}\cdot(\rho
w\textbf{v}_{f})=Q_{f}+Q_{fm}\ \text{on }\Gamma,$ (8)

where $w$ denotes the normal fracture aperture, $Q_{f}$ is a prescribed source
term within the fracture, and $Q_{fm}$ accounts for the exchange of fluid
between the fracture and the matrix111the flux interactions $Q_{fm}$ and
$Q_{mf}$ are modeled as in classical well models, following [35]. This ensures
the balance of mass between the fractures and matrix
$\int_{V}Q_{mf}dV=-\int_{\Gamma}Q_{fm}d\Gamma$.. In the above,
$\nabla_{\Gamma}$ indicates a gradient operator taken on the lower dimensional
manifold $\Gamma$. Boundary conditions similar to (6) and (7) can also be
applied.

Constitutive relationships are required to close the system and tie the
stresses to the displacements u and the fluid velocities to the pressures. In
particular, we adopt the basic assumptions of Biot’s theory of poroelasticity
[36], Darcy’s law for the flow in the matrix, and a lubrication theory
approximation for the flow in fractures. This gives rise to the following set
of constitutive relationships for the stress, porosity, and velocities:

$\boldsymbol{\sigma}=\mathbb{C}:\boldsymbol{\epsilon}(\textbf{u})-\alpha
p_{m}\mathbb{I}\ \text{on }\Omega\setminus\Gamma,$ (9)

$\dot{\phi}=\alpha\nabla\cdot\dot{\textbf{u}}+\dfrac{\dot{p_{m}}}{N}\ \text{on
}\Omega\setminus\Gamma,$ (10)

$\textbf{v}_{m}=-\dfrac{\kappa}{\mu}\nabla p_{m}\ \text{on
}\Omega\setminus\Gamma,$ (11)

$\textbf{v}_{f}=-\dfrac{w^{2}}{12\mu}\nabla_{\Gamma}p_{f}\ \text{on }\Gamma.$
(12)

In the above, $p_{m}$ is the pore-pressure, $\mathbb{C}$ is the fourth-order
isotropic tensor of drained elastic moduli,
$\boldsymbol{\epsilon}(\textbf{u})$ is the mechanical strain, $\alpha$ is the
Biot coefficient, and $\mathbb{I}$ is the second-order identity tensor. The
temporal evolution of the porosity is governed by the rate of dilatation and
time rate of change in the pore pressure, as modulated by the modulus $N$. The
fluid velocity in the matrix is related to the gradient of the pressure
through the ratio of the intrinsic permeability $\kappa$ to the viscosity
$\mu$.

The fluid is assumed to be linearly compressible. For both the fluid in the
fractures and the matrix, this implies that the density is updated from its
reference value $\rho_{ref}$ based on the change in pressure according to

$\rho=\rho_{ref}\left(1+\dfrac{p-p_{ref}}{K_{F}}\right)\ \text{on }\Omega,$
(13)

where $p_{ref}$ denotes a reference value for the pressure, and $K_{F}$ is the
fluid bulk modulus.

Finally, the initial conditions for the displacements and pressures are given
by

$\textbf{u}(\textbf{x},0)=\textbf{u}^{0}\ \text{on }\Omega\setminus\Gamma,$
(14)

$p_{m}(\textbf{x},0)=p^{0}_{m}\ \text{on }\Omega\setminus\Gamma,$ (15)

$p_{f}(\textbf{x},0)=p^{0}_{f}\ \text{on }\Gamma.$ (16)

For a given crack geometry, the combination of equations (1), (5) and (8),
with constitutive assumptions (9) - (13), boundary conditions (2)-(4),(6),(7)
and initial conditions (14)-(16) leads to a system of equations whose solution
can be approximated by many different numerical methods. What remains is a
model to describe the evolution of the crack geometry. In the next subsection,
we present the governing equations for a phase-field method for fracture, a
regularized approach for representing fractures and their evolution.

###  2.2 The phase-field method for fracture

The phase-field method for fracture started as an approximation [14] to the
variational approach for fracture by Francfort and Marigo [13] for the quasi-
static propagation of fracture in brittle materials. For the method adopted in
this work, we follow the work of Chukwudozie et al. [24] and associate the
following total energy to a crack configuration $\Gamma$ in a poroelastic
brittle solid $\Omega$:

$\mathcal{E}(\textbf{u},p_{m},p_{f},\Gamma)=\int\limits_{\Omega\setminus\Gamma}W(\boldsymbol{\epsilon}(\textbf{u}),p_{m})d\Omega-\int\limits_{\partial\Omega_{N}}\textbf{t}\cdot\textbf{u}ds-\int\limits_{\Omega\setminus\Gamma}\textbf{b}\cdot\textbf{u}d\Omega\\\
+\int\limits_{\Gamma}p_{f}\llbracket\textbf{u}\cdot\textbf{n}_{\Gamma}\rrbracket
ds+G_{c}\mathcal{H}^{n-1}(\Gamma),$ (17)

where the tractions applied to the boundary are denoted by t, the normal to
the crack is denoted by $\textbf{n}_{\Gamma}$ and $\mathcal{H}^{n-1}(\Gamma)$
is the $n-1$ dimensional Hausdorff measure of $\Gamma$.

The strain energy density $W(\boldsymbol{\epsilon}(\textbf{u}),p_{m})$ is
postulated as,

$W(\boldsymbol{\epsilon}(\textbf{u}),p_{m})=\dfrac{1}{2}\left(\boldsymbol{\epsilon}(\textbf{u})-\dfrac{\alpha}{nK}p_{m}\mathbb{I}\right):\mathbb{C}:\left(\boldsymbol{\epsilon}(\textbf{u})-\dfrac{\alpha}{nK}p_{m}\mathbb{I}\right),$
(18)

with $K$ denoting the bulk modulus and $n$ the system’s dimension (2 or 3).
The phase-field regularization, based on the Ambrosio-Tortorelli [37]
functional is then performed by the introduction of the damage parameter $d$
and the regularization length $\ell$,

$\mathcal{E}_{\ell}(\textbf{u},d,p_{m},p_{f})=\int\limits_{\Omega}\widetilde{W}(\boldsymbol{\epsilon}(\textbf{u}),p_{m},d)d\Omega-\int\limits_{\partial\Omega_{N}}\textbf{t}\cdot\textbf{u}ds-\int\limits_{\Omega}\textbf{b}\cdot\textbf{u}d\Omega+\int\limits_{\Omega}p_{f}\textbf{u}\cdot\nabla
dd\Omega\\\
+\dfrac{G_{c}}{c_{0}}\int\limits_{\Omega}\left(\dfrac{\zeta(d)}{\ell}+\ell\nabla
d\cdot\nabla d\right)d\Omega,$ (19)

where the function $\zeta(d)$ is the local dissipation function, which is
usually taken as $\zeta(d)=d$ or $d^{2}$. The constant $c_{0}$ is given by
$c_{0}=4\int_{0}^{1}\sqrt{\zeta(z)}dz$ and the regularized strain energy is
defined by,

$\widetilde{W}(\boldsymbol{\epsilon}(\textbf{u}),p_{m},d)=\dfrac{1}{2}\left((1-d)\boldsymbol{\epsilon}(\textbf{u})-\dfrac{\alpha}{nK}p_{m}\mathbb{I}\right):\mathbb{C}:\left((1-d)\boldsymbol{\epsilon}(\textbf{u})-\dfrac{\alpha}{nK}p_{m}\mathbb{I}\right),$
(20)

which is consistent with an assumption of damage arising in the sub pore scale
[24]. Finally, our regularized crack evolution problem is then stated as a
minimization principle for the functional
$\mathcal{E}_{\ell}(\textbf{u},d,p_{m},p_{f})$, with respect to the variables
u and $d$, with the added condition that the damage process is irreversible:

$\textbf{u},d=\underset{\textbf{u},d}{{\operatorname{argmin}}}\
\mathcal{E}_{\ell}(\textbf{u},d,p_{m},p_{f})\ \text{,\ \ subject to
}\dot{d}\geq 0.$ (21)

The following set of evolution equations can then be derived from the
Karush–Kuhn–Tucker (KKT) [38, 39] conditions:

$\nabla\cdot\left((1-d)^{2}\
\mathbb{C}:\boldsymbol{\epsilon}(\textbf{u})-(1-d)\alpha
p_{m}\mathbb{I}\right)+\textbf{b}=p_{f}\nabla d,$ (22)

$S_{d}\coloneqq\dfrac{2G_{c}\ell}{c_{0}}\Delta
d-\dfrac{G_{c}}{c_{0}\ell}\zeta^{\prime}(d)-g^{\prime}(d)W(\boldsymbol{\epsilon}(\boldsymbol{\textbf{u}}))+\alpha
p_{m}\nabla\cdot\textbf{u}+\nabla\cdot(p_{f}\textbf{u})\leq 0,$ (23)

$\dot{d}\geq 0,$ (24)

$S_{d}\dot{d}=0,$ (25)

with boundary conditions, $\boldsymbol{\sigma}\cdot\textbf{n}=\textbf{t}$ on
$\partial\Omega_{N}$ and $(2G_{c}\ell\nabla
d+c_{0}p_{f}\textbf{u})\cdot\textbf{n}\geq 0$ on $\partial\Omega$. The damaged
stress is defined as $\boldsymbol{\sigma}=(1-d)^{2}\
\mathbb{C}:\boldsymbol{\epsilon}(\textbf{u})-(1-d)\alpha p_{m}\mathbb{I}$.

##  3 Multi-resolution method

###  3.1 General overview

Figure 2: Global domain $\Omega$ on the left and local domain $\Omega_{L}$ on
the right. The surface corresponding to the local domain boundary
$\partial\Omega_{L}$ is indicated within the global domain with the dashed red
lines. The local domain is magnified to highlight the phase-field
representation of the global crack $\Gamma$.

We propose a multi-resolution method to approximate the solution of the
hydraulic fracture problem in porous media. It consists of coupling two
problems between a global domain and a local subdomain, as shown in Figure 2.
In the global domain problem, the governing equations (1)-(8) are discretized,
and the crack is represented with a sharp geometry. The crack geometry is
assumed to be fixed during a solution step in the global domain, and relevant
fields are calculated over the entire domain.

By contrast, the local subdomain concerns only a portion of the entire domain,
namely in the vicinity of crack tips. It is encapsulated within the global
domain as shown in Figure 2. In the local subdomain, a discretization of the
variational principle (21) is used to simulate crack evolution. During a
solution step in the local subdomain, the pressure fields within the fracture
and matrix are assumed to be fixed.

The two problems are coupled in the following manner. The displacement and
pressure fields are extracted from a solution step in the global domain and
passed to the local subproblem in different ways. In particular, the global
displacement fields are extracted along the surface $\partial\Omega_{L}$.
These fields are then applied as Dirichlet boundary conditions for the local
subproblem. The matrix pressure $p_{m}$ and pressure in the fracture $p_{f}$
are transfered as fields from the global domain to the local subdomain (see
Section 3.3.2 for details), and assumed to be fixed for the local subproblem.

The aforementioned operations provide everything that is needed for the local
subproblem from the global domain. Based on the imprint of the global crack
geometry on the local subdomain, a regularized fracture surface is created
(through an initial damage field) and then crack propagation is simulated in
the local subdomain. Once an extension of the crack in a scale that can be
represented in the global domain is identified, the sharp crack geometry in
the global problem is updated accordingly.

In principle, the aforementioned multi-resolution approach can be implemented
using a number of different discretization methods in the global domain and
the local subdomain. In the sections that follow, we describe the particular
choices used in this work as well as some important implementation details.
The method is described in a two-dimensional context, but many aspects can be
readily extended to three-dimensional problems.

###  3.2 Global problem discretization

The global problem in our multi-resolution approach encompasses the physics of
fluid flow in both the fractures and the pore structure of our domain, as well
as the deformation of the solid media. The set of governing equations consists
of (1)-(8), combined with constitutive assumptions (9)-(13) and appropriate
boundary conditions.

In this work, we use the discretization method proposed by Cusini et al. [31]
in the study of fluid flow in fractured porous media. The domain $\Omega$ is
partitioned with a mesh $\mathcal{T}$. Then, the intersection of the fracture
network $\Gamma$ with $\mathcal{T}$ defines the fracture triangulation
$\mathcal{F}$. These meshes are then used to define discrete counterparts
$\textbf{u}^{h}_{G}$, $p_{m}^{h}$ and $p_{f}^{h}$ of the unknown fields u,
$p_{m}$ and $p_{f}$, as well as discrete approximations of equations (1), (5)
and (8).

A finite element approximation is constructed for the displacement field and
employed in a standard Galerkin approximation to the global force balance (1).
The continuous part of the displacement field is approximated with a standard
space $\mathcal{U}$ of 4-node bilinear shape functions
$\\{{\boldsymbol{\eta}}_{a}\\}$. In the subset of elements that are “cut” by
the global crack geometry, a space $\mathcal{W}$ of discontinuous enrichment
functions $\\{\boldsymbol{\phi}_{b}\\}$ is constructed using the formulation
described in [40].

The full displacement field in the global problem is constructed using both
continuous and discontinuous parts as

$\textbf{u}^{h}_{G}=\underbrace{\sum_{a=1}^{n_{u}}u_{a}\boldsymbol{\eta}_{a}}_{\text{continuous
part}}+\underbrace{\sum_{b=1}^{n_{w}}w_{b}\boldsymbol{\phi}_{b}}_{\text{discontinuous
part}},$ (26)

where $\\{u_{a}\\},\\{w_{b}\\}$ are scalar degrees of freedom.

The flow equations (5) and (8) are discretized with a finite-volume method.
Piecewise-constant pressure fields are constructed for $p_{m}^{h}$ and
$p_{f}^{h}$ over the matrix mesh $\mathcal{T}$ and the fracture mesh
$\mathcal{F}$, respectively. Fluxes are computed using a two point flux
approximation. The interaction between the flow in the fracture and matrix is
effected via the embedded discrete fracture model (EDFM) [41, 35].

In terms of the temporal discretization, a backward Euler method is used
throughout. This gives rise to a fully-coupled system of nonlinear equations,
whose solution is obtained with a Newton method, in a monolithic fashion.

One limitation of the construction (26) is that the discontinuous enrichment
functions are not capable of representing a crack tip that terminates inside
of an element. As such, any new extension of the crack geometry has to
traverse from one side of a new element to another. The implementation of the
fracture flow solver requires all cells to have non-zero volumes. Accordingly,
new fracture cells are assigned a small aperture value ($w_{0}$). The total
discrete aperture of the cell is thus given by $w_{h}=w_{n}+w_{0}$, where
$w_{n}$ is the mechanical aperture that is consistent with the jump provided
by the displacement field. The minimum opening $w_{0}$ has been interpreted as
a representation of the roughness of the fracture surfaces, providing a
pathway for fluid flow even when the cracks are mechanically closed. See, for
example [31].

###  3.3 Local problem initialization

####  3.3.1 Construction of subdomain, submesh and damage-fixed nodes

A local subdomain of size $L\times L$ is constructed by simply centering it on
a global crack tip. The size $L$ is selected to be an even integer multiplier
of the global mesh spacing, leading to a square bounding box that conforms to
the background mesh. This choice is adopted for convenience in this work,
although other constructions are possible, such as in [27].

To obtain the submesh $\mathcal{T}_{L}$, we start by considering the
restriction of the global triangulation $\mathcal{T}$ to the subdomain
$\Omega_{L}$. The mesh for the local subdomain is constructed by uniformly
refining the set of elements in this restriction. This facilitates the
transfer of nodal data from the global to the local problem. The mesh size
$h_{local}$ in the local subdomain is chosen to be sufficiently small to
resolve the damage band, of size $\mathcal{O}(\ell)$. In this work, we use
$\ell/h_{local}\approx 4$.

Figure 3: Global domain $\Omega$ on the left and local domain $\Omega_{L}$ on
the right. The nodes in the subset $\mathcal{K}$, where $d$ is set to $1$ are
colored in yellow.

The crack is represented by prescribing $d=1$ in a set
$\mathcal{K}\in\mathcal{T}_{L}$ of nodes in the local subdomain, as shown in
Figure 3. This set is constructed in two steps. First, all elements of
$\mathcal{T}_{L}$ which are intersected by the global fracture triangulation
$\mathcal{F}$ are identified. Then, $\mathcal{K}$ is defined to be the set of
all nodes that belong to any of these intersected elements.

####  3.3.2 Transfer of global pressures to local mesh

Due to the different levels of resolution between the global and local
problems, we find it advantageous to transfer the pressure fields in a
particular manner. We note that in the global problem, the fracture pressures
$p_{f}^{h}$ are available at the cell centers of 1D finite volumes, and the
matrix pressures $p_{m}^{h}$ are available at the cell centers of 2D finite
volumes. These fields are transferred to quadrature points in the finite
element mesh for the local subproblem using the transfer operators
$\Pi^{\Omega_{L}}_{f}$ and $\Pi^{\Omega_{L}}_{m}$.

Figure 4: Illustration of global cells (in gray) in the neighborhood of a
local element, indicated in pink. The matrix pressure at a quadrature point in
the local element is calculated using an average of global pressures from
neighboring cells $\\{i\\}$, with weights corresponding to the distances
$r_{i}$ between the quadrature point and the cell center.

The operator $\Pi^{\Omega_{L}}_{f}$ that transfers the fracture pressure is
straightforward, and inspired by a similar operation described in Santillán et
al. [42]. Given any point ${\bf{x}}\in\Omega_{L}$, and a global fracture
pressure field $p_{f}^{h}$, $\Pi^{\Omega_{L}}_{f}p_{f}^{h}({\bf{x}})$ is
obtained by finding the closest cell of $\mathcal{F}$ to $x$ and taking the
value of $p_{f}^{h}$ at this cell. More precisely,

$\Pi^{\Omega_{L}}_{f}p_{f}^{h}({\bf{x}})=p_{f}^{h}(\underset{c\in\mathcal{F}}{\text{arg
min }}\text{dist}(c,{\bf{x}})).$ (27)

For the operator $\Pi^{\Omega_{L}}_{m}$ that transfers the matrix pressure
field, we use an averaging procedure. This has the effect of smoothing the
matrix pressure field at the resolution of the local mesh. At a quadrature
point ${\bf{x}}_{Q}\in\Omega_{L}$, the matrix pressure in the local domain is
obtained from a weighted average of the pressures in the global cells that
surround the point. Specifically,

$\Pi^{\Omega_{L}}_{m}p_{m}^{h}({\bf{x}}_{Q})=\dfrac{\sum_{i=1}^{s}r^{-1}_{i}p_{m}^{h}(c_{i})}{\sum_{i=1}^{s}r^{-1}_{i}},$
(28)

where $r_{i}$ denotes the distance between the quadrature point location and
the center of global cell $i$, as indicated in Figure 4. In the sum, all cells
that neighbor the global cell containing quadrature point ${\bf{x}}_{Q}$ are
used. In the particular case when a quadrature point happens to reside in the
center of the local element and $r_{1}=0$, the above sum is replaced with
$\Pi^{\Omega_{L}}_{m}p_{m}^{h}({\bf{x}}_{Q})=p^{h}_{m}(c_{1})$.

###  3.4 Local problem discretization

With the local subdomain properly identified and initialized, we now describe
the additional steps to discretize the displacement and damage fields in the
local subdomain and solve for their approximations. The governing equations
for the macro-force balance (22) and the damage evolution (23) are both
treated with the finite element method. The damage and displacement fields are
both approximated using four-node bilinear quadrilateral elements.

Let $\Omega_{L}$ denote the local domain. From the finite-element
approximation $\textbf{u}_{G}^{h}$ computed in the global problem, we extract
$\textbf{u}^{h}_{G}|_{\partial\Omega_{L}}$ and use it to constrain the
displacements on the boundary of the subproblem 222In the multi-resolution
method of Muixi et al. [28], the displacement boundary conditions near the
crack base were released. We did not find this to be necessary in our
approach.. As such, the trial space $\boldsymbol{\mathcal{U}}^{h}$ is given by

$\boldsymbol{\mathcal{U}}^{h}=\\{\textbf{u}_{L}^{h}\in
H^{1}(\Omega_{L})^{n}\mid\textbf{u}_{L}^{h}=\textbf{u}^{h}_{G}\text{ on
}\partial\Omega_{L}\\}.$ (29)

The function space $\mathcal{D}^{h}$ of admissible damage fields is given by

$\mathcal{D}^{h}=\\{d^{h}\in H^{1}(\Omega_{L})\mid d^{h}=1\text{ on
}\mathcal{K}\\},$ (30)

where $\mathcal{K}$ denotes the set of damage-fixed nodes that correspond to
the global crack at the beginning of the load step, as described in subsection
3.3.1.

The test spaces are given by
$\boldsymbol{\mathcal{W}}^{h}=H_{0}^{1}(\Omega_{L})^{n}$ for the displacements
and $\mathcal{C}^{h}=\\{c^{h}\in H^{1}(\Omega_{L})\mid c^{h}=0\text{ on
}\mathcal{K}\\}$ for the damage field. The Galerkin approximation to the
problem on the local subdomain is then:

Find $\textbf{u}_{L}^{h}\in\boldsymbol{\mathcal{U}}^{h}$ and
$d^{h}\in\mathcal{D}^{h}$, such that
$\forall\textbf{w}^{h}\in\boldsymbol{\mathcal{W}}^{h}$ and $\forall
c^{h}\in\mathcal{C}^{h}$,
$\displaystyle\left(\nabla\textbf{w}^{h},\boldsymbol{\sigma}_{L}^{h}\right)-\left(\textbf{w}^{h},\Pi^{\Omega_{L}}_{f}p_{f}^{h}\nabla
d^{h}\right)-\left(\textbf{w}^{h},\textbf{b}\right)-\left<\textbf{w}^{h},\textbf{t}\right>_{\partial_{t}\Omega}$
$\displaystyle=0,$ (31)
$\displaystyle\begin{split}\frac{2G_{c}^{h}\ell}{c_{0}}\left(\nabla
c^{h},\nabla
d^{h}\right)+\frac{G_{c}^{h}}{c_{0}\ell}\left(c^{h},\zeta^{\prime}(d^{h})\right)+\left(c^{h},2(d^{h}-1)\Psi^{<A>}(\textbf{u}_{L}^{h})\right)\\\
-\left(c^{h},\alpha\Pi^{\Omega_{L}}_{m}p_{m}^{h}\nabla\cdot\textbf{u}_{L}^{h}\right)+\left(\nabla
c^{h},\Pi^{\Omega_{L}}_{f}p_{f}^{h}\textbf{u}_{L}^{h}\right)&=0,\end{split}$
(32)

where we have used $(\cdot,\cdot)$ to denote the standard inner product in the
$L^{2}(\Omega)$ space, and $\left<\cdot,\cdot\right>$ to denote its
restriction on the boundary. In the above, $\boldsymbol{\sigma}^{h}_{L}$
denotes the discrete stress, given by

$\boldsymbol{\sigma}^{h}_{L}=\left((1-d^{h})^{2}\
\mathbb{C}:\boldsymbol{\epsilon}(\textbf{u}_{L}^{h})-(1-d^{h})\alpha
p_{m}^{h}\mathbb{I}\right).$ (33)

The discrete active strain energy $\Psi^{<A>}({\bf{u}}^{h}_{L})$ that serves
as a driver for the damage evolution can be constructed from the total strain
energy
$\Psi(\textbf{u})=(1/2)\boldsymbol{\epsilon}(\textbf{u}^{h}):\mathbb{C}:\boldsymbol{\epsilon}(\textbf{u}^{h})$
in a number of ways. Options include the spectral split [43] or the
volumetric-deviatoric split [44]. We employ the spectral split in this work,
unless otherwise indicated.

We employ an alternating minimization scheme to solve the coupled system of
equations (31)-(32). Convergence is measured with respect to the $L_{2}$-norms
of the change in the damage $d^{h}$ and displacement fields
$\textbf{u}_{L}^{h}$, using a relative tolerance of $10^{-4}$. With this
approach, each of the equations becomes linear with respect to its primary
variable, simplifying the solution process.

We note that (32) does not explicitly include the irreversibility constraint
(25). We have found this to be unnecessary in the multi-resolution scheme, as
the Dirichlet condition $d^{h}=1$ on the crack set $\mathcal{K}$ is sufficient
to prevent any healing of the fracture surface relative to the global crack.
In terms of the local dissipation function $\zeta(d)$, in this work we use
$\zeta(d)=d$ which corresponds to the AT-1 phase-field model of fracture. This
is selected due to the fact that it gives rise to a compactly supported damage
field and a fully elastic stage prior to damage initiation [45].

Finally, we note that phase-field models of fracture tend to give rise to a
mesh and regularization length dependent critical fracture energy that is
larger than $G_{c}$ [46]. To account for this, we use a discrete value of
$G_{c}^{h}$ that is obtained as a function of the local mesh spacing and
regularization length, as

$G_{c}^{h}=G_{c}\left(1+\dfrac{h_{local}}{c_{0}\ell}\right)^{-1},$ (34)

such that the effective critical fracture energy is very close to that of the
material.

###  3.5 Crack propagation: translating local damage updates into global crack
extensions

A key step of the algorithm concerns how changes to the damage in the local
subdomain are translated into updates to the crack geometry in the global
domain. We provide details of our algorithm for this procedure here. The
scheme is presented in the context of a single crack tip that is propagating
through the global domain, but the generalization to multiple cracks is
straightforward.

At each step in the solution algorithm, in the global problem, we keep track
of two geometric entities, namely the global crack tip and the global tip
element (Figure 5a). The global crack tip is defined by the interior endpoint
of the crack. The global tip element is the element that contains the tip on
one of its sides, but is not yet “cut” by the crack. In essence, it is the
global element just ahead of the propagating crack tip. In the local
subproblem, we identify a local version of the crack tip that is obtained from
the discrete damage field $d^{h}$ (Figure 5b and 5c). The algorithm to extend
the global crack geometry depends on the relative location of the global crack
tip, global tip element, and local crack tip, as described below.

(a)

(b)

(c)

Figure 5: (a): Illustration of the global crack tip and global tip element.
(b) Case when the local tip falls inside the tip element. (c) Case when the
local tip falls outside the tip element.

The enrichment strategy described in Section 3.2 requires global elements that
are completely cut by the crack geometry. As such, any extension of the crack
geometry at the global scale must correspond to the global tip element being
fractured. Accordingly, global crack propagation is triggered whenever the
local tip falls outside the tip element. In this case, the new global tip is
identified by connecting the current global tip and the local tip. Since the
local tip is outside the tip element, a new segment will intersect the
perimeter of the tip element exactly once (neglecting the obvious intersection
at the current global tip). This process is illustrated in Figure 6.

Any crack advance beyond this new tip location is neglected at this point, and
the algorithm returns to a new global solve with an updated crack geometry.

(a)

(b)

(c)

Figure 6: (a): Construction of the segment connecting global and local tips.
(b) New crack segment to be added. (c) Update of the global crack tip and
global tip element.

In general, the problem of extracting a sharp crack front from a diffuse
representation is not trivial. When the crack evolution involves complex
topological changes, it is particularly difficult[47]. In the numerical
examples studied in this manuscript, we take advantage of the relatively
simple fracture geometry and predicted crack patterns to simplify this
process. In particular, we first identify all elements in the local subdomain
with nodes whose damage values are all above a threshold $d_{tr}$. A similar
approach was proposed in [27]. The local crack tip is taken to be the center
of the element in this set that is farthest from the base of the crack in the
local subdomain. The threshold used for this process is taken to be
$d_{tr}=1-h_{local}/(2\ell)$, which is based on the estimate for a damage
field near a crack tip given in [25]. In essence, for a phase-field model of
fracture, this threshold identifies nodes that are expected to correspond to
the peaks of the discrete damage field. For a damage band resolved with a mesh
spacing of $\ell/h_{local}=4$, this gives rise to a threshold of
$d_{tr}=0.875$.

###  3.6 Algorithm summary

Having described the solution strategies for both the global and local
problems and the transfer of various quantities, we now detail the algorithm
(1) that couples the two problems together to simulate crack propagation.
Figure 7 provides an illustration of the algorithm. Within each time step
(outer loop), the algorithm employs an inner loop that allows the global
problem to be updated as soon as any large enough change in the crack geometry
is detected in the local subproblem. The inner loop is terminated when the
propagation step, described in subsection 3.5 does not identify any crack
advance. The construction with two nested loops can be viewed as an implicit
treatment of the fracture front position, which, according to Lecampion et al.
[2] tends to be more accurate and robust, permitting the use of larger time
steps.

Define initial and boundary conditions

$n=1$

$n_{F}=endStep$

while  _$n\leq n_{F}$_ do

$\mathcal{F}_{n}^{1}=\mathcal{F}_{n-1}$

for  _$1\leq k\leq maxIter$ 333The index $k$ is only a dummy variable for this
loop that searches for the correct fracture geometry at a given time step._ do

(1) Solve Global Problem and obtain
$\textbf{u}_{G}^{k,n},p_{m}^{k,n},p_{f}^{k,n}$

(2) Construct subdomain $\Omega_{L}$ and submesh $\mathcal{T}_{L}$. Identify
subset of cracked nodes $\mathcal{K}$.

(3) Prescribe local boundary conditions
$\textbf{u}^{h}_{L}|_{\partial\Omega_{L}}=\textbf{u}_{G}^{k,n}|_{\partial\Omega_{L}}$
and $d^{h}=1$ on $\mathcal{K}$.

(4) Construct local pressure fields $p_{m}=\Pi^{\Omega_{L}}_{m}(p_{m}^{k,n})$
and $p_{f}=\Pi^{\Omega_{L}}_{f}(p_{f}^{k,n})$

(5) Solve the Local Problem to obtain $d^{k}$.

(6) Use $d^{k}$ and the propagation step (subsection 3.5) to update discrete
fracture $\mathcal{F}_{n}^{k}$.

if  _$\mathcal{F}_{n}^{k+1}$ = $\mathcal{F}_{n}^{k}$_ then

n = n + 1

$\mathcal{F}_{n}=\mathcal{F}_{n}^{k}$

$\textbf{u}_{G}^{n}=\textbf{u}_{G}^{k,n}$

$p_{m}^{n}=p_{m}^{k,n}$

$p_{f}^{n}=p_{f}^{k,n}$

break

end if

end for

end while

Algorithm 1 Solution algorithm for multi-resolution hydraulic fracture Figure
7: Multi-resolution solution algorithm.

##  4 Numerical Results

We now present the results from various benchmark problems in fluid-driven
fracture propagation. The problems range from those in which flow is only
present within the cracks to fully coupled problems involving flow in both the
matrix and the evolving manifold that is the fracture geometry.

(a)

(b)

Figure 8: (a) Schematic of the KGD problem, inspired by [1]. (b) Computational
domain for the KGD problem (not to scale).

###  4.1 KGD problems

We begin with the well-known Khristianovic, Geertsma and De Klerk (KGD)
problems of hydraulic fracture [32, 33]. The problems concern the propagation
of a planar fracture in an infinite, impermeable domain, under the injection
of a viscous fluid at a constant rate (Figure 8).

The response of the system can be characterized by the dimensionless group
$\mathcal{K}$ [48]:

$\mathcal{K}=\dfrac{4G^{1/2}_{c}}{(6\pi^{2}E^{\prime}Q\mu)^{1/4}},$ (35)

where $E^{\prime}=E/(1-\nu^{2})$. Values of $\mathcal{K}>4$ correspond to the
“toughness dominated regime” in which crack growth is largely controlled by
the fracture toughness of the media. By contrast, when $\mathcal{K}<1$, the
fracture toughness is relatively small and the fluid viscosity is the main
factor controlling the speed of crack growth. In the following, we explore the
performance of the multi-resolution algorithm to simulate problems in both
regimes.

####  4.1.1 Toughness-dominated regime

The toughness dominated regime is characterized by the creation of new
fracture surfaces accounting for almost all of the energy dissipation. In this
scenario, the fluid can be considered to be inviscid, which leads to a
constant pressure distribution over the entire crack. The assumption of
inviscid flow reduces much of the complexity, allowing for the construction of
a simple analytical solution for this problem. Here, we consider a problem in
the toughness-dominated regime resulting from the material parameters and
settings given in Table 1. Using these values in (35), we obtain
$\mathcal{K}=6.54$.

Table 1: Material properties and problem parameters for the toughness-dominated KGD problem | Value | Unit
---|---|---
Young’s modulus ($E$) | 16.0 | GPa
Poisson’s ratio ($\nu$) | 0.18 | –
Fluid viscosity ($\mu$) | 1.0$\times 10^{-12}$ | GPa . s
Energy release rate ($G_{c}$) | 1.85$\times 10^{3}$ | N/m
Injection rate ($Q$) | 1.0$\times 10^{-3}$ | $\text{m}^{2}/\text{s}$
Initial crack size ($a_{0}$) | 4 | m

The analytical solution for this problem can be separated into two stages as a
function of time. In the first stage, the pressure builds linearly with time
and the crack does not propagate, as the pressure is below the critical
threshold $p_{cr}=(G_{c}E^{\prime}/\pi a_{0})^{1/2}$. The pressure then
reaches the critical value at $t=t_{cr}$, after which the crack begins to
propagate. In the second stage the propagation is stable, since the amount of
fluid injected is finite and crack propagation leads to a pressure drop as the
total space available for the fluid to occupy increases.

The solution for the crack length and pressure in the crack can be written as

$a(t)=\begin{cases}a_{0},&t\leq t_{cr},\\\
\left(\dfrac{E^{\prime}(Qt)^{2}}{\pi G_{c}}\right)^{1/3},&t\geq
t_{cr},\end{cases}$ (36)

$p(t)=\begin{cases}\dfrac{t}{t_{cr}}p_{cr},&\ t\leq t_{cr},\\\
\left(\dfrac{E^{\prime}G_{c}^{2}}{\pi Qt}\right)^{1/3},&t\geq
t_{cr},\end{cases}$ (37)

where $t_{cr}=(\pi G_{c}a_{0}^{3}/Q^{2}E^{\prime})^{1/2}$. A derivation of
this solution can be found in Yoshioka [25].

For simulations using the multi-resolution scheme, a computational domain of
size $W\times H=30\text{ m}\times 240\text{ m}$ is used (Figure 8(b)).
According to [49], this domain size should be sufficiently large to yield a
good approximation to an infinite plane. We first report results using a sub-
domain size of $L\times L=a_{0}\times a_{0}$. The sensitivity of the results
to the choice of sub-domain size will be discussed later in this subsection.

In what follows, we report results from a refinement study. In particular, we
report results for three different regularization lengths of decreasing value,
beginning with $\ell=a_{0}/15$. As the regularization length is decreased, we
maintain a mesh size in the subdomain of $h_{local}=\ell/{4}$. This allows the
finite-element approximation to sufficiently resolve the regularized fracture
band. We also note that the effective fracture toughness in the computational
problem depends on $h/\ell$ [25]. In the global domain, we maintain the mesh
spacing at the ratio of $h_{global}=3h_{local}$. For the coarsest global mesh,
this translates into roughly 20 elements over the span of the initial crack
geometry. This level of resolution is consistent with those found to be
sufficient for resolving pressurized cracks using the Embedded Finite Element
Method, in Cusini et al. [31].

The results obtained with the multi-resolution scheme are compared to the
analytical solutions 36 and 37 in Figures 9(a) and 9(b). We note the overall
good match between the simulation results and the analytical solution, as well
as the convergence of the results towards 36 and 37 when $\ell$ decreases.
Although this problem is relatively simple as the matrix is assumed to be
impermeable, it does test the coupling between the global and the local
domains and verifies that the phase-field method, even when used only in a
vicinity of the crack tip, can still provide accurate predictions of fracture
propagation.

(a)

(b)

Figure 9: Comparison of numerical results and analytical solution for the
toughness-dominated KGD problem: (a) crack length vs. time and (b) pressure
vs. time.

We now examine the sensitivity of the results to changes in the subdomain
size. This is accomplished by fixing the sizes of the global and local meshes
as well as the regularization length, and varying only $L$ in Figure 8(b). In
particular, we fix the regularization length to $\ell=0.13\text{ m}$ and vary
the sub-domain size between $15\ell$ and $45\ell$.

Figures 10(a) and 10(b) provide the results for the pressure and crack length
as a function of time, for the various choices of subdomain size. Table 2
shows the error in the computation of the crack length relative to the
analytical solution. The error is taken as an average over the time range,
starting at the beginning of propagation.

(a)

(b)

Figure 10: (a) Effects of subdomain size on predicted crack length, for reduced problem. (b) Effects of subdomain size on predicted pressure, for reduced problem.  Table 2: Effect of subdomain size on the relative error in the crack length.  Subdomain size | Relative error |
---|---|---
$L=a_{0}/2$ | 11.5% |
$L=a_{0}$ | 3.3% |
$L=3a_{0}/2$ | 1.1% |

The results indicate that the error reduces as the subdomain size is
increased. This is to be expected, as the subdomain mesh is at a higher
resolution than the global mesh. As a result, there is a trade-off between
accuracy and computational cost. In general, one should select a subdomain
size that strikes a balance between being large enough to capture the crack
evolution and small enough to render the overall calculation efficient.

####  4.1.2 Viscosity dominated KGD problem

We now consider a case in which the choice of parameters gives rise to
conditions wherein the energy dissipation is dominated by viscous dissipation.
In particular, we consider the material properties and problem parameters
provided in Table 3. Consistent with (35), these values give rise to
$\mathcal{K}=0.57$, a result that is clearly within the viscosity-dominated
regime.

Table 3: Material properties and problem parameters for the viscosity-dominated KGD problem | Value | Unit
---|---|---
Young’s modulus ($E$) | 0.17 | GPa
Poisson’s ratio ($\nu$) | 0.20 | –
Fluid viscosity ($\mu$) | 1.0$\times 10^{-10}$ | GPa . s
Energy release rate ($G_{c}$) | 21 | N/m
Injection rate ($Q$) | 1.0$\times 10^{-3}$ | $\text{m}^{2}/\text{s}$
Initial crack size ($a_{0}$) | 4.0 | m

In this scenario, the toughness of the medium can be neglected, and the
fracture evolution is effectively dictated by the motion of the fluid front.
In this case, the pressure varies with time and space over the length of the
fracture surface. Although this problem does not include a pressure field in
the matrix, it bears emphasis that it does require the transfer of the
pressure field in the fracture from the global scale to the phase-field
problem in the subdomain. As such, it serves as a test of that aspect of the
multi-resolution scheme.

The presence of dynamic terms in the fluid equation also implies the need for
a discretization with sufficiently accurate temporal resolution, and in what
follows we examine the sensitivity of the results to the choice of the time
step size.

Once again, consider a computational domain with dimensions $W\times
H=30\text{ m}\times 240\text{ m}$. The sub-domain size is taken as $L\times
L=a_{0}\times a_{0}$, and the regularization length is taken to be
$\ell=a_{0}/15$. In contrast to what we observed for the toughness regime, we
note that the results to the viscosity-dominated problem were found to be
largely insensitive to the choice of the regularization length, provided that
$\ell<a_{0}/10$. This is not surprising, given that the phase-field subproblem
only serves to propagate a crack that has effectively zero toughness.

In what follows, we report results using mesh spacings in the local and global
domains of $h_{local}=\ell/{4}$ and $h_{global}=3h_{local}$. As an initial
condition, the pressure field is prescribed to match the analytical solution.
We report on the temporal convergence of results obtained using fixed time
steps, beginning with $\text{dt}_{0}=0.5$s. As (discrete) negative pressures
arose in our simulations, the phase-field model was modified to account for a
tension-compression asymmetry. In particular, consistent with the model
described in Miehe et al. [43], the strain energy density was split into
active and inactive parts, and only the “tensile” part was degraded with the
damage.

(a)

(b)

Figure 11: Numerical results and analytical solution for the viscosity-
dominated KGD problem: (a) crack length vs. time; and (b) inlet pressure vs.
time.

Figure 11 compares the analytical solutions for the crack length and inlet
pressure to the model-based simulation results for a sequence of decreasing
time steps. The initial time ($t_{0}=8.5\text{ s}$), initial inlet pressure
($p_{0}=40\text{ kPa}$) and initial crack size ($a_{0}=4\text{ m}$) were used
to render the results dimensionless. The particular values of $t_{0}$ and
$p_{0}$ were chosen in order to make our initial state a snapshot of the
analytical solution. Overall, the good agreement between the results indicates
that the multi-resolution scheme is capable of handling a viscosity-dominated
case.

The results for the crack length show an excellent agreement with the
analytical solution as the time step is decreased. The results for the inlet
pressure (Figure 11b) appear to converge to a trajectory that is slightly
offset from the analytical solution as the time step is decreased. This
relatively small discrepancy can be explained by the presence of a minimum
aperture that is assigned to the newly-initialized finite volumes in the flow
solver. As explained in [7], this is necessary to preserve the stability of
the scheme when new crack segments are added. The recent work by Jin et al.
[50] proposes a method for removing this parameter by employing partially
fractured elements. The incorporation of similar modifications in the context
of the multi-resolution scheme is the subject of future work.

###  4.2 Poroelastic problem with coupled matrix-fracture flow

Figure 12: Geometry and notation for the coupled poroelastic-fracture problem.

We now examine a problem in which the matrix is permeable and there is a
coupling between the fluid flow in the fracture, fluid flow in the matrix, and
crack propagation. A schematic of the problem is shown in Figure 12. It
corresponds to a rectangular domain subject to the injection of a viscous
fluid at the mouth of an initial fracture. The fluid injection rate $Q$ is
assumed to be constant, and the displacement in the normal direction is fixed
on all sides of the domain.

Although this problem is relatively simple, to our knowledge an analytical
solution is not available. In what follows, we therefore compare our results
to those that are obtained using an existing phase-field method for hydraulic
fracture. In particular, we compare our results to those obtained using the
method described by Miehe and Mauthe [18, 19]. This is arguably one of the
simpler methods available for this class of problems, combining the equations
for phase-field fracture, Biot’s theory of linear poroelasticity in the
matrix, and lubrication theory for fluid flow within fractures.

The model assumes that a Darcy model of flow holds in the matrix, with
isotropic permeability. In the fractures, a lubrication flow model is assumed.
The transition between the two flow regimes is effected through the use of a
permeability tensor that varies as a function of the crack aperture and
orientation. In particular, the permeability tensor ${\boldsymbol{\kappa}}$ is
given by

$\boldsymbol{\kappa}=\boldsymbol{\kappa}_{0}+d^{\xi}\dfrac{w^{2}_{n}}{12}\left(\textbf{I}-\textbf{n}^{d}\otimes\textbf{n}^{d}\right),$
(38)

$w_{n}=h(\textbf{n}^{d}\cdot\bf{\epsilon}\cdot\textbf{n}^{d}),$ (39)

in which the normalized gradient $\textbf{n}^{d}=\nabla d/|\nabla d|$
approximates the normal to the crack plane. In the above,
$\boldsymbol{\kappa}_{0}$ is an isotropic part of the permeability that
accounts for the undamaged permeability of the matrix, $w_{n}$ is the crack’s
normal aperture, and $\xi$ is a weighting exponent. The weighting exponent
concentrates the effects of the second term in 38 to areas where $d\approx 1$.
Consistent with the work of Miehe and Mauthe [18, 19], we use $\xi=50$.

Table 4: Material properties for poroelastic problem | Value | Unit
---|---|---
Young’s modulus ($E$) | 16.0 | GPa
Poisson’s ratio ($\nu$) | 0.18 | –
Fluid viscosity ($\mu$) | 1.0$\times 10^{-12}$ | GPa . s
Energy release rate ($G_{c}$) | 3.67 | N/m
Biot coefficient ($\alpha$) | 0.79 | -
Rock permeability ($\boldsymbol{\kappa_{0}}$) | 1.0$\times 10^{-13}$ | $\text{m}^{2}$
Table 5: Problem settings | Value | Unit
---|---|---
Phase-Field reg. length ($\ell$) | 0.2 | m
Subdomain size ($L$) | 4.0 | m
Domain width ($W$) | 30 | m
Domain height ($H$) | 10 | m
Initial crack size ($a_{0}$) | 4.0 | m
Injection rate ($Q$) | 1.0$\times 10^{-3}$ | $\text{m}^{2}/\text{s}$

The material properties used in this problem are given in Table 4, while the
problem and model parameters are provided in Table 5. The results reported in
this section rely on the use of an eigen-decomposition of the strain energy,
as described in Miehe et al. [43]. The subdomain size for the multi-resolution
method was chosen to be proportional to the initial crack size, or
$a_{0}\times a_{0}$. The problem was discretized spatially using a uniform
mesh in the local domain with $h_{local}=\ell/4$, while
$h_{global}=3h_{local}$ was used in the global domain. For the temporal
discretization, a uniform time step of $\text{dt}=0.125\text{ s}$ was used.
This level of spatial and temporal discretization was found to yield results
that were sufficiently converged.

An important difference between the phase-field formulation used in Miehe and
Mauthe [18, 19] and that employed in the current multi-resolution method
concerns the use of pressure-dependent driving forces. We draw the reader’s
attention to the terms involving the pressure fields $p_{m}$ and $p_{f}$ in
the evolution equation (23) for the damage field. In our studies of the KGD
problem, we found these terms necessary to include in order to obtain
sufficiently accurate simulations. We refer to these terms as “driving
pressures” as the pressure fields contribute directly to the evolution of the
damage field. Importantly, the phase-field formulation of Miehe and Mauthe
[18, 19] does not include these terms in the damage evolution equation.
Accordingly, in what follows we find it useful to compare results from Miehe
and Mauthe [18, 19] to those obtained using our multi-resolution method with
and without the driving pressures.

(a)

(b)

Figure 13: Results from the fully coupled poroelastic-fracture problem for the
multi-resolution method (with and without driving pressures) and a standard
phase-field model: (a) crack size over time, (b) inlet pressure. The time for
crack propagation ($t_{ref}$) in the phase-field model and the corresponding
inlet pressure ($p_{ref}$) are used as references for the dimensionless charts
above.

(a)

(b)

Figure 14: Fields from the coupled poroelastic-fracture problem, taken at the
end of the simulation, contrasting results from the multi-resolution scheme
without the driving pressures(top row) to those obtained using a standard
phase-field model of fracture in poroelastic media (bottom row). (a) Contour
plots of the damage field. (b) Contour plots of the pressure field.

Results for the crack length and inlet pressure as a function of time are
provided in Figure 13. Contour plots of the damage and pressure fields for the
multi-resolution scheme and the full phase-field formulation are provided in
Figure 14. Despite the numerous differences in the models, the results
obtained with the multi-resolution compare very favorably to those obtained
using our implementation of the phase-field fracture model [19]. This is
particularly the case when the driving pressures are removed from the multi-
resolution method, such that the two phase-field models are as close as
possible. At early times the crack remains stationary, as the pressure at the
inlet increases. At some point the pressure near the crack tip reaches a
magnitude that is sufficient to give rise to crack propagation. After crack
propagation begins, the rate of pressure increase begins to decrease with
time, as crack growth allows for additional fluid to be accommodated within
the fracture. Clearly the phase-field subproblem with a frozen pressure field
that acts effectively as a body force is still capable of simulating crack
propagation, even when there is flow in the crack and the matrix. It bears
emphasis that the multi-resolution scheme does not need to rely on estimations
such as 38,39 to extract an approximation to the aperture from the regularized
crack geometry.

###  4.3 Propagation around an inclusion

Figure 15: Geometry for the inclusion problem.

We now consider a problem that gives rise to a non-planar crack evolution.
Specifically, we investigate the effects of a stiff inclusion on the
trajectory of a hydraulically-driven fracture in an impermeable medium. The
problem setup is shown in Figure 15. A circular inclusion of radius $r$ is
placed in a rectangular domain, at a distance $s$ from the left boundary, and
slightly offset from the axis of symmetry. The inclusion is assumed to have
properties that are identical to the matrix, with the exception of the Young’s
modulus. An initial crack of size $a_{0}$ is placed in the left boundary and a
fluid is injected at a constant rate $Q$. The complete set of material
properties and geometric parameters are listed in Tables 6 and 7.

Table 6: Material properties for inclusion problem | Value | Unit
---|---|---
Young’s modulus ($E$) | 9.0 | GPa
Poisson’s ratio ($\nu$) | 0.25 | –
Fluid viscosity ($\mu$) | 1.0$\times 10^{-12}$ | GPa . s
Energy release rate ($G_{c}$) | 2.5$\times 10^{6}$ | N/m
Table 7: Problem settings | Value | Unit
---|---|---
Phase-Field reg. length ($\ell$) | 0.4 | m
Subdomain size ($L$) | 4.0 | m
Domain width ($W$) | 20 | m
Domain height ($H$) | 40 | m
Initial crack size ($a_{0}$) | 4.0 | m
Injection rate ($Q$) | 1.2$\times 10^{-1}$ | $\text{m}^{2}/\text{s}$

Simulations using the multi-resolution scheme were conducted for three
different values of inclusion stiffness, with $h_{local}=\ell/6$ in the local
domain, and $h_{global}=3h_{local}$ in the global domain. The objective was to
test how the stiffness of the inclusion influenced the crack trajectory.
Consider the following two limiting cases. When the inclusion is just slightly
stiffer than the matrix, one would expect the crack trajectory to remain
straight. At the other extreme, when the inclusion is much stiffer, one would
expect the crack to propagate around it.

(a)

(b)

(c)

Figure 16: Crack paths for the inclusion problem: (a) 2X stiffer inclusion;
(b) 7X stiffer inclusion; and (c) 15X stiffer inclusion.

The crack paths predicted by the multi-resolution scheme are shown in Figure
16. The deformed meshes are shown with displacements exaggerated to highlight
the crack trajectories. As expected, for an inclusion with much higher
stiffness, the crack propagates around it. By contrast, when the stiffness of
the inclusion is relatively close to that of the matrix (2X case), the
trajectory remains nearly straight. Interestingly, for an intermediate value,
where the inclusion is 7 times stiffer than the matrix, the crack trajectory
gets deflected, but still goes through the inclusion.

The pressures at the injection point are plotted as a function of time in
Figure 17. Due to the very high toughness of the material, they are almost
uniform over the crack. For this problem, we did not observe much sensitivity
of the results to the choice of time step size. By analyzing these curves, we
can identify several different stages for this problem. At early times before
any crack propagation begins, we observe a linear increase in pressure with
time. Then, as we reach a certain pressure, which we denote as
$p_{ref}=60\text{ MPa}$, straight crack propagation starts and proceeds until
the tip approaches the inclusion. Due to the higher stiffness of the
inclusion, the crack arrests and pressure builds up again until it reaches a
level that is sufficient to allow propagation to continue. This is accompanied
by a pressure drop that is most pronounced for the higher-stiffness inclusion
cases.

This problem highlights several capabilities of the multi-resoltion scheme,
such as the simple handling of crack curving as well as changes in the crack
speed, including crack arrest. Although we are not aware of any analytical
solution for problems like this, the results we have obtained appear to make
sense qualitatively.

Figure 17: Inlet pressure as a function of time for the three inclusion
problems, corresponding to inclusions with different stiffnesses. The
reference values of the pressure $p_{ref}$ and time $t_{ref}$ were defined as
the minimum pressure that led to fracture growth and the associated time,
considering the simulated injection rate.

##  5 Summary and conclusions

This manuscript presents a new approach for developing model-based simulations
of hydraulic fracturing. The method relies on a global problem that captures
flow and deformation fields, and a local problem that captures crack growth
with the aid of a phase-field method. The two problems are coupled through the
transfer of displacement and pressure fields, as well as updates to the crack
geometry. This multi-resolution approach allows the phase-field method to be
used as a tool to update the crack geometry without the need to explicitly
reconstruct the crack aperture from the diffuse representation.

The accuracy of the multi-resolution approach is evaluated through several
numerical examples. These include the well-known KGD problem, for which a good
agreement with analytical solutions is obtained in both the toughness and
viscosity regimes. A problem with a curved crack trajectory around a stiff
inclusion also demonstrates the overall robustness of the approach. Several
areas for future work present themselves. These include an enhancement of the
algorithm that tracks the phase-field crack tip, as well as a generalization
to three-dimensional problems.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. This
research was supported by the Exascale Computing Project (17-SC-20-SC), a
collaborative effort of the U.S. Department of Energy Office of Science and
the National Nuclear Security Administration. The support is gratefully
acknowledged.

## References

  * [1] J. Adachi, E. Siebrits, A. Peirce, J. Desroches, Computer simulation of hydraulic fractures, Int. J. Rock Mech. Min. Sci. 44 (5) (2007) 739–757.
  * [2] B. Lecampion, A. Bunger, X. Zhang, Numerical methods for hydraulic fracture propagation: a review of recent trends, Journal of Natural Gas Science and Engineering 49 (2018) 66–83.
  * [3] M. W. McClure, C. A. Kang, A Three-Dimensional Reservoir, Wellbore, and Hydraulic Fracturing simulator that is Compositional and Thermal, Tracks Proppant and Water Solute Transport, Includes Non-Darcy and Non-Newtonian Flow, and Handles Fracture Closure, in: SPE Reservoir Simulation Conference, OnePetro, 2017.
  * [4] M. McClure, C. Kang, S. Medam, C. Hewson, ResFrac technical writeup, arXiv preprint arXiv:1804.02092 (2018).
  * [5] R. R. Settgast, S. M. Johnson, P. Fu, S. D. C. Walsh, F. Ryerson, Simulation of hydraulic fracture networks in three dimensions, Tech. rep., Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States) (2012).
  * [6] R. R. Settgast, S. M. Johnson, P. Fu, S. D. C. Walsh, Simulation of hydraulic fracture networks in three dimensions utilizing massively parallel computing resources, in: SPE/AAPG/SEG Unconventional Resources Technology Conference, OnePetro, 2014.
  * [7] R. R. Settgast, P. Fu, S. D. C. Walsh, J. A. White, C. Annavarapu, F. J. Ryerson, A fully coupled method for massively parallel simulation of hydraulically driven fractures in 3-dimensions, Int. J. Numer. Anal. Methods Geomech. 41 (5) (2017) 627–653.
  * [8] H. Zia, B. Lecampion, PyFrac: A planar 3D hydraulic fracture simulator, Comput. Phys. Commun. 255 (2020) 107368.
  * [9] P. Gupta, C. A. Duarte, Simulation of non-planar three-dimensional hydraulic fracture propagation, Int. J. Numer. Anal. Methods Geomech. 38 (13) (2014) 1397–1430.
  * [10] P. Gupta, C. A. Duarte, Coupled hydromechanical-fracture simulations of nonplanar three-dimensional hydraulic fracture propagation, Int. J. Numer. Anal. Methods Geomech. 42 (1) (2018) 143–180.
  * [11] N. Shauer, C. A. Duarte, A three-dimensional Generalized Finite Element Method for simultaneous propagation of multiple hydraulic fractures from a wellbore, Eng. Fract. Mech. 265 (2022) 108360.
  * [12] Y. Heider, A review on phase-field modeling of hydraulic fracturing, Eng. Fract. Mech. 253 (2021) 107881.
  * [13] G. A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (8) (1998) 1319–1342.
  * [14] B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids 48 (4) (2000) 797–826.
  * [15] A. Karma, D. A. Kessler, H. Levine, Phase-field model of mode III dynamic fracture, Phys. Rev. Lett. 87 (4) (2001) 045501.
  * [16] B. Bourdin, C. P. Chukwudozie, K. Yoshioka, A variational approach to the numerical simulation of hydraulic fracturing, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2012.
  * [17] M. F. Wheeler, T. Wick, W. Wollner, An augmented-lagrangian method for the phase-field approach for pressurized fractures, Comput. Methods Appl. Mech. Engrg. 271 (2014) 69–85.
  * [18] C. Miehe, S. Mauthe, S. Teichtmeister, 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 (2015) 186–217.
  * [19] C. Miehe, S. Mauthe, 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. Engrg. 304 (2016) 619–655.
  * [20] A. Mikelić, M. F. Wheeler, T. Wick, Phase-field modeling of a fluid-driven fracture in a poroelastic medium, Comput. Geosci. 19 (6) (2015) 1171–1195.
  * [21] A. Mikelic, M. F. Wheeler, T. Wick, A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium, Multiscale Modeling & Simulation 13 (1) (2015) 367–398.
  * [22] Z. A. Wilson, C. M. Landis, Phase-field modeling of hydraulic fracture, J. Mech. Phys. Solids 96 (2016) 264–290.
  * [23] H. C. Brinkman, A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles, Flow Turbul. Combust. 1 (1) (1949) 27–34.
  * [24] C. P. Chukwudozie, B. Bourdin, K. Yoshioka, A variational phase-field model for hydraulic fracturing in porous media, Comput. Methods Appl. Mech. Engrg. 347 (2019) 957–982.
  * [25] K. Yoshioka, D. Naumov, O. Kolditz, On crack opening computation in variational phase-field models for fracture, Comput. Methods Appl. Mech. Engrg. 369 (2020) 113210.
  * [26] R. J. M. Geelen, Y. Liu, J. E. Dolbow, A. Rodríguez-Ferran, An optimization-based phase-field method for continuous-discontinuous crack propagation, Internat. J. Numer. Methods Engrg. 116 (1) (2018) 1–20.
  * [27] B. Giovanardi, A. Scotti, L. Formaggia, A hybrid XFEM–Phase field (Xfield) method for crack propagation in brittle elastic materials, Comput. Methods Appl. Mech. Engrg. 320 (2017) 396–420.
  * [28] A. Muixí, O. Marco, A. Rodríguez-Ferran, S. Fernández-Méndez, A combined XFEM phase-field computational model for crack growth without remeshing, Comput. Mech. 67 (1) (2021) 231–249.
  * [29] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Internat. J. Numer. Methods Engrg. 46 (1) (1999) 131–150.
  * [30] Y. Sun, Z. Liu, X. Tang, A hybrid FEMM-Phase field method for fluid-driven fracture propagation in three dimension, Eng. Anal. Bound. Elem. 113 (2020) 40–54.
  * [31] M. Cusini, J. A. White, N. Castelletto, R. R. Settgast, Simulation of coupled multiphase flow and geomechanics in porous media with embedded discrete fractures, Int. J. Numer. Anal. Methods Geomech. 45 (5) (2021) 563–584.
  * [32] J. Geertsma, F. De Klerk, A rapid method of predicting width and extent of hydraulically induced fractures, J. Petrol. Technol. 21 (12) (1969) 1571–1581.
  * [33] A. K. Zheltov, Formation of Vertical Fractures by Means of Highly Viscous Liquid, in: 4th World Petroleum Congress, OnePetro, 1955.
  * [34] C. Landis, Z. Wilson, The Phase-field Modeling of Hydraulic Fracture, in: workshop: Variational Models of Fracture, 2016, http://www.birs.ca//workshops//2016/16w5090/files/Banff-Landis.pdf, visited on 2022-05-05.
  * [35] H. Hajibeygi, D. Karvounis, P. Jenny, A hierarchical fracture model for the iterative multiscale finite volume method, J. Comput. Phys. 230 (24) (2011) 8729–8743.
  * [36] M. A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (2) (1941) 155–164.
  * [37] L. Ambrosio, V. M. Tortorelli, Approximation of functional depending on jumps by elliptic functional via t-convergence, Comm. Pure Appl. Math. 43 (8) (1990) 999–1036.
  * [38] W. Karush, Minima of functions of several variables with inequalities as side constraints, M. Sc. Dissertation. Dept. of Mathematics, Univ. of Chicago (1939).
  * [39] H. W. Kuhn, A. W. Tucker, Nonlinear Programming, in: Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, University of California Press, 1951, pp. 481–492.
  * [40] C. Linder, F. Armero, Finite elements with embedded strong discontinuities for the modeling of failure in solids, Internat. J. Numer. Methods Engrg. 72 (12) (2007) 1391–1433.
  * [41] S. Lee, M. Lough, C. Jensen, Hierarchical modeling of flow in naturally fractured formations with multiple length scales, Water Resources Research 37 (3) (2001) 443–455.
  * [42] D. Santillán, R. Juanes, L. Cueto-Felgueroso, Phase field model of hydraulic fracturing in poroelastic media: Fracture propagation, arrest, and branching under fluid injection and extraction, Journal of Geophysical Research: Solid Earth 123 (3) (2018) 2127–2155.
  * [43] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Engrg. 199 (45-48) (2010) 2765–2778.
  * [44] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, J. Mech. Phys. Solids 57 (8) (2009) 1209–1229.
  * [45] K. Pham, H. Amor, J.-J. Marigo, C. Maurini, Gradient damage models and their use to approximate brittle fracture, Int. J. Damage Mech 20 (4) (2011) 618–652.
  * [46] B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to fracture, J. Elasticity 91 (1) (2008) 5–148.
  * [47] E. Tamayo-Mas, A. Rodríguez-Ferran, A medial-axis-based model for propagating cracks in a regularised bulk, International Journal for Numerical Methods in Engineering 101 (7) (2015) 489–520.
  * [48] E. Detournay, Mechanics of hydraulic fractures, Annu. Rev. Fluid Mech. 48 (2016) 311–339.
  * [49] M. Isida, Analysis of stress intensity factors for the tension of a centrally cracked strip with stiffened edges, Eng. Fract. Mech. 5 (3) (1973) 647–665.
  * [50] T. Jin, J. A. White, R. R. Settgast, A Robust Finite Element-Finite Volume Strategy for Viscosity-Dominated Hydraulic Fracture Propagation Using Asymptotic Tip Enrichment, Int. J. Numer. Anal. Methods Geomech. (2022 (under review)).