mixtures (Abichou et al.,, 2002). The curves’ parameters are:$\
{\mathrm{A}}_{\mathrm{1}}\mathrm{=4.2\ \&\
}{\mathrm{V}}^{\mathrm{*}}_{\mathrm{1,\ clay}}\mathrm{=0.25,\
}{\mathrm{A}}_{\mathrm{2}}\mathrm{=4.8\ \&\
}{\mathrm{V}}^{\mathrm{*}}_{\mathrm{2,\ clay}}\mathrm{=0.10,\
}{\mathrm{A}}_{\mathrm{3}}\mathrm{=7\ \&\
}{\mathrm{V}}^{\mathrm{*}}_{\mathrm{3,\ clay}}\mathrm{=0.12,\ and\
}{\mathrm{A}}_{\mathrm{4}}=7.4\mathrm{\ \&\
}{\mathrm{V}}^{\mathrm{*}}_{\mathrm{4,clay}}\mathrm{=0.18}$ with S = 1.38
$\times$ A.

Further analysis shows that the simulation predictions in Fig. 4.8A collapse
into a single relation between log(k/k0) and the clay matrix volume fraction
Vclay (i.e., the ratio of volume occupied by the clay matrix, including its
internal microporosity, to the total volume occupied by the clay and coarse
grains). As shown in Fig. 4.8B, the main preferential fluid flow path is
consistently closed at Vclay $\approx$ 25% in our simulations, in reasonable
agreement with previous predictions for spherical bead packs clogged either by
clay or by cementation (MacArt & Mueller,, 2016; Tuller & Or,, 2006). At clay
volume fractions below this percolation threshold, clay swelling influences
permeability by a relatively small (but still highly significant) two orders
of magnitude.

To compare the simulation predictions in Fig. 4.8B with experimental data on
the permeability of clayey media, we fitted the simulation predictions of
log(k/k0) vs. Vclay using an error function:

${\mathrm{log}\left(\frac{k}{k_{0}}\right)\
}=\frac{A}{2}\left({\mathrm{erf}\left(\ S\
\left(V_{clay}-V^{*}_{clay}\right)\right)\ }-1\right)$ (4.7)

In Eqn. 4.7, A is the overall magnitude of the decrease in log k, S describes
the sensitivity of permeability to $V_{clay}$ near the percolation threshold,
and $V^{*}_{clay}$ describes the location of the threshold. Equation 4.7
provides a reasonable fit to the simulation predictions in Fig. 4.8B with only
three parameters (though it underestimates the sharpness of the threshold). To
test the broader validity of Eqn. 4.7, we identified four experimental
datasets that reported the permeability of a series of porous media with
similar mineralogy differing only in their clay content and that expressed
clay content in the same manner as in Fig. 4.8B (Abichou et al.,, 2002; Gräfe
et al.,, 2017; Heap et al.,, 2017; Revil & Cathles,, 1999). As shown in Fig.
4.8C, Eqn. 4.7 provides accurate descriptions of all four datasets. In the
five parametric fits carried out with Eqn. 4.7 in Figs. 4.8B,C, the values of
S and A were consistently related (S/A = 1.38 $\pm$ 0.25). Fitted A values
ranged from 4.2 to 9.6, a range consistent with the reported permeability
values of different clays, which span 4 to 5 orders of magnitude and increase
in the order Na-smectite $\mathrm{<}$ Ca-smectite $\approx$ Illite
$\mathrm{<}$ kaolinite (Mesri & Olson,, 1971; Mondol et al.,, 2008). Fitted
values of $V^{*}_{clay}$ ranged from 0.10 to 0.25, in agreement with the
expectation that the location of the percolation threshold depends on the
porosity of the network of coarse grains (Aksu et al.,, 2015), i.e., it should
be lower for systems with greater cementation or with less uniform grain size
distributions. The simulation predictions shown in Fig. 4.8B, as expected, are
at the upper end of the range of both A values (because we used the lowest-
permeability clay, Na-smectite) and $V^{*}_{clay}$ values (because we used a
2-D slice of a 3-D bead pack with a uniform bead size and no cementation).

The simulations described in this section are designed to approximate systems
where the network of coarse grains is load-bearing and where clay influences
permeability only through swelling and shrinking in response to salinity
changes. Clay erosion, transport, and deposition are not included in our
model, while viscous deformation of the clay by the flowing fluid is minimized
by the use of a small fluid pressure gradient. Nevertheless, the good
agreement between the results obtained with smectite and Illite and between
the parametric fits to simulated and measured permeabilities suggest that the
computational framework used here may be applicable to clayey media with a
range of porosities, clay mineralogies, and configurations.

###  4.4 Model Application to Fluid Injection Into a Poroelastic Medium

As a final illustration of the versatility of the modeling framework developed
in Chapter 2, we briefly apply and verify our model for a simple poroelastic
system. More precisely, we simulate fluid injection into a symmetric
poroelastic material, a model system used to mimic the properties of filters,
biological tissues, membranes, and soils (Auton & MacMinn,, 2017; Barry et
al.,, 1995; Nagel & Kelly,, 2012). Barry et al., (1997) used Biot theory and
the assumption that material stress is directly proportional to the fluids’
frictional forces (i.e. drag) to obtain an analytical solution for the
pressure and deformation profiles brought about by a constant point fluid
source placed at different depths in the system shown in Fig. 4.9A. The system
studied by Barry (and replicated in our simulations) is defined by a no-slip
condition at the lower boundary (y = 0), a constant-pressure fluid point
source placed at height $y_{0}$ on the axis of symmetry (x = 0), and zero-
gradient boundary conditions at all other boundaries. Please refer to Table
4.1 and Barry et al., (1997) for the complete parametrization of the
simulation.

Figure 4.9: (a) Representation of the modeled system. (b) Quantitative
comparison of the maximum vertical deformation as a function of the point
fluid source position ${(y}_{0}/h)$. Our model correctly predicts the vertical
deformation of the material with a MAE equal to 0.01 times the overall
deformation. One potential source of error is the fact that the fluid source
in our numerical simulations was fixed in space and did not move upward with
solid deformation as it does in the analytical solution. (c-f) Qualitative
comparison of the deformation and pressure profiles obtained by our model on a
120 by 150 grid (c and e) and Barry’s analytical solution (d and f),
respectively. The red arrows in c and e represent the fluid source position in
the numerical model. The parameters for these simulations can be found in
Table 4.1.

The complete figure shows that our numerical framework replicates Barry’s
analytical solutions for the system’s maximum deformation in the y-direction
to a relatively high degree of accuracy for most values of $y_{0}$. As stated
in the original paper, the minimum at $y_{0}/h\approx 0.4$ results from two
phenomena: first, when the source is near $y/h=1\ $most of the flow leaves
through the top boundary, increasing drag near this boundary and maximizing
its deformation; second, the maximum vertical deformation is essentially the
sum of all local deformations along the axis of symmetry ($x=0$), meaning that
at low values of $y_{0}$, fluid drag is able to act over more of the solid,
thus increasing the overall accumulated vertical deformation. The sum of these
two effects is minimized at intermediate values of $y_{0},$ leading to the
observed minimum. The tendency of our simulations to underestimate surface
deformation at high values of $y_{0}$ likely reflects the fact that the
location of the fluid source in our simulations is not affected by the solid’s
deformation, contrary to the analytical solution, a difference most
significant at $y_{0}/h\approx 1$ because of the first phenomenon noted above.

One notable weakness of the simulations presented here is that the fluid-solid
interface is not represented as a sharp step-function, but rather as a
continuous interpolation of the stresses between both domains. This loss of
sharpness is particularly noticeable as the magnitude of the deformations
become large, as seen in Fig. 4.9C, perhaps reflecting the breakdown of the
validity of the elastic momentum equation in systems with large deformations.
This loss of sharpness is much less apparent in the case of the swelling
porous media modeled in the previous section, where the swelling pressure acts
in favor of a sharp interface. Errors emanating from the lack of sharpness in
elastic porous media could be addressed by the introduction of an artificial
compression term (as used in VOF models) (Rusche,, 2002), higher resolution
meshes, or the addition a suitable swelling pressure model.

Table 4.1: Parameter values were obtained from Bourg et al., (2006), Diamond &
Kinter, (1956), Liu, (2013), Mesri & Olson, (1971), and Spearman, (2017) and
are consistent across the simulations unless specified otherwise.

###  4.5 Conclusions

We have implemented, tested, and verified a multi-scale Hydro-Chemo-
Mechanically coupled modeling framework for porous media containing both
macropores and deformable microporous regions. We demonstrated the framework’s
accuracy and versatility by modeling HCM-coupled viscoplastic and poroelastic
dynamics.

The model was implemented and parametrized in particular to simulate the
swelling behavior and fluid-induced deformation of water-saturated
viscoplastic Na-smectite clay. Predicted couplings between hydrology,
mechanics, and salinity were successfully validated against experimental data
and numerical models by Di Maio, (1996) and Hewitt et al., (2016) describing a
broad range of simulated conditions. The model was then used to predict the
permeability of a spherical bead pack as a function of clay content and
salinity. For this simple model, a master parametric equation for permeability
as a function of clay matrix volume fraction was extracted. This parametric
equation was found to be consistent with experimental datasets on the
permeability of smectite-coated glass bead packs and of different types of
siliciclastic sedimentary rocks. Finally, the model’s versatility was
demonstrated by qualitatively predicting water-induced formation damage in a
propped fracture in clayey rock and, also, by quantitatively predicting the
pressure fields and deformation profiles resulting from fluid injection in an
elastically-deformable axisymmetric porous medium, for which analytical
solutions were developed by Barry et al., (1997).

Although the proposed framework has proven relatively accurate in the
conditions examined, it comes with several limitations, all of which will be
discussed in Chapter 9. [75mm] I have something good… but it’s weird …
Francisco J. Carrillo – Cyprien Soulaine

##  Chapter 5 Applications to Multiphase-Phase Flow in Static Multiscale
Porous Media

In this chapter we focus exclusively on applying and verifying the multiphase
aspect of our model. The objective is to show that the multiscale solver
converges effectively towards its two asymptotic limits in static porous
media, namely the two-phase Darcy model at the continuum scale and the Volume-
Of-Fluid formulation at the pore scale. Additionally, we present two exemplary
cases that show the potential and wide-range applicability of our model, where
we simulate wave propagation in porous coastal barriers and
drainage/imbibition in micro-fractures. The work shown in this chapter
represents the second step toward verifying and showcasing the Multiphase DBB
framework. This chapter is adapted from Carrillo et al., (2020).

###  5.1 Darcy Scale Validation

The model’s ability to predict multiphase flow at the Darcy scale is validated
against three well-known analytical and semi-analytical solutions. Together,
these assessments test for the correct implementation of the relative
permeability, gravity, and capillary terms derived in section 2.3. This
validation follows the steps outlined in Horgue et al., (2015) for the
development and validation of their own multiphase Darcy scale solver:
impesFoam. A complete list of parameters used is provided in Tables 5.2 and
5.2.

Property | Value
---|---
Water Density | $1000\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$
Water Viscosity | $1\text{\times}{10}^{-3}\text{\,}\mathrm{P}\mathrm{a}\,\mathrm{s}$
Air Density | $1\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$
Air Viscosity | $1.76\text{\times}{10}^{-5}\text{\,}\mathrm{P}\mathrm{a}\,\mathrm{s}$
Oil Density | $800\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$
Oil Viscosity | $0.1\text{\,}\mathrm{P}\mathrm{a}\,\mathrm{s}$
Gravity | $9.81\text{\,}\mathrm{m}\,\mathrm{s}^{-2}$
Table 5.1: Table of Fluid Properties

Model Parameter | Value
---|---
$p_{c,0}$ | 100 Pa
m (Van Genuchten) | 0.5
m (Brooks-Corey) | 3
$\beta$ (Brooks-Corey) | 0.5
Table 5.2: Table of Model Parameters

####  5.1.1 Buckley-Leverett

We first consider the well-established Buckley-Leverett semi-analytical
solution for two-phase flow in a horizontal one-dimensional system with no
capillary effects (4 m long, 2000 cells, $\phi_{f}=0.5,k^{-1}_{0}=1\times
10^{11}$ $\mathrm{m}^{-2}$). In this case, water is injected into an air-
saturated reservoir at a constant flow rate with the following boundary
conditions: $\textbf{U}_{water}$ =
$1\text{\times}{10}^{-5}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$,
$\frac{\partial p_{inlet}}{\partial x}=0\
$\mathrm{P}\mathrm{a}\,\mathrm{m}^{-1}$$, and $p_{outlet}=0\
$\mathrm{P}\mathrm{a}$$. As water flows into the reservoir, it creates a
saturation profile that is characterized by a water shock at its front, an
effective shock velocity, and a saturation gradient behind said front. Figure
5.1 shows a good agreement between our numerical solutions and the semi-
analytical solutions presented in Leverett, (1940) for all three features
regardless of the chosen relative permeability model.

Figure 5.1: Comparison of the time-dependent saturation profiles calculated
from our numerical framework and Buckley-Leverett’s semi-analytical solution
for water injection into air-saturated (B) and oil-saturated (C) reservoirs.
Figure A is a visual representation of the water saturation in the reservoir
over time. Figures B and C show the semi-analytical (lines) and numerical
(symbols) solutions of the system when using the Brooks-Corey and Van
Genuchten relative permeability models, respectively.

####  5.1.2 Gravity dominated Buckley-Leverett

We then tested the exact same air-saturated system, but this time with the
addition of gravity in the same direction of the water injection velocity (see
Figure 5.2). Under these conditions, gravity becomes the dominating driving
force and the following equation can be used to calculate the water saturation
at the front (Horgue et al.,, 2015):

$\overline{\boldsymbol{v}_{l}}^{l}-\frac{k_{0}k_{r,l}(\alpha_{l}^{front})}{\mu_{l}}\rho_{l}\boldsymbol{g}=0,$
(5.1)

where the symbols are consistent with the ones presented in previous sections.
Given the parameters presented in Tables 5.2 and 5.2, Eq. 5.1 is solved
iteratively to obtain $\alpha_{l}^{front}=0.467$ and
$\alpha_{l}^{front}=0.753$ when using the Brooks-Corey and Van Genuchten
relative permeability $(k_{r,l})$ models, respectively (Appendix A). Figure
5.2 shows that our numerical solutions agree with the semi-analytical
solutions.

Figure 5.2: Comparison of the time-dependent saturation profiles calculated
from our numerical framework and the semi-analytical solution presented in
section 4.1.2. Figure A is a visual representation of water saturation in the
reservoir over time. Figures B and C show the semi-analytical (lines)
numerical (symbols) solutions of the systems parameterized through the Brooks
& Corey and Van-Genuchten relative permeability models, respectively.

