# A Darcy–Cahn–Hilliard model of multiphase fluid-driven fracture

Alexandre Guével  Yue Meng  Christian Peco  Ruben Juanes  John E. Dolbow

###### Abstract

A Darcy–Cahn–Hilliard model coupled with damage is developed to describe
multiphase-flow and fluid-driven fracturing in porous media. The model is
motivated by recent experimental observations in Hele–Shaw cells of the fluid-
driven fracturing of a synthetic porous medium with tunable fracture
resistance. The model is derived from continuum thermodynamics and employs
several simplifying assumptions, such as linear poroelasticity and viscous-
dominated flow. Two distinct phase fields are used to regularize the interface
between an invading and a defending fluid, as well as the ensuing damage. The
damage model is a cohesive version of a phase-field model for fracture, in
which model parameters allow for control over both nucleation and crack
growth. Model-based simulations with finite elements are then performed to
calibrate the model against recent experimental results. In particular, an
experimentally-inferred phase diagram differentiating two flow regimes of
porous invasion and fracturing is recovered. Finally, the model is employed to
explore the parameter space beyond experimental capabilities, giving rise to
the construction of an expanded phase diagram that suggests a new flow regime.

###### keywords:

multiphase flow , hydraulic fracturing , phase-field , porous media

††journal: Journal of the Mechanics and Physics of Solids

[inst1]organization=Department of Mechanical Engineering and Materials, Duke
University, city=Durham, postcode=27708, state=NC, country=USA

[inst2]organization=Department of Geosciences, Princeton University,
city=Princeton, postcode=08544, state=NJ, country=USA

[inst3]organization=Department of Engineering Science and Mechanics, Penn
State University, city=State College, postcode=16801, state=PA, country=USA

[inst4]organization=Department of Civil and Environmental Engineering,
Massachusetts Institute of Technology, city=Cambridge, postcode=02139,
state=MA, country=USA

##  1 Introduction

The past few decades have witnessed an increased focus on obtaining a better
understanding of multiphase flow in deformable porous media. This interest is
driven by the prevalence of such phenomena in both natural and engineered
systems, including enhanced geothermal energy, geological carbon storage, and
geo-hazards [1]. Models that account for various aspects of the process can be
extremely useful, as bench-scale experiments that mimic a wide range of
subsurface conditions have proven to be elusive. From a modeling standpoint,
this class of problems is nonetheless challenging, with difficulties stemming
from a three-way coupling between poromechanics, fracture mechanics, and
multiphase fluid flow. This work introduces a novel double phase-field model
in which both the fluid-fluid interface and the fracture surfaces are
regularized and distinct. This permits the use of the model to explore a wide
range of processes, including hydraulic fracturing, viscous fingering, and
their combination. The model is developed within a simplified poroelastic
setting in order to limit the number of undetermined parameters and keep the
model’s response physically tractable. This paves the way for direct
comparison with recent experiments examining multiphase-flow driven fracturing
[2, 3, 4]. Beyond its ability to explain these recent experimental
observations, the broad potential of the model is utilized to build phase
diagrams indicating both existing and new flow regimes.

The current model builds on a considerable amount of previous theoretical
developments in which aspects of the underlying physics were examined in
isolation or with some degree of coupling. This includes work on multiphase
flow in rigid porous media, regularized models of fracture, and coupled
poromechanics and fracture. In what follows, important theoretical works and
recent modeling approaches that informed the current approach are discussed.

To begin with, the topic of multiphase flow in rigid porous media is one with
a rich history. Here, emphasis is placed on the seminal work of Saffman and
Taylor [5], who studied the fluid-fluid displacement instability that occurs
when an invading fluid is less viscous than the defending fluid. This
phenomena has since become known as the viscous fingering instability.
Subsequently, Lenormand [6] identified some of the conditions for viscous
fingering to occur, in addition to stable invasion and capillary fingering.
More recently, Zhao et al. [7] have augmented the parameter space of
consideration by taking into account the role of wettability in multiphase
flow. For a discussion of other important efforts in modeling multiphase flow
through rigid porous media, the reader is referred to the recent paper by
Juanes et al. [1].

The modeling of multiphase flow often gravitates around a Cahn–Hilliard [8]
description of the interface between the fluids. Such a regularization
circumvents a shortcoming of sharp-interface approaches when it comes to
modeling large interface deformation, and, in particular, interface
instabilities. The associated general description of multiphase flow led to
the Navier–Stokes–Cahn–Hilliard model, introduced by Lowengrub and
Truskinovsky [9], which led to the Hele–Shaw–Cahn–Hilliard model [10] in the
case of viscous flow in Hele–Shaw cells. The latter was proven to
asymptotically converge to the corresponding Hele–Shaw sharp interface
problem. The Navier–Stokes–Cahn–Hilliard model was later formally derived from
a continuum thermodynamic framework by Abels et al. [11]. For that, they took
inspiration from Gurtin et al.’s work [12], formally deriving from the laws of
thermodynamics the previously ad-hoc Navier–Stokes–Cahn–Hilliard models. In
Lowengrub and Truskinovsky’s derivation [9], the binary fluid velocity was
mass-averaged and not divergence-free, whereas in Abels et al.’s derivation
[11], the velocity was volume-averaged and, advantageously, divergence-free.
Dede et al.[13] obtained the Hele–Shaw counterpart to the model of Abels et
al. [11], along with the proof of sharp interface convergence. Importantly,
Papatzacos [14], Cueto-Felgueroso and Juanes [15, 16], and others, have
applied such Cahn–Hilliard approaches to multiphase flow through porous media.

The potential for porous media to fracture complicates modeling efforts
considerably due to the coupling between fluid flow and fracturing process. In
essence, fluid pressure can drive the formation of new fracture surfaces, and
the fracture openings in turn provide preferential paths for fluid flow. Most
of the existing modeling efforts in this area have focused on hydraulic
fracture and have been limited to instances in which the media is saturated
with only one fluid, as opposed to multiphase flow. We refer the reader to the
paper of Chen et al. [17] for a recent review of models for hydraulic
fracture.

More recently, advances in fracture mechanics have driven new developments in
hydraulic fracture modeling. We focus in particular on models based on the
variational treatment of fracture and accompanying phase-field regularizations
[18, 19]. Early efforts to extend the basic phase-field for fracture approach
to coupling with fluid flow and hydraulic fracture were spearheaded by Miehe
and collaborators [20, 21]. Their approach was based on Biot’s [22] and
Coussy’s [23] poromechanics framework cast in a variational formulation. This
model opened avenues for embedding the well-established phase-field fracture
approach into poromechanics. In particular, it proposed modifying the
permeability to account for the flow enhancement due to fracturing. An
alternative approach was derived from the Theory of Porous Media [24, 25, 26],
which is the counterpart of Biot’s empirical theory [22] embedded in continuum
thermodynamics and the theory of mixtures [27, 28]. We refer in particular to
the models of Ehlers and Luo [29] and Wilson and Landis [30]. The main
advances consist of the use of a thermodynamically-consistent framework for
the model derivation and constitutive assumptions that permit a transition
from Darcy-type flow in the pores to Stokes-type flow within open fractures.
For further references on phase-field models of hydraulic fracture, the reader
is referred to a recent review by Heider [31].

The full three-way coupling of multiphase flow, poromechanics, and fracture
has been investigated to a lesser extent, with some notable recent advances.
For instance, Holtzman et al. [32] experimentally identified, for non-cohesive
granular media, the fracturing, capillary fingering, and viscous fingering
regimes in a phase diagram representing fracturing number versus capillary
number. On the modeling side, Lee et al. [33] introduced a hydraulic
fracturing model for two-phase flows using fracture phase-field modeling and
lubrication theory. In a similar vein, Heider and Sun [34] recently proposed a
multiphase hydraulic fracturing model but based on the Theory of Porous Media,
for both capillary and viscous flow. Importantly, neither the work of Lee et
al. [33] nor Heider and Sun [34] regularized the fluid-fluid interfaces. As a
result, these models are limited in their ability to account for aspects of
the interfacial fluid dynamics.

The closest analog to the model proposed in this work is that of Carrillo and
Bourg [35], who adopted a fluid mechanics approach. In particular, they
treated soft porous media as an intermediate state between solid and fluid,
and which was effectively modeled through visco-plastic rheology. The model of
Carrillo and Bourg [35] allows for a clear distinction between the pore and
continuum scales due to volume-averaging upscaling, as well as between viscous
and capillary effects. While the distinction between viscous and capillary
effects was made in the non-fracturing regime in the experimental study of
Holtzman et al. [32], the work of Carrillo and Bourg [35] focused on this
distinction in the fracturing regime. The latter also identified a regime of
non-invasive fracturing, where invasion only occurs in the cracks and not in
the porous medium. Importantly, the model of Carrillo and Bourg [35] relies on
a local yield stress criterion to govern fracture and is limited in its
ability to examine the potentially important role of fracture mechanics in
governing the response of the system. Such a local, stress-based criterion for
fracture is also likely to suffer from the kinds of spurious strain
localization effects that plague local damage models (see, e.g., Pijaudier-
Cabot and Bazant [36]). Finally, we note that much like the aforementioned
three-way coupling models, the fluid-fluid interfaces are not regularized in
the work of Carrillo and Bourg [35].

Lastly, attention is called to the recent work of Paulin et al. [37] who
developed a model for non-cohesive soft porous media which accounts for
regularized damage and finite-deformation kinematics. The novelty of this work
lies in envisioning the solid, liquid, and gas as three distinct phases, whose
interactions are governed by interfacial energies of Cahn–Hilliard type.
Although the model was developed for non-cohesive media, it employed a phase-
field for fracture regularization that prevents the transmission of traction
across fully damaged surfaces. Distinct from many of the aforementioned works,
the model of Paulin et al.[37] does not insist that damage be irreversible,
such that fracture healing can occur. As the model was developed within the
framework of continuum thermodynamics, it is thermodynamically consistent.

The approach described in this manuscript incorporates aspects that are
similar to some of the aforementioned literature, with several important
distinctions that are motivated by an interest in explaining recent
experiments of fluid-driven fracturing of cohesive media [4]. The
configuration consists of a monolayer of beads that are lightly cemented
together and placed into a Hele–Shaw cell, and saturated with a viscous
defending fluid. The injection of a far less viscous invading fluid into this
system gives rise to a clear coupling between fracturing and fluid invasion in
which the fracture resistance of the skeleton plays an important role in
delineating various regimes of the response. In the configuration examined in
this work, one can expect new features such as a non-trivial injection
pressure curve and the possibility of viscous fingering instabilities. The
former feature can be explained by the fact that both the solid’s cement and
the defending fluid’s viscosity act as a resistance to fracture and invasion
of the injected fluid.

Accordingly, for the fracture model, this work adopts a cohesive variation on
the well-established phase-field regularization for fracture. In particular,
the model described in [38] is used. Importantly, it incorporates an energetic
threshold for damage that is insensitive to the choice of regularization
length. To accommodate for fluid-fluid interface instabilities, a second phase
field is employed, representing the invading fluid saturation. This allows
merging the two-pronged effort described above, and, in particular, expanding
the space of potential flow regimes. The aforementioned works have identified
six different types of flow regimes: uniform invasion, capillary and viscous
fingering, invasive capillary and viscous fracturing, and non-invasive
fracturing. A supplementary decomposition could be made between cohesive and
non-cohesive fracturing. In this contribution, a new regime where viscous
fracturing and viscous fingering are combined is revealed. For that, and in
line with the experiments, attention is confined to viscous flows (large
capillary numbers), leaving aside the capillary regimes.

This manuscript is divided into two main sections, a theoretical derivation
followed by a numerical study in relation to the aforementioned Hele–Shaw
experiments. The modeling strategy is first described, yielding a
Darcy–Cahn–Hilliard model of fluid-driven fracture with multiphase flow. The
composite model is assembled through a continuum thermodynamic derivation to
ensure a consistent coupling of the different components. Namely, these
components are inspired from Gurtin et al. [12] and Cogswell and Szulczewski
[39] for the multiphase flow, Borja [40] for poromechanics, and Ehlers and Luo
[29] for hydraulic fracturing. The novelty in the fracture modeling will be to
employ a cohesive fracture model [38, 41, 42] which allows for controlling
both the onset of cracks and their propagation, but also accounts for rate
effects through a damage viscosity. The resulting system of coupled partial
differential equations is discretized using the finite-element method. Model-
based simulations are then employed to calibrate various parameters against
the experimental observations, and to reproduce a phase diagram representing
modified capillary number versus cement volume ratio. Finally, the parameter
space to build a new phase diagram that extends beyond the experimental
observations is explored, yielding a new dimensionless group of interest and a
new type of flow regime mixing fracturing and viscous fingering.

##  2 Derivation

###  2.1 Overall framework and approach

The main novelty in the model lies in the use of two distinct phase fields to
regularize geometric features of interest (see Figure 1). The first phase
field, denoted by $S$, corresponds to the invading fluid saturation, so that
$S=1$ in the invading fluid and $S=0$ in the defending fluid. Intermediate
values of $S$ correspond to the diffuse interface, of thickness $l_{S}$, which
regularizes the otherwise sharp interface and ensures a smooth transition from
0 to 1. The second phase field $d$ can be similarly described in that it
discriminates between fully cracked regions, where $d=1$, and fully intact
regions, where $d=0$. The thickness of the damage band is characterized by the
regularization length $l_{d}$.

Although the resulting phase-field equations will end up looking similar, the
saturation phase-field stems from a gradient energy regularization, following
the pioneering work of Cahn [8], whereas the damage phase-field theory stems
from a variational regularization of Griffith’s crack surface energy, as
introduced by Francfort and Marigo [18]. Both phase-field approaches are
reconciled by the configurational theory of Fried and Gurtin [43, 44, 45],
whereby phase-field equations are derived from a fundamental law, the micro-
force balance. The latter approach is followed in this work to derive the
evolution equations governing both phase fields.

Figure 1: Schematic of the modeling of the fluid-fluid interface and crack
surfaces through the phase fields $S$ (water saturation) and $d$ (damage),
respectively.

###  2.2 Notation and terminology

Consider a porous medium consisting of three components: one solid and two
immiscible pore fluids. In keeping with Coussy’s theory and terminology [23],
a porous medium is modeled as the superimposition of a solid skeleton
continuum (formed by the solid matrix and the connected pore space emptied of
fluid) and of a pore-filling (binary) fluid continuum. Both subscripts and
superscripts will be used to denote fields associated with one of the
components. In particular, $s$ is used to denote the solid component and
$\alpha=\\{1,2\\}$ is used to denote the fluid components, with $\alpha=1$
representing the invading fluid and $\alpha=2$ the defending fluid.