####  5.1.3 Gravity-capillarity equilibrium

Lastly, we tested the validity of the capillary pressure term derived in Eqs.
2.97, 2.94, and 2.95 by solving for the steady state saturation profile of a
one-dimensional porous column filled with water and air (1 m tall, 1500 cells,
$\phi_{f}=0.5,\ k^{-1}_{0}=1\times 10^{11}$ $\mathrm{m}^{-2}$). Here, the
initial water saturation of the column is set far from its thermodynamic
equilibrium in a step-wise fashion: the lower half is partially saturated with
water ($\alpha_{water}=0.5$) while the upper half is initially dry as shown in
Figure 5.3A. To ensure proper equilibriation, both fluids are allowed to flow
freely through the column’s top boundary, but not through its lower one:
$\frac{\partial\textbf{v}_{top}}{\partial y}=0\
$\mathrm{m}\,\mathrm{s}^{-1}\,\mathrm{m}^{-1}$,\frac{\partial
p_{top}}{\partial y}=0\ $\mathrm{P}\mathrm{a}\,\mathrm{m}^{-1}$,\
\textbf{U}_{bottom}=0\ $\mathrm{m}\,\mathrm{s}^{-1}$,\ p_{bottom}=0\
$\mathrm{P}\mathrm{a}$$. For this simplified case, the theoretical steady-
state can be described as the balance between capillary and gravitational
forces, where gravity pulls the heavier fluid (water) downwards while
capillarity pulls it upwards. This behaviour can be described by the following
equation (Horgue et al.,, 2015):

$\frac{\partial p_{c}}{\partial y}=(\rho_{g}-\rho_{l})g_{y},$ (5.2)

which can be rearranged to yield:

$\frac{\partial\alpha_{l}}{\partial
y}=\frac{(\rho_{g}-\rho_{l})g_{y}}{\frac{\partial
p_{c}}{\partial\alpha_{l}}}.$ (5.3)

This last expression allows for the explicit calculation of the equilibrium
water saturation gradient by using the closed-form Brooks-Corey or Van
Genuchten capillary pressure models to obtain $\frac{\partial
p_{c}}{\partial\alpha_{l}}$ (Appendix A). Figure 5.3 shows that our numerical
model accurately replicates the results obtained from Eq. 5.3 regardless of
the chosen capillary pressure model.

Figure 5.3: Comparison of the steady state water saturation profiles
calculated from our numerical framework and the analytical solution shown in
equation 5.3. Figure A is a visual representation of the initial and final
water saturation profiles in the reservoir. Figures B and C show the steady
state saturation profiles and the resulting equilibrium saturation gradients
for both implemented capillary pressure models, respectively.

####  5.1.4 Darcy scale application: Oil drainage in a heterogeneous reservoir

As an illustration of the applicability of our model to more complex systems
at the Darcy scale, we simulate water injection into a heterogeneous oil-
saturated porous medium (1 by 0.4 m, 2000 by 800 grid, water injection
velocity =
$1\text{\times}{10}^{-4}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$,
$p_{outlet}=0\ $\mathrm{P}\mathrm{a}$$). Oil drainage is commonly used in the
energy sector, particularly as a form of enhanced oil recovery (Alvarado &
Manrique,, 2010). Although analytical solutions such as the ones presented
above are useful approximations for simple systems, they become greatly
inaccurate when modeling complex multi-dimensional systems with spatially
heterogeneous permeability. To illustrate this effect, we initialize our
reservoir’s permeability field as a grid of 0.25 by 0.1 m blocks with $k_{0}$
values ranging from $1\text{\times}{10}^{-13}\text{\,}\mathrm{t}$o
$4\text{\times}{10}^{-13}\text{\,}\mathrm{m}^{2}$ (see Figure 5.4). The
relative permeabilities within the reservoir are modeled through the Van
Genuchten model with negligible capillary effects (Table 5.2). We note that
this case was originally presented in Horgue et al., (2015) for the
development of impesFoam, a solver that uses the Implicit Pressure Explicit
Saturation (IMPES) algorithm to solve the two-phase Darcy model, making it a
convenient benchmark for comparison with hybridPorousInterFoam.

Under the aforementioned parametric conditions and with equivalent numerical
setups (i.e. same grid, time-stepping, and solver tolerances), Figure 5.5
shows that the simulations performed with hybridPorousInterFoam and impesFoam
develop very similar, yet not perfectly equivalent, saturation profiles. Of
particular interest is the development of fingering instabilities that form
due to the viscosity difference between the two fluids (Saffman & Taylor,,
1958; Chen & Wilkinson,, 1985). These instabilities are know to greatly reduce
the efficiency of enhanced oil recovery, as they essentially trap residual oil
behind the main water saturation front (Figure 5.5). Previous numerical
studies have shown that the evolution of viscous fingering is highly dependent
on the model’s hyper-parameters, grid refinement, and/or solution algorithms
(Ferrari & Lunati,, 2013; Riaz & Tchelepi,, 2006; Horgue et al.,, 2015; Chen &
Meiburg,, 1998; Holzbecher,, 2009). This characteristic explains why
hybridPorousInterFoam and impesFoam develop slightly different viscous
fingering instabilities despite having virtually perfect agreement with the
previously-presented analytical solutions: the two solvers rely on entirely
distinct sets of governing equations, boundary conditions, discretization
schemes, and pressure-solving algorithms (PISO vs IMPES). Nevertheless, this
example application shows that our solver can readily simulate complex porous
systems that have traditionally been modeled using conventional single-scale
Darcy solvers.

Figure 5.4: Simulation setup for oil drainage within a heterogeneous
reservoir. The different colored blocks represent the spatially variable
permeability field.  Figure 5.5: Oil drainage in a heterogeneous porous medium
solved at the continuum scale using hybridPorousInterFoam or impesFoam. The
white rectangular grid represent the blocks with $k_{0}$ values ranging from
$1\text{\times}{10}^{-13}\text{\,}\mathrm{t}$o
$4\text{\times}{10}^{-13}\text{\,}\mathrm{m}^{2}$ as shown in the previous
figure.

###  5.2 Pore Scale Validation

Having validated all aspects of the model within the porous domain, we now
test our model’s ability to recover known multiphase Navier-Stokes solutions
within a non-porous domain. This validation follows the steps used in previous
validations of multi-phase CFD solvers by Horgue et al., (2014), Xu et al.,
(2017), and Maes & Soulaine, (2019) and involves testing the implementation of
the imposed contact angle boundary condition against several well-known
numerical and analytical cases. Some of the simulation results obtained with
our multi-scale solver are compared with simulations performed using
_interFoam_ , the algebraic VOF solver of OpenFOAM®. In the following
simulations, we implement a static contact angle as an approximate description
of multiphase behaviour at solid interfaces, while noting the existence of
more sophisticated formulations including dynamic contact angles with viscous
bending or surface roughness (Wenzel,, 1936; Cassie & Baxter,, 1944; Voinov,,
1976; Cox,, 1986; Whyman et al.,, 2008; Meakin & Tartakovsky,, 2009).

####  5.2.1 Contact angle on a flat plate

We first test the implementation of the penalized contact angle within
hybridPorousInterFoam by initializing several “square” water droplets on a 2-D
flat porous plate with negligible permeability (6 by 2.4 mm, 480 by 192 cells,
$k_{0}^{-1}=$ $1\text{\times}{10}^{20}\text{\,}\mathrm{m}^{-2}$) and allowing
them to reach equilibrium for different prescribed contact angles
($\theta_{water}$ = $None$, $None$, $None$). These tests are compared against
equivalent droplets initialized on conventional non-porous boundaries and
solved through interFoam. Figure 5.6A shows excellent agreement between the
numerical simulations and the target equilibrium contact angle
$\theta_{water}$. The lack of a perfectly sharp interface (an intrinsic
feature of the VOF method) makes it difficult to accurately measure the
contact angle at the solid interface. However, we can confidently state that
all our results are within $None$ of the target equilibrium contact angle.
These tests are virtually identical to the ones shown in Horgue et al., (2014)
and are consistent with their results.

####  5.2.2 Capillary rise

As a second classic test for the correct implementation of multiphase flow at
the pore-scale, we model water capillary rise in an air-filled tube (1 by 20
mm, 40 by 400 cells, $\theta_{water}$ = $None$) and measure the steady-state
position of the water meniscus. To ensure a proper numerical setup, the tube’s
lower boundary is modeled as an infinite water reservoir and its upper
boundary as open to the atmosphere. To prevent initialization bias, the
meniscus is initialized about 2 mm lower than the theoretical equilibrium
height of 10 mm, which is given by the following equation (Jurin,, 1719):

$h_{eq.}=\frac{\sigma cos(\theta)}{R\rho_{l}g_{y}},$ (5.4)

where $R$ is the tube’s radius. We then numerically simulate the system with
hybridPorousInterFoam and interFoam, using impermeable porous boundaries with
the former ($k_{0}^{-1}=$ $1\text{\times}{10}^{20}\text{\,}\mathrm{m}^{-2}$)
and conventional sharp boundaries with the latter. Figure 5.6B shows the
steady state configurations of both cases, which have a meniscus height of 8.8
mm. According to Eq. 5.4, this height is equivalent to an imposed contact
angle of $None$, a small yet significant difference to the imposed contact
angle of $None$. We are not the first to note that interFoam (the standard
pore scale multiphase flow solver in OpenFOAM®) presents minor inaccuracies in
its ability to impose a prescribed contact angle (Horgue et al.,, 2014;
Gründing et al.,, 2019). The comparisons presented here show that our solver’s
accuracy in this regard is similar to that of interFoam.

Figure 5.6: Compilation of all test cases performed for the verification of
the solver within the Navier-Stokes domain. Parts A, B, and C refer to the
experiments described in Sections 5.2.1, 5.2.2, and 5.2.3, respectively. When
present, the shaded walls show the porous boundaries used in
hybridPorousInterFoam, as opposed to the standard boundary (no-slip boundary
condition at an impermeable wall) using in interFoam. For reference and easy
comparison, the white lines in Part A show the input equilibrium contact
angle.

####  5.2.3 Taylor film

We now model the drainage of ethanol ($\mu_{eth.}$
=$1.2\text{\times}{10}^{-3}\text{\,}\mathrm{P}\mathrm{a}\,\mathrm{s}$,
$\rho_{eth.}$ = $789\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$) by air in
a 2-D micro-channel (800 by 100 $\mu$m, 280 by 116 cells, $\theta_{eth.}$ =
$None$, injection velocity $U$ = 0.4 m/s, $p_{outlet}=0\
$\mathrm{P}\mathrm{a}$$). Under these circumstances, a film forms at the
channel’s boundaries as a result of competing viscous and capillary forces at
the solid interface (see Figure 5.6C). The height of this film is given by the
following analytical solution, which we use as a benchmark to verify our
numerical simulations (Aussillous et al.,, 2000),

$\frac{h_{film}}{R}=\frac{1.34Ca^{2/3}}{1+3.35Ca^{2/3}},$ (5.5)

where $Ca$ is the capillary number defined as $Ca=\frac{\mu_{eth.}U}{\sigma}$.
We can solve Eq. 5.5 with the given simulation parameters to obtain a film
thickness of 4.35 $\mu$m. Simulations of this system performed using
hybridPorousInterFoam with impermeable porous boundaries ($k_{0}^{-1}=$
$1\text{\times}{10}^{20}\text{\,}\mathrm{m}^{-2}$) and interFoam with
conventional boundaries yield a value of 4.50 $\mu$m, representing a relative
error of about 3% or 0.15 $\mu$m. These tests and their results are consistent
with numerical simulations reported by Graveleau et al., (2017) and Maes &
Soulaine, (2019) using interFoam.

####  5.2.4 Pore scale application: Oil drainage in a complex pore network

As we did at the end of the Darcy scale verification section, we now
illustrate our model’s applicability to more complex systems by presenting a
simulation of oil drainage, this time at the pore scale. The relevance of the
simulated system follows from our previous illustrative problem, as this is
simply its un-averaged equivalent at a smaller scale. The complexity of the
simulated system (1.7 by 0.76 mm, 1700 by 760 cells, water injection velocity
= 0.1 m/s, $\theta_{oil}$ = $None$, $p_{outlet}=0\ $\mathrm{P}\mathrm{a}$$)
stems from the initialization of a heterogeneous porosity field as a
representation of a cross-section of an oil-wet rock. Here, the porosity is
set to one in the fluid-occupied space and close to zero in the rock-occupied
space (See Fig. 5.7A). This allows for the solid grains to act as virtually
impermeable surfaces ($k_{0}^{-1}=$
$1\text{\times}{10}^{20}\text{\,}\mathrm{m}^{-2}$) with wettability boundary
condition (Horgue et al.,, 2014). To verify the accuracy of our solver, we
solved an equivalent system with interFoam by removing the rock-occupied cells
from the mesh and imposing the same contact angle at these new boundaries
through conventional methods.

Figure 5.7 shows that the results of the two simulations are practically
identical, down to the creation of the same preferential fluid paths and the
same droplet snap-off at 5 ms. Nevertheless, there are minor differences in
the results, where some interfaces are displaced at slightly different rates
than their counterparts (see the upper right corner at 10 ms). We attribute
these slight differences to the differing implementations of the contact angle
at the solid boundaries. We invite interested readers to find this case in the
accompanying toolbox and to refer to the extensive literature on this topic
for further discussion on numerical and experimental studies of drainage and
imbibition (Lenormand et al.,, 1988; Ferrari & Lunati,, 2013; Datta et al.,,
2014; Roman et al.,, 2016; Zacharoudiou et al.,, 2018; Liu et al.,, 2019;
Singh,, 2019) .

Figure 5.7: Oil drainage in a complex porous medium solved at the pore scale
using hybridPorousInterFoam and interFoam. The shaded sections represent solid
grains (modeled using $\phi_{f}$ = 0.001 and $k_{0}^{-1}=$
$1\text{\times}{10}^{20}\text{\,}\mathrm{m}^{-2}$ in hybridPorousInterFoam)
and the blue and red colors represent oil and water, respectively.

###  5.3 Hybrid Scale Applications

The complete body of work presented in the previous two sections verifies the
capability of our model to perform simulations of multiphase flow in complex
porous media at the pore and continuum scales. We now show how
hybridPorousInterFoam makes the simulation of hybrid scale multiphase systems
a fairly straightforward endeavor (a task that is challenging to perform with
conventional methods). The main challenge when modeling these systems can be
summarized by the following question: How can we rigorously model the porous
interface between coupled Navier-Stokes and Darcy scale domains? Although this
is still an open question, we attempt to approximate an answer by guaranteeing
three of its necessary components in the present micro-continuum framework:
first, mass conservation across the interface; second, continuity of stresses
across the interface and; third, a wettability formulation at the interface.
The first two components are intrinsic features of the solver which have been
proven necessary and sufficient to accurately model single phase flow in
hybrid scale simulations (Neale & Nader,, 1974) and have been used as closure
conditions in previous multiphase models (Lacis et al.,, 2017; Zampogna et
al.,, 2019). The latter, as explained in the pore scale validation section, is
roughly approximated through a constant contact angle boundary condition. We
recognize that these components represent an approximation to the complete
description of the boundary. Nevertheless, to the best of our knowledge, there
does not exist a better way to model this interfacial behavior, a testament to
the novelty and potential of the proposed modeling framework.

The following illustrative cases are used to show our model’s capability to
simulate multiphase systems at the hybrid scale. They are also included as
tutorial cases in the accompanying toolbox.

####  5.3.1 Wave propagation in a coastal barrier

Coastal barriers are used throughout the world to prevent flooding, regulate
water levels, and protect against inclement weather (Morton,, 2002). Accurate
simulation of water interaction with these barriers is challenging as it
requires predicting the behavior of open water at large scales (Navier-Stokes)
while also resolving small-scale multiphase effects within the barrier itself
(Darcy flow).

Figure 5.8: Coastal barrier simulation at different times. The thin black line
represents the boundary of the porous solid: $\phi_{f}$ and $k_{0}^{-1}$ are
set to 0.5 and $5\text{\times}{10}^{7}\text{\,}\mathrm{m}^{-2}$ below said
line and to 1 and 0 above it, respectively. The 2-D representation shows a
plane that cuts through the middle of the 3-D simulation.

We created a three-dimensional coastal barrier (8.3 by 2.7 m by 0.25 m, 1660
by 540 by 50 cells) by initializing a heterogeneous porosity field in the
shape of a barrier ($k_{0}^{-1}=$
$5\text{\times}{10}^{7}\text{\,}\mathrm{m}^{-2}$, $\phi_{f,barrier}$ = 0.5)
and by setting the water level such that it partially covers the barrier (see
Figure 5.8). In this particular case, we chose not to impose a contact angle
at the barrier-water interface as its effects would be minimal when compared
to macroscopic gravitational effects (Bond Number
$=\frac{\Delta\rho(\textnormal{Length Scale})^{2}g_{y}}{\sigma}>>1$). To
ensure proper initialization, we allowed the water saturation profile on the
above-water section of coastal barrier to reach its capillary-induced steady
state (similarly to the capillary rise simulation presented in Section 5.2.2).
This process was modeled using the Van Genuchten relative permeability and
capillary pressure models (m = 0.8, $p_{c,0}$ = 1000 Pa). After equilibration,
we started the simulation by initializing a wave as a rectangular water
extrusion above the water surface. To ensure proper wave propagation behavior,
we tuned the simulation’s numerical parameters (discretization schemes, linear
solvers, time stepping strategy) according to the guidelines established in
Larsen et al., (2019).

The results from this simulation show that we can simultaneously model coupled
wave and Darcy dynamics in three dimensions. The snapshots shown in Figure 5.8
illustrate how water saturation within the porous domain is controlled by the
crashing of waves, gravity, and capillary effects. The associated wave
absorption and dissipation cycle brought about by the porous structure is
repeated every few seconds with lowering intensity until the initial
configuration is eventually recovered. To the best of our knowledge, Figure
5.8 shows the first existing numerical simulation coupling multiphase flow
with real capillary effects at two different scales without the use of
different meshes, solvers, or complex interfacial conditions. Other models
such as olaFlow have been developed to simulate similar wave dynamics with
coastal barriers (Higuera et al.,, 2013). Many of these models rely on the
assumption that the pores within the coastal barrier are large ($>$10 cm),
meaning that they can reasonably ignore capillary effects within the porous
medium. Contrastingly, our model makes no such assumption, meaning it can also
be used to model coastal barriers with arbitrarily small pores (such as in
sand or gravel structures) and also should be applicable to other types of
groundwater-surface water interaction (Maxwell et al.,, 2014).

####  5.3.2 Drainage and imbibition in a fractured microporous matrix

A second conceptually similar, yet completely different hybrid scale
application of hybridPorousInterFoam involves the injection of fluids into
fractured porous materials. Accurately capturing the fluid behavior in these
systems is especially challenging due to the fact that it requires accounting
for multiphase effects simultaneously within the fracture (Navier-Stokes), in
the surrounding microporous matrix (Darcy), and at the porous boundary (the
contact angle implementation).

Here, we model drainage and imbibition in a water-wet fracture system, where
we inject air into a 90% water-saturated microfracture in the former and we
inject water into a 90 % air-saturated microfracture in the latter (1.2 by 0.5
mm, 1200 by 500 cells, $\theta_{water}$ = $None$, fluid injection velocity =
0.1 m s-1, $p_{outlet}=0\ $\mathrm{P}\mathrm{a}$$). The relative
permeabilities and capillary pressures in the heterogeneously-initiated porous
domain ($\phi_{f}=0.5,$ $k_{0}^{-1}=$
$4\text{\times}{10}^{12}\text{\,}\mathrm{m}^{-2}$) are modeled through the
Brooks-Corey model with m = 3, $p_{c,0}=100\ $\mathrm{P}\mathrm{a}$$, and
$\beta=0.5$.

Figure 5.9 presents the results of these simulations and illustrates how
strongly multi-scale wettability effects can influence simulations results. In
both cases, the injected fluid is able to invade the microporous matrix, but
the mechanism through which it does is completely different. In the case of
water injection (imbibition), the wetting contact angle boundary condition
encourages complete water saturation of the whole fracture such that air is
completely displaced by time = 125 ms. Furthermore, throughout the whole
process, the microporous capillary pressure acts as an additional driving
force for water invasion into the surrounding microporous matrix, leading to
the almost complete saturation of the whole system by time = 500 ms. This
process has some strong analogies to the flow of water in into hygroscopic
materials(Zhou et al.,, 2019).

The drainage case is slightly less intuitive, yet conceptually more
interesting. Here, the contact angle and microporous capillary pressures act
against the invasion of air into the fracture and into the surrounding porous
material, respectively. The result is that the air cannot effectively displace
water from the fracture, leading to the trapping of water droplets in fracture
ridges. Initially, these droplets act as barriers that prevent air entry into
the porous matrix (see time = 125 ms). However, as the flow-induced pressure
gradient pushes air into the porous matrix, the water saturation in the pores
surrounding the droplets decreases. The system then responds by increasing the
capillary pressure at the porous interface, which eventually leads the water
droplets to imbibe into the matrix. Lastly, we highlight the clear time scale
separation between the imbibition and drainage cases, as the invading
interface progresses about three times more slowly within the microporous
matrix in the latter case.

Figure 5.9: Drainage and imbibition in a microporous fracture. Shades of blue
and red represent of the degree of air and water saturation, respectively. The
thin white line shows the fracture outline (i.e. the fluid-solid interface),
which separates the open fracture ($\phi_{f}=1,k_{0}^{-1}=0$) from the porous
fracture walls ($\phi_{f}=0.5,$ $k_{0}^{-1}=$
$4\text{\times}{10}^{12}\text{\,}\mathrm{m}^{-2}$) located above and below it.

Several similar dual porosity models have been proposed to model the types of
effects illustrated in Figure 5.9, but never in this way or to this degree of
detail (Douglas et al.,, 1991; Di Donato et al.,, 2003). Many of these models
rely on a description of fractures as single-dimensional features with high
porosity and permeability values within a pure Darcy scale simulation (Nandlal
& Weijermars,, 2019; Yan et al.,, 2016). Although very useful, many of these
simulations ignore the geometric capillary effects and non-linear couplings
presented above. Our approach can therefore be seen as the missing link
between pore-scale modeling and Discrete Fracture Networks (Karimi-Fard &
Durlofsky,, 2016) and as a useful tool for refining the transfer functions
used in these large scale models.

###  5.4 Conclusions

We have successfully tested and verified a multiscale model for two-phase flow
in porous media. This modeling framework and its open-source implementation
hybridPorousInterFoam can be used to simultaneously model multiphase flow at
two different length scales: a Darcy Scale where sub-voxel fluid-fluid
interactions within a porous medium are modeled through relative permeability
and capillary pressure constitutive models and a pore scale (or Navier-Stokes
scale) where the solid material is non-porous and fluid-fluid interactions are
modeled through a continuum representation of the Young-Laplace equation.
Furthermore, our model is able to do this through the use of a single momentum
conservation equation without the need to define different meshes, separate
solvers/domains, or complex interfacial conditions. The proposed framework is
an accurate and straightforward way to introduce the physics of two-phase flow
in porous media in CFD softwares.

Through this study, we showed that our model can successfully simulate
multiphase Darcy and Navier-Stokes flow to the same standard (and with the
same assumptions and limitations) as conventional single-scale solvers
impesFoam and interFoam. The coupling between the two scales at porous
interfaces is handled by ensuring mass conservation and continuity of stresses
at said boundary, as well as by implementing a constant contact angle
wettability condition. We then leveraged all these features to show that our
model can be used to model hybrid scale systems such as wave interaction with
a porous coastal barrier and drainage and imbibition in a fractured porous
matrix.

Although the proposed formulation represents a significant advance in the
simulation of multiscale multiphase systems, we note that further study is
required in particular to properly and rigorously model the multi-scale porous
interface. The implemented interface, as it stands, has been shown to
accurately predict single phase flow into porous media (Neale & Nader,, 1974),
impose static contact angles over porous boundaries (Horgue et al.,, 2014),
and approximate multiphase flow in porous media (Lacis et al.,, 2017).
However, its accuracy when modeling multiphase flow at rough porous interfaces
is still an open question, as there does not currently exist a rigorous
formulation to model such behaviour. The derivation, implementation, and
verification of such a boundary condition and the inclusion of erosion,
chemical reactions (Soulaine et al.,, 2017, 2018), and solid mechanics
(Carrillo & Bourg,, 2019) into this framework will be the focus of subsequent
chapters and investigations. [75mm] All models are wrong, but some are useful
George E. P. Box

##  Chapter 6 Application to Multiphase Flow in Multiscale Deformable Porous
Media

Most of the underlying components of the Multiphase DBB approach described in
Chapter 2 have been already tested and verified. Chapter 4 validated the
momentum exchange terms as an effective coupling mechanism between a single
fluid phase and a deformable plastic or elastic porous medium. The effects of
confining and swelling pressures on porous media were also examined in said
study. Then, Chapter 5 extensively validated the extension of the Darcy-
Brinkman equation into multiphase flow within and around static porous media
by comparison with reference test cases in a wide range of flow, permeability,
capillarity, and wettability conditions. Therefore, the only thing left to
validate is the ability of the Multiphase DBB model to accurately predict the
behavior of multiscale systems that exhibit coupling effects between multiple
fluids and a deformable porous matrix. This chapter focuses on addressing this
last point by using our model to replicate known solutions to coupled
problems, with a particular emphasis on fracture mechanics. Additionally, it
presents two exemplary cases that show the potential and wide-range
applicability of our coupled model, where we simulate wave propagation in
poroelastic coastal barriers and surface cracking/subsidence as a result of
subsurface hydraulic fracturing. The work presented here is adapted from
Carrillo & Bourg, 2021b  and represents the third and last step in the journey
of verifying and showcasing the Multiphase DBB model.

###  6.1 Model Validation

We begin with two validation cases relating to multiphase poroelasticity and
the coupling between solid deformation and fluid pressure. Then, we proceed
with two poroplastic cases that validate this framework for multiscale plastic
systems. Finally, we conclude with two additional cases that verify the
implementation of the capillary force interaction terms. These cases can be
found in the accompanying CFD simulation package.

####  6.1.1 Terzaghi Consolidation Problem

The Terzaghi uniaxial compaction test has been extensively used as a benchmark
for the validation of numerical codes relating to poroelasticity (Terzaghi et
al.,, 1996). Its main utility is to test the accuracy of the solid-fluid
couplings that relate fluid pressure to solid deformation and vice versa. The
problem consists of a constrained saturated elastic porous medium that is
abruptly compressed from its upper boundary by a constant uniaxial load
(Figure 6.1). This creates a sudden increase in pore pressure, which is then
dissipated by flow through the upper boundary (all other boundaries have
impermeable boundary conditions). In the case of a one-dimensional porous
medium, the resulting temporal and spatial evolution in fluid pressure can be
described by the following simplified analytical solution (Verruijt,, 2013).

$\displaystyle\frac{p}{p_{max}}={\mathrm{erf}\left(\frac{h-z}{2\sqrt{c_{v}t}}\right)\
}\ \ \ \mathrm{for\ }\ \ \frac{c_{v}t}{h^{2}}\ll 1$ (6.1)

where ${c}_{v}=({k}_{0}E\ (\nu-1))/(\eta(2{\nu}^{2}+\nu-1))$ is the
consolidation coefficient, ${k}_{0}$ is permeability, $E$ is Young’s modulus,
$\nu$ is Poisson’s ratio, $\eta$ is the fluid’s unit weight, $h$ is the column
height, and $z$ is the vertical coordinate. Our equivalent numerical setup is
shown in Figure 6.1. The values of the relevant parameters in our simulations
are $h=10\ \mathrm{m}$, ${k}_{0}=5\times{10}^{-11}\ \mathrm{m^{2}}$, $E=2\
\mathrm{MPa}$, and $\nu=0.25$. To show the accuracy of our model across
different conditions, the loading pressure was varied from $10$ to $200\
\mathrm{kPa}$ (Figure 6.1B) and the porosity from $0.25$ to 0.75 (Figure
6.1C). Lastly, the column was partially saturated $({\alpha}_{w}=0.5)$ with
fluids with equal densities $({\rho}_{f}=1000\ \mathrm{kg/m^{3}})$,
viscosities $(\mu_{f}=1\ \mathrm{cp})$, and negligible capillary effects. This
last point allowed for testing the validity of the fluid-solid couplings
irrespective of the simulated phases without violating any of the assumptions
present in the analytical solution. Our numerical results show excellent
agreement with Equation 6.1 for all tested conditions.

Figure 6.1: One-dimensional Terzaghi consolidation problem. (A) Simulation
setup. (B) Analytical (solid lines) and numerical (symbols) pressure profiles
at $t=100\ \mathrm{s}$ for different loading pressure values. (C) Time-
dependent pressure profiles for different column porosity values (From top to
bottom: ${\phi}_{s}=0.75,\ 0.5,\ 0.25$).

####  6.1.2 Pressure Oscillation in Poroelastic Core