The volume fractions of the components are denoted by $\phi^{s}$, $\phi^{1}$,
and $\phi^{2}$, and satisfy $\phi^{s}+\phi^{1}+\phi^{2}=1$. Similarly, the
respective partial densities (at the macroscale) are denoted by $\rho^{s}$ and
$\rho^{\alpha}$, and the partial Cauchy stresses by $\boldsymbol{\sigma}^{s}$
and $\boldsymbol{\sigma}^{\alpha}$. The corresponding intrinsic quantities (at
the mesoscale) are denoted with subscripts. The volume fraction of each fluid
in the pores is described by the saturations $S^{\alpha}$, whereby
$S^{1}+S^{2}=1$.

Given the above notation, some key relationships can be summarized as:

$\displaystyle\phi^{s}$ $\displaystyle=1-\phi,$ $\displaystyle\phi^{\alpha}$
$\displaystyle=\phi S^{\alpha},$ $\displaystyle\phi^{f}$
$\displaystyle=\phi^{1}+\phi^{2}=\phi,$ (1) $\displaystyle\rho^{s}$
$\displaystyle=(1-\phi)\rho_{s},$ $\displaystyle\rho^{\alpha}$
$\displaystyle=\phi S^{\alpha}\rho_{\alpha},$ $\displaystyle\rho^{f}$
$\displaystyle=\rho^{1}+\rho^{2},$ $\displaystyle\boldsymbol{\sigma}^{s}$
$\displaystyle=(1-\phi)\boldsymbol{\sigma}_{s},$
$\displaystyle\boldsymbol{\sigma}^{\alpha}$ $\displaystyle=\phi
S^{\alpha}\boldsymbol{\sigma}_{\alpha},$
$\displaystyle\boldsymbol{\sigma}^{f}$
$\displaystyle=\boldsymbol{\sigma}^{1}+\boldsymbol{\sigma}^{2},$

where tensorial quantities are written in bold font.

###  2.3 Main assumptions

Consider the isothermal motion of a binary fluid through a deformable porous
skeleton, and assume that the binary fluid is incompressible, as well as its
two individual components, 1 and 2. For simplicity, further assume that the
difference in (constant) densities $\rho_{1}$ and $\rho_{2}$ is negligible.
The reader is referred to Abels et al. [11] for a derivation when this
assumption does not hold. Following Gurtin et al. [12], also assume that “the
momenta and kinetic energies of the constituents are negligible when computed
relative to the gross motion of the fluid.” This means in particular that the
binary fluid is tracked in time through the material time derivative
associated with the average velocity
$\boldsymbol{v}_{f}=S\boldsymbol{v}_{1}+(1-S)\boldsymbol{v}_{2}$. Furthermore,
assume that the capillary number is sufficiently high such that the pore
pressure $p$ is identical in both phases 1 and 2, and also that the flow is
sufficiently slow to neglect the fluid’s kinetic energy. In other words,
consider a viscous flow.

With regard to the porous skeleton, small deformations and deformation
gradients are assumed, so that the formulation falls in the framework of
linear poroelasticity. Unlike the fluid components, the solid matrix is
assumed to be compressible, so that its density is not constant. The resulting
continuum is approached through the Theory of Porous Media (see, e.g., Ehlers
and Luo [29] and references therein), that is, it stems from the
superimposition of the continua formed by the different solid and fluid
components.

###  2.4 Mass balance

To derive the mass balance equations for the three phases, the derivation here
follows Coussy [23] and Borja [40], combined with Gurtin et al. [12].

The partial mass fluxes associated with the densities $\rho^{1}$ and
$\rho^{2}$ are denoted by $\tilde{\boldsymbol{h}}^{1}$ and
$\tilde{\boldsymbol{h}}^{2}$, respectively. As these are internal to the fluid
mixture, they satisfy
$\tilde{\boldsymbol{h}}^{1}+\tilde{\boldsymbol{h}}^{2}=\boldsymbol{0}$. In
what follows, a superposed dot is used to denote the material time derivative
following the solid skeleton motion. The solid velocity $\boldsymbol{v}_{s}$
is therefore given by $\boldsymbol{v}_{s}=\dot{\boldsymbol{u}}$, where
$\boldsymbol{u}$ denotes the solid displacement.

With these definitions in hand, the balances for the solid skeleton and the
two fluid phases read

$\displaystyle\dot{\rho}^{s}+\rho^{s}\nabla\cdot\boldsymbol{v}_{s}=0,$ (2)
$\displaystyle\dot{\rho}^{\alpha}+\rho^{\alpha}\nabla\cdot\boldsymbol{v}_{s}+\nabla\cdot\rho_{\alpha}\boldsymbol{w}_{\alpha}=-\nabla\cdot\tilde{\boldsymbol{h}}^{\alpha},$

where $\nabla=\partial/\partial\boldsymbol{x}$ denotes the spatial gradient,
and
$\boldsymbol{w}_{\alpha}=\phi^{\alpha}(\boldsymbol{v}_{\alpha}-\boldsymbol{v}_{s})$
is the Darcy velocity of phase $\alpha$ with respect to the solid skeleton.

Following Borja’s formulation [40], the compressibility of the solid is
assumed to depend only on the pore pressure $p$. More precisely, at the matrix
scale (mesoscale), the intrinsic density $\rho_{s}$ depends only on $p$, and
at the skeleton scale (macroscale), the partial density
$(1-\phi)\rho_{s}=\rho^{s}$ depends only on the partial pressure $(1-\phi)p$.
The functional relationships for the matrix and the skeleton follow as,
respectively:

$\dot{\rho_{s}}=\frac{\rho_{s}}{K_{s}}\dot{p},\quad(1-\phi)\dot{p}=-K\nabla\cdot\boldsymbol{v}_{s},$
(3)

where $K_{s}=\rho_{s}p^{\prime}(\rho_{s})$ is the matrix bulk modulus (i.e.
the bulk modulus of the grains in the experiments studied in this work) and
$K$ is the skeleton (drained) bulk modulus, and $p^{\prime}$ denotes the
derivative of $p$. The solid mass balance (2)1 can then be written in the
equivalent following forms:

$\displaystyle\dot{\overline{1-\phi}}+(1-\phi)\left(\nabla\cdot\boldsymbol{v}_{s}+\frac{\dot{p}}{K_{s}}\right)=0,$
(4)
$\displaystyle\dot{\overline{1-\phi}}+(b-\phi)\nabla\cdot\boldsymbol{v}_{s}=0,$

where $b=1-K/K_{s}$ denotes Biot’s coefficient.

As for the mass balances of the fluid constituents, dividing (2)2 by the
constant densities yields

$\displaystyle\dot{\overline{\phi S_{1}}}+\phi
S_{1}\nabla\cdot\boldsymbol{v}_{s}+\nabla\cdot\boldsymbol{w}_{1}=-\nabla\cdot\boldsymbol{h}^{1},$
(5) $\displaystyle\dot{\overline{\phi S_{2}}}+\phi
S_{2}\nabla\cdot\boldsymbol{v}_{s}+\nabla\cdot\boldsymbol{w}_{2}=-\nabla\cdot\boldsymbol{h}^{2},$

where $\boldsymbol{h}^{\alpha}=\tilde{\boldsymbol{h}}^{\alpha}/\rho_{\alpha}$.
Summing up the two previous equations, the mass balance of the binary fluid
can be written as

$\dot{\phi}+\phi\nabla\cdot\boldsymbol{v}_{s}+\nabla\cdot\boldsymbol{w}=0,$
(6)

where
$\boldsymbol{w}=\boldsymbol{w}_{1}+\boldsymbol{w}_{2}=\phi\boldsymbol{v}_{f/s}$
is the Darcy velocity of the binary fluid and
$\boldsymbol{v}_{f/s}=S(\boldsymbol{v}_{1}-\boldsymbol{v}_{s})+(1-S)(\boldsymbol{v}_{2}-\boldsymbol{v}_{s})$,
the velocity of the binary fluid with respect to the solid, is defined here
for later use. Then, summing up the mass balances of the solid
(LABEL:eq:equivalent_solid_mass_balance)2 and of the binary fluid (6) yields
the mass balance of the mixture:

$b\nabla\cdot\boldsymbol{v}_{s}+\nabla\cdot\boldsymbol{w}=0.$ (7)

Finally, for later use, it will be helpful to have at hand some alternative
expressions of the mass balances. A superposed circle (e.g.,
$\overset{\circ}{a}$), is used to denote the material time derivative
following the binary fluid. The mass balances of each fluid component and of
the resulting binary fluid can thus be rewritten as

$\displaystyle\overset{\circ}{\overline{\phi
S_{\alpha}}}=-\nabla\cdot\boldsymbol{h}^{\alpha},$ (8)
$\displaystyle\overset{\circ}{\phi}=0,$

where the incompressiblity of each fluid component
($\nabla\cdot\boldsymbol{v}_{\alpha}=0$) and of the binary fluid
($\nabla\cdot(S\boldsymbol{v}_{1}+(1-S)\boldsymbol{v}_{2})=0$) have been used.

###  2.5 Force balances

The macro-force balances of the solid and binary fluid read

$\displaystyle\nabla\cdot\boldsymbol{\sigma}^{s}+\boldsymbol{f}^{fs}=0,$ (9)
$\displaystyle\nabla\cdot\boldsymbol{\sigma}^{f}+\boldsymbol{f}^{sf}=0,$

where $\boldsymbol{\sigma}^{s}$ and $\boldsymbol{\sigma}^{f}$ denote the
partial solid and fluid Cauchy stresses, respectively, while
$\boldsymbol{f}^{fs}$ and $\boldsymbol{f}^{sf}$ are the forces exerted by the
fluid upon the solid and vice-versa. The effects of external body forces like
gravity have been neglected. The forces satisfy the reciprocity relationship
$\boldsymbol{f}^{fs}+\boldsymbol{f}^{sf}=0$. They are also referred to as
momentum productions in Ehlers and Luo’s work [29]. They result from the
superimposition of the continua formed by the solid and the binary fluid.

The above form of the force balance for the binary fluid stems from neglecting
the relative momenta of the constituents, as in Gurtin et al.’s derivation
[12]. In particular, recall that the pore pressure is assumed to be the same
in both phases. As a result, the force balance for the mixture, described by
the total stress
$\boldsymbol{\sigma}=\boldsymbol{\sigma}^{s}+\boldsymbol{\sigma}^{f}$ reads

$\nabla\cdot\boldsymbol{\sigma}=0.$ (10)

A derivation similar to that in [12] would also show that the symmetry of
$\boldsymbol{\sigma}^{s}$ and $\boldsymbol{\sigma}^{f}$ follows from the
corresponding balances of angular momentum. The fluid stress tensor can be
written in terms of its shear component $\boldsymbol{\tau}$ and the pore
pressure $p$ as
$\boldsymbol{\sigma}^{f}=\phi(\boldsymbol{\tau}-p\boldsymbol{I})$. It bears
emphasis that, unlike in Gurtin et al.’s derivation [12], due to the
poromechanical coupling, the pore pressure $p$ is not indeterminate.

Finally, following the configurational theory of Fried and Gurtin [43, 44],
the phase fields $S$ and $d$ are governed by the following microforce
balances:

$\displaystyle\nabla\cdot\boldsymbol{\xi}_{d}+\pi_{d}=0,$ (11)
$\displaystyle\nabla\cdot\boldsymbol{\xi}_{S}+\pi_{S}=0,$

where $\boldsymbol{\xi}$ denotes the microstress that is energetically
conjugate to the gradient of the respective phase field, and $\pi$ the
microforce that is conjugate to the phase field.

###  2.6 Dissipation inequality

Within the framework of continuum thermodynamics, constitutive restrictions
for the governing equations are now invoked. The partial Helmholtz free
energies per unit volume for the solid and binary fluid are denoted by
$\psi^{s}$ and $\psi^{f}$, respectively. The chemical potential of constituent
$\alpha$ associated with the flux $\boldsymbol{h}^{\alpha}$ is denoted by
$\mu_{\alpha}$. Since the mass fluxes are complementary, in this section
attention is restricted to constituent 1 and $\boldsymbol{h}$ is used in place
of $\boldsymbol{h}^{1}$ and $\mu$ is used to denote the difference in chemical
potentials, i.e. $\mu=\mu_{1}-\mu_{2}$. Similarly, $S$ is used in place of
$S_{1}$.

The second law expressed in the form of the dissipation inequality requires
that, for a given control volume $\Omega_{0}$ in the reference configuration,
the rate at which the energy of $\Omega_{0}$ increases (following Coussy’s
formulation [23] and neglecting the kinetic energy),

$\int\limits_{\Omega_{0}}{\left((\dot{\psi}^{s}+\overset{\circ}{\psi^{f}}\right)}\
\text{d}V,$ (12)

is limited by the external work on $\Omega_{0}$,

$\int\limits_{\partial\Omega_{0}}{\boldsymbol{\sigma}^{s}\boldsymbol{n}\cdot\dot{\boldsymbol{u}}}\
\text{d}A+\int\limits_{\partial\Omega_{0}}{\dot{d}\boldsymbol{\xi}_{d}\cdot\boldsymbol{n}}\
\text{d}A+\int\limits_{\partial\Omega_{0}}{\boldsymbol{\sigma}^{f}\boldsymbol{n}\cdot\boldsymbol{v}_{f}}\
\text{d}A+\int\limits_{\partial\Omega_{0}}{\overset{\circ}{S}\boldsymbol{\xi}_{S}\cdot\boldsymbol{n}}\
\text{d}A,$ (13)

plus the rate at which energy is transported to $\Omega_{0}$ by diffusion
[12],

$-\int\limits_{\partial\Omega_{0}}{\mu_{\alpha}\boldsymbol{h}^{\alpha}\cdot
n}\ \text{d}A=-\int\limits_{\partial\Omega_{0}}{\mu\boldsymbol{h}\cdot n}\
\text{d}A.$ (14)

The local form of the dissipation inequality expressed above is then given by

$\dot{\psi}^{s}+\overset{\circ}{\psi^{f}}\leq{\boldsymbol{\sigma}^{s}:\nabla{\boldsymbol{v}_{s}}}+{\boldsymbol{v}_{s}\cdot\nabla\cdot\boldsymbol{\sigma}^{s}}+{\boldsymbol{\xi}_{d}\cdot\nabla\dot{d}}+{\dot{d}\nabla\cdot\boldsymbol{\xi}_{d}}+{\boldsymbol{\sigma}^{f}:\nabla\boldsymbol{v}_{f}}+{\boldsymbol{v}_{f}\cdot\nabla\cdot\boldsymbol{\sigma}^{f}}+\boldsymbol{\xi}_{S}\cdot\nabla{\overset{\circ}{S}}+\overset{\circ}{S}\nabla\cdot\boldsymbol{\xi}_{S}-\mu\nabla\cdot\boldsymbol{h}-\boldsymbol{h}\cdot\nabla\mu,$
(15)

where the linear poroelasticity approximation
$\partial/\partial\boldsymbol{X}\approx\partial/\partial\boldsymbol{x}$ has
been invoked to retain the spatial gradients. Using the force balances and the
mass balance of the invading fluid, and noting in particular that
$\overset{\circ}{\phi}=0$, one obtains

$\dot{\psi}^{s}+\overset{\circ}{\psi^{f}}\leq{\boldsymbol{\sigma}^{s}:\nabla{v_{s}}}-{\boldsymbol{v}_{s}\cdot\boldsymbol{f}^{fs}}+{\boldsymbol{\xi}_{d}\cdot\nabla\dot{d}}-{\pi_{d}\dot{d}}+{\boldsymbol{\sigma}^{f}:\nabla\boldsymbol{v}_{f}}-{\boldsymbol{v}_{f}\cdot\boldsymbol{f}^{sf}}+\boldsymbol{\xi}_{S}\cdot\nabla{\overset{\circ}{S}}-\pi_{S}{\overset{\circ}{S}}+\mu\phi\overset{\circ}{S}-\boldsymbol{h}\cdot\nabla\mu.$
(16)

The familiar poromechanical concept of total effective stress can be obtained
by adding the null term
$(\boldsymbol{\sigma}-\boldsymbol{\sigma}^{s}-\boldsymbol{\sigma}^{f}):\nabla\boldsymbol{v}_{s}$
to (16), yielding

$\dot{\psi}^{s}+\overset{\circ}{\psi^{f}}\leq{\boldsymbol{\sigma}:\nabla\boldsymbol{v}_{s}}+{\boldsymbol{\xi}_{d}\cdot\nabla\dot{d}}-{\pi_{d}\dot{d}}+{\boldsymbol{\sigma}^{f}:\nabla\boldsymbol{v}_{f/s}}-{\boldsymbol{v}_{f/s}\cdot\boldsymbol{f}^{sf}}+\boldsymbol{\xi}_{S}\cdot\nabla{\overset{\circ}{S}}+(\phi\mu-\pi_{S}){\overset{\circ}{S}}-\boldsymbol{h}\cdot\nabla\mu.$
(17)

Following the derivation by Ehlers and Luo [29], the pore pressure is then
employed as a Lagrange multiplier, multiplying (7) by $p$ and adding it to the
previous inequality to obtain:

$\dot{\psi}^{s}+\overset{\circ}{\psi^{f}}\leq{\boldsymbol{\sigma}^{\prime}:\nabla\boldsymbol{v}_{s}}+{\boldsymbol{\xi}_{d}\cdot\nabla\dot{d}}-{\pi_{d}\dot{d}}+{\phi\boldsymbol{\tau}:\nabla\boldsymbol{v}_{f/s}}+{p\boldsymbol{v}_{f/s}\cdot\nabla\phi}-{\boldsymbol{v}_{f/s}\cdot\boldsymbol{f}^{sf}}+\boldsymbol{\xi}_{S}\cdot\nabla\overset{\circ}{S}+(\phi\mu-\pi_{S}){\overset{\circ}{S}}-\boldsymbol{h}\cdot\nabla\mu,$
(18)

where $\boldsymbol{\sigma}^{\prime}=\boldsymbol{\sigma}+bp\boldsymbol{I}$
denotes the total effective stress. The last step consists in identifying the
energy-conjugates, that is, the products of forces and rates of state
variables. For that, one must carefully commute gradients and time
derivatives, especially when the former are spatial quantities whereas the
latter are material quantities, as here, in which case one must employ the
following commutator identity (see e.g. [12]):

$\nabla\left(\frac{\text{d}^{\alpha}\varphi}{\text{d}t}\right)=\frac{\text{d}^{\alpha}}{\text{d}t}\left(\nabla\varphi\right)+(\nabla\boldsymbol{v}_{\alpha})^{T}\nabla\varphi.$
(19)

Thereupon,

$\displaystyle\nabla\dot{d}=\dot{\overline{\nabla{d}}}+(\nabla\boldsymbol{v}_{s})^{T}\nabla{d}\approx\dot{\overline{\nabla{d}}},$
(20)
$\displaystyle\nabla\dot{\boldsymbol{u}}=\dot{\overline{\nabla\boldsymbol{u}}}+(\nabla\boldsymbol{v}_{s})^{T}\nabla\boldsymbol{u}\approx\dot{\overline{\nabla\boldsymbol{u}}},$
$\displaystyle\nabla\overset{\circ}{S}=\overset{\circ}{\overline{\nabla{S}}}+(\nabla\boldsymbol{v}_{f})^{T}\nabla{S},$

where the approximations in the two first equations stem from the linear
poroelasticity assumption that second-order terms in deformation are
negligible (see e.g. section 3 in MacMinn et al. [46]). In particular, assume
that $d$ varies like $\phi$ since their evolution is tightly coupled, as the
model-based simulations in this work will indicate. The local dissipation
inequality (18) then reads:

$\dot{\psi}^{s}+\overset{\circ}{\psi^{f}}\leq{(\boldsymbol{\sigma}^{\prime}+\nabla{S}\otimes\boldsymbol{\xi}_{S}):\dot{\overline{\nabla\boldsymbol{u}}}}+\boldsymbol{\xi}_{d}\cdot\dot{\overline{\nabla{d}}}-{\pi_{d}\dot{d}}+{(\phi\boldsymbol{\tau}+\nabla{S}\otimes\boldsymbol{\xi}_{S}):\nabla\boldsymbol{v}_{f/s}}+{(p\nabla\phi-\boldsymbol{f}^{sf})\cdot\boldsymbol{v}_{f/s}}+\boldsymbol{\xi}_{S}\cdot\overset{\circ}{\overline{\nabla{S}}}+(\phi\mu-\pi_{S})\overset{\circ}{S}-\boldsymbol{h}\cdot\nabla\mu.$
(21)

The change of variable $\hat{\mu}=\phi\mu$ is applied and the hat notation is
dropped thereafter for clarity. Identifying as state variables
$\\{\nabla\boldsymbol{u},d,\nabla{d},\dot{d},\boldsymbol{v}_{f/s},\boldsymbol{D},S,\nabla{S},\mu,\nabla\mu\\}$,
where
$\boldsymbol{D}=(\nabla\boldsymbol{v}_{f/s}+{\nabla\boldsymbol{v}_{f/s}}^{T})/2$
denotes the strain rate tensor for the binary fluid [12], and invoking Ehlers’
principle of phase separation [29], consider the constitutive dependencies:

$\displaystyle{\psi^{s}=\hat{\psi^{s}}(\nabla\boldsymbol{u},d,\nabla{d},\dot{d})=\hat{\psi^{s}}(\dots)_{s}},$
(22)
$\displaystyle{\psi^{f}=\hat{\psi^{f}}(\boldsymbol{v}_{f/s},\boldsymbol{D},S,\nabla{S},\mu,\nabla\mu)=\hat{\psi^{f}}(\dots)_{f}},$
$\displaystyle{\boldsymbol{\sigma}^{\prime}=\hat{\boldsymbol{\sigma}}^{\prime}(\dots)_{s}},$
$\displaystyle{\boldsymbol{\xi}_{d}=\hat{\boldsymbol{\xi}}_{d}(\dots)_{s}},$
$\displaystyle{\pi_{d}=\hat{\pi}_{d}(\dots)_{s}},$
$\displaystyle{\boldsymbol{\tau}=\hat{\boldsymbol{\tau}}(\dots)_{f}},$
$\displaystyle{\boldsymbol{\xi}_{S}=\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}},$
$\displaystyle{\pi_{S}=\hat{\pi}_{S}(\dots)_{f}},$
$\displaystyle{\boldsymbol{\tau}=\hat{\boldsymbol{\tau}}(\dots)_{f}},$
$\displaystyle{p=\hat{p}(\dots)_{f}},$
$\displaystyle{\boldsymbol{h}=\hat{\boldsymbol{h}}(\dots)_{f}}.$

This hereby-defined constitutive model satisfies the dissipation inequality
(21) if and only if

${\left(\hat{\boldsymbol{\sigma}}^{\prime}(\dots)_{s}+\nabla{S}\otimes\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}-\frac{\partial\hat{\psi}^{s}(\dots)_{s}}{\partial{\nabla\boldsymbol{u}}}\right):\dot{\overline{\nabla\boldsymbol{u}}}}+{\left(\hat{\boldsymbol{\xi}}_{d}(\dots)_{s}-\frac{\partial\hat{\psi}^{s}(\dots)_{s}}{\partial{\nabla{d}}}\right)\cdot\dot{\overline{\nabla{d}}}}-{\left(\hat{\pi}_{d}(\dots)_{s}+\frac{\partial\hat{\psi}^{s}(\dots)_{s}}{\partial{d}}\right)\dot{d}}-\frac{\partial\psi^{s}}{\partial\dot{d}}\ddot{d}+{\left(\phi\hat{\boldsymbol{\tau}}(\dots)_{f}+\nabla{S}\otimes\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}\right):\boldsymbol{D}}+{\mathcal{A}\left(\nabla{S}\otimes\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}\right):\boldsymbol{A}}+{\left(\hat{p}(\dots)_{f}\nabla\phi-\boldsymbol{f}^{sf}\right)\cdot\boldsymbol{v}_{f/s}}-\frac{\partial\psi^{f}}{\partial\boldsymbol{v}_{f/s}}\overset{\circ}{\boldsymbol{v}_{f/s}}-\frac{\partial\psi^{f}}{\partial\boldsymbol{D}}\overset{\circ}{\boldsymbol{D}}+{\left(\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}-\frac{\partial\hat{\psi}^{f}(\dots)_{f}}{\partial{\nabla{S}}}\right)\cdot\overset{\circ}{\overline{\nabla{S}}}}+\left(\mu-\hat{\pi}_{S}(\dots)_{f}-\frac{\partial\hat{\psi}^{f}(\dots)_{f}}{\partial{S}}\right)\overset{\circ}{S}-\hat{\boldsymbol{h}}(\dots)_{f}\cdot\nabla\mu-\frac{\partial\psi^{f}}{\partial\mu}\overset{\circ}{\mu}-\frac{\partial\psi^{f}}{\partial\nabla\mu}\overset{\circ}{\overline{\nabla\mu}}\geq
0,$ (23)

where the tensor $\nabla\boldsymbol{v}_{f/s}$ has been decomposed into its
symmetric $\boldsymbol{D}$ and antisymmetric part $\boldsymbol{A}$, and the
symmetry of $\boldsymbol{\tau}$ has been employed.

Now, since $\dot{\overline{\nabla\boldsymbol{u}}}$,
$\dot{\overline{\nabla{d}}}$, $\ddot{d}$, $\boldsymbol{A}$,
$\overset{\circ}{\boldsymbol{v}_{f/s}}$, $\overset{\circ}{\boldsymbol{D}}$,
$\overset{\circ}{\overline{\nabla{S}}}$, $\overset{\circ}{S}$,
$\overset{\circ}{\mu}$ and $\overset{\circ}{\overline{\nabla\mu}}$ appear
linearly in (23), employing the Coleman–Noll procedure [47] yields the
following constitutive restrictions:

$\displaystyle\frac{\partial\psi^{s}}{\partial\dot{d}}=0,\quad\frac{\partial\psi^{f}}{\partial\boldsymbol{v}_{f/s}}=0,\quad\frac{\partial\psi^{f}}{\partial\boldsymbol{D}}=0,\quad\frac{\partial\psi^{f}}{\partial\mu}=0,\quad\frac{\partial\psi^{f}}{\partial\nabla\mu}=0,$
(24)
$\displaystyle\mathcal{A}\left(\nabla{S}\otimes\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}\right)=0,$
$\displaystyle\hat{\boldsymbol{\sigma}}^{\prime}(\dots)_{s}+\nabla{S}\otimes\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}=\frac{\partial\hat{\psi}^{s}(\dots)_{s}}{\partial{\nabla\boldsymbol{u}}},$
$\displaystyle\hat{\boldsymbol{\xi}}_{d}(\dots)_{s}=\frac{\partial\hat{\psi}^{s}(\dots)_{s}}{\partial{\nabla{d}}},$
$\displaystyle\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}=\frac{\partial\hat{\psi}^{f}(\dots)_{f}}{\partial{\nabla{S}}},$
$\displaystyle\hat{\pi}_{S}(\dots)_{f}=\mu+\frac{\partial\hat{\psi}^{f}(\dots)_{f}}{\partial{S}},$

leaving the following reduced dissipation inequality:

$-{\left(\hat{\pi}_{d}(\dots)_{s}+\frac{\partial\hat{\psi}^{s}(\dots)_{s}}{\partial{d}}\right)\dot{d}}+{\left(\phi\hat{\boldsymbol{\tau}}(\dots)_{f}+\nabla{S}\otimes\hat{\boldsymbol{\xi}}_{S}(\dots)_{f}\right):\boldsymbol{D}}+{\left(\hat{p}(\dots)_{s}\nabla\phi-\boldsymbol{f}^{sf}\right)\cdot\boldsymbol{v}_{f/s}}-\hat{\boldsymbol{h}}(\dots)_{f}\cdot\nabla\mu\geq
0.$ (25)

From (LABEL:eq:constitutive_restrictions)2, one can deduce that
$\nabla{S}\otimes\boldsymbol{\xi}_{S}=\boldsymbol{\xi}_{S}\otimes\nabla{S}$,
which implies that
$\left(\nabla{S}\otimes\boldsymbol{\xi}_{S}\right)\boldsymbol{\xi}_{S}=\left(\boldsymbol{\xi}_{S}\otimes\nabla{S}\right)\boldsymbol{\xi}_{S}$,
i.e.
$\left(\nabla{S}\cdot\boldsymbol{\xi}_{S}\right)\boldsymbol{\xi}_{S}=\left(\boldsymbol{\xi}_{S}\cdot\boldsymbol{\xi}_{S}\right)\nabla{S}$:
the vectors $\boldsymbol{\xi}_{S}$ and $\nabla{S}$ are collinear. The
corresponding collinearity coefficient $\kappa$ is introduced and chosen
constant for simplicity, so that

$\boldsymbol{\xi}_{S}=\kappa\nabla{S}=\frac{\partial\psi^{f}}{\partial\nabla{S}}.$
(26)