This verification quantifies the effects of the seismic stimulation of a
poroelastic core saturated with water and trichloroethene (TCE). Our
simulations follow the experimental and numerical set up described in Lo et
al., (2012), where a horizontal one-dimensional sand core ($0.3\ \mathrm{m}$
long, $30\times 1$ grid cells, ${\phi}_{f}=0.5$, ${\alpha}_{w}=0.9$,
${k}_{0}=1.1\times{10}^{-11}\ \mathrm{m^{2}}$) is subjected to constant
uniaxial compression and oscillatory pore pressure variations imposed by time-
dependent boundary conditions (Figure 6.2). In this case, flow is allowed
through both boundaries, which results in a system that continuously undergoes
a relaxation-compression cycle. The ensuing cyclical change in the core’s
fluid content as a function of time can be described by a semi-analytical
solution first derived in Lo et al., (2012) and reproduced in Appendix E.

Figure 6.2: Change in fluid content of an oscillating poroelastic core. (A)
Simulation setup. (B) Semi-analytical (solid lines) and numerical solutions
(symbols) for the percent change in the core’s fluid volume as a function of
time.

For our matching simulations, the porous structure’s Young’s modulus was set
to $E=53\ \mathrm{MPa}$ and its Poisson ratio to $\nu=0.32$. Here, water
density was ${\rho}_{w}=1000\ \mathrm{kg/m^{3}}$, water viscosity was $1\
\mathrm{cp}$, TCE density was ${\rho}_{TCE}=1480\ \mathrm{kg/m^{3}}$, and TCE
viscosity was ${\mu}_{TCE}=0.57\ \mathrm{cp}$. Furthermore, the pressure at
the left boundary was held at $p=1\ \mathrm{kPa}$ while the pressure at the
right boundary was set by $p={p}_{0}{sin\left(2\pi ft\right)}$, with
${p}_{0}=1$ \- $2\ \mathrm{MPa}$ and $f=35$ \- $70\ \mathrm{Hz}$. Lastly, the
core was uniaxially compressed through a constant stress of $1\ \mathrm{kPa}$
applied at both boundaries. A comparison between our numerical solutions and
Lo’s semi-analytical solution is presented in Figure 6.2, yielding excellent
agreement for all tested cases.

We note that the Multiphase DBB formulation should be able to describe ‘Slow”
Biot pressure waves caused by the relative motion of the solid and fluid
phases which occurs at much higher frequencies than the ones simulated here
(i.e. 10 MHz). However, capturing these effects and modelling
“Fast/Compressional” pressure waves would require the implementation of a
pressure-velocity coupling algorithm that allows for compressible flow (Lo et
al.,, 2012). Such an endeavour is outside the scope of this dissertation.

####  6.1.3 Capillary Pressure Effects in a Poroelastic Column

Having verified the two-way coupling between solid deformation and fluid
pressure, we now verify the implementation of the capillary pressure terms
within the solid mechanics equation. To do so, we simulate a poroelastic
column ($1$ m tall, $1500$ cells, $\phi_{f}=0.5$) bounded by two non-wetting
fluid reservoirs at its upper and lower boundaries. The column is initialized
with a linear saturation profile spanning from $\alpha_{w}=0$ to $1$ (see Fig.
6.3). Fluid saturation is kept fixed by not solving Equation 2.89, and the
mobilities of both fluids are set to very high values ($M_{i}=1\times 10^{10}\
\mathrm{m^{3}/kg.s}$) to minimize drag-related effects. Under these
conditions, the solid’s effective stress is exclusively controlled by
capillary effects and is described by the following analytical solution:

$\mathrm{Effective\ Stress}=\phi_{s}\times\alpha_{w}\times p_{c}$ (6.2)

We used the Van Genutchen capillary pressure model with $m=0.6$ or $0.8$ and
$p_{c,0}=50$ to $2000$ Pa to calculate the solutions to said problem. The
resulting agreement between the numerical and analytical solutions, shown in
Fig. 6.3, confirms the accuracy of the fluid-solid capillary pressure coupling
implemented in our model. Furthermore, the transitional behaviour of the
effective stress at the macroscopic solid-fluid interface confirms the
applicability of the interfacial condition described in Section 2.11: as
expected, solid stresses are dictated by standard elasticity theory in the
porous region and become negligible in solid-free regions.

Given that the fluid-solid couplings in a poroelastic solid are now verified,
we proceed to verify said terms for poroplastic materials.

Figure 6.3: Capillary effects in a poroelastic column. (A) Simulation setup.
(B $\&$ C) Analytical (solid lines) and numerical (symbols) effective stress
profiles for different capillary pressure values ($p_{c,0}=50$ to $2000$ Pa)
and Van Genuchten coefficients ($m=$ $0.6$ and $0.8$).

####  6.1.4 Fluid Invasion and Fracturing in a Hele-Shaw Cell

The third verification case (and the first poroplastic case) consists in the
qualitative replication of a set of fracturing experiments that examined the
injection of aqueous glycerin into dry sand within a 30 by 30 by 2.5 cm Hele-
Shaw cell (Huang et al., 2012b, ; Huang et al., 2012a, ). These experiments
are inherently multiscale, in that the characteristic length scale of
fractures in this system ($\sim\mathrm{cm}$) is orders of magnitude larger
than that of pores within the microporous matrix ($\sim 100\ \mathrm{\mu m})$.
They are also multiphysics, as they clearly exemplify the drag-controlled
transition from Darcy flow within the porous medium to Stokes flow in the open
fractures and the coupling between the hydrodynamics of fluid flow and the
mechanics of fracture propagation (Fig. 6.4).

The experimental setup involved the injection of aqueous glycerin at various
flow rates $q$ between $5$ and $50$ ml/min while also varying the fluid’s
viscosity ${\mu}_{gly}$ between $5$ and $176\ \mathrm{cp}$ for different
experiments. Our numerical simulations were parameterized using measured
values of the glycerin-air surface tension ($\gamma=0.063\
\mathrm{kg/s^{2}}$), the density of pure glycerin (${\rho}_{gly}=1250\
\mathrm{kg/m^{3}}$), the density of air (${\rho}_{air}=1\ \mathrm{kg/m^{3}}$),
the viscosity of air (${\mu}_{air}=0.017\ \mathrm{cp}$), and the average
radius and density of sand grains ($100\ \mathrm{\mu m}$ and $2650\
\mathrm{kg/m^{3}}$, respectively). To mimic the sand’s experimental
configuration and permeability, the simulated solid fraction field was set to
a random initial normal distribution such that ${\phi}_{s}=0.64\ \pm 0.05$ and
the permeability was modelled as a function of the solid fraction through the
Kozeny-Carman relation with ${k}_{0}=6.7\times{10}^{-12}\ \mathrm{m^{2}}$.
Relative permeabilities were calculated through the Van Genutchen model with
the Van Genuchten coefficient $m$ set to $0.99$ (see Appendix A), while
capillary pressures were deemed negligible (as
$2\gamma{r}^{-1}\ll\mu{k}^{-1}{U}_{f}L)$. Finally, the porous medium was
modeled as a continuous Hershel-Bulkley-Quemada plastic (Appendix B) with
kinematic yield stress ${\tau}_{0}=16.02\ \mathrm{m^{2}/s^{2}}$ (Quemada,,
1977). Plasticity was used as the preferred mode of solid rheology due to its
ability to account for the compressive and irreversible effects caused by
fracturing within these experiments (Abou-Sayed et al.,, 2004; van Dam et
al.,, 2002).

Numerically speaking, the simulations were carried out in a 30 by 30 cm 2-D
grid (500 by 500 cells) with constant velocity and zero-gradient pressure
boundary conditions at the inlet, zero-gradient velocity and zero pressure
boundary conditions at the boundary walls, and a solid velocity tangential
slip condition at all boundaries (i.e. the solid cannot flow across the
boundaries, but the fluids can). Lastly, to enable a closer comparison between
our 2D simulation and the 3D experiment we added an additional drag term to
the fluid momentum equation equal to $12\mu{a}^{-2}{\boldsymbol{U}}_{f}$,
which accounts for viscous dissipation through friction with the walls in a
Hele-Shaw cell with aperture $a$ (Ferrari et al.,, 2015).

Figure 6.4: Comparison of experimental (A, B, C, G, H, I) and simulated (D, E,
F, J, K, L) fracturing in a Hele-Shaw cell. The color bar represents the solid
fraction within the simulations (where red implies a pure solid and blue pure
fluids) and the black lines represent the advancing glycerin saturation front.
The experiments shown here are part of the results presented in Huang et al.,
2012b .

As shown in Figure 6.4, a dramatic transition in the mode of fluid invasion is
observed with increasing fluid injection velocity and viscosity. At low flow
rates and low viscosity $\left(q=5\ \mathrm{ml/min},\ \mu=5\
\mathrm{cp}\right)$, there is no discernible solid deformation and the main
mode of fluid flow is through uniform invasion of the porous medium (Figure
6.4A). At intermediate flow rates and low viscosity ($q=25\ \mathrm{ml/min}$
to $30\ \mathrm{ml/min},\ \mu=5\ \mathrm{cp}$), we still observe a uniform
invasion front, but small fractures begin to appear (Figure 6.4B, C). At high
viscosity $\left(\mu=176\ \mathrm{cp}\right)$, we see clear fracturing
patterns preceded by a non-uniform fluid invasion front (Figure 6.4H, I).

Figure 6.4 shows that our simulation predictions are qualitatively consistent
with the experiments presented in Huang et al., 2012b  with regard to both the
stability of the capillary displacement front and the observed fracturing
transition behavior. As suggested above, accurate prediction of this
transition requires not only proper handling of fluid-fluid interactions
(surface tension and relative permeability effects), but also accurate
descriptions of their relationship with solid mechanics (drag) and the proper
implementation of a solid rheological model that can replicate irreversible
and unstable fracturing processes. We note that in our simulations, fracture
initialization and propagation are predicted based on continuum-scale
equations for the rheology and mechanics of the bulk microporous solid, with
no specific treatment of grain-scale mechanics. Grid-level instabilities are
brought about by the normally distributed porosity and permeability fields, as
shown in Appendix C. The microstructural differences between the experiments
and our simulations (most clear in Figure 6.4C, F, and H, K) likely arise at
least in part from the fact that the solid is modelled as a continuum rather
than a granular material.

This section demonstrates that the multiphase DBB model can be used to
replicate and predict the main mode of fluid flow and solid deformation within
fracturing systems. A comprehensive study of the controlling parameters for
multiphase fracturing in the presence of both viscous and capillary stresses
will be the focus of the next chapter (Chapter 7).

####  6.1.5 Modeling Fracturing Wellbore Pressure

Having shown that our model can qualitatively predict fracturing behavior, we
now aim to determine whether it can do so in a quantitative matter. As
depicted in Figure 6.5, fluid-induced fracturing of low-permeability rocks
proceeds through the following well-established series of stages: First, fluid
pressure increases linearly as fracturing fluid is injected into the wellbore.
Second, as wellbore pressure increases and approaches the leak-off pressure, a
small amount of pressure is propagated by fluid leakage into the rock. Third,
fluid pressure continues to increase until it reaches the breakdown pressure,
at which point it is high enough to fracture the rock. Fourth, a fracture is
initiated and propagates; the wellbore pressure slowly decreases. Fifth,
injection stops, fracture propagation stops, and wellbore pressure rapidly
dissipates (Abass et al.,, 1996; Abou-Sayed et al.,, 2004; Huang et al.,
2012b, ; Papanastasiou,, 2000; Santillán et al.,, 2017).

Figure 6.5: Conceptual representation of wellbore pressure evolution during
fluid-induced fracturing of low permeability rocks. In this section, we are
interested in modeling the behavior between $t_{frac}$ and $t_{stop}$.

In this section we aim to numerically replicate the time-dependent fracturing
wellbore pressure during fracture propagation (i.e., the fourth stage outlined
above) as described by an analytical solution presented in Barros-Galvis et
al., (2017).

$p_{well}=p_{0}-\frac{\mu q}{4\pi
k_{0}h}\left[ln\left(\frac{tk_{0}\tau_{0}}{\phi_{f}\mu
r_{well}^{2}}\right)+0.81\right]$ (6.3)

where $t$ is the time elapsed since fracture initialization, $q$ is the fluid
injection rate, $p_{well}$ is the wellbore pressure, $p_{0}$ is the minimum
pressure required for starting a fracture (a function of the solid’s yield
stress $\tau_{0}$), $h$ is the formation thickness, and $r_{well}$ is the
wellbore radius. The remaining variables follow the same definitions described
earlier.

The general numerical setup is almost identical to the one presented in the
previous section. The key difference is that we now inject aqueous glycerin
into a strongly-non wetting (and thus almost impermeable) porous material.
This is done to ensure an accurate replication of the analytical solution and
its related assumptions, where fracturing is the main mode of fluid flow and
there is virtually no fluid invasion into the porous matrix. The exact
simulation parameters are $q=46$ to $110\ \mathrm{ml/min}$, ${\tau}_{0}=0.2$
or $2\ \mathrm{m^{2}/s^{2}}$, ${k}_{0}=6.7\times{10}^{-11}$ or
$6.7\times{10}^{-12}\ \mathrm{m^{2}}$, ${\mu}_{gly}=5\ \mathrm{cp}$, and
$m=0.05$. Note that low values of $m$ indicate that the porous formation is
strongly non-wetting to the injected fluid. All other parameters are as in the
previous section.

Lastly, as hinted at before, a notable characteristic of our model is that
different normally-distributed solid fraction field initializations give
different fracturing results (Appendix C). For this reason, we performed four
simulations for each parameter set. In Figure 6.6, we present the average
predicted wellbore pressure evolution with errors bar representing the 95%
confidence interval.

Figure 6.6: Wellbore pressure as a function of injection rate and time. (A)
The initial simulation setup showing the initial wellbore radius
${r}_{well}=1.3\ \mathrm{cm}$, as well as the normally distributed solid
fraction field. (B) The fractured system, where the thin black line represents
the position of the advancing glycerin saturation front. C and D show the
wellbore pressure as a function of time for different flow rates and different
combinations of solid yield stress and permeability. Solid curves represent
analytical solutions, while symbols represent simulation predictions. The
color scheme in A and B is the same as in Figure 6.4, and ${p}_{max}$ is the
maximum analytically-predicted pressure in each simulation.

Figure 6.6 shows that our model can accurately and reliably predict the
pressure and deformation behavior of a variety of fracturing systems, as all
curves exhibit excellent agreement with their respective analytical solution.
Note that the length of each curve relates inversely to the injection speed.
This is because fractures at higher injection rates consistently reach the
system’s boundary faster than their counterparts, at which point there is a
sharp decrease in pressure and the analytical solution no longer applies.
Therefore, each curve’s cutoff point represents the time at which the fracture
effectively becomes an open channel between the wellbore and the outer
boundary, normalized to the average value of that time for the slowest-moving
fracture (i.e. $t=t_{max}$).