The reasoning above is similar to that invoked by Abels et al. [11].
Furthermore, Gurtin’s invariance argument [44] implies that the constitutive
functions can depend on $\nabla\boldsymbol{u}$ only through the strain tensor
of the solid
$\boldsymbol{\epsilon}=(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{T})/2$, so
that $\nabla\boldsymbol{u}$ can be replaced here above by
$\boldsymbol{\epsilon}$. Finally, a sufficient condition to satisfy (25) is to
require that each term be non-negative:

$\displaystyle\beta\dot{d}=\nabla\cdot\frac{\partial\psi^{s}}{\partial{\nabla{d}}}-\frac{\partial\psi^{s}}{\partial{d}},$
(27)
$\displaystyle\phi\boldsymbol{\tau}+\kappa\nabla{S}\otimes\nabla{S}=C_{St}\boldsymbol{D},$
$\displaystyle p\nabla\phi-\boldsymbol{f}^{sf}=C_{Da}\boldsymbol{v}_{f/s},$
$\displaystyle\boldsymbol{h}=-M\nabla\mu,$

where $\beta\geq 0$ is the damage viscosity, $C_{St}\geq 0$ is a coefficient
proportional to the binary fluid viscosity coefficient (in keeping with
Newtonian fluids’ rheological law [11, 12]), $C_{Da}\geq 0$ is a coefficient
inversely proportional to the binary fluid mobility, and $M\geq 0$ is the
Cahn-Hilliard mobility. These quantities need not be constant, and their
dependencies will be specified later.

###  2.7 Specification of the free energy

Following the previous thermodynamic analysis, the solid free energy is of the
form $\hat{\psi}^{s}(\boldsymbol{\epsilon},d,\nabla{d})$ and is expressed in
the following well-established form (see, e.g., [38, 41, 42]):

$\hat{\psi}^{s}(\boldsymbol{\epsilon},d,\nabla{d})=g(d)\psi_{e}^{A}(\boldsymbol{\epsilon})+\psi_{e}^{I}(\boldsymbol{\epsilon})+G_{c}h(d,\nabla{d}),$
(28)

where $g(d)$ is a stiffness degradation function, with $g(0)=1$ and $g(1)=0$,
and $\psi_{e}^{A}$ and $\psi_{e}^{I}$ denote the active (tensile) and inactive
(compressive) parts of the elastic energy, respectively. The critical fracture
energy is denoted by $G_{c}$, and $h$ is the crack density function defined by

$h(d,\nabla{d})=\frac{1}{c_{0}l_{d}}\left(w(d)+l_{d}^{2}\lvert\nabla{d}\rvert^{2}\right),$
(29)

where $c_{0}>0$ is a normalization constant, $w(d)$ is the local dissipation
function [48], and $l_{d}$ is the regularization length for the phase field
$d$.

The fluid free energy is of the form $\hat{\psi}^{f}(S,\nabla{S})$ and is
expressed following Cogswell and Szulczewski’s formulation [39] as:

$\hat{\psi}^{f}(S,\nabla{S})=\frac{12\gamma}{l_{S}}\left(f(S)+\frac{l_{S}^{2}}{16}\lvert\nabla{S}\rvert^{2}\right),$
(30)

where $\gamma$ is the surface tension between the two fluid constituents,
$l_{S}$ is the regularization length for the phase field $S$, and
$f(S)=S^{2}(1-S)^{2}$ is chosen as a double-well potential. Thereby, the
gradient coefficient $\kappa$ discovered in the thermodynamic analysis (26) is
given by $\kappa=3\gamma l_{S}/2$.

###  2.8 Fracture phase-field model

Different choices of the functions $g(d)$ and $w(d)$ lead to different
fracture phase-field models. Here, we adopt a cohesive fracture model [38, 41,
42] so as to control both crack nucleation and propagation. The former is
controlled by a nucleation energy $\psi_{c}$, acting as a damage threshold,
and the latter by the critical fracture energy $G_{c}$. Importantly, such a
model allows the regularization length $l_{d}$ to be prescribed independently
of the fracture properties, provided it satisfies an upper bound. Following
[42], the functions $g(d)$ and $w(d)$ are taken to be

$\displaystyle g(d)=\frac{(1-d)^{2}}{(1-d)^{2}+md(1+d)},$ (31) $\displaystyle
w(d)=d,$

where $m=G_{c}/c_{0}l_{d}\psi_{c}$ and
$c_{0}=4\int_{0}^{1}{\sqrt{w(d)}\text{d}d}=8/3$. To ensure the aforementioned
advantageous properties of the cohesive model, the condition $m\geq 1$ should
be satisfied [38], which yields the following upper bound for the
regularization length:

$l_{d}\leq\frac{3G_{c}}{8\psi_{c}}.$ (32)

The fracture phase-field evolution equation is then obtained from
(LABEL:eq:solution_dissipation_inequality)1, (28) and (32):

$\displaystyle\beta\dot{d}$
$\displaystyle=\frac{G_{c}}{c_{0}l_{d}}\left(2l_{d}^{2}\triangle{d}-w^{\prime}(d)\right)-g^{\prime}(d,\psi_{c})\psi_{e}^{A},$
(33)
$\displaystyle=\frac{3G_{c}}{8l_{d}}\left(2l_{d}^{2}\triangle{d}-1\right)-g^{\prime}(d,\psi_{c})\psi_{e}^{A},$

where the dependence of the degradation function $g$ on the the nucleation
energy $\psi_{c}$ has been emphasized. The term $\beta\dot{d}$ accounts for
the overall rate-dependency of fracture. A source of this rate-dependency is
the viscous effect experimentally observed during fracturing (see Figure 5b),
due to the latent cement that slows down the separation of the beads. A second
potential source is the frictional contact between the beads and the two
plates of the Hele–Shaw cell.

To complete the phase-field for the fracture part of the model, the
irreversibility condition $\dot{d}\geq 0$ is enforced. Finally, the spectral
decomposition of the elastic energy, as described in the work of Miehe et al.
[49], is invoked:

$\displaystyle\psi_{e}^{A}(\boldsymbol{\epsilon})=\frac{\lambda_{s}}{2}\langle\epsilon_{v}\rangle^{2}+\mu_{s}\sum_{i=1}^{3}\langle\epsilon_{i}\rangle^{2},$
(34)
$\displaystyle\psi_{e}^{I}(\boldsymbol{\epsilon})=-\frac{\lambda_{s}}{2}\langle-\epsilon_{v}\rangle^{2}-\mu_{s}\sum_{i=1}^{3}\langle-\epsilon_{i}\rangle^{2},$

where $\langle\cdot\rangle$ denotes the Macaulay brackets and $\epsilon_{i}$
the principal strains.

###  2.9 Binary fluid evolution equations

In light of the experimental conditions examined in this work, the flow in the
cracks is assumed to take the form of a Poiseuille flow between the Hele–Shaw
cell plates, which are maintained at a constant separation. More complex
models for the fluid flow in the cracks can be found in the works of Wilson
and Landis [30] or Ehlers and Luo [29], for example. Consistent with the
simplifying assumption employed in this work, $C_{St}$ is assumed to vanish
and only the terms describing a Darcy flow are maintained, such that
(LABEL:eq:solution_dissipation_inequality)2,3 can be injected into the fluid
force balance (9)2 to obtain:

$C_{Da}\boldsymbol{v}_{f/s}=-\phi\nabla{p}-\nabla\cdot(\kappa\nabla{S}\otimes\nabla{S}).$
(35)

The fluid mobilities are noted $\lambda_{\alpha}$. Upon choosing
$C_{Da}=\phi^{2}/(\lambda_{1}+\lambda_{2})$ and recalling that
$\boldsymbol{w}=\phi\boldsymbol{v}_{f/s}=\boldsymbol{w}_{1}+\boldsymbol{w}_{2}$,
a sufficient condition to describe the flow of the two constituents within the
current thermodynamic framework is

$\boldsymbol{w}_{\alpha}=-\lambda_{\alpha}\left(\nabla{p}+\phi^{-1}\nabla\cdot(\kappa\nabla{S}\otimes\nabla{S})\right),$
(36)

which is Darcy’s law augmented by a capillary (or Korteweg) stress [12]. More
specifically, the fluid mobilities are given by

$\lambda_{\alpha}=\frac{k(\phi)k_{r\alpha}(S_{\alpha})}{\eta_{\alpha}},$ (37)

where $k$ denotes the absolute permeability, $k_{r\alpha}$ the relative
permeability of phase $\alpha$, and $\eta_{\alpha}$ the dynamic viscosity of
phase $\alpha$. The absolute permeability is expressed as a function of $\phi$
following Kozeny–Carman’s formula (see, e.g., [23]), through

$k(\phi)=\frac{\phi^{3}}{(1-\phi)^{2}}\bar{k},$ (38)

where $\bar{k}$ is the intrinsic permeability. Importantly for the coupling
between poromechanics and fracturing, the porosity is expected to reach
$\phi=1$ in the cracks, whereby the permeability $k$ is expected to diverge.
This potential issue is circumvented by imposing an upper bound
$k_{max}=h^{2}/12$, where $h$ denotes the height of the Hele–Shaw cell, in
line with the aforementioned assumption of Poiseuille flow within the crack.
Thereby, the combination of Darcy flow in the intact porous skeleton and
Poiseuille flow in the cracks is obtained through imposing
$k=\text{min}(k(\phi),k_{max})$.

As for the relative permeabilities, the model of Fourar and Lenormand [50]
(also used in Cueto–Felgueroso and Juanes’s model [51]) derived for viscous
coupling in fractures is employed in this work:

$\displaystyle k_{r1}=S_{1}^{3}+\frac{3}{2\text{M}}S_{1}(1-S_{1})(1+S_{1}),$
(39) $\displaystyle k_{r2}=\frac{1}{2}S_{2}^{2}(3-S_{2}),$

where constituent 1 is chosen to be the non-wetting (invading) phase and
constituent 2 is chosen to be the wetting (defending) phase, and
$\text{M}=\eta_{2}/\eta_{1}$ is the viscosity ratio.

The Cahn–Hilliard mobility $M$ in (LABEL:eq:solution_dissipation_inequality)4
is then specified following Cogswell and Szulczewski’s model [39]:

$M=\frac{k}{\overline{\eta}},$ (40)

where $\overline{\eta}=2\eta_{1}\eta_{2}/(\eta_{1}+\eta_{2})$ is the harmonic
average of the viscosities of the two fluid constituents.

###  2.10 Model summary

To conclude this theoretical section, a summary of the evolution equations in
the model that will be discretized is provided below
(LABEL:eq:summary_equations). Since the two constituents are directly related
through $S_{1}+S_{2}=1$ and
$\boldsymbol{h}^{1}+\boldsymbol{h}^{2}=\boldsymbol{0}$, only the evolution of
constituent 1 (the invading fluid) is described; recall the notation
$S=S_{1}$. The model therefore consists of the mass balance for the invading
fluid, the total mass balance for the binary fluid, the solid mass balance,
the total macro-force balance, and the damage evolution (through the damage
micro-force balance), which are respectively given by

$\boxed{\begin{aligned} &\dot{\overline{\phi
S}}+\nabla\cdot\boldsymbol{w}_{1}(p,\phi)+\phi S\dot{\epsilon}_{v}=\nabla\cdot
M(\phi)\nabla\mu(S),\\\
&\dot{\phi}+\nabla\cdot\boldsymbol{w}(p,\phi)+\phi\dot{\epsilon}_{v}=0,\\\
&\phi=\phi_{0}+(1-\phi_{0})\epsilon_{v}+\frac{1-\phi_{0}}{K_{s}}p,\\\
&\nabla\cdot\left(\boldsymbol{\sigma}^{\prime}(\boldsymbol{\epsilon},S)-bp\boldsymbol{I}\right)=0,\\\
&\beta\dot{d}=\frac{3G_{c}}{8l_{d}}\left(2l_{d}^{2}\triangle{d}-1\right)-g^{\prime}(d,\psi_{c})\psi_{e}^{A},\quad\dot{d}\geq
0.\end{aligned}}$ (41)

The porosity evolution equation (LABEL:eq:summary_equations)3 derives from
integrating the solid mass balance (LABEL:eq:equivalent_solid_mass_balance)1
(also see equation (3.34) in Borja [40]) and linearizing the resulting
exponential function around the reference state
$\\{\epsilon_{v}=0,p/K_{s}=0\\}$. This system of equations
(LABEL:eq:summary_equations) will be solved for $\mu$, $p$, $\phi$,
$\boldsymbol{u}$, and $d$, respectively. It is closed with initial and
boundary conditions, and the following system of constitutive equations, given
by (36), (LABEL:eq:constitutive_restrictions)6, (11)2, (30),
(LABEL:eq:constitutive_restrictions)3, and (28), as:

$\displaystyle\boldsymbol{w}_{\alpha}=-\lambda_{\alpha}\left(\nabla{p}+\phi^{-1}\nabla\cdot(\kappa\nabla{S}\otimes\nabla{S})\right),$
(42)
$\displaystyle\mu=\frac{12\gamma}{l_{S}}\left(f^{\prime}(S)-\frac{l_{S}^{2}}{8}\triangle
S\right),$
$\displaystyle\boldsymbol{\sigma}^{\prime}=g(d)\frac{\partial\psi_{e}^{A}}{\partial\boldsymbol{\epsilon}}+\frac{\partial\psi_{e}^{I}}{\partial\boldsymbol{\epsilon}}-\kappa\nabla{S}\otimes\nabla{S}.$

It bears emphasis that the third terms in (LABEL:eq:summary_equations)1 and
(LABEL:eq:summary_equations)2 are of second order in terms of the deformation,
which may appear at odds with the assumption of linear poroelasticity. Through
extensive testing, the retention of these higher order terms was found to
significantly improve the calibration of model-based simulations to
experimental observations. In particular, these terms permit the model to
better represent the apparent coupling between the evolution of the saturation
front and fracture propagation.

###  2.11 Dimensionless form

The governing equations can be rendered dimensionless with the assistance of a
characteristic length $l_{0}$, velocity $v_{0}$, fluid viscosity
$\eta_{0}=\eta_{2}$, permeability $k_{0}$, and surface tension
$\gamma_{0}=\gamma$. It is also useful to define the characteristic pressure
$p_{0}=v_{0}\eta_{0}/l_{0}$, fluid mobility $\lambda_{0}=k_{0}/\eta_{0}$,
Cahn–Hilliard mobility $M_{0}=k_{0}/\eta_{0}=\lambda_{0}$, and chemical
potential $\mu_{0}=\gamma_{0}/l_{0}$. In the following, dimensionless
quantities are noted with an asterisk. The previous system
(LABEL:eq:summary_equations) then becomes

$\displaystyle\frac{v_{0}}{l_{0}}\dot{\overline{\phi
S}}^{*}+\frac{v_{0}}{l_{0}}\nabla^{*}\cdot\boldsymbol{w}_{1}^{*}+\frac{v_{0}}{l_{0}}\phi
S_{1}\dot{\epsilon}_{v}^{*}=\frac{M_{0}\gamma_{0}}{l_{0}^{3}}\nabla^{*}\cdot
M^{*}\nabla^{*}\mu^{*},$ (43)
$\displaystyle\frac{v_{0}}{l_{0}}\dot{\phi}+\frac{v_{0}}{l_{0}}\nabla^{*}\cdot\boldsymbol{w}^{*}+\frac{v_{0}}{l_{0}}\phi\dot{\epsilon}_{v}^{*}=0,$
$\displaystyle\phi=\phi_{0}+(1-\phi_{0})\epsilon_{v}+\frac{1-\phi_{0}}{p_{0}K_{s}^{*}}p_{0}p^{*},$
$\displaystyle\frac{p_{0}}{l_{0}}\nabla\cdot\left({\boldsymbol{\sigma}^{\prime}}^{*}-bp^{*}\boldsymbol{I}\right)=0,$
$\displaystyle\frac{v_{0}\beta}{l_{0}}\dot{d}^{*}=\frac{3G_{c}}{8l_{d}}\left(2\frac{l_{d}^{2}}{l_{0}^{2}}\triangle^{*}{d}-1\right)-g^{\prime}(d,\psi_{c})\psi_{e}^{A},\quad\dot{d}^{*}\geq
0.$

Similarly, the dimensionless forms of the secondary equations
(LABEL:eq:summary_equations_secondary) read:

$\displaystyle\boldsymbol{w}_{\alpha}^{*}=-\lambda_{\alpha}^{*}\left(\frac{\lambda_{0}p_{0}}{v_{0}l_{0}}\nabla^{*}{p^{*}}+\frac{\lambda_{0}\kappa}{v_{0}l_{0}^{3}}\phi^{-1}\nabla\cdot(\nabla{S}\otimes\nabla{S})\right),$
(44)
$\displaystyle\frac{\gamma_{0}}{l_{0}}\mu^{*}=\frac{12\gamma_{0}}{l_{S}}\left(f^{\prime}(S)-\frac{l_{S}^{2}}{8l_{0}^{2}}\triangle
S\right),$ $\displaystyle
p_{0}{\boldsymbol{\sigma}^{\prime}}^{*}=p_{0}g(d)\frac{\partial{\psi_{e}^{A}}^{*}}{\partial\boldsymbol{\epsilon}}+p_{0}\frac{\partial{\psi_{e}^{I}}^{*}}{\partial\boldsymbol{\epsilon}}-\frac{\kappa}{l_{0}^{2}}\nabla^{*}{S}\otimes\nabla^{*}{S}.$

The final system of equations in dimensionless form is then given by:

$\displaystyle\dot{\overline{\phi S}}+\nabla\cdot\boldsymbol{w}_{1}+\phi
S_{1}\dot{\epsilon}_{v}=\text{Pe}^{-1}\nabla\cdot M\nabla\mu,$ (45)
$\displaystyle\dot{\phi}+\nabla\cdot\boldsymbol{w}+\phi\dot{\epsilon}_{v}=0,$
$\displaystyle\phi=\phi_{0}+(1-\phi_{0})\epsilon_{v}+\frac{1-\phi_{0}}{K_{s}}p,$
$\displaystyle\nabla\cdot\left({\boldsymbol{\sigma}^{\prime}}-bp\boldsymbol{I}\right)=0,$
$\displaystyle\text{B}\dot{d}=2\text{L\textsubscript{d}}^{2}\triangle{d}-1-g^{\prime}(d,\psi_{c})\text{D\textsubscript{f}},\quad\dot{d}\geq
0,$

along with

$\displaystyle\boldsymbol{w}_{\alpha}=-\lambda_{\alpha}\text{Da}\left(\nabla{p}+\frac{3}{4}\text{Ca}^{-1}\text{L\textsubscript{S}}\phi^{-1}\nabla\cdot(\nabla{S}\otimes\nabla{S})\right),$
(46)
$\displaystyle\mu=\frac{12}{\text{L\textsubscript{S}}}\left(f^{\prime}(S)-\frac{\text{L\textsubscript{S}}^{2}}{8}\triangle
S\right),$
$\displaystyle{\boldsymbol{\sigma}^{\prime}}=g(d)\frac{\partial{\psi_{e}^{A}}}{\partial\boldsymbol{\epsilon}}+\frac{\partial{\psi_{e}^{I}}}{\partial\boldsymbol{\epsilon}}-\text{K}\nabla{S}\otimes\nabla{S},$

where we have dropped the asterisk notation for clarity. The dimensionless
groups employed in the above are given by

Pe $\displaystyle=\frac{v_{0}l_{0}^{2}}{M_{0}\gamma_{0}},$ Ca
$\displaystyle=\frac{\eta_{0}v_{0}}{\gamma_{0}},$ Da
$\displaystyle=\frac{k_{0}}{l_{0}^{2}},$ (47) K
$\displaystyle=\frac{\kappa}{p_{0}l_{0}^{2}},$ M
$\displaystyle=\frac{\eta_{2}}{\eta_{1}},$ LS
$\displaystyle=\frac{l_{S}}{l_{0}},$ B
$\displaystyle=\frac{8l_{d}v_{0}\beta}{3l_{0}G_{c}},$ Df
$\displaystyle=\frac{8l_{d}\psi_{e}^{A}}{3G_{c}},$ Ld
$\displaystyle=\frac{l_{d}}{l_{0}},$

denoting the Péclet number, the capillary number, the Darcy number, the
Korteweg number, the viscosity ratio (appearing indirectly through
(LABEL:eq:relative_permeabilities)), the normalized saturation regularization
length, the damage viscosity number, the fracture driving force, and the
normalized damage regularization length, respectively. The two first rows
represent groups associated with fluid fields, whereas the third row
represents the groups associated with solid fields. Note that there are only
five independent fluid groups since $\text{Pe}=4\text{Da}^{-1}\text{Ca}/3$.
The expression for the Péclet number matches the one obtained by Abels et
al.[11] for instance. Recall that $\kappa=3\gamma l_{S}/4$.

In practice, it is also useful to define a modified capillary number
$\text{Ca}^{*}$ introduced by Holtzman et al. [32] and employed in the
experiments [4] considered in this work. This number is the ratio of the
driving force (i.e., the viscous pressure drop $\delta{p}_{vis}$) to the
capillary pressure force $\delta{p}_{cap}$. The former can be estimated as

$\delta{p}_{vis}\sim\frac{\eta_{o}w_{i}l_{0}}{k_{0}}\sim\frac{\eta_{o}ql_{0}}{h\delta
k_{0}},$ (48)

where $w_{i}$ is the injection Darcy velocity of water, $q$ is the
corresponding volumetric flow rate, which scales as $q\sim w_{i}h\delta$, $h$
is the height of the Hele–Shaw cell, $\delta$ is a characteristic grain size,
which will be the diameter of the beads used in the experiments studied in
this work, and the characteristic length $l_{0}$ is taken as the radius
$r_{o}$ of the Hele–Shaw cell. The capillary pressure drop is estimated as

$\delta{p}_{cap}\sim\frac{\gamma}{\delta}.$ (49)

Thereby, the modified capillary number reads

$\text{Ca}^{*}\sim\frac{q\eta_{o}l_{0}}{\gamma
hk_{0}}\sim\text{Ca}\frac{l_{0}\delta}{k_{0}}.$ (50)

Given the assumption of viscous-dominated flow, $\text{Ca}^{-1}\gg 1$, and
therefore $\text{Pe}^{-1}\ll 1$, but also $\text{K}\ll 1$, so that one can
anticipate that the terms containing these numbers should not play a
significant role. That said, the Cahn–Hilliard term on the right side of
(LABEL:eq:dimensionless_equations) plays an important interface regularization
role, as mentioned by Cogswell and Szulczewski [39].

##  3 Numerical study

In this section, the ability of the model to reproduce the salient aspects of
recent experimental observations, and in particular, to recover a phase
diagram distinguishing two flow regimes is examined. Model-based simulations
are then employed to make predictions for system responses corresponding to
regions of parameter space that are beyond recent experimental observations.

###  3.1 Numerical implementation

The system of equations (LABEL:eq:summary_equations),
(LABEL:eq:summary_equations_secondary) is discretized through the finite-
element method and implemented within the Multiphysics Object Oriented
Simulation Environment (MOOSE) [52]. For the coupling with phase-field
fracture, the RACCOON extension [53] is used. The resulting system of
equations is solved in a staggered fashion, whereby the poromechanical
equations (LABEL:eq:summary_equations)1-4 and the damage equation
(LABEL:eq:summary_equations)5 are solved iteratively at each time step until
convergence is reached. For the damage equation, the irreversibility
constraint $\dot{d}\geq 0$ is enforced with a primal-dual active set strategy
[41]. Each of the two sub-systems of nonlinear equations are solved with the
Newton–Raphson method. Therein, the rate terms are integrated in time using an
implicit backward Euler approximation. Finally, the matrix inversions are
obtained through an LU decomposition, using the preconditioner MUMPS
(Multifrontal Massively Parallel sparse direct Solver) [54] available within
the PETSc library.

###  3.2 Description of recent experiments

Recent experiments conducted at MIT concerning the fluid-driven fracturing of
a cemented pack of beads in a Hele–Shaw cell [55, 4] are briefly described in
this subsection. Various detailed observations from these experiments are
useful to both calibrate various parameters in the current model and
demonstrate some of its capabilities and limitations.

The experimental configuration consists of a monolayer of spherical beads that
are cemented together. Both the beads and cement are made from polyurethane
rubber, which enables the use of photoporoelasticity, a technique developed by
Li et al. [56] to visualize the effective stress field. The amount of cement
can be controlled, effectively modulating the fracture resistance of the
skeleton. The skeleton is placed between the two plates of the Hele–Shaw cell,
and the plates are held at a constant distance from each other with four outer
clamps. The system is then saturated with oil and water is injected at the
center. As such, the entire system effectively approximates a deformable
porous media interacting with multiphase flow. Aspects of the experimental
setup are very similar to those conducted for non-cohesive granular media, as
described in [3].

Importantly, in the most recent experiments [4], the cement renders the
granular media cohesive, with the degree of cohesion effectively delineating
two regimes in the response of the system. For sufficiently high cement volume
fractions, the solid skeleton deforms but does not fracture as the water
invades the media and expels the defending oil. Conversely, as the amount of
cement is decreased, the fracture resistance of the solid skeleton also
decreases, and at some point the forces on the skeleton are sufficiently large
to break the cement bonds and permit fracture patterns to emerge.

The values of the main parameters are listed in Table 1. Note in particular
the oil viscosity $\eta_{o}$, which is five orders of magnitude larger than
the water viscosity $\eta_{w}$. This gives rise to a relatively large viscous
pressure drop $\delta{p}_{vis}$ (see (48)), acting as a driving force
responsible for fracturing the cemented beads. The high contrast in fluid
viscosities was necessary to compensate for the limited injection pressure
delivered by the laboratory pump.

Another particularity of the type of experiment focused on here is the fact
that water leaks off the fracture: the invading water flows both in the
opening fractures and in the pore spaces of the undamaged domain. Following
the description by Carrillo and Bourg [35], this type of response is referred
to as “invasive fracturing”. A last feature that bears emphasis is that the
injection pressure keeps increasing as the cracks form, almost until the water
is expelled out of the cell. This is expected to be the case as long as cracks
are continuing to grow (see, e.g., Figure 7 in [57]), which is the case in the
experiments and in the corresponding model-based simulations.

It bears emphasis that the mechanical and fracture response of the cement-bead
network was experimentally characterized independently in dry conditions. The
experimentally measured bulk moduli are listed in Table 1. As for the fracture
toughness, the critical mode-I stress intensity factor $K_{1c}$ was measured
via a dog-bone tension test. As a function of the cement ratio $C$ (expressed
in %), it was experimentally found that
$K_{1c}=$0.15C+0.59\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{m}^{1/2}$$,
so that

$G_{c}=\frac{K_{1c}(C)^{2}}{E/(1-\nu^{2})}.$ (51)

This assumes brittle failure, which is confirmed by the stress-strain output
of the tensile test showing a linear increase of the stress until failure,
followed by a sharp decrease; that is, there is negligible plastic deformation
and ductile failure can be neglected.

###  3.3 Boundary-value problem and initial conditions

Consider the two-dimensional annular domain defined by an inner and outer
radius, as shown in Figure 2. The boundary conditions are prescribed on the
inner and outer surfaces in a manner meant to approximate the aforementioned
experimental setup. The values of all model parameters contained within the
boundary conditions are provided in Table 1.

Figure 2: Schematic of the boundary-value problem along with boundary and
initial conditions. For the boundary conditions, $\boldsymbol{n}_{i}$ and
$\boldsymbol{n}_{o}$ denote the inward normal and outward normal at the inner
and outer boundaries, respectively.

For the mechanical fields, boundary conditions that approximate fluid
injection at a constant flow rate at the inner boundary and constrained
outflow conditions at the outer boundary are prescribed. In particular, at the
inner boundary, the Darcy velocity of the invading fluid is fixed at a value
$w_{i}$ chosen to match the experimental injection pressure curve. This is
akin to a traction boundary condition on the displacement field [58], whereby
the injection pressure is applied to the total stress. At the outer boundary,
the pressure is fixed to 0 (the air pressure in the laboratory, used as a
reference pressure) and the displacement field is held fixed, consistently
with the mechanical constraint on the cemented bead network in this region.
Finally, zero Neumann boundary conditions on the damage field are imposed at
both boundaries, as is standard in phase-field models of fracture.

For the saturation field $S$, Dirichlet conditions of $S=1$ and $S=0$ are
prescribed at the inner and outer boundaries, respectively. These conditions
reflect the fact that the invading fluid enters at the inner boundary, while
the defending fluid exits the domain at the outer boundary. The latter
boundary condition obviously limits the applicability of the model to the
point in time when the leading edge of the invading fluid just reaches the
outer boundary.

The boundary conditions on the chemical potential are inspired by the work of
Dong [59]. At the inner boundary, the chemical potential is fixed to a small
value $\mu_{i}>0$ to facilitate the incoming flow of fluid. On the outer
boundary, Dong’s zero-flux condition is replaced by a Dirichlet boundary
condition $\mu_{o}>0$. In the calibration of the model-based simulations
against the experimental observations, the point in time when the injection
pressure peaks is observed to be sensitive to $\mu_{o}$. Accordingly, the
magnitude of $\mu_{o}$ is adjusted to obtain the best match between the
simulated and experimental injection curves.