The successful replication of the analytical pressure profiles in this section
verifies the model components pertaining to the pressure-velocity-deformation
coupling and the two-way momentum transfer between the fluid and solid phases
(drag). Therefore, the only model component left to verify is the
implementation of the capillary force terms during fracturing of a plastic
solid.

####  6.1.6 Capillary Effects on Fracturing Wellbore Pressure

Our fifth verification systematically varies the capillary entry pressure
within non-wetting fracturing systems to quantify its effects on wellbore
pressure. For this, we consider two different complementary cases: one where
capillary forces are comparable to their viscous counterparts, and another
where they are significantly larger than them. All parameters are the same as
in the previous experiments (Section 6.1.5) unless otherwise specified.

The first set of experiments expands the previous section’s analysis into
strongly non-wetting systems with the addition of a constant capillary
pressure jump at the fracture interface imposed by a flat capillary pressure
curve ($p_{c}=p_{c,0}=0$ to $2$ kPa, $\tau_{0}=2\ \mathrm{m^{2}/s^{2}},\
{k}_{0}=6.7\times{10}^{-12}\ \mathrm{m^{2}},\ m=0.05$, and $q=78\
\mathrm{ml/min}$ ). In this case, all the assumptions present in the
fracturing analytical solution (Eqn. 6.3) are satisfied. However, said
solution still does not account for capillarity. For constant flow in non-
wetting systems, the addition of a constant capillary entry pressure jump at
the fluid-solid interface would increase the calculated propagation pressure
in Eqn. 6.3 by said value such that ${p}^{new}_{well}={p}_{well}+{p}_{c}$.
This effect is exemplified in Figure 6.7A, where we present the updated
analytical results in conjunction with our equivalent numerical results,
demonstrating excellent agreement between them. Note that the predicted linear
relationship between wellbore pressure and capillary entry pressure is not
explicitly imposed in the numerical model. On the contrary, it arises
naturally from the balance of viscous, capillary, and structural forces in
Eqns. 2.90-2.92.

Figure 6.7: Effect of capillary entry pressure on fracturing wellbore
pressure. (A-B) Wellbore pressure as a function of time and entry pressure for
low and high permeability systems, respectively. In B, curves at increasingly
high pressures were cut off for illustrative purposes and the solid line
represents a fitted reference logarithmic pressure descent curve. (C-H) Time
evolution of fractured system with a 1 kPa capillary entry pressure and high
permeability. (C) Initial fluid invasion ($t/t_{max}<0$): at early times the
wellbore pressure rises rapidly and becomes larger than the entry capillary
pressure. The fluid invades the porous formation symmetrically. (D) Fracture
initiation ($t/t_{max}=0$): The wellbore pressure continues to rise until it
is larger than the breakdown pressure, at which point small fractures start to
form. Fluid invasion continues. (E-F) Fracture propagation ($t/t_{max}>0\mid$
$p_{well}>p_{c,0}$): the wellbore pressure drops as fractures propagate. Fluid
invasion continues asymmetrically around said fractures. (G) Fluid invasion
stops ($t/t_{max}>0\mid$ $p_{well}\sim p_{c,0}$): As the wellbore pressure
keeps dropping, the entry capillary pressure condition at the porous interface
ensures that that wellbore pressure never goes below $p_{c,0}$, at which point
fluid invasion stops. (H) Fracture reaches the simulation boundary
($t/t_{max}=1$). The color convention in Figures C-H is the same as in Figure
6.4.

The second set of experiments modifies the previous experiments by making the
porous medium significantly more permeable, while still maintaining a constant
capillary pressure jump at the fracture interface ($p_{c}=p_{c,0}=1$ to $3$
kPa, ${\tau}_{0}=0.2\ \mathrm{m^{2}/s^{2}},\ k_{0}=6.7\times{10}^{-11}\
\mathrm{m^{2}},\ m=0.99$, and $q=78\ \mathrm{ml/min}$). This results in a set
of cases where the wellbore pressure is increasingly controlled by the
capillary pressure drop rather than by the viscous pressure drop across the
fracture and porous formation.

Figure 6.7 demonstrates precisely this effect. Our simulations show that the
wellbore pressure always decays towards the capillary entry pressure once
viscous effects are dissipated by fracture growth, i.e., we observe a
transition between viscous- and capillary-dominated regimes. At low values of
${p}_{c,0}$ $\left(<2500\ \mathrm{Pa}\right)$ the entry pressure is not high
enough to prevent fluid flow into the surrounding porous matrix during
fracturing (Figure 6.7B-H). The resulting pressure drop cannot be modeled by
the previously presented analytical solution (as it violates the no leak-off
assumption), but still follows a logarithm-type curve that is characteristic
of flow in fracturing systems. With increasing fracture propagation, the
viscous pressure drop decreases until the wellbore pressure equals the entry
pressure, which is, by definition, the minimum pressure drop required for
fluid flow in highly permeable non-wetting systems. Finally, we note that in
cases where capillary entry pressure is high relative to the pressure required
to fracture the solid (i.e., at ${p}_{c,0}>2.250\ \mathrm{Pa}$ in the
conditions simulated in Fig. 6.7b), fracturing begins before the wellbore
pressure can exceed ${p}_{c,0}$. This prevents essentially all flow into the
porous formation, and the wellbore pressure is immediately stabilized at
$\sim{p}_{c,0}$. For all cases, fractures continue to propagate until they
reach the system boundary, at which point the pressure drops rapidly as noted
in Section 6.1.5.

In this section we reduced the inherent complexity of the model’s capillary
force terms $F_{c,i}$ (Eqns. 2.94-2.95) into a simple set of intuitive
verifications. The quantitative agreement between these two analytical cases
and their corresponding numerical simulations validate the implementation of
the impact of capillary pressure effects on the mechanics of a ductile porous
solid within our model.

###  6.2 Illustrative Applications

Having verified and tested the model, we now proceed with two illustrations
that demonstrate how hybridBiotInterFoam enables the simulation of relatively
complex coupled multiphase multiscale systems. The following cases serve as
illustrative examples of our model’s features and capabilities as well as
tutorial cases within the accompanying toolbox.

####  6.2.1 Elastic Failure in Coastal Barriers

Coastal barriers are ubiquitous features in coastal infrastructure
development. When designed appropriately, these structures can be very
effective in regulating water levels and protecting against inclement weather
(Morton,, 2002). However, accurate prediction of the coupled fluid-solid
mechanics of these structures (which can lead to barrier failure) is
inherently challenging as it requires modeling large-scale features (waves)
while also considering small-scale viscous and capillary interactions within
the barrier.

The following case represents the continuation of the three-dimensional
coastal barrier illustration presented in Chapter 5 with the addition of
linear-elastic poromechanics. As such, the simulation was created by
initializing a heterogeneous porosity field (with $k_{0}=2\times{10}^{-8}\
\mathrm{m^{2}}$ and ${\phi}_{f}=0.5$) in the shape of a barrier within a 8.3
by 2.7 by 0.25 m rectangular grid (1600 by 540 by 50 cells). The relevant
solid mechanics parameters were $E=5\ \mathrm{MPa}$, $\nu=0.45$, and
${\rho}_{s}=2350\ \mathrm{kg/m^{3}}$. Relative permeabilities and capillary
pressures were evaluated through the Van Genuchten model with $m=0.8$ and
$p_{c,0}=1\ \mathrm{kPa}$. Before the start of the simulation, the water level
was set to partially cover the barrier and then allowed to equilibrate. A
single wave was then initialized at $t=0$. This results in a simulation that
exhibits a clear wave absorption cycle that gradually dissipates in time, as
seen in Figure 6.8. Detailed discussion on the fluid mechanics of this problem
can be found in Chapter 5 and Carrillo et al., (2020).

Figure 6.8: Waves crashing against a poroelastic coastal barrier. Here, the
thin black line represents the water-air interface
$\left({\alpha}_{w}=0.5\right)$ and red-blue colors outside the coastal
barrier represent water and air, respectively. Colored contours within the
barrier are the calculated Von Mises stresses and are shown in 5 kPa
increments in the general downwards direction. Note that the largest stresses
are seen during the initial wave crash and increase towards the base of the
barrier due to gravitational effects.

Here, however, we are interested in evaluating the barrier’s propensity to
failure. We do this by applying the Von Mises yield criterion, which is
commonly used to predict material failure in elastic systems. It states that
if the second invariant of the solid’s deviatoric stress (the Von Mises
stress) is greater than a critical value (the yield strength) the material
will begin to deform non-elastically (Von Mises,, 1913). Although we do not
specify said critical value within our simulations, we can map the time-
evolution of Von Misses stresses within the coastal barrier as a result of a
wave absorption cycle (Figure 6.8). Our results illustrate the potential
utility of our simulation framework in predicting the location and time-of-
formation of stress induced defects within coastal barrier as a function of
wave characteristics, permeability, and barrier geometry.

####  6.2.2 Flow-Induced Surface Deformation

Surface deformation due to subsurface fluid flow is a common geological
phenomenon occurring in strongly coupled systems and has clear implications in
studies related to induced seismicity (Shapiro & Dinske,, 2009),
${\mathrm{CO}}_{2}$ injection in the subsurface (Morris et al.,, 2011), land
subsidence (Booker & Carter,, 1986), and the formation of dykes and volcanoes
(Abdelmalak et al.,, 2012; Mathieu et al.,, 2008). In order to properly model
these systems, it is necessary to be able to capture the time-evolution of
surface uplift, cracks, and hydraulic fractures, as well as the effects that
these features have on the overall flow field. Here, we use the terms
hydraulic fracture vs. crack to refer to solid failure at vs. away from the
injected fluid, respectively.

This illustrative case was inspired by the experiments reported by Abdelmalak
et al., (2012), where the authors injected a highly viscous fluid into a dry
silica powder in a Hele-Shaw cell in order to study the impact of hydraulic
fractures on surface deformation, e.g., during the creation of volcanic
structures. The system also bears some analogy to situations involving the
injection of fluids into subsurface reservoirs, e.g., during geologic CO2
sequestration (Rutqvist,, 2012). The base case of our simulations consists of
an impermeable rectangular container (50 by 30 cm, 500 by 300 cells) that is
open to the atmosphere, is partially filled with a dry porous medium
(${\phi}_{s}=0.6\ \pm 0.05$, ${\rho}_{s}=2650\ \mathrm{kg/m^{3}},\
k_{0}=5\times 10^{-11}\ \mathrm{m^{2}})$, and has an injection well at its
lower boundary that injects water at $q=6.5\ \mathrm{ml/s}$ (Fig. 6.9). To
account for irreversible solid deformation, the porous medium is modeled as a
plastic with a kinematic yield stress $\tau_{0}=0.22\ \mathrm{m^{2}/s^{2}}.$
The solid is represented as impermeable to the invading fluid through the use
of the Van Genuchten model with $m=0.05$ and $p_{c}=0$. Then, using this base
case as a standard, we individually varied each of the main parameters ($q,\
k_{0},\ \tau_{0},\ m,\ {\phi}_{s},\ {\mu}_{water})$ over several simulations
in order to model the resulting solid deformation processes: fracturing,
cracking, surface uplift, and subsidence (Figure 6.9).

Figure 6.9: Study of the effects of hydraulic fracturing and cracking on
surface deformation. (A-I) Representative cases showing the effects of
changing permeability $k_{0}$ (purple), solid yield stress $\tau_{0}$ (green),
injection rate $q$ (brown), and injected fluid viscosity $\mu$ (red) on
surface deformation. The blue and yellow subsections contain the results of
increasing or decreasing the controlling parameters, respectively. (J-L) Time
evolution of the fracturing base case. (M) Surface subsidence example. The
difference between the base case (E) and all other simulations is shown in
each case’s legend. Dotted white lines represent the surface height of the
initial solid fraction configuration. Note that the color scheme in all
simulations is the same as in Figure 6.4.

The resulting cases demonstrate that cracking (solid failure away from the
injected fluid) is strictly dependent on the number and orientation of
existing hydraulic fractures, as it only occurs when there is more than one
fracture branching off from the main injection point (Figure 6.9B, C, D, H,
and I). This is likely because in cases presenting a single vertical fracture
solid displacement is almost exclusively perpendicular to the fracturing
direction, leading to virtually no surface deformation or cracking (Figure
6.9A, E and M). Contrastingly, the creation of inclined fractures exerts
vertical forces on the solid, resulting in surface uplift and crack formation.
The above diagram strongly suggests that deformation is controlled by the
balance between viscous and structural forces: larger fractures occur within
softer solids with higher momentum transfer, and smaller fractures occur in
tougher solids with lower momentum transfer. As stated above, a comprehensive
examination of the parameters that control solid fracturing will be the focus
of Chapter 7.

In addition to the surface uplift presented above, subsurface subsidence is
observed in the simulated system in conditions where the porous solid is
rendered permeable to the invading fluid (i.e., $m\gg 0.05$). This phenomenon
is not primarily controlled by momentum transfer, but rather by a
gravitational effect whereby the displacement of air by water within the
porous medium around the advancing hydraulic fracture renders the solid-fluid
mixture heavier. Once it is heavy enough to overcome the plastic yield stress,
the solid sinks and compresses around the fluid source (Figure 6.9M).

###  6.3 Conclusions

Throughout this chapter and its predecessors, we have shown that our modeling
framework is flexible and readily applicable to a large variety systems: from
single-phase flow in static porous media, to elastic systems under
compression, to viscosity- or capillarity-dominated fracturing systems, all
the way up to multiscale wave propagation in poroelastic coastal barriers. We
invite the interested reader to tune, adapt, and expand the present
illustrative simulations, which are included in the accompanying CFD toolbox.

Lastly, we would like to note that this framework cannot be applied to every
multiscale system, at it is based on strong assumptions of lengths-scale
separation, sub-volume homogeneity, Eulerian descriptions, and relatively
small elastic deformations. These assumptions come with their own limitations.
A thorough discussion on the model’s constraints and future work can be found
in Chapter 9.

[75mm] The most exciting phrase to hear in science, the one that heralds the
most discoveries, is not ‘Eureka!’ but ‘That’s funny…’Isaac Asimov

##  Chapter 7 Capillary and Viscous Fracturing During Drainage in Porous Media

In this chapter, we use the Multiphase DBB framework to study the transition
from uniform fluid invasion to capillary and viscous fracturing during
drainage in porous media. To do so we examine multiphase flow in deformable
porous media in a broad range of flow, wettability, and solid rheology
conditions. We then demonstrate the existence of three distinct fracturing
regimes controlled by two non-dimensional numbers that quantify the balance of
viscous, capillary, and structural forces in a porous medium. We then use
these parameters to establish a first-of-its-kind phase diagram for material
failure caused by multiphase flow in poroplastic media. Lastly, we examine the
effects of compaction on said dimensional numbers and the system’s propensity
to fracture. This chapter is adapted from Carrillo & Bourg, 2021a .