The initial conditions, as indicated in Figure 2, consist in setting all
variables to 0, except the porosity field. In this work, the influence of the
microstructure stemming from the arrangement of the beads is not examined in
detail. Instead, a continuum perspective and focus on the influence of the
hydro-mechanical parameters is adopted. However, some heterogeneity in the
initial conditions is required to break the inherent symmetry of the problem
and facilitate localization. To effect this, the initial porosity field is
assumed to be uniformly random. More precisely, a porosity value, varying
between 0.2 and 0.6, is randomly assigned to each mesh element (see Figure
3a). This ensures that the mean porosity value is 0.4, which corresponds to
the measured experimental value. Although the effect of this porosity range is
not studied extensively, in limited testing, changes to the upper and lower
bounds in the porosity range are found to yield mostly quantitative and not
qualitative differences in the model-based simulation results.

In addition, to establish spatial convergence with the discretization, a
different porosity field that does not vary with mesh refinement is
considered. To construct such a field, a constant porosity field is used with
three small circular regions of high porosity near the inner boundary, giving
rise to the initial porosity field shown in Figure 3b.

\begin{overpic}[width=195.12767pt]{poro_t=0}
\put(-8.0,65.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$a)$} \end{overpic}

\begin{overpic}[width=195.12767pt]{phi_mesh_conv}
\put(-8.0,65.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$b)$} \end{overpic}

Figure 3: Initial porosity fields used in this work. a) Random porosity field
with mean value $\bar{\phi}=0.4$. b) Constant porosity field with existing
porosity pockets near the inner boundary employed in the mesh convergence
study.

###  3.4 Spatial convergence study

The ability of the model and accompanying discretization to obtain spatial
convergence in selected quantities of interest is first established. A series
of simulations with increasing mesh refinement for the initial boundary-value
problem described in Figure 2 using the initial porosity field shown in Figure
3b are performed. Figure 4 shows plots of the evolution of the injection
pressure over time and of the maximum hoop stress for a reference mesh with
uniform spacing $h=$1.0\text{\,}\mathrm{mm}$$ (20,672 elements) and two
coarser meshes. The results indicate that the reference mesh satisfies the
requirement of mesh convergence. Accordingly, subsequent studies reported in
this manuscript use unstructured meshes of triangular elements with mesh
spacing $h=$1\text{\,}\mathrm{mm}$$. Importantly, this level of spatial
resolution is sufficient to capture the regularization lengths indicated in
Table 1. More precisely, the phase field interfaces for $S$ and $d$ span
approximately four elements.

\begin{overpic}[width=208.13574pt]{mesh_conv_pressure}
\end{overpic}\begin{overpic}[width=208.13574pt]{mesh_conv_hoop_stress}
\end{overpic}

Figure 4: Graphs of the injection pressure over time and of the maximum hoop
stress for different mesh refinements, with reference element size
$h=$1.0\text{\,}\mathrm{mm}$$.

###  3.5 Reproduction of a given experiment

The extent to which the model can reproduce qualitative aspects of a given
Hele–Shaw experiment is now examined. In particular, consider an experiment in
which the cement fraction is 1.2%, the oil viscosity is
$300\text{\,}\mathrm{k}\mathrm{c}\mathrm{S}\mathrm{t}$ (kilo centistokes) or
$0.29\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$, and the injection
flow rate is
$100\text{\,}\mathrm{m}\mathrm{L}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$. In
line with the provided experimental parameters, the model-based simulation
parameters for this problem are summarized in Table 1 in the Appendix. The
random field shown in Figure 3a is used as an initial condition for the
porosity.

The Hele–Shaw cell is represented in the simulations by an annulus of inner
radius $r_{i}=$5\text{\,}\mathrm{mm}$$ and outer radius
$r_{o}=$53\text{\,}\mathrm{mm}$$, as indicated in Figure 2. The inner radius
is estimated by measuring the average radius of the inner cavity delimited by
the beads in the experiment. Within the experimental estimated range of values
for the initial average permeability, a value toward the lower end of this
range is found to work best in the simulations, namely
$0.0015\text{\,}\mathrm{m}\mathrm{m}^{2}$. The elastic moduli of the porous
skeleton are estimated from experimental measurements [55].

The simulation results are obtained with fracture parameters that scales with
the cement volume ratio $C$ as follows. For the fracture toughness, the
experimentally-determined relationship (51) is adopted. The nucleation energy
$\psi_{c}$ governs damage initiation. Assuming that crack nucleation occurs at
a critical value $p_{c}$ of the injection pressure, a simple estimate of the
critical fracture energy is given by (see equation 24 in [38])

$\psi_{c}=\frac{p_{c}(C)^{2}}{2E}.$ (52)

From the experimental measurements, the critical injection pressure $p_{c}$
that corresponded to crack initiation appears to scale linearly with the
cement volume ratio, such that $p_{c}=\alpha C$. A magnitude of
$\alpha=$1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$$ is found to provide a good
match with the experimental observations for crack initiation as a function of
cement volume ratio. Finally, the damage viscosity $\beta$ in
(LABEL:eq:summary_equations)5 is also set to be an increasing function of the
cement ratio. The precise relationship is given in the next section, when the
simulation results are calibrated against experimental observations.

Figure 5 provides snapshots of the experimental results at the start of the
experiment and at $t=$2.0\text{\,}\mathrm{s}$$, along with simulation results
of the damage and porosity fields at the corresponding times. All simulation
results are shown in the reference configuration. The time
$t=$2.0\text{\,}\mathrm{s}$$ corresponds to the moment in the experiment when
the injection pressure reaches a plateau of approximately
$470\text{\,}\mathrm{kPa}$ (see Figure 7). To facilitate comparisons between
the experimental and simulation results at $t=$2.0\text{\,}\mathrm{s}$$, the
water front is indicated with white contour lines. For the experiment, this
contour was extracted manually based on the color of the fluid (blue water vs.
gray oil). The fractured surfaces are also sketched on the experimental image
as blue contours. In the simulation results, the white contour corresponds to
the saturation level $S=0.5$. For reference, the full saturation field at the
start of the simulation and at $t=$2.0\text{\,}\mathrm{s}$$ is shown in Figure
6.

\begin{overpic}[width=169.11153pt]{exp_t=0}
\put(-8.0,80.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$a)$} \end{overpic}

\begin{overpic}[width=169.11153pt]{exp_t=2}
\put(-8.0,80.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$b)$} \end{overpic}

\begin{overpic}[width=208.13574pt]{poro_t=0}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$c)$} \end{overpic}

\begin{overpic}[width=208.13574pt]{poro_t=2}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$d)$} \end{overpic}

\begin{overpic}[width=208.13574pt]{damage_t=0}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$e)$} \end{overpic}

\begin{overpic}[width=208.13574pt]{damage_t=2}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$f)$} \end{overpic}

Figure 5: Comparison of the experimental results (a,b) with the simulation
results (c, d, e, f). The left and right columns correspond to $t=0$ and
$t=$2.0\text{\,}\mathrm{s}$$, respectively. The second row shows the porosity
field, while the third row shows the damage field. In the images on the right,
the white contours indicate the leading edge of the invading fluid front.

\begin{overpic}[width=195.12767pt]{saturation_t=0}
\put(-8.0,80.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$a)$} \end{overpic}

\begin{overpic}[width=195.12767pt]{saturation_t=2}
\put(-8.0,80.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$b)$} \end{overpic}

Figure 6: Saturation field corresponding to the white contour shown in Figure
5 at $t=0$ (a) and $t=$2.0\text{\,}\mathrm{s}$$ (b).

An inspection of the results in Figure 5 indicates that the model is able to
replicate the main salient features of the water invasion. First, there are
three main cracks, corresponding to regions of porosity equal to 1 in the
simulations, preceded by partially damaged zones. At
$t=$2.0\text{\,}\mathrm{s}$$, in the reference configuration, the cracks
obtained in the simulations measure around $11\text{\,}\mathrm{m}\mathrm{m}$
on average, similar to the experiment. Secondly, at
$t=$2.0\text{\,}\mathrm{s}$$, the saturation front from the simulation lies
ahead of the cracks at a distance of $45\text{\,}\mathrm{m}\mathrm{m}$ from
the center, calculated from an average of the distance of the three main
invasion branches. The analogous distance in the experiment is found to be
approximately $44\text{\,}\mathrm{m}\mathrm{m}$. Both the experimental and
simulation results indicate three types of fluid invasion: 1) into fully
intact regions; 2) into partially open regions; and 3) into fully open regions
(where the porosity $\phi=1.0$). Importantly, a fully open region (where
$\phi=1$) does not necessarily coincide with a fully damaged region (where
$d=1$). This is due to the damage viscosity term $\beta\dot{d}$. In other
words, fully open regions can still bear some tensile loading. This is
justified by the presence of residual strands of glue observed in the
experiment (see Figure 5b), which are assumed sufficiently small so as not to
impede flow but large enough to still allow the skeleton to transmit
tractions. As a result, particularly in the present setting with viscous
damage, it is important to interpret hydraulic cracks through the assessment
of the porosity field and not the damage field.

Figure 7 provides a comparison of the injection pressure as measured in the
experiment to that extracted from the model-based simulation. The results
compare favorably, with a first stage of pressure linearly increasing with
time, followed by a plateau starting around $t=$1.7\text{\,}\mathrm{s}$$. In
both the experiment and the simulation, the injection pressure peaks at
approximately $470\text{\,}\mathrm{kPa}$. The decrease of the rate of
injection pressure coincides with the decrease of the rate of the average
damage in the domain, as shown in Figure 15 in Appendix B. The delay between
the damage plateau and the injection pressure plateau may be attributable to
the damage viscosity.

Figure 7: Injection pressure curve obtained from the experiment and matched in
the simulation.

Figure 8 compares the volumetric strain response between the simulation and
the experiment at $t=$2.0\text{\,}\mathrm{s}$$. In both cases, three regions
of tensile volumetric strain are observed to develop out of the injection
site. For both the simulation and the experimental results, these regions
correspond to zones in which fracturing has occurred. On one hand, the high
tensile regions in the experiment are somewhat larger and more diffuse than
those in the simulation. On the other, the transition from regions of tensile
volumetric strain to compressive volumetric strain is more abrupt in the
experiment, at least in the areas between the three main branches. Finally,
the outermost boundary is seen to be in a state of compressive volumetric
strain with a spatial variation that is more uniform in the experiment than
the simulation. This discrepancy may be due to the fact that the four outer
clamps in the experiment are not explicitly modeled in the simulation.
Alternatively, it could also be due to the fact that in the experiments the
outer boundary containing the beads within the domain is not strictly rigid,
and therefore allows some release of elastic energy.

\begin{overpic}[width=164.77771pt,trim=0.0pt -312.9803pt 0.0pt
0.0pt,clip]{vol_strain_t=2.png} \end{overpic}

\begin{overpic}[width=216.81pt]{exp_vol_strain_t=2.png} \end{overpic}

Figure 8: Comparison of the volumetric strain field between the model-based
simulation (left) and the experiment (right). A fluid saturation of $S=0.5$ is
indicated in the simulation result as a white contour.

Finally, Figure 9 shows the simulation results for the pressure, hoop stress,
permeability, and Darcy velocity of water at $t=$2.0\text{\,}\mathrm{s}$$. To
facilitate the interpretation of the results, a white contour line is once
again used to denote the $S=0.5$ saturation front. The peak values of the
pressure field can be seen to be concentrated within the invading fluid phase,
as expected. The hoop stress is observed to be somewhat similar to the
volumetric strain. The permeability is largest within the highly-damaged
regions, until it reaches an upper bound of
($1.92^{2}/12=$0.31\text{\,}\mathrm{m}\mathrm{m}^{2}$$). Recall that this
upper bound follows from assuming a Poiseuille flow confined by the two plates
of the Hele–Shaw cell. As anticipated, the Darcy velocity is maximal in the
damaged areas and at the water front.

\begin{overpic}[width=173.44534pt]{pressure_t=2}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$a)$} \end{overpic}

\begin{overpic}[width=173.44534pt]{hoop_stress_t=2}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$b)$} \end{overpic}

\begin{overpic}[width=173.44534pt]{perm_t=2}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$c)$} \end{overpic}

\begin{overpic}[width=177.78578pt]{darcy_velocity_t=2}
\put(-8.0,70.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
$d)$} \end{overpic}

Figure 9: Simulation results for: a) the pressure; b) the hoop stress; c) the
permeability; and d) the magnitude of the Darcy velocity of water.

###  3.6 Reproduction of an experimental phase diagram

Experiments were conducted using monolayers of cemented beads with varying
degrees of cement and a range of injection rates [4]. As expected, larger
cement volume fractions gave rise to layers with increased fracture
resistance. Figure 10 shows a phase diagram delineating two regimes of
response as a function of the modified capillary number $\text{Ca}^{*}$ and
the cement volume ratio $C$. The symbols in the diagram correspond to
individual experiments or simulations, with the symbol type and color
indicating the observed flow regime and source, respectively.

Results obtained from simulations using the model developed in this work
appear to replicate the phase diagram reasonably well, as indicated in Figure
10. The phase diagram can be interpreted as follows. For a given cement ratio,
fluid injection at relatively low capillary numbers gives rise to porous
invasion without any fracturing of the monolayer. As the capillary number is
increased, the mechanical loading on the monolayer becomes large enough to
give rise to fracturing. Note that the transition from porous invasion to
fracturing regime is progressive, in the sense that a simulation result
labeled as porous invasion may display some small cracks, albeit with little
influence on the invasion pattern. Throughout this work, these two regimes are
differentiated by measuring the length of the longest crack. If this length is
larger (resp. smaller) than an arbitrarily small length, then the simulation
result is labeled as belonging to the fracturing regime (resp. porous invasion
regime). Following the discussion in Section 3.5 regarding the importance of
interpreting hydraulic cracks through the porosity field, crack lengths are
calculated based on the distance from the center of the domain to the end of
regions with porosity $\phi=1$. The outcome of both experiments and
simulations is that the threshold separating porous invasion from fracturing
was found to scale exponentially with the cement volume ratio, as indicated by
the solid black line in Figure 10.

Figure 10: Phase diagram representing the modified capillary number
$\text{Ca}^{*}$ vs. the cement volume ratio $C$. An exponential separatrix
(black line) was found to separate a porous invasion regime from a fracturing
regime for both the experimental and numerical results.