###  7.1 Introduction

Multiphase flow in deformable porous media is a ubiquitous phenomenon in
natural and engineered systems that underlies key processes in water and
energy resource engineering and materials science, including membrane
filtration, soil wetting/drying, unconventional hydrocarbon recovery, and
geologic carbon sequestration (Bächer & Gekle,, 2019; Räss et al.,, 2018;
Towner,, 1987). A key obstacle to more accurate representations of this
phenomenon is our limited understanding of the transition from fluid invasion
to flow-induced fracturing, i.e., material failure caused by multiphase flow.
In large part, this limitation is caused by a lack of computational approaches
capable of representing multiphase flow in fractured deformable porous media.

Previous work on multiphase flow within static porous media is extensive and
includes detailed examinations of the influence of wettability, viscosity, and
flow rate on flow in unsaturated porous media at multiple scales. In
particular, existing studies have demonstrated how capillary forces give rise
to differences between drainage and imbibition (Lenormand,, 1986); how the
ratio of fluid viscosities controls the stability of the invading fluid front
(Måløy et al.,, 1985; Saffman & Taylor,, 1958; Stokes et al.,, 1986); and how
the magnitude of the capillary number delineates distinct flow regimes (Ferer
et al.,, 2004; Yortsos et al.,, 1997; Datta & Weitz,, 2013; Lu et al.,, 2020).
Each of the aforementioned controls is highly dependent on the system of
interest. This complicates efforts to develop general relative permeability
and capillary pressure models that apply to most systems of interest (Picchi &
Battiato,, 2018, 2019; Brooks & Corey,, 1964; van Genuchten,, 1980).

Flow of a single fluid phase through deformable porous media also has been
studied in depth. Numerical modeling studies are largely based on the work of
Biot and Terzaghi (Biot,, 1941; Terzaghi,, 1943) and have been used to
reproduce the behavior of arteries, boreholes, swelling clays, and gels (Auton
& MacMinn,, 2017; Bertrand et al.,, 2016; Carrillo & Bourg,, 2019; MacMinn et
al.,, 2015). In the last decade, fundamental studies have generated detailed
information on the dynamics that arise from fluid-solid couplings beyond the
ideal poroelastic regime, including fracturing, granular fingering, and
frictional fingering (Campbell et al.,, 2017; Sandnes et al.,, 2011; Zhang et
al.,, 2013). In particular, these studies have shown that the main parameters
controlling the deformation of a porous solid by single phase flow are the
material softness and the magnitude of the fluid-solid momentum transfer
(Sandnes et al.,, 2011).

The study of multiphase flow in a deformable porous medium is inherently more
complex than the problems outlined above, as it requires simultaneous
consideration of capillarity, wetting dynamics, fluid rheology, and solid
deformation. Deformation modes associated with material failure (i.e.,
multiphase fracturing) are particularly challenging as they require
simultaneous representation of multiphase flow in fractures and in the
surrounding porous matrix. The existing detailed examinations of this
phenomenon have focused exclusively on granular systems. Notably, Holtzman &
Juanes (Holtzman & Juanes,, 2010; Holtzman et al.,, 2012) used experiments and
discrete element models to demonstrate that the transitions between capillary
fingering, viscous fingering, and fracturing during multiphase flow in
granular media reflect two non-dimensional numbers: a fracturing number (ratio
of fluid driving force to solid cohesive force) and a modified capillary
number (the ratio between viscous and capillary pressure drops). Other
discrete element approaches have shown that fracturing is highly dependent on
the invading fluid’s capillary entry pressure (Jain & Juanes,, 2009; Meng et
al.,, 2020). However, it is not clear how these conclusions translate to
continuous non-granular systems.

To the best of our knowledge, no experimental or numerical investigation has
simultaneously explored the effects of flow rate, wettability, and
deformability during multiphase flow in deformable porous media at the
continuum scale and identified the controlling parameters that relate all
three properties within a single phase diagram. Here, we use simulations
carried out with our new Multiphase Darcy-Brinkman-Biot (DBB) framework
(Carrillo & Bourg, 2021b, ) to fill this knowledge gap and identify non-
dimensional parameters that govern viscously-stable fluid drainage and
fracturing in deformable porous media. We also find that the fracturing
dynamics predicted by our continuum-scale framework is consistent with those
observed or predicted for granular systems. In other words, in systems with a
large length scale separation between pores and fractures, volume-averaged
properties are sufficient to capture the onset and propagation of fractures at
the continuum scale.

###  7.2 Numerical Simulations

####  7.2.1 Crossover from Imbibition to Fracturing in a Hele-Shaw Cell

In addition to the derivation and extensive quantitative validation of the
modeling framework, Chapter 6 included a qualitative validation of the ability
of the Multiphase DBB model to predict the transition from invasion to
fracturing during multiphase flow. Briefly, this validation replicated
experiments by Huang et al., 2012b  involving the injection of aqueous
glycerin into dry sand at incremental flow rates within a 30 by 30 by 2.5 cm
Hele-Shaw cell. As shown in Fig. 7.1, these experiments are inherently
multiphysics as fluid flow is governed by Stokes flow in the fracture
(aperture $\sim$cm) and by multiphase Biot Theory in the porous sand (pore
width $\sim 100\mu$m).

Figure 7.1:  Continuous transition from fluid imbibition to fracturing in a
Hele-Shaw cell. Experimental images (A, B, C) were taken from Huang et al.,
2012b  and numerically replicated using equivalent conditions (D, E, F). Black
lines represent the advancing saturation front. Additional cases can be found
in Chapter 6.

As discussed in Chapter 6, the similarities between our model and the
experimental results are evident: as the viscous forces imposed on the solid
increase, so does the system’s propensity to exhibit fracturing as the primary
flow mechanism (as opposed to imbibition). Minor microstructural differences
between our simulations and the experiments reflect the manner in which the
implemented continuum-scale rheology model approximates the solid’s granular
nature. It is clear, however, that both systems are controlled by the balance
between viscous forces and solid rheology at the scale of interest (Carrillo &
Bourg, 2021b, ). As such, these experiments present an ideal starting point
for our investigation.

####  7.2.2 Creation of Fracturing Phase Diagrams

Here, we use the same simulation methodology developed in (Carrillo & Bourg,
2021b, ) and illustrated in Figure 7.1 to identify the general non-dimensional
parameters that control the observed transitional behavior between invasion
and fracturing in a plastic porous medium. To do so, we systematically vary
the solid’s porosity ($\phi_{f}$ = 0.4 to 0.8), density-normalized plastic
yield stress ($\tau_{yield}$ = $1.5$ to ${24\ \mathrm{m^{2}/s^{2}}}$),
capillary entry pressure ($p_{c,0}$ = $100$ to $50,000\ \mathrm{Pa}$), and
permeability ($k$ = $1\times 10^{-13}$ and $5\times 10^{-9}\ \mathrm{m^{2}}$)
as well as the invading fluid’s viscosity ($\mu_{n}$ = $0.5$ to $50\
\mathrm{cP}$) and injection rate ($\boldsymbol{U}_{f}$ = $1\times 10^{-4}$ to
$8\times 10^{-2}\ \mathrm{m/s}$). As in our previous work, the solid’s
porosity was initialized as a normally-distributed field, the deformable solid
was modeled as a Hershel-Bulkley-Quemada plastic (Spearman,, 2017; Quemada,,
1977), the porosity-dependence of permeability was modeled through the Kozeny-
Carman relation, and relative permeabilities where calculated through the van-
Genuchten model (van Genuchten,, 1980). Further details regarding the base
numerical implementation of this model can be found in Carrillo & Bourg, 2021b
and the accompanying code (Carrillo, 2020a, ). The only major differences
relative to our previous simulations are that we now include capillary effects
and represent viscously-stable drainage as opposed to imbibition (i.e., the
injected glycerin is now non-wetting to the porous medium). A representative
sample of the more than 400 resulting simulations is presented in the phase
diagrams shown in Fig. 7.2.

Figure 7.2: Phase diagrams describing the effects of varying permeability,
plastic yield stress, fluid injection rate, and capillary entry pressure on
the transition from fluid drainage to fracturing. All cases are at
${\phi}_{s}=0.60\pm 0.05$ and ${\mu}_{n}=5\ \mathrm{cP}$. The remaining
parameters are case-specific and can be found in each figure’s upper legend.
The areas separated by thin blue lines highlight and label the four
deformation regimes described in Section 7.3. The vertical red lines represent
where these diagrams intersect in 3-dimensional space. The color scheme is the
same as in Fig. 7.1.

Overall, the results make intuitive sense. Figure 7.2A shows that, ceteris-
paribus, less permeable solids are more prone to fracturing. This is due to
the fact that, given a constant flow rate, lower permeability solids
experience greater drag forces. Our results also show that solids with lower
plastic yield stresses fracture more readily, as their solid structure cannot
withstand the effects of relatively large viscous or capillary forces. The
y-axis behavior of Fig. 7.2B further shows that systems with higher entry
pressures are more likely to fracture, i.e., the capillary stresses are more
likely to overwhelm the solid’s yield stress, in agreement with grain scale
simulations (Jain & Juanes,, 2009). Finally, Fig. 7.2B also shows that higher
injection rates lead to more fracturing, as these increase viscous drag on the
solid structure.

###  7.3 Characterization of Fracturing Mechanisms

The deformation regimes observed in the previous experiments can be delineated
by defining two simple non-dimensional parameters that quantify the balance
between viscous pressure drop, solid softness, and capillary entry pressure.

$N_{vF}=\frac{\Delta
p}{{\tau}_{yield}{\rho}_{s}}=\frac{{\mu}U{r}_{in}}{k{\tau}_{yield}{\rho}_{s}}{\mathrm{ln}\left(\frac{r_{out}}{r_{in}}\right)}\
\ \ \ $ (7.1)

$N_{cF}=\frac{p_{c,0}}{{\tau}_{yield}{\rho}_{s}}=\frac{2\gamma}{r_{pore}{\tau}_{yield}{\rho}_{s}}$
(7.2)

Here, the viscous fracturing number ($N_{vF}$) represents the ratio between
the viscous pressure drop and the solid’s structural forces. It embodies the
question: Does fluid flow generate sufficient friction to induce fracturing?
As shown in Fig. 7.3, the answer is no if $N_{vF}<1$ and yes if $N_{vF}>1.$
This number is the continuum scale analog to the fracturing number presented
by Holtzman et al., (2012) for granular solids. It also explains the
experimental finding by Zhou et al., 2010a  that fracture initiation is only a
function of the resulting fluid pressure drop, irrespective of the injection
rate or fluid viscosity used to create it. Furthermore, it illustrates why
increasing the injection rate and decreasing the permeability have similar
effects in Fig. 7.2.

Complementarily, the capillary fracturing number ($N_{cF}$) represents the
ratio between the capillary entry pressure and the solid’s structural forces;
it embodies the question: Does multiphase flow generate sufficient capillary
stresses to fracture the solid? Figure 7.3 shows that when $N_{cF}<1$ drainage
is the preferential flow mechanism and when $N_{cF}>1$ fracturing becomes the
dominant phenomenon.

Figure 7.3:  Fluid invasion and fracturing in plastic porous media as a
function of the viscous fracturing number $N_{vF}$ and the capillary
fracturing number $N_{cF}$. Green triangles denote uniform invasion, red
diamonds denote the transitional fracturing regimes, and blue circles denote
non-invasive fracturing. The four image insets are representative samples of
each fracturing regime.

This analysis yields the rudimentary conclusion that fracturing should occur
if either of the fracturing numbers is greater than unity, as confirmed by our
simulations. However, our simulations further demonstrate the existence of
three distinct fracturing regimes (Figs. 7.2-7.3). The first regime, referred
here as non-invasive fracturing ($N_{vF}>1$ and $N_{cF}>1$) is characterized
by fracturing of the porous solid with minimal fluid invasion, where fractures
precede any invasion front. In the second regime, referred to here as the
viscous fracturing transition ($N_{vF}>1$ and $N_{cF}<1$), only the viscous
stresses are sufficiently large to fracture the solid. This leads to the
formation of relatively wide fractures enveloped and preceded by a non-uniform
invasion front. Finally, in the third regime, referred to here as the
capillary fracturing transition ($N_{vF}<1$ and $N_{cF}>1$), only the
capillary stresses are sufficiently large to fracture the solid. Given a
constant injection rate, this leads to the formation of fractures preceded by
an invasion front, as in the viscous fracturing transition regime, but with a
more uniform saturation front (due to lower viscous stresses) and less solid
compaction (hence narrower fractures). We note that the crossover between each
of the four regimes is continuous, meaning that systems with $N_{vF}$ or
$N_{cF}\sim 1$ can share elements of neighboring regimes.

Although $N_{vF}$ and $N_{cF}$ are fairly intuitive numbers, their impacts on
fracture propagation mechanisms are not. For this reason, we also studied the
dynamics of fracture nucleation and growth and the evolution of the solid’s
strain for all three fracturing regimes. As seen in Fig. 7.4, fracturing in
the two transition zones is characterized by the initial formation of non-
flow-bearing failure zones (hereafter referred to as cracks), which function
as nucleation sites for propagating flow-bearing fractures. These cracks are
formed by the simultaneous movement of large contiguous sections of the porous
medium in different directions, a process induced by uniform fluid invasion
into the porous medium. However, the similarities between both transition
zones end here. In the viscous fracturing transition regime, fractures quickly
become the dominant deformation mechanism, localizing the majority of the
stresses and solid compaction around the advancing fracture tip. Conversely,
in the capillary fracturing transition regime, fluid-invasion continues to
serve as the main flow mechanism and source of deformation, where fractures
and cracks are slowly expanded due to the more evenly-distributed capillarity-
induced stresses localized at the advancing invasion front. Finally, non-
invasive fracturing follows a different process, where there is little-to-no
crack formation and fracture propagation is the main source of deformation and
flow. Here, the co-advancing fracture and saturation fronts uniformly compress
the solid around and in front of them until this deformation reaches the outer
boundary of the simulated system (see the “jet” like-structures at fracture
tips in Fig. 7.4C.)

Figure 7.4:  Dynamic fracture formation mechanisms. Each row represents the
time-dependant fracture formation process for each fracturing type, where time
advances from left to right. Here, the red-blue color scheme represents the
log-normalized strain-rate magnitude specific to each simulated case,
fractures are shown in white and the advancing fluid-fluid interface is shown
as a thin black line.

###  7.4 Influence of Localized and Uniform Deformation