It bears emphasis that a key parameter in the model was adjusted to yield the
best match between the experimental and simulation results in this phase
diagram. As mentioned in Section 3.5, the damage viscosity $\beta$ is assumed
to be an increasing function of the cement volume ratio $C$. A relatively
simple relationship that appears to provide a good match with the experimental
phase diagram is $\beta=c\sqrt{C}$, with
$c=$65.7\text{\,}\mathrm{k}\mathrm{Pa}\,\mathrm{s}$$. Having sufficiently
calibrated the model, in the remainder of this paper, simulations are employed
to explore regions of the parameter space well outside the reach of the
experimental setup.

###  3.7 Effect of the fracture parameters $\psi_{c}$, $G_{c}$, and $\beta$

Before examining the influence of the main hydraulic parameters (permeability
and viscosity), results illustrating the sensitivity of the model to the
fracture parameters are reported. The fracturing response of the system is
governed by 3 parameters in the model: the nucleation energy $\psi_{c}$ (crack
initiation), the fracture toughness $G_{c}$ (crack propagation), and the
damage viscosity $\beta$ (rate effects). The sensitivity of the simulation
results to variations in each of these parameters is illustrated in Figure 11.
The middle column in the Figure corresponds to the base state that was used to
make the comparisons in Section 3.5. For the sake of comparison, all
simulation results are shown when the saturation front is at a similar
distance from the outer boundary.

Figure 11: Effect of the fracture parameters $\psi_{c}$ (nucleation energy),
$G_{c}$ (fracture toughness), $\beta$ (damage viscosity). The damage field is
shown along with the white contour of the saturation front. The middle column
corresponds to a reference simulation common to the three rows, obtained with
the fracture parameters values stemming from the experimental study. The left
and right columns are obtained by dividing and multiplying these reference
values by 10, respectively. The simulations results are shown when the
saturation front is at a same certain distance from the outer boundary.

As expected, increasing either the fracture toughness $G_{c}$ or the damage
viscosity $\beta$ is observed to inhibit crack propagation. Less expected is
the effect of the nucleation energy $\psi_{c}$: while increasing its value
delays crack initiation, the final cracks are observed to progress further in
the domain. This may be explained by the accumulation of strain energy in the
system while crack initiation is prevented until the energy reaches the higher
threshold.

###  3.8 Phase diagram representing fracturing number vs. branching number

A wide range of the parameter space is now explored to infer the types of
possible flow regimes predicted by the model. As described in the
introduction, the model is expected to indicate various types of viscous
flows, namely porous invasion, invasive and non-invasive fracturing, in
possible combination with viscous fingering. The primary model parameters
determining the type of flow are found to be, unsurprisingly, the 3 main fluid
parameters $\eta_{o}$, $\eta_{w}$, and $k_{0}$, which are therefore widely
varied in this section. Note that the characteristic permeability $k_{0}$ is
taken as the initial average permeability in this section.

First, to determine whether the flow is of fracturing type or not, the
fracturing number, denoted by $\text{N}_{f}$, is employed. Taking inspiration
from Holtzman et al. [32], the fracturing number is defined as the ratio of
the driving force to the resisting force. The driving force is taken to be the
viscous pressure drop $\delta p_{vis}$ introduced in (48). The resisting force
is the critical injection pressure $\delta p_{cri}$ at which the flow regime
transitions from porous invasion to fracturing. The fracture number thus
reads:

$\text{N}_{f}=\frac{\delta p_{vis}}{\delta p_{cri}}.$ (53)

While the driving force characterizes the invasion and contains only flow
parameters, the resisting force characterizes the fracture resistance of the
medium and depends on pertinent strength quantities such as the cement
fraction. In this section, the fracture parameters are held fixed and,
therefore, the value of $\delta p_{cri}$ is common for all simulations.

Second, to further characterize the type of flow within the fracturing and
porous invasion regimes, a new dimensionless number is constructed, called the
branching number and denoted by $\text{N}_{b}$, given by

$\text{N}_{b}=\frac{\eta_{w}l_{0}^{2}}{\eta_{o}k_{0}}.$ (54)

Note that the branching number is related to the Darcy number and the mobility
ratio defined in (47) through $\text{N}_{b}=\text{M}^{-1}\text{Da}^{-1}$. For
a fixed fracturing number $\text{N}_{f}$, varying the branching number
$\text{N}_{b}$ was found to give rise to two types of flow. In the porous
invasion regime ($\text{N}_{f}<1$), low values of $\text{N}_{b}$ indicate a
uniform porous invasion, whereas large values of $\text{N}_{b}$ indicate a
viscous fingering instability. The number of branches of this instability was
found to increase with $\text{N}_{b}$. In the fracturing regime
($\text{N}_{f}>1$), low values of $\text{N}_{b}$ indicate invasive fracturing,
whereas large values of $\text{N}_{b}$ indicate non-invasive fracturing. The
number of “crack branches” was found to increase with $\text{N}_{b}$.

To illustrate these observations, Figure 12 shows the individual influence of
the main parameters of $\text{N}_{b}$ on the type of flow. A first intuitive
result is that decreasing the permeability tends to result in the invading
fluid flow being more confined to the cracks, thereby indicating a transition
from invasive to non-invasive fracturing. Non-invasive fracturing is
characterized by the invading fluid front being confined to the open channels
that form. As for the influence of the fluid viscosities, similar trends are
observed with increasing $\eta_{w}$ and $\eta_{o}$, as shown in the figure.
Note that while $k_{0}$ and $\eta_{w}$ play a role that is consistent with the
aforementioned rationale behind the branching number, it is not the case for
$\eta_{o}$, at least individually. However, as a group, the branching number
does play the role explained above, for fixed values of the fracturing number.
Finally, the left simulation result of the third row of Figure 12 displays a
new type of flow, which is the combination of the fracturing regime and
viscous fingering. However, the branching number alone is not sufficient to
identify the parameter space where this regime dominates; this will be done in
combination with the fracturing number later in this section.

Figure 12: Plots of the porosity field and saturation front (white contour)
illustrating the sensitivity to changes in various flow parameters. Influence
on the branching pattern of permeability (from left to right:
$\bar{k}=$5\times 10^{-2}\text{\,}\mathrm{m}\mathrm{m}^{2}$$,
$\bar{k}=$10^{-2}\text{\,}\mathrm{m}\mathrm{m}^{2}$$,
$\bar{k}=$10^{-3}\text{\,}\mathrm{m}\mathrm{m}^{2}$$, at fixed
$\eta_{o}=$1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$ and
$\eta_{w}=$10^{-3}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$), of
invading fluid viscosity $\eta_{w}$ (from left to right:
$\eta_{w}=$10^{-6}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$,
$\eta_{w}=$10^{-4}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$,
$\eta_{w}=$10^{-3}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$, at
fixed $\eta_{o}=$1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$ and
$\bar{k}=$10^{-3}\text{\,}\mathrm{m}\mathrm{m}^{2}$$), and of defending fluid
viscosity (from left to right:
$\eta_{o}=$10^{-3}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$,
$\eta_{o}=$10^{-2}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$,
$\eta_{o}=$1\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$, at fixed
$\eta_{w}=$10^{-3}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$ and
$\bar{k}=$10^{-3}\text{\,}\mathrm{m}\mathrm{m}^{2}$$).

Following the previous qualitative study of the role of the main fluid
parameters on the type of flow, a quantitative study is now presented, with
the phase diagram shown in Figure 13.

Figure 13: Phase diagram representing the effect of the three main flow
parameters on the type of flow regime.

Note in particular that although the invading fluid viscosity $\eta_{w}$ has
an effect on the branching pattern, it has a very limited influence on the
fracture onset, as experimentally observed by Zhou et al. [60], and
numerically confirmed by Carrillo and Bourg [35]. The model-based simulations
appear to recover this effect.

To build a phase diagram of the fracturing number $\text{N}_{f}$ vs. the
branching number $\text{N}_{b}$ representing all possible flow regimes, a
large number of simulations are conducted over wide ranges of the parameters
$\eta_{o}$, $\eta_{w}$, $k_{0}$. The resulting phase diagram is shown in
Figure 14, with each symbol corresponding to an individual simulation and
positioned at the corresponding values of $\text{N}_{b}$ and $\text{N}_{f}$.
For the resulting values of $\text{N}_{b}$, the characteristic length was
chosen as the inner radius of the Hele–Shaw cell (i.e.,
$l_{0}=r_{i}=$5\text{\,}\mathrm{mm}$$).

The particular flow regime indicated by each simulation is tagged as either a
blue circle (uniform invasion), a blue star (viscous fingering), an orange
triangle (invasive fracturing), an orange square (non-invasive fracturing), or
an orange cross (fracturing/fingering). The resulting regions where each type
of flow dominates are also colored to help indicate threshold values for
$\text{N}_{f}$ or $\text{N}_{b}$ that separate various regimes.

Note that $\text{N}_{b}$ and $\text{N}_{f}$ are inherently different since,
unlike the latter, the former is not a conventional dimensionless number since
it is not the ratio of two competing forces. Nevertheless, the phase diagram
indicates that both the dimensionless numbers $\text{N}_{b}$ and
$\text{N}_{f}$ do, in fact, play the role of regime delimiters, with values of
$\text{N}_{f}=1$ and $\text{N}_{b}=1$ roughly separating fracturing vs. non-
fracturing for the former, and uniform invasion vs. viscous fingering and
invasive fracturing vs. non-invasive fracturing for the latter. A new flow
regime, mixing fracturing and fingering is found between the regions of non-
invasive fracturing and viscous fingering, for $\text{N}_{b}>1$ and a certain
range of $\text{N}_{f}$. Interestingly, in the region $\text{N}_{b}>1$,
sufficiently large values of the fracturing number effectively prohibits
viscous fingering, apparently favoring flow within the fractures over porous
invasion. Finally, the red rectangle shown in Figure 14 represents the region
of the diagram explored by the experiments. The latter was obviously able to
showcase the two regimes of invasive fracturing and uniform invasion, while
our simulations indicate the possibility of three more regimes that are beyond
the range of conditions studied experimentally.

Figure 14: Phase diagram of the fracturing number $N_{f}$ vs the branching
number $N_{b}$, mapping the five types of flow regimes. The red rectangle
indicates the region of diagram explored by the experiments. The circular
inserts are simulation results representing the porosity field along with the
saturation front in white.

For practical purposes, this phase diagram can provide insight on the type of
flow regime that can be expected at the field scale. For instance, it is well
known that viscous fingering decreases the efficiency of fluid-fluid
displacement in industrial applications (e.g., enhanced oil recovery [61]). It
is expected that industrial hydraulic fracturing may also suffer from this
efficiency decrease, insofar as our results suggest the existence of a flow
regime combining hydraulic fracturing and viscous fingering, provided that the
branching number $\text{N}_{b}$ is sufficiently large. To estimate the latter,
consider three types of field applications. First, in the context of oil
recovery, if water
($\eta_{w}\sim$10^{-6}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$)
is used to displace oil
($\eta_{o}\sim$10^{-3}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$$),
through a reservoir of permeability
$k\sim$10^{-9}\text{\,}\mathrm{m}\mathrm{m}^{2}$$, then $\text{N}_{b}\sim
10^{10}$. The injection borehole is assumed to have a radius of
$l_{0}\sim$100\text{\,}\mathrm{m}\mathrm{m}$$. Second, in the context of CO2
sequestration, where brine is displaced for water co-production, one can
estimate $\text{N}_{b}\sim 10^{12}$. The invading fluid, supercritical CO2, is
assumed to have a viscosity of
$10^{-7}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$, whereas the
defending fluid, brine, is assumed to have the same viscosity as water, and
the permeability is the same as in the first application. Third, in the
context of geothermal reinjection, even though the invading fluid, water, is
the same as the defending fluid, the viscosity of the former can be lower than
that of the latter when its temperature is larger [62]. As in Macdowell et al.
[62], when the viscosity of the invading and defending fluid is
$10^{-4}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$ and
$10^{-3}\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}\,\mathrm{s}$, respectively,
and the permeability is $10^{-7}\text{\,}\mathrm{m}\mathrm{m}^{2}$, then
$\text{N}_{b}\sim 10^{10}$. In all, the branching numbers for field
applications are much larger than the range studied in Figure 14, so that the
type of flow cannot be ascertained for the moment. It also bears emphasis that
in real reservoirs, the flow regimes are not expected to be as distinguishable
as in the simulations shown here. In particular, invasive fracturing, also
called leak-off, can occur through existing natural cracks, even for tight
geological formations [17]. This means that invasive fracturing could occur
for high values of $\text{N}_{b}$.

##  4 Conclusion

In this work, a double phase-field approach regularizing both cracks and
fluid-fluid interfaces is introduced. Derived within the frameworks of
continuum thermodynamics and linear poroelasticity, the model behavior is
characterized by a tight three-way coupling between multiphase fluid flow,
poromechanics, and fracturing. Through finite-element discretization and
numerical simulations, the model was validated against Hele–Shaw experiments,
both directly and via the reconstruction of a phase diagram discriminating
porous invasion from hydraulic fracturing. The parameter space was then
explored beyond experimental capabilities to discover a variety of flow
regimes, including a combination of fracturing and viscous fingering.

The model employs many simplifying assumptions that permit relatively low
computational cost while satisfyingly reproducing relatively complex
experiments. However, in future works, it will be of interest to explore the
influence of higher-fidelity modeling choices. First, the small strain
assumption underpinning linear poroelasticity is at the limit of validity in
the present setting and may be replaced by finite-strain kinematics, as in
Paulin et al. [37]. Second, interpreting the packing of the experimental beads
as a continuum is at the limit of validity as well. This could be addressed by
resorting to explicitly modeling the microstructure on one hand or by
experiments where the microstructural scale is much smaller than the domain
size, on the other. Third, the binary fluid flow description can be enriched
by replacing our Darcy–Poisseuille flow model by a Darcy–Stokes flow model, as
introduced by Ehlers and Luo [29], and Wilson and Landis [30]. Fourth,
provided light modifications such as the introduction of a capillary pressure,
our model can accommodate for capillary flows in addition to viscous flows.

Ultimately, the model presented in this contribution could be used to predict
and control the type of flows in field applications, as it bridges laboratory
and field scales.

## Acknowledgements

This work was the result of a collaboration between researchers at Duke and
MIT, largely supported by a set of collaborative NSF research grants, in
particular NSF grant CMMI-1933367 to Duke University, NSF grant CMMI-1933416
to MIT, and NSF grant CMMI-1826221 to Penn State. The support is gratefully
acknowledged. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of the National Science Foundation.

##  Appendix A Summary table of the experiment’s and model’s parameters

Table 1: Summary of the values chosen for the different model parameters. Name | Symbol | Unit | Exp. value | Num. value
---|---|---|---|---
| height of
---
Hele–Shaw cell
$h$ | $\text{\,}\mathrm{mm}$ | 1.96 | /
| outer radius of
---
Hele–Shaw cell
$r_{o}$ | $\text{\,}\mathrm{mm}$ | 53 | 53
beads diameter | $d$ | $\text{\,}\mathrm{mm}$ | 2 | /
packing fraction | $1-\phi_{0}$ | - | 0.60 | 0.60
cement volume ratio | $C$ | % | $[0,3]$ | /
oil viscosity | $\eta_{o}$ | $\text{\,}\mathrm{k}\mathrm{Pa}\,\mathrm{s}$ | $[0.029,0.29]$ | $[10^{-4},1]$
water viscosity | $\eta_{w}$ | $\text{\,}\mathrm{k}\mathrm{Pa}\,\mathrm{s}$ | $10^{-6}$ | $[10^{-6},10^{-3}]$
permeability | $k$ | $\text{\,}\mathrm{m}\mathrm{m}^{2}$ | $[(0.02d)^{2},(0.06d)^{2}]$ | /
permeability coeff | $\bar{k}$ | $\text{\,}\mathrm{m}\mathrm{m}^{2}$ | / | $[10^{-4},10^{-1}]$
surface tension | $\gamma$ | $\text{\,}\mathrm{kPa}\text{\,}\mathrm{mm}$ | 0.03 | 0.03
| Young’s modulus
---
of granular pack
$E$ | $\text{\,}\mathrm{kPa}$ | 331 | 331
| drained bulk modulus
---
of granular pack
$K$ | $\text{\,}\mathrm{kPa}$ | 184 | 184
| Poisson ratio
---
of granular pack
$\nu$ | $\text{\,}\mathrm{kPa}$ | 0.2 | 0.2
| bulk modulus
---
of solid grains
$K_{s}$ | $\text{\,}\mathrm{kPa}$ | 1600 | 1600
nucleation energy | $\psi_{c}$ | $\text{\,}\mathrm{kPa}$ | / | $0.0015C^{2}$
fracture toughness | $G_{c}$ | $\text{\,}\mathrm{kPa}\text{\,}\mathrm{mm}$ | $(0.26C+1.00)^{2}$ | $(0.26C+1.00)^{2}$
damage viscosity | $\beta$ | $\text{\,}\mathrm{k}\mathrm{Pa}\,\mathrm{s}$ | / | $65.7\sqrt{C}$
| regularization length
---
for $S$
$l_{S}$ | $\text{\,}\mathrm{mm}$ | / | 4.0
| regularization length
---
for $d$
$l_{d}$ | $\text{\,}\mathrm{mm}$ | / | 4.0
| Inner boundary
---
condition for $\boldsymbol{w}_{w}$
$w_{i}$ | $\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$ | / | 50
| Inner boundary
---
condition for $\mu$
$\mu_{i}$ | $\text{\,}\mathrm{kPa}$ | / | 0.006
| Outer boundary
---
condition for $\mu$
$\mu_{o}$ | $\text{\,}\mathrm{kPa}$ | / | 0.02

##  Appendix B Comparison injection pressure and average damage

Figure 15: Comparison of the injection pressure curve with the average damage
curve in the simulation reproducing the experiment.

## References

  * [1] R. Juanes, Y. Meng, B. K. Primkulov, Multiphase flow and granular mechanics, Phys. Rev. Fluids 5 (2020) 110516.
  * [2] Y. Meng, B. K. Primkulov, Z. Yang, C. Y. Kwok, R. Juanes, Jamming transition and emergence of fracturing in wet granular media, Physical Review Research 2 (2) (2020) 022012.
  * [3] Y. Meng, W. Li, R. Juanes, Fracturing in wet granular media illuminated by photoporomechanics, Physical Review Applied 18 (6) (2022) 064081.
  * [4] Y. Meng, W. Li, R. Juanes, Crossover from viscous fingering to fracturing in cohesive wet granular media: a photoporomechanics study, submitted for publication (2023).  arXiv:arXiv:2305.11037.
  * [5] P. G. Saffman, G. I. Taylor, The penetration of a fluid into a porous medium or Hele–Shaw cell containing a more viscous liquid, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 245 (1242) (1958) 312–329.
  * [6] R. Lenormand, Pattern growth and fluid displacements through porous media, Physica A: Statistical Mechanics and its Applications 140 (1-2) (1986) 114–123.
  * [7] B. Zhao, C. W. MacMinn, R. Juanes, Wettability control on multiphase flow in patterned microfluidics, Proceedings of the National Academy of Sciences 113 (37) (2016) 10251–10256.
  * [8] J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. i. interfacial free energy, The Journal of chemical physics 28 (2) (1958) 258–267.
  * [9] J. Lowengrub, L. Truskinovsky, Quasi-incompressible Cahn–Hilliard fluids and topological transitions, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454 (1978) (1998) 2617–2654.
  * [10] H.-G. Lee, J. Lowengrub, J. Goodman, Modeling pinchoff and reconnection in a Hele–Shaw cell. i. the models and their calibration, Physics of Fluids 14 (2) (2002) 492–513.
  * [11] H. Abels, H. Garcke, G. Grün, Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities, Mathematical Models and Methods in Applied Sciences 22 (03) (2012) 1150013.
  * [12] M. E. Gurtin, D. A. Polignone, J. Viñals, Two-phase binary fluids and immiscible fluids described by an order parameter, Mathematical Models and Methods in Applied Sciences 6 (1995) 815–831.
  * [13] L. Dede, H. Garcke, K. F. Lam, A Hele–Shaw–Cahn–Hilliard model for incompressible two-phase flows with different densities, Journal of Mathematical Fluid Mechanics 20 (2) (2018) 531–567.
  * [14] P. Papatzacos, Macroscopic two-phase flow in porous media assuming the diffuse-interface model at pore level, Transport in porous media 49 (2) (2002) 139–174.
  * [15] L. Cueto-Felgueroso, R. Juanes, Nonlocal interface dynamics and pattern formation in gravity-driven unsaturated flow through porous media, Physical Review Letters 101 (24) (2008) 244504.
  * [16] L. Cueto-Felgueroso, R. Juanes, A phase field model of unsaturated flow, Water resources research 45 (10) (2009).
  * [17] B. Chen, B. R. Barboza, Y. Sun, J. Bai, H. R. Thomas, M. Dutko, M. Cottrell, C. Li, A review of hydraulic fracturing simulation, Archives of Computational Methods in Engineering (2021) 1–58.
  * [18] G. A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (8) (1998) 1319–1342.
  * [19] B. Bourdin, G. Francfort, J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (2000) 797–826.
  * [20] 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, Journal of the Mechanics and Physics of Solids 82 (2015) 186–217.  doi:https://doi.org/10.1016/j.jmps.2015.04.006.
  * [21] 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, Computer Methods in Applied Mechanics and Engineering 304 (2016) 619–655.  doi:https://doi.org/10.1016/j.cma.2015.09.021.
  * [22] M. A. Biot, General theory of three-dimensional consolidation, Journal of applied physics 12 (2) (1941) 155–164.
  * [23] O. Coussy, Poromechanics, John Wiley & Sons, 2004.
  * [24] R. M. Bowen, Incompressible porous media models by use of the theory of mixtures, International Journal of Engineering Science 18 (9) (1980) 1129–1148.
  * [25] R. De Boer, Theory of porous media: highlights in historical development and current state, Springer Science & Business Media, 2012.
  * [26] W. Ehlers, Foundations of multiphasic and porous materials, in: Porous media, Springer, 2002, pp. 3–86.
  * [27] C. Truesdell, R. Toupin, The classical field theories, in: Principles of classical mechanics and field theory/Prinzipien der Klassischen Mechanik und Feldtheorie, Springer, 1960, pp. 226–858.
  * [28] R. M. Bowen, Part i - theory of mixtures, in: A. C. ERINGEN (Ed.), Continuum Physics, Academic Press, 1976, pp. 1–127.  doi:https://doi.org/10.1016/B978-0-12-240803-8.50017-7.
  * [29] W. Ehlers, C. Luo, A phase-field approach embedded in the theory of porous media for the description of dynamic hydraulic fracturing, Computer Methods in Applied Mechanics and Engineering 315 (2017) 348–368.
  * [30] Z. A. Wilson, C. M. Landis, Phase-field modeling of hydraulic fracture, Journal of the Mechanics and Physics of Solids 96 (2016) 264–290.
  * [31] Y. Heider, A review on phase-field modeling of hydraulic fracturing, Engineering Fracture Mechanics 253 (2021) 107881.
  * [32] R. Holtzman, M. L. Szulczewski, R. Juanes, Capillary fracturing in granular media, Physical review letters 108 (26) (2012) 264504.
  * [33] S. Lee, A. Mikelic, M. F. Wheeler, T. Wick, Phase-field modeling of two phase fluid filled fractures in a poroelastic medium, Multiscale Modeling & Simulation 16 (4) (2018) 1542–1580.
  * [34] Y. Heider, W. Sun, A phase field framework for capillary-induced fracture in unsaturated porous media: Drying-induced vs. hydraulic cracking, Computer Methods in Applied Mechanics and Engineering 359 (2020) 112647.
  * [35] F. J. Carrillo, I. C. Bourg, Modeling multiphase flow within and around deformable porous materials: a darcy-brinkman-biot approach, Water Resources Research 57 (2) (2021) e2020WR028734.
  * [36] G. Pijaudier-Cabot, Z. P. Bažant, Nonlocal damage theory, Journal of engineering mechanics 113 (10) (1987) 1512–1533.
  * [37] O. W. Paulin, L. C. Morrow, M. G. Hennessy, C. W. MacMinn, Fluid–fluid phase separation in a soft porous medium, Journal of the Mechanics and Physics of Solids 164 (2022) 104892.
  * [38] R. Geelen, Y. Liu, T. Hu, M. Tupek, J. E. Dolbow, A phase-field formulation for dynamic cohesive fracture, Computer Methods in Applied Mechanics and Engineering (2019).
  * [39] D. A. Cogswell, M. L. Szulczewski, Simulation of incompressible two-phase flow in porous media with large timesteps, Journal of Computational Physics 345 (2017) 856–865.
  * [40] R. I. Borja, On the mechanical energy and effective stress in saturated and unsaturated porous continua, International Journal of Solids and Structures 43 (2006) 1764–1786.
  * [41] T. Hu, J. Guilleminot, J. E. Dolbow, A phase-field model of fracture with frictionless contact and random fracture properties: Application to thin-film fracture and soil desiccation, Computer Methods in Applied Mechanics and Engineering 368 (2020) 113106.
  * [42] T. Hu, B. Talamini, A. J. Stershic, M. Tupek, J. E. Dolbow, A variational phase-field model for ductile fracture with coalescence dissipation, Computational Mechanics (2021) 1–25.
  * [43] E. Fried, M. E. Gurtin, Dynamic solid-solid transitions with phase characterized by an order parameter, Physica D: Nonlinear Phenomena 72 (4) (1994) 287–308.  doi:10.1016/0167-2789(94)90234-8.
  * [44] M. E. Gurtin, Generalized Ginzburg–Landau and Cahn–Hilliard equations based on a microforce balance, Physica D: Nonlinear Phenomena 92 (1996) 178–192.
  * [45] M. N. da Silva Jr, F. P. Duda, E. Fried, Sharp-crack limit of a phase-field model for brittle fracture, Journal of the Mechanics and Physics of Solids 61 (11) (2013) 2178–2195.
  * [46] C. W. MacMinn, E. R. Dufresne, J. S. Wettlaufer, Large deformations of a soft porous material, Physical Review Applied 5 (4) (2016) 044020.
  * [47] B. D. Coleman, W. Noll, The thermodynamics of elastic materials with heat conduction and viscosity, Archive for Rational Mechanics and Analysis 13 (1) (1963) 167–178.  doi:10.1007/bf01262690.
  * [48] J.-J. Marigo, C. Maurini, K. Pham, An overview of the modelling of fracture by gradient damage models, Meccanica 51 (2016) 3107–3128.
  * [49] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations, International journal for numerical methods in engineering 83 (10) (2010) 1273–1311.
  * [50] M. Fourar, R. Lenormand, A viscous coupling model for relative permeabilities in fractures, in: Spe annual technical conference and exhibition, OnePetro, 1998\.
  * [51] L. Cueto-Felgueroso, R. Juanes, A phase-field model of two-phase Hele–Shaw flow, Journal of Fluid Mechanics 758 (2014) 522–552.  doi:10.1017/jfm.2014.512.
  * [52] C. J. Permann, D. R. Gaston, D. Andrš, R. W. Carlsen, F. Kong, A. D. Lindsay, J. M. Miller, J. W. Peterson, A. E. Slaughter, R. H. Stogner, R. C. Martineau, MOOSE: Enabling massively parallel multiphysics simulation, SoftwareX 11 (2020) 100430.  doi:10.1016/j.softx.2020.100430.
  * [53] https://github.com/hugary1995/raccoon.
  * [54] https://mooseframework.inl.gov/modules/porous_flow/solvers.html.
  * [55] Y. Meng, W. Li, R. Juanes, personal communication.
  * [56] W. Li, Y. Meng, B. K. Primkulov, R. Juanes, Photoporomechanics: an experimental method to visualize the effective stress field in fluid-filled granular media, Physical Review Applied 16 (2) (2021) 024043.
  * [57] 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.
  * [58] G. Rehbinder, Analytical solutions of stationary coupled thermo-hydro-mechanical problems, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 32 (5) (1995) 453–463.  doi:https://doi.org/10.1016/0148-9062(95)00035-F.
  * [59] S. Dong, An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach, Journal of Computational Physics 266 (2014) 47–73.  doi:https://doi.org/10.1016/j.jcp.2014.02.011.
  * [60] Experimental Study of the Impact of Shear Dilation and Fracture Behavior during Polymer Injection for Heavy Oil Recovery in Unconsolidated Reservoirs, Vol. All Days of SPE Canada Unconventional Resources Conference, sPE-137656-MS.  doi:10.2118/137656-MS.
  * [61] A. Pinilla, M. Asuaje, N. Ratkovich, Experimental and computational advances on the study of viscous fingering: An umbrella review, Heliyon 7 (7) (2021) e07614.
  * [62] A. Mcdowell, S. J. Zarrouk, R. Clarke, Modelling viscous fingering during reinjection in geothermal reservoirs, Geothermics 64 (2016) 220–234.  doi:https://doi.org/10.1016/j.geothermics.2016.06.003.