So far we have explored how independently changing $k$, $p_{c}$, and
$\tau_{yield}$ (among others) can affect the fracturing of plastic materials.
However, our results also have implications for situations in which these
variables are all varied simultaneously, such as during the compaction of
soils, sediments, or viscoplastic sedimentary rocks (i.e. mudstones or clay-
shales). In such situations, with increasing compaction, $k^{-1}$, $p_{c}$,
and $\tau_{yield}$ should all increase, although at different rates. As such,
we now study the effects of local and uniform deformation on the outlined
fracturing regimes.

####  7.4.1 Localized Deformation

The simulations presented above were carried out using the simplifying
assumption that $p_{c}$ is invariant with $\phi_{f}$ (whereas $k$ and
$\tau_{yield}$ are not). To evaluate the impact of this simplification on the
results shown in Figs. 7.2-7.3, we carried out additional simulations for all
four regimes with a deformation-dependent capillary entry pressure based on a
simplified form of the Leverett J-function where
$p_{c,0}=p_{c,0}^{*}(\phi_{s}/\phi_{s}^{avg})^{n}$, $p_{c,0}^{*}$ is the
capillary pressure at $\phi_{s}=\phi_{s}^{avg}$, and $n>0$ is a sensitivity
parameter (Leverett,, 1941; Li & Benson,, 2015). The results show that non-
zero values of $n$ promote the creation of finger-like instabilities and the
nucleation of cracks at the fluid invasion front, particularly in the
capillary fracturing transition regime. Simulation predictions with different
n values are shown in Fig. 7.5.

Despite the additional complexity of the resulting fluid invasion and
fracturing patterns, results with $n>0$ conform to the overall phase diagram
presented in Fig. 7.3. The results at $n=0$ are therefore highlighted in the
previous sections due to the greater simplicity of their fluid and solid
distribution patterns.

Figure 7.5:  Influence of the $\phi_{f}$-dependence of $p_{c}$ on fluid
invasion and fracturing patterns for all three fracturing regimes. Here, $n$
represents the sensitivity parameter in the Leverett J-function analogue
presented above. The color scheme is the same as in Fig. 7.1.

####  7.4.2 Uniform Deformation

Having verified that the applicability of the fracturing numbers holds for
systems were $k$, $\tau_{yield}$, and $p_{c}$ all vary with $\phi_{f}$, we now
examine the effects of uniform compaction on said numbers. A direct analysis
using the widely-used porosity-parameter relationships implemented above (the
Kozeny-Carman relation for $k$, Leverett J-Function for $p_{c}$, and Quemada
model for $\tau_{yield}$ (Leverett,, 1941; Quemada,, 1977; Spearman,, 2017))
yields the following fracturing number - porosity dependence:

$N_{vF}\propto\frac{(1-\phi_{f})^{2-D}(1-\phi_{f,min}/\phi_{f})}{\phi_{f}^{2}}$
(7.3)

$N_{cF}\propto(1-\phi_{f})^{2-D}(1-\phi_{f,min}/\phi_{f})$ (7.4)

where $D$ is a rheological parameter based on the solid’s fractal dimension
(common values range for 1.7-2.9 for different clayey sediments (Spearman,,
2017)) and $\phi_{f,min}$ is the maximum possible degree of compaction.
Through these relations, we can see that uniform compaction (or expansion) has
a highly non-linear effect on fracturing. Equations 7.3-7.4 indicate that
whereas $N_{cF}$ tends to consistently decrease with increasing compaction,
$N_{vF}$ is considerably more susceptible to changes in $\phi_{f}$ and
exhibits multiple changes in the sign of its first derivative when $D>2$, non-
intuitively suggesting that fracturing can be either induced or suppressed
through uniform compression. Plots of $N_{vF}$ and $N_{cF}$ as a function of
solid fraction are shown in Fig. 7.6.

Figure 7.6:  Fracturing numbers’ dependence on Porosity. A) For $N_{vF}$ B)
For $N_{cF}$. Colored arrows represent the overall trends in the fracturing
number behaviours when changing the fractal parameter $D$ and the minimum
porosity parameter $\phi_{f,min}$ in Eqns. 7.3-7.4. Note the appearance of
“valleys” at high values of D in sub-figure A.

###  7.5 Conclusions

In this chapter, we used the Multiphase DBB modeling framework to create a
phase diagram that identifies two non-dimensional parameters that categorize
the crossover between viscously-stable fluid drainage and fracturing as a
function of wettability, solid deformability, and hydrodynamics. To the best
of our knowledge, our results are the first to relate all three of these
properties to characterize multiphase flow in viscoplastic porous media. As
expected intuitively, we observe that fracturing occurs if the viscous and/or
capillary stresses are sufficient to overcome the solid’s structural forces.
Thus, when it comes to systems with multiple fluids, it is necessary to
consider the effects of surface tension, wettability, and pore size on the
fluids’ propensity to fracture or invade the permeable solid. Furthermore, we
found that the two non-dimensional fracturing numbers described above
delineate the existence of three fracturing regimes with distinct fracture
propagation mechanisms. Lastly, we examined how uniform compression or
expansion affect said non-dimensional numbers and a system’s propensity to
fracture . [75mm] Science without religion is lame, religion without science
is blind Albert Einstein

##  Chapter 8 Simulation and Prediction of Stochastic Clogging Processes in
Porous Media

In this chapter, we take a small detour from modeling multiphase flow in
deformable porous media to focus on a conceptually simple, yet physically
complex phenomenon: Clogging mechanics in porous media. This work was done as
part of Princeton’s Center for Statistics and Machine Learning (CSML) graduate
certificate program. Although this chapter has not been published yet, we hope
to do so in the near future.

###  8.1 Introduction

The erosion, transport, and eventual deposition of fluid-suspended particles
is ubiquitous within both natural and engineered porous systems. These
processes control the evolution of sedimentary formations, the distribution of
contaminants in the environment, and the clogging of porous media, pipes, and
arteries (Molnar et al.,, 2015; Phenrat et al.,, 2009; Robert De Saint Vincent
et al.,, 2016; Zuriguel et al.,, 2014). Although material transport in porous
media is relatively well understood, the topics of particle deposition,
accumulation, and clogging are still relatively new in the field.

The difficulty of predicting clogging in porous media arises from the
complexity, heterogeneity, and stochastic nature of the relevant systems.
Clogging in porous media is notoriously hard to probe and characterize, as it
is necessary to account for each system’s 3D geometry, grain size
distribution, pore size distribution, rock/fluid chemistry, particle size, and
even surface charge (Ding et al.,, 2015; Sahimi & Imdakm,, 1991; Liu et al.,,
1995; Gerber et al.,, 2018; Mirabolghasemi et al.,, 2015; Pham &
Papavassiliou,, 2017). Added to these spatial challenges is the fact that
clogging is also a temporal process. Clogging mechanisms often don’t reach a
steady state: Clogs can form, redirect fluid fields, change pressure
gradients, break, and form again downstream (Bizmark et al.,, 2020).

Experimentalists have taken two separate routes in addressing these
challenges. The first approach attempts to characterize particle deposition
and clogging by reducing the complexity of the studied systems through the
creation of simplified micro-models (Wyss et al.,, 2006; Agbangla et al.,,
2012; Auset & Keller,, 2006; Gerber et al.,, 2019). This allows a measure of
control of a handful of relevant variables (fluid flow rate, particle size,
pore size, flow time) while significantly constraining others (flow path
complexity and material heterogeneity). These studies have yielded noteworthy
results. In particular, they have shown that clogging scales as a function of
the pore to particle size ratio and that, under certain conditions, it is
independent of particle injection rate and system porosity (Wyss et al.,,
2006). However, it is not clear if these conclusions still hold for more
complex natural systems.

The second approach relies on using advanced imaging techniques to probe real
porous systems (or close approximations to them). These often rely on computed
X-ray micro-tomography (XCT) and/or confocal microscopy. The former has taken
a key role in characterizing particle aggregate formation in natural rock
samples and identifying their dependence on fluid chemistry and rock geometry
(Liu et al.,, 1995; Li et al., 2006b, ; Li et al., 2006a, ; Chen et al.,,
2009). However, although highly informative, this technique suffers from not
having large enough fields-of-view or high enough temporal resolutions.
Confocal microscopy takes advantage of the fact that it is possible to match
the refractive index of an artificial porous medium with the index of its
permeating fluid phase in order to “see” through the whole solid-fluid
mixture. This technique has allowed scientists to characterize particle
deposition and clogging as a function of time while having a high level of
control of the relevant flow variables and material characteristics (Bizmark
et al.,, 2020; Mays et al.,, 2011).

However, all experimental studies suffer from the same practical problem: they
cannot probe the full parametric phase-space required to properly generalize
their results to naturally-occurring porous media due to the difficulty of
creating and running an experiment with more than 3-5 independent variables.
The necessity to investigate a larger range of systems is evidenced by the
fact that experimental studies often reach conflicting conclusions. For
example, some studies maintain that particle transport suppresses fluid flow
by reducing the permeability of porous media (Civan,, 2010; Liu et al.,, 1995;
Wiesner et al.,, 1996), while others state that particles can actually help
enhance flow (Weber et al.,, 2009; Kersting et al.,, 1999; Ryan & Elimelech,,
1996; Schneider et al.,, 2021). Even so, assuming that we could obtain a
large-enough data set, there is no guarantee that we could derive general
relationships between variables: some conclusions may only hold at low flow
rates and/or low grain heterogeneity, while others may only hold at low fluid
viscosities and/or high particle concentrations. One way to address this
problem is through the application of Machine Learning (ML) algorithms to
analyze such large data sets. These approaches have been shown to learn and
predict complex relations in a wide variety of fields, from predicting fluid
turbulence in jet engines (Sirignano et al.,, 2020) to enabling image
recognition in autonomous vehicles (Kocić et al.,, 2019). These models,
however, require highly extensive data sets to become effective predictors.

Fortunately, due to the rise of large-scale parallelized cloud computing, it
is now possible to run thousands of particle transport simulations in order to
systematically study 10-20 dimensional parametric phase spaces. However, to
the best of our knowledge, there have not been any studies that have used
numerical models for this purpose. Current popular approaches rely on
continuum-level (i.e. volume averaging) approximations such as filtration
theory to study particle deposition (Molnar et al.,, 2015; Messina et al.,,
2016; Boccardo et al.,, 2018) and erosion (Hilpert & Johnson,, 2018; Jäger et
al.,, 2017). Contrastingly, coupled CFD-DEM approaches (Zhao & Shan,, 2013;
Mirabolghasemi et al.,, 2015; Zuriguel et al.,, 2014; Natsui et al.,, 2012)
attempt to model particle deposition directly through careful consideration of
coupled particle-fluid mechanics and direct numerical simulations (DNS).
However, due to their complexity and relatively-high computational cost, very
few studies have attempted to model particle mechanics through these means.

In this chapter, we leverage Princeton’s University’s large computational
resources in order to explore the underlying principles that control clogging
mechanics. In particular, we explore the feasibility of using ML approaches to
create a computational tool that can predict clogging a-priori and that can be
used to improve and streamline the design process of engineered porous
systems. To do so, we ran 2000 CFD-DEM simulations in randomly-generated
porous geometries while systematically varying 12 of the system’s design
parameters (pore size, pore size heterogeneity, grain size heterogeneity,
particle size, particle-particle attraction, particle-wall attraction,
particle flux, porosity, system size, fluid velocity, and fluid viscosity).
The resulting cases where labeled and analyzed in order to train and evaluate
several different ML classifiers and regressors, which we then used to predict
clogging in equivalent systems. Then, through model optimization, we
identified which training features were most indicative of clogging in
heterogeneous porous media. To the best of our knowledge, this is the first
time ML approaches have been successfully used to reliably predict and control
clogging processes.

###  8.2 Methods

####  8.2.1 CFD-DEM

Our numerical simulations where performed on the CFDEM® computational modeling
framework, which couples Computational Fluid Dynamics (CFD) simulations
performed in OpenFOAM® with Discrete Element Models (DEM) performed in
LIGGGHTS®. Both C++ libraries are free-to-use, open-source, and parallelizable
platforms. OpenFOAM® uses the Finite Volume Method to discretize and solve
partial differential equations to describe fluid motion in complex 3-D grids.
In turn, LIGGGHTS® uses Lagragian approaches and discrete particle models to
describe the motion and interaction of a large number of individual particles
within a grid-free environment. The object-oriented structure of both codes
and multitude of supporting libraries allows the user to easily customize each
simulation’s setup with different numerical discretization schemes, time-
stepping procedures, matrix-solution algorithms, and supporting physical
models.

Figure 8.1:  Conceptual representation of the models we used for clogging
simulation and prediction.

CFDEM® successfully couples both models through an iterative algorithm that
super-imposes their physical domains. Here, an additional particle-drag term
is added to the fluid’s momentum conservation equation to account for particle
movement, and an additional fluid-drag term is added to each particle’s
equation of motion to account for fluid-induced motion (i.e. flow and
buoyancy). Each model can probe its counterpart at any given moment in time,
at which point the CFDEM® algorithm iterates over both solver’s solutions
until the desired convergence metric is reached. Just as in its two base
models, CFDEM® counts with several customizable options regarding the coupling
frequency, numerical schemes, and coupling-physical models. Our specific
numerical implementation will be discussed in Section 8.4.2.

####  8.2.2 Machine Learning

The statistical and ML analysis performed in this study was carried out
through Scikit-learn, a free-to-use, python-based library. Sckikit-learn was
selected for its ability to allow users to easily train, optimize, and
interpret a large number of supervised and unsupervised machine learning
algorithms. More complex and customizable ML libraries such as TensorFlow and
pyTorch would have been equally adequate. In our case, since we worked
exclusively with labeled data, we only focused on training and testing a
representative sample of ML classifiers and regressors (see Tables 8.1 and
8.3).

###  8.3 Workflow

Having established the numerical simulation and ML basis of our investigation,
we now describe our project workflow: First, we created randomized porous
systems. In this step, we generated the geometry and computational mesh on
which we performed our numerical simulations. A combination of randomization
and systematic-variation of the relevant geometric variables (e.g. pore size)
was necessary for the eventual generalization of our results. Second, we ran
CFDEM simulations. During this step, we ran thousands of clogging processes
over a large parameter space by varying the specified fluid and particle
variables (e.g. fluid flow rate and particle size). Third, we identified
clogged systems by developing a standardized metric for labeling and
characterizing our completed numerical simulations. Fourth, we identified and
engineered features for training. In this step, we converted the varied system
parameters into a set of standardized non-dimensionalized variables. This was
done for two reasons: it reduces redundancy and creates more predictive
features (e.g. the ratio between particle size and pore size is a better
predictor for clogging than either variable by itself) while also allowing us
to apply the resulting ML model to any potential system that can be similarly
non-dimensionalized. Fifth, we trained the ML models. Here, we trained,
optimized, and identified the best ML classifier and regressor for predicting
our simulated clogging processes. Finally, we interpreted and analysed the
models. In other words we identified the most predictive features for particle
clogging in porous media. In addition, we also evaluated the efficacy of
different models, their shortcomings, and their potential improvements.

The rest of the chapter will focus on developing and discussing each of these
individual steps.

###  8.4 Experimental Setup

####  8.4.1 Creating Randomized Porous Systems

The randomized porous geometries used in this study were generated through a
custom Matlab® script which allowed us to specify the following criteria:
average grain size, grain size standard deviation, average pore size (i.e.
shortest distance between two given grains), pore size standard deviation, and
porosity (see Figure 8.2). Note that, for simplicity, all grains where assumed
to be cylindrical in shape and the domain’s length and width where kept
constant at 50 cm. Given a fixed domain space, the porosity and pore size
distribution become coupled variables, which is why only one of these could be
specified for a given geometry (the other one was measured at the end of
geometry generation). Furthermore, to allow for a measurable particle flux
into the geometry, a funnel-type feature was added to the top of the generated
geometry during the numerical simulation setup (see Figure 8.2C). The result
is an algorithm that can generate well-characterized, yet randomized porous
systems on which we can perform our numerical simulations.

Figure 8.2:  A,B) Examples of randomly-generated porous media configurations
with different grain sizes and porosities. C) Example of the numerical mesh
and simulation setup. Here, the areas within the cylinders (i.e. grains) are
removed during meshing to making them inaccessible to fluids and particles.
Also note the direction of flow and the fact that the only inlets and outlets
are at the upper and lower boundaries of the geometry, respectively.

####  8.4.2 CFDEM Simulations

We now present our numerical simulation setup. As discussed earlier, the main
advantage that numerical simulations have over conventional experiments is
their ability to probe large parameter spaces with ease. As such, we
systematically varied the following 12 variables over the following ranges in
2000 CFDEM simulations:

  1. 1.

Porosity ($\phi$): $0.26$ to $0.98$

  2. 2.

Average pore size ($D$): ${10}^{-4}$ to ${10}^{-2}$ m

  3. 3.

Pore size standard deviation ($D_{std}$): $5\times{10}^{-4}$ to
$5\times{10}^{-3}$ m

  4. 4.

Average grain size ($G$): ${10}^{-3}$ to $5\times{10}^{-2}$ m

  5. 5.

Grain size standard deviation ($G_{std}$): ${10}^{-3}$ to ${10}^{-2}$ m

  6. 6.

Geometric thickness ($T$): ${10}^{-4}$ to ${10}^{-2}$ m

  7. 7.

Particle diameter ($d$): ${10}^{-4}$ to ${10}^{-2}$ m

  8. 8.

Particle flux ($F$): ${10}^{2}$ to $1.2\times{10}^{4}$ particles/s

  9. 9.

Particle-particle attraction ($PP$): $0$ to $5\times{10}^{6}$
J/$\mathrm{m^{3}}$

  10. 10.

Particle-wall attraction ($PW$): $0$ to $5\times{10}^{6}$ J/$\mathrm{m^{3}}$

  11. 11.

Fluid velocity ($U$): $0$ to $1.5$ m/s

  12. 12.

Fluid kinematic viscosity ($\nu$): ${10}^{-2}$ to ${10}^{-8}$
$\mathrm{m^{2}}$/s.

All other parameters were kept constant. In particular, fluid density was 1000
$\mathrm{kg}/\mathrm{m^{3}}$, particle density was 1200
$\mathrm{kg}/\mathrm{m^{3}}$, and gravity ($g=9.8m/s^{2}$) was set to be
constant in the direction of flow. Particle-Particle/Wall attraction was
modeled through the SJKR model, Particle-Particle/Wall collisions where
captured through the Granular Hertz model with a Poisson Ratio of 0.22, and
the fluid-particle drag coupling was calculated through the DiFelice Drag
model every 0.001 seconds (Barthel,, 2008). These models imply the following:
A) the attraction force is normal to the inter-particle contact area and is
only activated once particles are in direct contact with each other. B) Inter-
particle collisions are modeled by defining both a normal force (spring +
damping forces) and a tangential force (shear + damping forces). This allows
particles to both bounce and roll around obstacles and other particles. C)
Fluid-particle interactions are governed by a drag force that is proportional
to the relative velocity between a given particle and the fluid, the solid
volume fraction in the specified control volume, and the square of the
Reynolds number (Zhou et al., 2010b, ).

Figure 8.3:  Examples of different clogging simulations. Note the wide range
of particle sizes (A and D), porosities (B and C), pore size distributions (C
and D), and grain size distributions (B and C).

In order to avoid modeling trivial clogging cases dictated by size exclusion,
we exclusively considered systems where the particles’ diameter was less than
the pore size and geometric thickness (See Figure 8.3). Each simulation was
set to run for 5 hydraulic residence times (0.5 m / fluid velocity), which
took on average about 5.5 hours of computing time on a 28-core Broadwell Xeon
node.

###  8.5 Data Analysis and Model Training

####  8.5.1 Identifying Clogged Systems

Accurate characterisation of our final configurations was crucial for the
training and testing of our ML algorithms. This process was complicated by the
fact that clogging in a heterogeneous porous medium is not necessarily a
discrete state. As shown Figure 8.4, clogging might only occur at certain pore
throats, might not occur at all, or might occur throughout the porous medium
Bizmark et al., (2020). As such, in order to properly label our clogging
geometries for ML classification (or obtain a continuous prediction variable
for ML regression) we had to define an objective metric that could
characterize the degree of clogging in each system. We called this variable
the “Clogging Number” (CN).

Figure 8.4:  Examples of different clogging levels. Each panel in this figure
shows the final state of four separate clogging simulations with identical
geometries and flow conditions. Their only difference is the size of the
simulated particles (shown here as the particle size to pore size ratio, d/D,
which increases from left to right). A) Free geometry. B, C) Semi-clogged
geometries D) Fully-clogged geometry.

The Clogging Number (Eqn. 8.1) is the product of two independent factors
obtained at the end of each simulation : 1) The average relative distance
between particles (DP). 2) The symmetry of the particle’s final velocity
distribution (SYMM). The significance of the first factor is fairly intuitive;
systems with particles that end up closer together tend to have a higher
degree of clogging. The second factor describes the relative behaviour of all
the particles throughout the porous medium, where relatively high symmetry
indicates homogeneous particle behaviour (i.e. the system is either fully
clogged or fully unclogged) and low symmetry indicates heterogeneous behaviour
(i.e. some particles are moving and some are stationary).

$\mathrm{CN}=\mathrm{DP}*\mathrm{SYMM}$ (8.1)

Figure 8.5:  Clogging sample distribution as a function of the average final
distance between particles (DP) and symmetry of the particles’ final velocity
distribution. Here, red and blue colors represent clogged and unclogged
samples, respectively. Note that point size is directly proportional to DP and
that all cases where plotted from lowest DP to highest DP values for
visualization purposes. A) Result of using CN ¡1 to separate 300 manually-
labelled fully clogged and unclogged samples. B) Result of using CN ¡1 to
label all data samples. The color contours in both graphs represent the phase-
space prediction of the optimized MLP classifier trained with their respective
samples and discussed in Section 8.5.3.

Therefore, multiplying both factors together into the CN allows us to
characterize each simulation’s degree of clogging, where DP tells us if the
system tends to be clogged or unclogged, and SYMM tells us the degree to which
they are one or the other. After testing several different averaging
procedures and symmetry measures, the following two metrics yielded the best
performance when used to create a CN that could separate fully clogged from
fully unclogged systems (see Figure 8.5A). DP was calculated by averaging the
average distance between each particle and its closest 50 neighbors, dividing
by the particle diameter, and subtracting 1 from said value (to avoid double
counting). In turn, SYMM was determined by calculating the ratio between the
particles’ median velocity and their mean velocity. This ratio describes the
symmetry of the distribution by using the fact that the median of a sample is
independent of changes in its extreme values, while the mean of a sample is.
Therefore, whenever a group of particles form a blockage in an otherwise
unclogged geometry, the median of the velocity will change at a much slower
rate than the mean velocity and their ratio will decrease. Conversely, if we
have a fully clogged or fully unclogged system, the velocity distribution will
have a higher degree of symmetry, and this ratio will be close to one.

Figure 8.5 shows that CN now provides us with a suitable way of characterizing
the final state our clogging simulations. It also gives us a straightforward
way to form discrete labels for classification out of continuous variables,
where all samples with CN ¡ $1$ will be labeled as “Clogged” ($1$) and all
others as “Unclogged” ($0$). This threshold was chosen by identifying the best
linear classifier that could separate a sample of 300 manually-labeled fully
clogged and fully unclogged systems (without considering any semi-clogged
samples). Applying this threshold to the complete data set yields 1012 samples
classified as Clogged and 988 classified as Unclogged.

####  8.5.2 Feature Standardization and Non-dimensionalization

In order to generalize our simulations and improve our models’ prediction
power, we now non-dimensionalize and/or re-scale our training features. The
objective is to create a standardized set of features that can be readily
quantified from any system involving the flow of solid bodies through a set of
static obstacles (such as cars through a road, sheep through a gate,
contaminants through soil, ext…). Nondimensionalization also allows us to add
information to the model by explicitly dictating important relationships
between variables. The following features are all system characteristics than
can be measured a-priori, meaning that, if trained correctly, a ML model will
be able to use these features to predict clogging without actually having to
run any simulations or experiments.

  1. 1.

Particle - pore size ratio ($d/D$): $d/D$

  2. 2.

Particle - geometric thickness ratio ($d/T$): $d/T$

  3. 3.

Particle - grain size ratio ($d/G$): $d/G$

  4. 4.

Geometric thickness - pore size ratio ($T/D$): $T/D$

  5. 5.

Grain size - pore size ratio ($G_{nd}$): $G/D$

  6. 6.

Non-dimensional grain size standard deviation ($G_{std,nd}$): $G_{std}/D$

  7. 7.

Non-dimensional pore size standard deviation ($D_{std,nd}$): $D_{std}/D$

  8. 8.

Non-dimensional particle-wall attraction ($PW_{nd}$):
$\frac{PW\times\pi\times(D/2)^{2}}{(\mathrm{ParticleMass}\times g)}$

  9. 9.

Non-dimensional particle-particle attraction ($PP_{nd}$):
$\frac{PP\times\pi\times(D/2)^{2}}{(\mathrm{ParticleMass}\times\textbf{g})}$

  10. 10.

Standardized particle flux ($F_{nd}$): $(F\times D)/\phi$

  11. 11.

Standardized Fluid velocity ($U_{std}$): $U/D$

  12. 12.

Logarithm of kinematic viscosity ($log(\nu)$): $log_{10}(\nu)$.

The resulting 11 variables where used to train the classifiers and the
regressors showcased in the following sections.

####  8.5.3 Training and Identifying Best Classifiers

The ML classifiers shown in Table 8.1 were trained and tested by using k-folds
cross-validation (k=5) on all 2000 cases, resulting in training sample of 1600
and a testing sample of 400. These classifiers were chosen by virtue of their
differing approaches, as we had little-to-no intuition of which one would work
best for our purposes (no one has ever used ML to predict clogging
mechanisms). Classifier performance was evaluated by calculating their testing
accuracy, precision, the area under their Receiving Operating Characteristic
(ROC) curve, and quantifying the total number of true positives (TP), false
positives (FP), true negatives (TN), and false negatives (FP). In this case, a
true positive pertains to correctly labeling a case as “Clogged”. These
metrics where chosen due to their ability to quantify overall classifier
performance (accuracy), while also testing for its ability to correctly
predict positive clogging labels (precision + ROC). Furthermore, we
implemented a grid-search algorithm in order to optimize our choice of model
hyperparameters and thus improve the predictive power of each type of
classifier. Table 8.1 shows the performance metrics of each classifier as a
result of average precision grid-search optimization.

Table 8.1:  Clogging prediction performance of the chosen ML classifiers. The
highlighted sections represent the classifiers with the best performance.

The results in Table 8.1 are very encouraging, showing that the best two
classifiers, Extra Trees and the Multi-Layer Perceptron (MLP), reliably
achieve a classifying precision of 0.94 and an accuracy of 0.96, meaning their
false positive/negative rates are practically indistinguishable. Further
analysis into the misclassified samples show that all of them are partially-
clogged samples that lie within the classification decision boundary shown in
Figure 8.5B. Please refer to Section 8.6 for a description of the optimized
models.

####  8.5.4 Classification Feature Analysis

Although the testing metrics of our models are relatively high, it might not
be practical for a potential user to have to obtain all 11 training variables
for every system of interest that he/she wishes to test. It might be more
efficient to reduce a model’s accuracy in favor of reducing its complexity and
number of features. As such, we now aim to identify the most predictive
features of the two most succesful models identified above.

To do so, we implemented a permutation feature importance algorithm, where we
quantified the decrease in the models’ precision as a result of shuffling the
values of a single feature. This effectively disassociates the feature values
from each sample’s labels. The relative magnitude of the drop in the model’s
precision is proportional to the relative importance of the shuffled feature.
This drop is then normalized with respect to all other shuffled features, and
a feature importance ranking is created. In our case, we chose to shuffle and
measure the score effect of each feature 50 times before continuing on to the
next one.

Table 8.2:  Tabulated relative feature importance for the top 6 classifier
features. The rightmost column identifies the accuracy of a model trained and
tested exclusively with the top “n” features.

Table 8.2 shows the top 6 features for our optimized Extra Trees and MLP
classifiers. Not surprisingly, both cases have the ratio between particle
diameter and pore size ($d/D$) as their top feature, followed by the
standardized particle flux ($F_{nd}$) and the ratio between the particle
diameter and the geometric thickness ($d/T$). These are then followed by the
non-dimensionalized attractive forces and the fluid viscosity. What is
surprising however, is the relatively small predictive power that $d/D$ has,
specially in the Extra Trees classifier, where a model trained and tested with
only said feature cannot predict clogging any better than a coin flip. The
addition of the next two features increases prediction accuracy to about 0.88
for both classifiers, meaning that the bulk of the clogging cases can be
explained/predicted by just 2 complementary physical processes: particle size
exclusion and particle flux. The remaining cases rely on more complex physics,
where it is important to consider particle-particle-wall attractions and
viscous flow-particle couplings.

Porosity, pore size and grain heterogeneity, and the normalized fluid velocity
do not seem to have as much as an influence on clogging. This might be
explained by two different factors: 1) We did not probe the sections of the
11-dimensional parameter space for which these features become significant. 2)
These features really do not affect clogging in a significant way. Since fluid
viscosity (i.e. drag) did appear to have a quantifiable influence on clogging
mechanics, we also believe that fluid velocity should play a part as well. For
this reason we are inclined to believe the first explanation is more likely,