Dr. Ian C. Bourg

Doctor of Philosophy Chemical and Biological Engineering 2021 September
Chemical and Biological Engineering

B.S. University of Texas at Austin 2016

M.A. Princeton University 2018

# Modeling Multiphase Flow Through and Around Multiscale Deformable Porous
Materials

Francisco J. Carrillo

###### Contents

  1. 1 Introduction
  2. 2 Derivation of the PDEs for Multiphase Flow in Deformable Porous Media
    1. 2.1 Volume Averaging
    2. 2.2 Solid Mass Conservation Equation
    3. 2.3 Fluid Mass Conservation Equation
    4. 2.4 Saturation Convection Equation
    5. 2.5 Fluid Momentum Conservation Equation
    6. 2.6 Closure of the Fluids’ Multiscale Interaction Terms
      1. 2.6.1 Matching the Volume-of-Fluid Model in Solid-Free Regions
      2. 2.6.2 Matching the Standard Two-Phase Darcy Model in Porous Regions
    7. 2.7 Momentum Conservation Equation for a Linear Elastic Solid
    8. 2.8 Momentum Conservation Equation for a Viscoplastic Solid
    9. 2.9 General Solid Momentum Conservation Equation
    10. 2.10 Closure of the Solid’s Multiscale Interaction Terms
    11. 2.11 Interfacial Conditions between Solid-Free Regions and Porous Regions
    12. 2.12 Advection-Diffusion Equation for Tracer Particles
    13. 2.13 Summary of Equations and Multiscale Parameters
    14. 2.14 Conclusion
  3. 3 Numerical Implementation
    1. 3.1 Numerical Platform
    2. 3.2 Solution Algorithm
    3. 3.3 Open-Source Implementation
  4. 4 Applications to Single-Phase Flow in Multiscale Deformable Porous Media
    1. 4.1 Introduction
    2. 4.2 Parametrization to Porous Media: Na-Smectite Clay
      1. 4.2.1 Parametrization of the Permeability Function
      2. 4.2.2 Parametrization of the Swelling Pressure Term
      3. 4.2.3 Parametrization of the Diffusion Coefficient of Dissolved Salts
      4. 4.2.4 Parametrization of Clay Rheology
    3. 4.3 Model Validation and Application to Plastically Deformable Porous Media (Clay)
      1. 4.3.1 Model Verification by Comparison to the Fluid-Driven Deformation of a Clay Filter Cake
      2. 4.3.2 Model Verification by Comparison to Oedometric Clay Swelling Experiments
      3. 4.3.3 Model Prediction of the Permeability of a Bead Pack Containing Na-Smectite as a Function of Salinity and Clay Content
    4. 4.4 Model Application to Fluid Injection Into a Poroelastic Medium
    5. 4.5 Conclusions
  5. 5 Applications to Multiphase-Phase Flow in Static Multiscale Porous Media
    1. 5.1 Darcy Scale Validation
      1. 5.1.1 Buckley-Leverett
      2. 5.1.2 Gravity dominated Buckley-Leverett
      3. 5.1.3 Gravity-capillarity equilibrium
      4. 5.1.4 Darcy scale application: Oil drainage in a heterogeneous reservoir
    2. 5.2 Pore Scale Validation
      1. 5.2.1 Contact angle on a flat plate
      2. 5.2.2 Capillary rise
      3. 5.2.3 Taylor film
      4. 5.2.4 Pore scale application: Oil drainage in a complex pore network
    3. 5.3 Hybrid Scale Applications
      1. 5.3.1 Wave propagation in a coastal barrier
      2. 5.3.2 Drainage and imbibition in a fractured microporous matrix
    4. 5.4 Conclusions
  6. 6 Application to Multiphase Flow in Multiscale Deformable Porous Media
    1. 6.1 Model Validation
      1. 6.1.1 Terzaghi Consolidation Problem
      2. 6.1.2 Pressure Oscillation in Poroelastic Core
      3. 6.1.3 Capillary Pressure Effects in a Poroelastic Column
      4. 6.1.4 Fluid Invasion and Fracturing in a Hele-Shaw Cell
      5. 6.1.5 Modeling Fracturing Wellbore Pressure
      6. 6.1.6 Capillary Effects on Fracturing Wellbore Pressure
    2. 6.2 Illustrative Applications
      1. 6.2.1 Elastic Failure in Coastal Barriers
      2. 6.2.2 Flow-Induced Surface Deformation
    3. 6.3 Conclusions
  7. 7 Capillary and Viscous Fracturing During Drainage in Porous Media
    1. 7.1 Introduction
    2. 7.2 Numerical Simulations
      1. 7.2.1 Crossover from Imbibition to Fracturing in a Hele-Shaw Cell
      2. 7.2.2 Creation of Fracturing Phase Diagrams
    3. 7.3 Characterization of Fracturing Mechanisms
    4. 7.4 Influence of Localized and Uniform Deformation
      1. 7.4.1 Localized Deformation
      2. 7.4.2 Uniform Deformation
    5. 7.5 Conclusions
  8. 8 Simulation and Prediction of Stochastic Clogging Processes in Porous Media
    1. 8.1 Introduction
    2. 8.2 Methods
      1. 8.2.1 CFD-DEM
      2. 8.2.2 Machine Learning
    3. 8.3 Workflow
    4. 8.4 Experimental Setup
      1. 8.4.1 Creating Randomized Porous Systems
      2. 8.4.2 CFDEM Simulations
    5. 8.5 Data Analysis and Model Training
      1. 8.5.1 Identifying Clogged Systems
      2. 8.5.2 Feature Standardization and Non-dimensionalization
      3. 8.5.3 Training and Identifying Best Classifiers
      4. 8.5.4 Classification Feature Analysis
      5. 8.5.5 Training and Identifying Best Regressors
      6. 8.5.6 Regression Feature Analysis
    6. 8.6 Conclusions
  9. 9 Conclusions
  10. A Relative Permeability and Capillary Pressure Models
    1. A.1 Relative Permeability Models
    2. A.2 Capillary Pressure Models
  11. B Solid Rheology Models
    1. B.1 Hershel-Bulkley Plasticity
    2. B.2 Quemada Rheology Model
    3. B.3 Linear Elasticity
  12. C Fracturing Instabilities
  13. D Additional Derivation Steps
    1. D.1 Skipped Steps in the Derivation of the Fluid Continuity Equation
    2. D.2 Derivation of the Single-Field Viscous Stress Tensor
    3. D.3 Derivation of the Single-Field Material Derivative
    4. D.4 Recovery of Biot Theory from the Fluid and Solid Momentum Equations
    5. D.5 Integral Geometric Relation for Vectors
  14. E Semi-Analytic Solution for the Seismic Stimulation of a Poroelastic Core
  15. References

###### List of Figures

  1. 1.1 Conceptual representation of a multiscale deformable porous medium and its related processes
  2. 1.2 Schematic representations of a porous medium with two characteristic pore sizes
  3. 2.1 Conceptual representation of the Multiphase DBB modelling framework
  4. 2.2 Distribution of the fluid phases in a continuous physical domain and a discrete Eulerian grid
  5. 4.1 Conceptual representation of the two key challenges posed by fluid flow in clay-rich sedimentary rocks
  6. 4.2 Illustration of the three main theoretical frameworks combined in this study to model fluid flow in porous media with two characteristic length scales and a deformable microporous matrix
  7. 4.3 Parameterization of clay models (Na-smectite)
  8. 4.4 Fluid-driven compaction of a clay filter cake
  9. 4.5 Macroscopic swelling of NaCl-saturated smectite clay
  10. 4.6 Effect of clay mass fraction and NaCl concentration on the spatial distribution of smectite clay and available macropore space in a 2-D cross section of an ideal 3-D spherical packed bed
  11. 4.7 Steady-state fluid velocity profiles with a 10 % clay mass fraction at different salinity values
  12. 4.8 Permeability of a simplified sedimentary rock model as a function NaCl concentration for a range of smectite clay mass fractions
  13. 4.9 Fluid-driven deformation of a poroelastic membrane
  14. 5.1 Comparison of predicted saturation profiles against Buckley-Leverett’s solution
  15. 5.2 Comparison of predicted saturation profiles against Buckley-Leverett’s solution with gravity
  16. 5.3 Gravity-capillary balance in an unsaturated column
  17. 5.4 Simulation setup for oil drainage within a heterogeneous reservoir.
  18. 5.5 Oil drainage in a heterogeneous porous medium solved at the continuum scale
  19. 5.6 Compilation of all test cases performed for the verification of the solver within the Navier-Stokes domain
  20. 5.7 Oil drainage in a complex porous medium solved at the pore scale
  21. 5.8 Wave dissipation in a porous coastal barrier
  22. 5.9 Drainage and imbibition in a microporous fracture
  23. 6.1 One-dimensional Terzaghi consolidation problem
  24. 6.2 Pressure oscillation in poroelastic core
  25. 6.3 Capillary effects in a poroelastic column
  26. 6.4 Experimental and simulated fracturing in a Hele-Shaw cell
  27. 6.5 Conceptual representation of wellbore pressure evolution during fluid-induced fracturing of low permeability rocks
  28. 6.6 Wellbore pressure as a function of injection rate and time
  29. 6.7 Effect of capillary entry pressure on fracturing wellbore pressure
  30. 6.8 Wave dissipation in a poroelastic coastal barrier
  31. 6.9 Study of the effects of hydraulic fracturing and cracking on surface deformation
  32. 7.1 Continuous transition from fluid imbibition to fracturing in a Hele-Shaw cell
  33. 7.2 Phase diagrams describing the impacts of permeability, plastic yield stress, fluid injection rate, and capillary entry pressure on the transition from fluid drainage to fracturing
  34. 7.3 Fluid invasion and fracturing in plastic porous media as a function of the viscous and capillary fracturing numbers
  35. 7.4 Dynamic fracture formation mechanisms
  36. 7.5 Influence of the porosity-dependence of capillary pressure on fluid invasion and fracturing patterns
  37. 7.6 Fracturing numbers’ dependence on porosity
  38. 8.1 Conceptual representation of the models used for clogging simulations and predictions
  39. 8.2 Examples of geometries used for clogging simulations.
  40. 8.3 Examples of different clogging simulations
  41. 8.4 Examples of different clogging levels
  42. 8.5 Clogging as a function of the average final distance between particles and the symmetry of the particles’ final velocity distribution.
  43. 8.6 Optimized Random Forest regression results
  44. C.1 Effects of the solid fraction field’s standard deviation on fracturing
  45. C.2 Effects of different solid fraction field initializations on fracturing

###### Acknowledgements.

[75mm] Research is what I’m doing when I don’t know what I’m doing. Wernher
von Braun

##  Chapter 1 Introduction

Fluid flow in deformable porous media is a ubiquitous phenomenon with
important implications in many energy and environmental technologies including
geologic ${\mathrm{CO}}_{2}$ sequestration, soil bioremediation, water
treatment, enhanced biochemical production, nuclear waste disposal, and fuel
cell design (Bächer & Gekle,, 2019; Bock et al.,, 2010; Cunningham et al.,,
2003; Räss et al.,, 2018; Towner,, 1987). It also underlies iconic geophysical
features at many scales, from coastal, riparian, and volcanic landforms to
fractures in subsurface reservoirs, cracks in clay soils, and bubbles in soft
sediments.

Whereas single-phase flow in porous media is relatively well understood from
atomistic to continuum scales, the dynamics of systems containing multiple
phases remain challenging to describe at all scales (Gray et al.,, 2015; Li et
al.,, 2018). Multiphase flow involves strong feedback between inertial,
viscous, capillary, and interfacial forces (Meakin & Tartakovsky,, 2009; Datta
et al.,, 2014). This coupling is intrinsically multiscale, as inertial and
viscous forces dominate in large pores or fractures while capillary forces and
interfacial energetics dominate within smaller porous or microporous
structures. The complex linkage between microscopic geometric heterogeneities
and macroscopic processes makes it necessary to consider multiple scales
across porous media in order to create truly predictive models, from the scale
of microscopic interfaces ($\sim\mu$m), to pore networks and lab columns
($\sim$cm), all the way up to the field scale ($\sim$km).

An important and largely unresolved challenge in the areas outlined above is
the difficulty of describing situations where multiple fluids interacts with a
deformable porous material. For example, when modeling flow through biofilms
or membranes it is imperative to understand how fluid flow behaves inside the
microporous medium (in pores with length scales of $\sim 10^{-6}\ \mathrm{m}$)
while simultaneously understanding how the deformation of this medium affects
the overall flow field (often controlled by much larger flow paths with length
scales on the order of $\sim 10^{-2}\ \mathrm{m}$) (Bottero et al.,, 2010).
Similarly, the propagation of flow-driven fractures in porous materials and
the propagation of waves in coastal barriers involve feedbacks between flow
and mechanics in systems with characteristics pore widths that differ by three
or more orders of magnitude. In this dissertation, we develop a framework
capable of representing multiphase flow and solid mechanics in systems with
two characteristic pore length scales, as required to simulate many of the
aforementioned phenomena (see Figs. 1.1 and 1.2).

Figure 1.1: Conceptual representation of a multiscale deformable porous medium
and its related processes. The porous domain is shown in the lower half
(orange), the free-fluid domain is shown in the upper half (blue), two
immiscible fluids (left and right) are shown in different shades of blue and
are separated by an interface (black). $\theta$ is the contact angle.

The starting point for our study is based on the present ample understanding
of multiphase flow dynamics within and around static porous materials, from
viscous and capillary fingering (Ferer et al.,, 2004; Lenormand & Zarcone,,
1989; Lenormand et al.,, 1988) to temperature and surface tension driven flows
(Shih & Megaridis,, 1996), all the way to turbulent multiphase flows (Colombo
& Fairweather,, 2015; Soulaine & Quintard,, 2014). This knowledge, in
conjunction with numerical techniques such as the Lattice Boltzmann Method,
the Finite Volume Method, Homogenization Theory, and Averaging Theory, forms
the basis of fast and accurate models that are routinely applied to help
design and improve hydrocarbon production (Burrus et al.,, 1991; Mehmani &
Tchelepi,, 2019), ${\mathrm{CO}}_{2}$ sequestration (Hassan & Jiang,, 2012),
and even nuclear reactors (Tentner et al.,, 2008). However, the study of
multiphase flow across different scales remains limited as shown by the
absence of well-established approaches to describe how bubbles or waves
propagate into an unsaturated porous medium or how a multiphase fluid mixture
is pushed out of a porous medium into open space.

A similar situation pertains with regard to computational models that couple
fluid flow and solid mechanics. Theoretical and numerical approaches based on
Biot’s Theory of poroelasticity (Biot,, 1941), Terzaghi’s effective stress
principle (Terzaghi,, 1943), and Mixture Theory (Siddique et al.,, 2017) have
been successful at modeling systems with flow in deformable porous media
including arteries, biofilms, boreholes, hydrocarbon reservoirs, seismic
systems, membranes, soils, swelling clays, and fractures (Auton & MacMinn,,
2017; Barry et al.,, 1997; Jha & Juanes,, 2014; Lo et al.,, 2005, 2002;
MacMinn et al.,, 2016; Mathias et al.,, 2017; Santillán et al.,, 2017).
However, as mentioned above, we still have very little understanding of how
flow-induced deformation of these solid materials affects the macroscopic flow
around them (and thus their boundary conditions) or how fluid-fluid interfaces
behave when pushed against a soft porous medium and vice-versa.

Three major approaches have been proposed to resolve the challenge posed by
fluid flow in porous media containing both solid-free regions and microporous
domains (hereafter referred to as multiscale systems). The most
straightforward of these involves performing direct numerical simulations
(DNS) throughout the entire multiscale domain, both within and outside the
porous medium (Breugem & Boersma,, 2005; Hahn et al.,, 2002; Krafczyk et al.,,
2015). Although rigorous, this technique is impractical in situations with a
large difference in length scales between the largest and smallest pores,
where it requires exceedingly fine grids and tremendous computational
resources.

To save time and resources, other studies have relied on hybrid DNS-Darcy
approaches, where fluid and solid mechanics within a porous medium are modeled
as averaged quantities through Darcy’s law, pore-network models, or Biot
theory (Weishaupt et al.,, 2019; Ehrhardt,, 2010). One such approach relies on
the use of the Beavers-Joseph (BJ) boundary condition to couple fluid flow in
solid-free domains (simulated using the Navier-Stokes equations) and in
microporous domains (simulated using Darcy’s law) for single phase flow and
static porous media (Beavers & Joseph, 1967a, ; Fetzer et al.,, 2016). Recent
studies have extended this BJ approach to allow multiphase flow in the solid-
free domain (Baber et al.,, 2016) or to include the effects of poroelasticity
within the porous medium (Lacis et al.,, 2017; Zampogna et al.,, 2019).
However, to the best of our knowledge, no BJ based technique has yet been
developed to couple solid mechanics with multiphase flow simultaneously within
the solid-free and porous domains.

The Darcy-Brinkman (DB) approach –also referred to as Darcy-Brinkman-Stokes
(DBS) equations– (Brinkman,, 1947) presents a well-known alternative to the BJ
interface matching technique. These equations arise from volume averaging the
Stokes (or Navier-Stokes) equations in a control volume that contains both
fluids and solids (Vafai & Tien,, 1981; Hsu & Cheng,, 1990; Bousquet-Melou et
al.,, 2002; Goyeau et al.,, 2003). It consists in a Stokes-like momentum
equation that is weighted by porosity and contains an additional drag force
term that describes the mutual friction between the fluids and solids within
said control volume. Unlike standard continuum scale equations for flow and
transport in porous media such as Darcy’s law, the DB equation remains valid
in solid-free regions (see Figure 1.2A) where the drag force term vanishes and
the DB equation turns into the Stokes (or Navier-Stokes) equation. In porous
regions (see Figure 1.2C), in contrast, viscous dissipation effects are
negligible compared with the drag force exerted onto the pore walls and the DB
momentum equation tends asymptotically towards Darcy’s law (Tam,, 1969;
Whitaker, 1986a, ; Auriault,, 2009). Therefore, the “micro-continuum” DB
equation has the ability to simultaneously solve flow problems through porous
regions and solid-free regions (Neale & Nader,, 1974), paving the path to
hybrid scale modeling (see Figure 1.2B). In the case of single phase flow, it
is known to be analogous (in fact, formally equivalent) to the previously
mentioned and well-established Beavers-Joseph boundary conditions (Beavers &
Joseph, 1967b, ; Neale & Nader,, 1974).

Figure 1.2: Schematic representations of a porous medium with two
characteristic pore sizes depending on the scale of resolution: (a) full pore
scale (Navier-Stokes), (b) intermediate or hybrid scale, and (c) full
continuum scale (Darcy). Our objective is to derive a framework that can
describe multiphase flow at all three scales described in the figure based on
a single set of equations resolved throughout the entire system.

The resulting so-called “micro-continuum” approach has been extensively used
to solve single phase flow through static multiscale porous media. A prime
example is the modelling of fluid flow in three-dimensional images of rock
samples that contain unresolved sub-voxel porosity (Knackstedt et al.,, 2006;
Apourvari & Arns,, 2014; Scheibe et al.,, 2015; Soulaine et al.,, 2016; Kang
et al.,, 2019; Singh,, 2019). It also has been used to simulate dissolution
wormholing during acid stimulation in cores by updating the weighting porosity
field through geochemical reactions (Liu et al.,, 1997; Golfier et al.,, 2002;
Soulaine & Tchelepi,, 2016; Tomin & Voskov,, 2018). Moreover, it has been
shown that whenever low-porosity low-permeability porous regions are present,
the velocity within these regions drops to near zero, such that the micro-
continuum DB framework can be used as a penalized approach to map a solid
phase onto a Cartesian grid with a no-slip boundary at the solid surface
(Angot et al.,, 1999; Khadra et al.,, 2000; Soulaine & Tchelepi,, 2016).
Therefore, this approach can be used to move fluid-solid interfaces
efficiently in a Cartesian grid without a re-meshing strategy. For example,
Soulaine et al., (2017) used a micro-continuum framework to predict the
dissolution kinetics of a calcite crystal and successfully benchmarked their
model against state-of-the-art pore scale dissolution solvers with evolving
fluid-solid interfaces (Molins et al.,, 2019).

In the present thesis, we propose the Multiphase Darcy-Brinkman-Biot (DBB)
approach, a fully coupled multiscale model for two-phase flow in deformable
porous media based on the micro-continuum approach, rooted in elementary
physical principles and rigorously derived using the method of the volume
averaging (S. Whitaker,, 1999). We show that there exists a single set of
partial differential equations that can be applied in pore, continuum, and
hybrid scale representation of multiphase flow in porous media. Particular
attention is paid to the derivation of gravity and capillary effects in the
porous domain for both fluid and solid mechanics. The resulting two-phase
micro-continuum framework is validated using an extensive series of test cases
where reference solutions exist. We verify that the multiscale solver
converges to the standard Darcy/Biot scale solutions (Buckley-Leverett,
capillary-gravity equilibrium, drainage in a heterogeneous reservoir, Terzaghi
consolidation tests, clay swelling experiments) when used at the continuum
scale in porous media and to the two-phase Navier-Stokes solutions (droplet on
a flat surface, capillary rise, drainage with film deposition, two-phase flow
in a complex porous structure) when used at the pore scale. The fully
implemented numerical model, along with the aforementioned verification and
tutorial cases, is provided as an accompanying open-source solver:
hybridBiotInterFoam.

This thesis is organized as follows. Chapter 2 introduces the concept of
volume averaging and describes the derivation of the governing equations for
coupled fluid and solid mechanics. Chapter 3 explains the numerical
implementation and algorithm development for the coupled mass and momentum
equations and introduces the resulting open-source solver. Then, Chapters 4 to
6 focus on verifying and showcasing the model’s ability to capture coupled
fluid-solid mechanics in multiscale porous media. This was done in an
incremental manner, where Chapter 4 focuses on single phase flow through
deformable porous media, Chapter 5 on multiphase flow through static porous
media, and Chapter 6 on multiphase flow through deformable porous media. After
that, Chapter 7 uses the developed model to identify the governing non-
dimensional parameters that predict capillary and viscous fracturing in porous
media. Lastly, Chapter 8 presents a slight detour through the use of Machine
Learning approaches to predict and simulate stochastic clogging process in
heterogeneous porous media. Chapter 9 concludes with a summary of this thesis
and a discussion on future work.

[75mm] The difference between screwing around and science is writing it down.
Adam Savage

##  Chapter 2 Derivation of the PDEs for Multiphase Flow in Deformable Porous
Media

In this chapter, we present the complete derivation of the micro-continuum
equations for multiphase flow in static and deformable porous media. This
Computational Fluid Dynamics (CFD) approach allows for simulation of two-phase
flow in hybrid systems containing solid-free regions and porous matrices, as
illustrated schematically in Fig. 2.1. Our approach consists of a unique set
of five volume-averaged partial differential equations that asymptotically
approach the Navier-Stokes Volume-of-Fluid equations in solid-free-regions and
multiphase Biot Theory in porous regions. This set of equations consists of 1)
a solid mass conservation equation, 2) a fluid mass conservation equation, 3)
a fluid saturation conservation equation, 4) a fluid momentum conservation
equation, and 5) a solid momentum conservation equation. This work is adapted
from Carrillo & Bourg, (2019), Carrillo et al., (2020), and Carrillo & Bourg,
2021b .

Figure 2.1:  Conceptual representation of the Multiphase DBB modelling
framework. The porous domain is shown in the lower half (orange), the free-
fluid domain is shown in the upper half (blue), the two immiscible fluids
(left and right) are shown in different shades of blue and are separated by an
interface (black), $\theta$ is the contact angle, and ${\phi}_{f}$ is the
porosity. REV is the “Representative Elementary Volume” over which all
conservation equations are averaged. Note that the stated relation between the
averaging volume’s length scale $L_{V}$ and the porous length scale $L_{P}$ is
required for the creation of a REV.

###  2.1 Volume Averaging

We start by introducing the concept of volume averaging. This technique forms
the basis of the micro-continuum equations, as it allows the classical mass
and momentum conservation equations to account for the coexistence of solids
($s$), wetting fluids ($w$), and non-wetting fluids ($n$) within a given
control volume. This method is well suited for use in conjunction with the
Finite Volume Method (Patankar,, 1980), as the latter method’s numerical grid
provides an intuitive and straightforward numerical interpretation of what we
will define as the averaging volume ($V$). We start by defining the volume
averaging operator:

${\overline{\beta}}_{i}=\frac{1}{V}\int_{V_{i}}{{\beta}_{i}dV}$ (2.1)

where ${\beta}_{i}$ is a function defined in each phase’s respective volume
$V_{i}\left(i=w,n,s\right).$ We also define the phase averaging operator:

${\overline{\beta}}^{i}_{i}=\frac{1}{V_{i}}\int_{V_{i}}{{\beta}_{i}dV}$ (2.2)

Each fluid’s ($f$) phase averaged variables are thus intrinsically related by
the porosity ${\phi}_{f}=V_{f}/V=(V_{w}+V_{n})/V$ and saturation
${\alpha}_{i}=V_{i}/V_{f}$ fields (where $V_{f}\equiv V_{w}+V_{n}$) such that
${\overline{\beta}}_{i}={\phi}_{f}{\alpha}_{i}{\overline{\beta}}^{i}_{i}\
\left(i=w,n\right)$. For phase-averaged solid variables, the equivalent
relationship only involves the solid fraction ${\phi}_{s}$, such that:
${\overline{\beta}}_{s}=\ {\phi}_{s}{\overline{\beta}}^{s}_{s}$. We note that
${\phi}_{f}+{\phi}_{s}=1$ and ${\alpha}_{1}+{\alpha}_{2}=1$; thus knowledge of
one of the $\phi$ or $\alpha$ variables always leads to knowledge of the
other. Volume averaging then allows for the definition of several regions
within a hybrid-scale numerical simulation:

$\phi_{f}=\begin{cases}1,&\textnormal{in solid-free regions}\\\
\left]0;1\right[,&\textnormal{in porous regions}\end{cases}$ (2.3)

$\alpha_{w}=\begin{cases}0,&\textnormal{in regions saturated with non wetting
fluid}\\\ \left]0;1\right[,&\textnormal{in unsaturated regions}\\\
1,&\textnormal{in regions saturated with wetting fluid}\\\ \end{cases}$ (2.4)

The application of an averaging transformation to mass and momentum
conservation equations will result in variables and equations that are
weighted differently in each of the regions identified in Eqns. 2.3-2.4.
However, the averaging of differential equations is not straightforward, which
is why we introduce the following spatial averaging theorems for volumes
containing three distinct phases (Howes & Whitaker,, 1985; S. Whitaker,,
1999):

$\overline{\frac{\partial{\beta}_{i}}{\partial t}}=\
\frac{\partial\overline{{\beta}_{i}}}{\partial
t}-\frac{1}{V}\int_{A_{i,j}}{{\beta}_{i}{\boldsymbol{v}}_{i,j}\cdot{\boldsymbol{n}}_{i,j}dA}-\frac{1}{V}\int_{A_{i,k}}{{\beta}_{i}{\boldsymbol{v}}_{i,k}\cdot{\boldsymbol{n}}_{i,k}dA}$
(2.5)
$\overline{\nabla{\beta}_{i}}=\nabla\overline{{\beta}_{i}}+\frac{1}{V}\int_{A_{i,j}}{{\beta}_{i}{\boldsymbol{n}}_{i,j}dA}+\frac{1}{V}\int_{A_{i,k}}{{\beta}_{i}{\boldsymbol{n}}_{i,k}dA}\
\ \ $ (2.6)
$\overline{\nabla\cdot{\boldsymbol{\beta}}_{i}}=\nabla\overline{\boldsymbol{\beta}_{i}}+\frac{1}{V}\int_{A_{i,j}}{{\boldsymbol{\beta}}_{i}\cdot{\boldsymbol{n}}_{i,j}dA}+\frac{1}{V}\int_{A_{i,k}}{{\boldsymbol{\beta}}_{i}{\cdot\boldsymbol{n}}_{i,k}dA}\
\ \ $ (2.7)

where $A_{i,j}$ represents the interfacial area between phase $i$ and $j$,
$\boldsymbol{n}_{i,j}$ is a vector normal to the interface and oriented toward
phase $j$, and ${\boldsymbol{v}}_{i,j}$ is the velocity of the interface. The
notation for symbols with subscript pair $i,k$ is equivalent, and the symbols
$i,j,k$ represent any combination of the solid, wetting, and non-wetting
phases. The surface integrals in Eqns. 2.5 \- 2.7 are crucial components of
the following derivations as they convert the interfacial conditions at the
fluid-fluid and fluid-solid interfaces into body forces within the averaged
partial differential equations.

The following properties will also become useful, which follow directly from
the basic averaging theorem in a system with a single fluid and solid phase
(Whitaker, 1986a, ).

$\ \frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{n}}_{f,s}dA}\ =\ -\nabla{\phi}_{f}\
$ (2.8)

However, if the integral is over a fluid-solid surface in a multiphase system,
the previous equation needs to be modified to account for the fact that fluid
phase $i$ only partially covers the full solid surface,

$\ \frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{n}}_{i,s}dA}\ =\
-{\alpha}_{i}\nabla{\phi}_{f}\ $ (2.9)

The original property can then be easily recovered by integrating over the
solid surface spanned by both fluid phases:

$\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{n}}_{i,s}dA}+\
\frac{1}{V}\int_{A_{j,s}}{{\boldsymbol{n}}_{j,s}dA}\ =\
-{\alpha}_{i}\nabla{\phi}_{f}-{\alpha}_{j}\nabla{\phi}_{f}=\
-\nabla{\phi}_{f}\ $ (2.10)

Figure 2.2: Distribution of the fluid phases in (a) the continuous physical
domain, (b) a discrete Eulerian grid. In this figure, the subscripts $g$ and
$l$ represent gaseous and liquid phases, respectively.

###  2.2 Solid Mass Conservation Equation

We start the derivation of our solid mass conservation equation with an
elementary mass conservation equation:

$\frac{\partial{\rho}_{s}}{\partial
t}+\nabla\cdot\left({\rho}_{s}{\boldsymbol{U}}_{s}\right)=0$ (2.11)

where ${\rho}_{s}$ and $\boldsymbol{U}_{s}$ are the solid’s density and
velocity, respectively. We then apply the averaging operator within a volume
containing the solid ($s$) and two immiscible fluid phases ($w,n$) such that:

$\overline{\frac{\partial{\rho}_{s}}{\partial
t}}+\overline{\nabla\cdot\left({\rho}_{s}{\boldsymbol{U}}_{s}\right)}=\\\
\frac{\partial\overline{{\rho}_{s}}}{\partial
t}-\frac{1}{V}\int_{A_{s,w}}{{\rho}_{s}{\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,n}}{{\rho}_{s}{\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,n}dA}\\\
+\nabla\cdot\left(\overline{{\rho}_{s}{\boldsymbol{U}}_{s}}\right)+\frac{1}{V}\int_{A_{s,w}}{{\rho}_{s}{\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,w}dA}+\frac{1}{V}\int_{A_{s,n}}{{\rho}_{s}{\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,n}dA}$
(2.12)

It then becomes clear that the integral operators that arise from the
transient term cancel with the ones that arise from the convection term such
that:

$\frac{\partial\overline{{\rho}_{s}}}{\partial
t}+\nabla\cdot\left(\overline{{\rho}_{s}{\boldsymbol{U}}_{s}}\right)=0$ (2.13)

We can now decompose the equation in terms of the intrinsic phase averages,

$\frac{\partial{\phi}_{s}\overline{{\rho}^{s}_{s}}}{\partial
t}+\nabla\cdot\left({\phi}_{s}{\overline{\rho}}^{s}_{s}\
\overline{{\boldsymbol{U}}^{s}_{s}}\right)=0$ (2.14)

Assuming the solid density is constant and given that intrinsic phase averages
of constants are equal to themselves, we can state the following:

$\frac{\partial{\phi}_{s}}{\partial
t}+\nabla\cdot\left({\phi}_{s}\overline{{\boldsymbol{U}}^{s}_{s}}\right)=0$
(2.15)

This finalizes the derivation of the solid mass conservation equation.

###  2.3 Fluid Mass Conservation Equation

Just as before, we start with a mass conservation equation for a given fluid
$"i"$ $\left(i=w,n\right)$,

$\frac{\partial{\rho}_{i}}{\partial
t}+\nabla\cdot\left({\rho}_{i}{\boldsymbol{U}}_{i}\right)=0$ (2.16)

where ${\rho}_{i}$ and $\boldsymbol{U}_{i}$ are the fluid’s density and
velocity, respectively. We then follow the same procedure outlined in the
Section 2.2 (also shown in Appendix D.1) in order to obtain the following
equation:

$\frac{\partial\overline{{\rho}_{i}}}{\partial
t}+\nabla\cdot\left(\overline{{\rho}_{i}{\boldsymbol{U}}_{i}}\right)=0$ (2.17)

Just as before, we want to decompose the averaged terms into their
intrinsically averaged components. This time, however, we do so by noting that
${\overline{\beta}}_{i}={\phi}_{f}{\alpha}_{i}{\overline{\beta}}^{i}_{i}$.

$\frac{\partial{\phi}_{f}{\alpha}_{i}\overline{{\rho}^{i}_{i}}}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{i}\overline{{\rho}^{i}_{i}}\overline{{\boldsymbol{U}}^{i}_{i}}\right)=0$
(2.18)

Assuming the fluid density is constant and given that intrinsic phase averages
of constants are equal to themselves, we can show that:

$\frac{\partial{\phi}_{f}{\alpha}_{i}}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{i}\overline{{\boldsymbol{U}}^{i}_{i}}\right)=0$
(2.19)

This step completes the derivation of the mass conservation equation for a
given fluid “$i$”. We now sum together the mass conservation equations for
both wetting ($w)$ and non-wetting ($n)$ fluids, such that:

$\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\frac{\partial{\phi}_{f}{\alpha}_{n}}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}\overline{{\boldsymbol{U}}^{w}_{w}}\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{n}\overline{{\boldsymbol{U}}^{n}_{n}}\right)=$
(2.20)

$\frac{\partial{\phi}_{f}({\alpha}_{w}+{\alpha}_{n})}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}\overline{{\boldsymbol{U}}^{w}_{w}}+{\phi}_{f}{\alpha}_{n}\overline{{\boldsymbol{U}}^{n}_{n}}\right)=0$
(2.21)

The result is the single-field mass conservation equation in terms of the
porosity and the single-field fluid velocity
${\boldsymbol{U}}_{f}={\phi}_{f}\left[{\alpha}_{w}{\overline{\boldsymbol{U}}}^{w}_{w}\
{+\ \alpha}_{n}{\overline{\boldsymbol{U}}}^{n}_{n}\right]$,

$\frac{\partial{\phi}_{f}}{\partial t}+\nabla\cdot{\boldsymbol{U}}_{f}=0$
(2.22)

where the phrase “single-field” refers to averaged variables that depend on
the properties of both fluids. This concludes the derivation of the volume
averaged fluid mass conservation equation.

###  2.4 Saturation Convection Equation

We start this derivation with the volume-averaged mass conservation equation
for the wetting fluid phase (Equation 2.19 in Section 2.3):

$\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}\overline{{\boldsymbol{U}}^{w}_{w}}\right)=0$
(2.23)

This equation can be used to advect the saturation function ${\alpha}_{w}.$
However, for practical reasons that will be made clear later, we wish to
express this as a function of the single-field velocity ${\boldsymbol{U}}_{f}$
(rather than $\overline{{\boldsymbol{U}}^{w}_{w}}$). We start by multiplying
by 1 and by noting that ${1=\alpha}_{w}+{\alpha}_{n}$,

$\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}({\alpha}_{w}+{\alpha}_{n})\overline{{\boldsymbol{U}}^{w}_{w}}\right)=0$
(2.24) $\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{w}\overline{{\boldsymbol{U}}^{w}_{w}}\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{n}\overline{{\boldsymbol{U}}^{w}_{w}}\right)=0$
(2.25)

We then add
${\phi}_{f}{\alpha}_{w}{\alpha}_{n}\overline{{\boldsymbol{U}}^{n}_{n}}$ to one
side and subtract it from the other (essentially adding zero to the equation),

$\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\alpha}_{w}\left({\phi}_{f}{\alpha}_{w}\overline{{\boldsymbol{U}}^{w}_{w}}{+\phi}_{f}{\alpha}_{n}\overline{{\boldsymbol{U}}^{n}_{n}}\right)\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{n}\left(\overline{{\boldsymbol{U}}^{w}_{w}}-\overline{{\boldsymbol{U}}^{n}_{n}}\right)\right)=0$
(2.26) $\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\alpha}_{w}{\boldsymbol{U}}_{f}\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{n}\left(\overline{{\boldsymbol{U}}^{w}_{w}}-\overline{{\boldsymbol{U}}^{n}_{n}}\right)\right)=0$
(2.27)

Finally, we define the relative velocity
${\boldsymbol{U}}_{r}=\overline{{\boldsymbol{U}}^{w}_{w}}-\overline{{\boldsymbol{U}}^{n}_{n}}$
in order to obtain the following saturation equation:

$\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\alpha}_{w}{\boldsymbol{U}}_{f}\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{n}{\boldsymbol{U}}_{r}\right)=0$
(2.28)

The definition of the relative velocity is domain dependent (it is defined
differently in porous and solid-free regions, (Carrillo et al.,, 2020)). Its
full analytical derivation and description is provided in Section 2.6.2.

###  2.5 Fluid Momentum Conservation Equation

We start with the classical form of the momentum conservation equation for a
single-phase fluid $i$,

$\frac{\partial{\rho}_{i}{\boldsymbol{U}}_{i}}{\partial
t}+\nabla\cdot\left({\rho}_{i}{\boldsymbol{U}}_{i}{\boldsymbol{U}}_{i}\right)=-\nabla
p_{i}+{\rho}_{i}\boldsymbol{g}+\nabla\cdot{\boldsymbol{S}}_{i}$ (2.29)

where ${\rho}_{i}$ is the density, ${\boldsymbol{U}}_{i}$ is the velocity,
$p_{i}$ is the fluid pressure, $\boldsymbol{g}$ is gravity, and
${\boldsymbol{S}}_{i}=\
{\mu}_{i}\left(\nabla{\boldsymbol{U}}_{i}+{\left(\nabla{\boldsymbol{U}}_{i}\right)}^{T}\right)$
is the viscous stress tensor for a Newtonian fluid. For simplicity, we neglect
the inertial and convective terms for the time being:

$0=-\overline{\nabla
p_{i}}+\overline{{\rho}_{i}\boldsymbol{g}}+\overline{\nabla\cdot{\boldsymbol{S}}_{i}}$
(2.30)

Applying the averaging theorems in a volume containing an additional
immiscible fluid ($j$) and a solid ($s$), we obtain:

$0=-\nabla{\overline{p}}_{i}-\frac{1}{V}\int_{A_{i,j}}{p_{i}{\boldsymbol{n}}_{i,j}dA}-\frac{1}{V}\int_{A_{i,s}}{p_{i}{\boldsymbol{n}}_{i,k}dA}+\
\overline{{\rho}_{i}\boldsymbol{g}}+\nabla\cdot\overline{{\boldsymbol{S}}_{i}}+\\\
\frac{1}{V}\int_{A_{i,j}}{{\boldsymbol{S}}_{i}\cdot{\boldsymbol{n}}_{i,j}dA}+\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{S}}_{i}{\cdot\boldsymbol{n}}_{i,s}dA}$
(2.31)

For improved readability, we now group the terms integrating over the fluid-
fluid interface into a single term ${\boldsymbol{D}}_{i,j}=\
\frac{1}{V}\int_{A_{i,j}}{{\boldsymbol{n}}_{i,j}\cdot(-\boldsymbol{I}p_{i}+{\boldsymbol{S}}_{i})dA}$.
Furthermore, we now convert the volume averages into their intrinsic phase
averages such that:

$0=-\nabla{{({\phi}_{f}\alpha}_{i}p}^{i}_{i})+{\phi}_{f}{\alpha}_{i}{\rho}_{i}\boldsymbol{g}+\nabla\cdot\overline{{\boldsymbol{S}}_{i}}-\frac{1}{V}\int_{A_{i,s}}{p_{i}{\boldsymbol{n}}_{i,k}dA}+\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{S}}_{i}{\cdot\boldsymbol{n}}_{i,s}dA}+{\boldsymbol{D}}_{i,j}$
(2.32)

Before we continue, we note that we can separate the pressure into its
intrinsic phase average and deviation terms (Whitaker,, 2013), such that
$p_{i}=\ {\overline{p}}^{i}_{i}+{\tilde{p}}_{i}$,

$0=-\nabla{{({\phi}_{f}\alpha}_{i}\overline{p}}^{i}_{i})+{\phi}_{f}{\alpha}_{i}{\rho}_{i}\boldsymbol{g}+\nabla\cdot\overline{{\boldsymbol{S}}_{i}}-\frac{1}{V}\int_{A_{i,s}}{{\hat{p}}^{i}_{i}{\boldsymbol{n}}_{i,k}dA}\
\\\ -\frac{1}{V}\int_{A_{i,s}}{{\hat{p}}_{i}{\boldsymbol{n}}_{i,k}dA}\
+\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{S}}_{i}{\cdot\boldsymbol{n}}_{i,s}dA}+{\boldsymbol{D}}_{i,j}$
(2.33)

Assuming that the intrinsic phase averages are constant along the integration
surface and by expanding the terms inside the gradient operators we get the
following,

$0=-{\phi}_{f}\nabla{{(\alpha}_{i}\overline{p}}^{i}_{i})-{{\alpha}_{i}\overline{p}}^{i}_{i}\nabla{\phi}_{f}+{\phi}_{f}{\alpha}_{i}{\rho}_{i}\boldsymbol{g}+\nabla\cdot\overline{{\boldsymbol{S}}_{i}}-\\\
{\overline{p}}^{i}_{i}\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{n}}_{i,k}dA}+\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{n}}_{i,j}\cdot\left(-\boldsymbol{I}{\tilde{p}}_{i}+{\boldsymbol{S}}_{i}\right)dA}+{\boldsymbol{D}}_{i,j}$
(2.34)

We now apply the averaging theorem that relates the surface integral of the
unit normal to the porosity gradient (Equation 2.10),

$0=-{\phi}_{f}\nabla{{(\alpha}_{i}\overline{p}}^{i}_{i})-{{\alpha}_{i}\overline{p}}^{i}_{i}\nabla{\phi}_{f}+{\phi}_{f}{\alpha}_{i}{\rho}_{i}\boldsymbol{g}+\nabla\cdot\overline{{\boldsymbol{S}}_{i}}\\\
+{{\alpha}_{i}\overline{p}}^{i}_{i}\nabla{\phi}_{f}\
+\frac{1}{V}\int_{A_{i,s}}{{\boldsymbol{n}}_{i,j}\cdot\left(-\boldsymbol{I}{\tilde{p}}_{i}+{\boldsymbol{S}}_{i}\right)dA}+{\boldsymbol{D}}_{i,j}$
(2.35)

We then cancel out like-terms and group the remaining integrals into a single
term that represents the total momentum exchange across the solid-fluid
interface
${\boldsymbol{D}}_{i,s}=\frac{1}{V}\int_{A_{i,s}}{\boldsymbol{n}}_{i,j}\cdot$
$\left(-\boldsymbol{I}{\tilde{p}}_{i}+{\boldsymbol{S}}_{i}\right)dA$. The
result is the averaged fluid momentum equation over a volume containing
another immiscible fluid and a solid,

$0=-{\phi}_{f}\nabla{{(\alpha}_{i}\overline{p}}^{i}_{i})+{\phi}_{f}{\alpha}_{i}{\rho}_{i}\boldsymbol{g}+\nabla\cdot\overline{{\boldsymbol{S}}_{i}}\
+{\boldsymbol{D}}_{i,s}+{\boldsymbol{D}}_{i,j}$ (2.36)

Adding the momentum equations for a wetting ($w)$ and a non-wetting ($n$)
fluid we obtain the single-field averaged momentum conservation equation in
terms of ${\boldsymbol{U}}_{f}$ and the single-field pressure
$p=\alpha_{w}\overline{p}^{w}_{w}+\alpha_{n}\overline{p}^{n}_{n}$.

$0=-{\phi}_{f}\nabla
p+{\phi}_{f}{\rho}_{f}\boldsymbol{g}+\nabla\cdot\overline{\boldsymbol{S}}\
+{\boldsymbol{D}}_{w,s}+{\boldsymbol{D}}_{n,s}+{\boldsymbol{D}}_{w,n}+{\boldsymbol{D}}_{n,w}$
(2.37)

where
$\overline{\boldsymbol{S}}={\mu}_{f}(\nabla{\boldsymbol{U}}_{f}+{(\nabla{\boldsymbol{U}}_{f})}^{T})$
is the averaged single-field viscous stress tensor (Appendix D.2), ${\mu}_{f}$
is the arithmetic average of each fluid’s viscosity
${\mu}_{f}={\alpha}_{w}{\mu}_{w}+{\alpha}_{n}{\mu}_{n}$, and ${\rho}_{f}$ is
the arithmetic average of each fluid’s density
${\rho}_{f}={\alpha}_{w}{\rho}_{w}+{\alpha}_{n}{\rho}_{n}$. If necessary, we
can now add the averaged single-field material derivative (Appendix D.3) to
obtain,

$\frac{\partial{\rho}_{f}{\boldsymbol{U}}_{f}}{\partial
t}+\nabla\cdot\left(\frac{{\rho}_{f}}{{\phi}_{f}}{{\boldsymbol{U}}_{f}\boldsymbol{U}}_{f}\right)=-{\phi}_{f}\nabla
p+{\phi}_{f}{\rho}_{f}\boldsymbol{g}+\nabla\cdot\overline{\boldsymbol{S}}+\\\
{\boldsymbol{D}}_{w,s}+{\boldsymbol{D}}_{n,s}+{\boldsymbol{D}}_{w,n}+{\boldsymbol{D}}_{n,w}$
(2.38)

Finally, for this equation to be valid in both a porous medium and within a
solid-free domain, we now need closure for the ${\boldsymbol{D}}_{i,j}$
interaction terms. This was done in Carrillo et al., (2020) and described in
the following section. The results is shown here:

${\boldsymbol{D}}_{w,s}+{\boldsymbol{D}}_{n,s}+{\boldsymbol{D}}_{w,n}+{\boldsymbol{D}}_{n,w}=-{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)+{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{f}{\boldsymbol{F}}_{c,2}$
(2.39)

where $\mu k^{-1}$ is the single-field mobility of the fluids (i.e. a
representation for drag) and ${\boldsymbol{F}}_{c}$ represents the capillary
forces arising from fluid-fluid interfacial dynamics. The analytical
representations of these terms are derived in the next section (Section 2.6)
and summarized in Section 2.13. The result is the complete single-filed fluid
momentum equation for two-phase flow in deformable porous media (i.e. a re-
formulated VOF approach for flow in dynamic porous media):

$\frac{\partial{\rho}_{f}{\boldsymbol{U}}_{f}}{\partial
t}+\nabla\cdot\left(\frac{{\rho}_{f}}{{\phi}_{f}}{\boldsymbol{U}}_{\boldsymbol{f}}{\boldsymbol{U}}_{\boldsymbol{f}}\right)=-{\phi}_{f}\nabla
p+{\phi}_{f}{\rho}_{f}\boldsymbol{g}+\nabla\cdot\overline{\boldsymbol{S}}-\\\
{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)+{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{f}{\boldsymbol{F}}_{c,2}$
(2.40)

The equation presented above tend towards the standard VOF approach in solid-
free regions (where the drag term becomes negligible) and towards the
multiphase Darcy equations in porous regions. The latter can be explained by
the fact that the viscous stress tensor $\nabla\cdot\overline{\boldsymbol{S}}$
becomes negligible under the scale-separation assumption, inertial terms
become negligible under the assumption of low Reynold’s number flow in the
porous medium, and the $\boldsymbol{F}_{c}$ terms are set to fit standard
multiphase Darcy’s law (Whitaker, 1986a, ; Carrillo et al.,, 2020). Therefore,
Eqn. 2.40 $\approx$

$\begin{cases}\frac{\partial{\rho}_{f}{\boldsymbol{U}}_{f}}{\partial
t}+\nabla\cdot\left({\rho}_{f}{{\boldsymbol{U}}_{f}\boldsymbol{U}}_{f}\right)=-\nabla
p+\nabla\cdot\overline{\boldsymbol{S}}+{\rho}_{f}\boldsymbol{g}+{\boldsymbol{F}}_{c,1}&\textnormal{in
solid-free regions}\\\
\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)=-\frac{k}{\mu}\left(\nabla
p-{\rho}_{f}\boldsymbol{g}-{\boldsymbol{F}}_{c,1}-{\boldsymbol{F}}_{c,2}\right)&\textnormal{in
porous regions}\end{cases}$ (2.41)

###  2.6 Closure of the Fluids’ Multiscale Interaction Terms

To form the multiscale momentum equation, we express the sum of the average
shear stress at the fluid-solid and fluid-fluid interfaces as the sum of two
independent terms, an Eulerian drag force $-{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)$ and a
capillary force $\boldsymbol{F}_{c}$:

${\boldsymbol{D}}_{w,s}+{\boldsymbol{D}}_{n,s}+{\boldsymbol{D}}_{w,n}+{\boldsymbol{D}}_{n,w}=-{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)+{\phi}_{f}{\boldsymbol{F}}_{c}$
(2.42)

these so called “multi-scale parameters” reflect sub-grid scale variables and
processes such as the location of the fluid-fluid interface and the
hydrodynamic impact of the porous micro-structure. They have different
meanings and different formulations depending on whether their averaging
volume contains solid material or not. Therefore we will present distinct
derivations of these terms with the purpose of matching the solution of the
two-phase Navier-Stokes equations in solid-free regions and of the two-phase
Darcy equations in porous regions.

####  2.6.1 Matching the Volume-of-Fluid Model in Solid-Free Regions

In CFD, the Volume of Fluid (VOF) method (Hirt & Nichols,, 1981) is a standard
approach to track the interface movement of two immiscible fluids in a fixed
Eulerian grid. In it, a phase indicator representing the volume of fluid in
each grid block is used to track the distribution of the fluid phases in a
computational domain as illustrated in the upper part of Figure 2.2B. This
phase indicator has the same form as the saturation field $\alpha_{w}$ defined
before, where $\alpha_{w}=1$ in cells saturated by the wetting phase and
$\alpha_{w}=0$ in cells that contain the non-wetting phase only. Note that
$0<\alpha_{w}<1$ in cells that contain the immiscible interface between both
fluids.

In standard VOF approaches, the computational cells do not contain solid
($\phi_{f}=1$). In these cases, the relative velocity $\boldsymbol{U}_{r}$ is
used as a compression term to force the fluid-fluid interface to be as sharp
as possible (Rusche,, 2002). This compression velocity acts in the direction
normal to the interface. In the VOF framework, the normal to the fluid-fluid
interface is computed using the gradient of the saturation. Rusche, (2002)
proposes a relative velocity oriented in the direction normal to interface
with a value based on the maximum magnitude of $\boldsymbol{U}_{f}$:

$\boldsymbol{U}_{r}=C_{\alpha}\max\left(\left|\boldsymbol{U}_{f}\right|\right)\overline{\boldsymbol{n}_{w,n}},$
(2.43)

where $C_{\alpha}$ is a model parameter used to control the compression of the
interface and $\overline{\boldsymbol{n}_{w,n}}$ is the mean normal vector. For
low values of $C_{\alpha}$, the interface diffuses. For higher values, the
interface is sharper, but excessive values are known to introduce parasitic
velocities and lead to unphysical solutions. In practice, $C_{\alpha}$ is
often chosen between 0 and 4. The mean normal vector $\boldsymbol{n}_{w,n}$ is
computed by using the gradient of the phase indicator function $\alpha_{w}$.
The relation between these two vectors can be obtained by applying Eq. 2.6 to
the liquid phase indicator function $\mathds{1}_{l}$ (a function equal to 1 in
$V_{w}$ and 0 elsewhere) in solid-free regions such that,

$\nabla\alpha_{w}=-\frac{1}{V}\int_{A_{w,n}}\mathds{1}_{l}\boldsymbol{n}_{w,n}dA.$
(2.44)

Therefore,

$\overline{\boldsymbol{n}_{w,n}}=-\frac{\nabla\alpha_{w}}{\left|\nabla\alpha_{w}\right|},$
(2.45)

is a unit vector defined at the cell centers that describes the mean normal to
the fluid-fluid interface in a control volume.

Another consequence of the absence of solid in the VOF equations is that the
forces describing the shear stresses of the fluids onto the solid surface are
null, hence $\boldsymbol{D}_{w,s}=\boldsymbol{D}_{n,s}=0$. Therefore, the
Darcy term in the momentum equation vanishes:

${\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)=0.$
(2.46)

The integration of the shear boundary condition at the fluid-fluid interface,
yields a relationship between the mutual shear between the two fluids and the
surface integral of the surface tension effects:

$\boldsymbol{D}_{w,n}+\boldsymbol{D}_{n,w}=\phi_{f}\boldsymbol{F}_{c}=\frac{1}{V}\int_{A_{lg}}\boldsymbol{n}_{lg}\cdot\gamma\kappa
dA.$ (2.47)

This equation cannot be used directly, because the terms under the volume
integral require knowledge of the location and curvature of the fluid-fluid
interface within a grid block. This information is unknown in a grid-based
formulation for which all the physical variables and forces are averaged on
control volumes. In the VOF method, the curvature of the interface $\kappa$ is
approximated by a mean interface curvature $\overline{\kappa}$. Brackbill et
al., (1992) assumes that the mean curvature of the interface can be
approximated by calculating the divergence of the mean normal vector,
$\overline{\kappa}=\nabla\cdot\overline{\boldsymbol{n}_{w,n}}$. Because
$\overline{\kappa}$ and $\gamma$ are constant within a control volume, they
can be extracted from the integral in Eq. 2.47 to obtain (after applying Eq.
2.44) the so-called Continuum Surface Force (CSF) formulation (Brackbill et
al.,, 1992):

$\boldsymbol{F}_{c}=\phi_{f}^{-1}\gamma\nabla\cdot\left(\frac{\nabla\alpha_{w}}{\left|\nabla\alpha_{w}\right|}\right)\nabla\alpha_{w}.$
(2.48)

####  2.6.2 Matching the Standard Two-Phase Darcy Model in Porous Regions

In this subsection, we recall the formulation of the standard two-phase Darcy
model that is classically used to describe two-phase flow in porous media at
the continuum scale. The model can be derived by applying the volume averaging
operators on a REV of the porous structure (Whitaker, 1986b, ). Unlike the
present micro-continuum model, the two-phase Darcy model is a two-field model,
meaning that instead of one velocity field describing the flow, there are two
velocities ($\overline{\boldsymbol{U}_{i}}$ with $i=w,n$) with separate
pressure fields ($\overline{p_{i}}$ with $i=w,n$).

The incompressible, immiscible two-phase Darcy model in static porous media
consists of a saturation equation for the wetting phase,

$\frac{\partial\phi\alpha_{w}}{\partial
t}+\nabla\cdot\overline{\boldsymbol{U}_{w}}=0,$ (2.49)

a mass balance equation,

$\nabla\cdot\overline{\boldsymbol{U}_{i}}=0,$ (2.50)

and two momentum balance equations, one for each phase,

$\displaystyle\overline{\boldsymbol{U}_{i}}=\phi_{f}\alpha_{i}\overline{\boldsymbol{U}_{i}^{i}}$
$\displaystyle=-\frac{k_{0}k_{r,i}}{\mu_{i}}\left(\nabla\overline{p_{i}}^{i}-\rho_{i}\boldsymbol{g}\right),\quad
i=w,n,$
$\displaystyle=-M_{i}\left(\nabla\overline{p_{i}}^{i}-\rho_{i}\boldsymbol{g}\right),\quad
i=w,n.$ (2.51)

These can also be written as,

$0=-\nabla\overline{p_{i}}^{i}+\rho_{i}\boldsymbol{g}-M_{i}^{-1}\overline{\boldsymbol{U}_{i}},\quad
i=w,n,$ (2.52)

where $k_{0}$ is the absolute permeability of the porous structure, $k_{r,w}$
and $k_{r,n}$ are the relative permeabilities with respect to each fluid
(classically represented here as functions of wetting fluid saturation; more
complex formulations exist that account for viscous coupling between the two
fluids or for the Klinkenberg effect in the non-wetting (often gaseous) phase
(Picchi & Battiato,, 2018)), and $M_{i}=\frac{k_{0}k_{r,i}}{\mu_{i}}$ are the
fluid mobilities. Equations 2.52 arise from further simplification of the
volume averaged Stokes equations, Eq. 2.36, where the drag forces are combined
and described as a single Darcy term. Moreover, Whitaker, 1986a  showed that
the viscous dissipative term, $\nabla\cdot\overline{\boldsymbol{S}}_{i}$ is
negligible in comparison to the drag forces whenever the system’s macroscopic
length scale is significantly larger than the length scale of that averaging
volume. This feature is a fundamental aspect of the multiscale framework
described above because it means that even though the viscous dissipative term
is retained in the single-field momentum equation, it naturally vanishes when
the computational cells contain solid content. This allows the continuity of
stresses between porous and solid-free domains (Neale & Nader,, 1974).

Because it involves four equations and five unknown variables, the two-phase
Darcy model is complemented by the definition of macroscopic capillary
pressure $p_{c}$, which provides an additional relationship between the two
averaged pressure fields:

$p_{c}\left(\alpha_{w}\right)=\left(\overline{p_{g}}^{g}-\overline{p_{l}}^{l}\right).$
(2.53)

This equation has been theoretically derived through homogenization techniques
and is assumed to be independent of sub-volume fluid properties. (Whitaker,
1986b, ; Torres,, 1987).

As the two-phase Darcy model explicitly represents the two phase-averaged
velocities, it can be used to derive an expression for the relative velocity
$\boldsymbol{U}_{r}$ in the porous region. Before going through the
derivation, we note that the application of the gradient operator to the
definition of the single-field pressure $p$, along with the definition of
capillary pressure, Eq. 2.53, results in:

$\displaystyle\nabla\overline{p_{w}}^{w}$ $\displaystyle=\nabla
p-\nabla\left(\alpha_{n}p_{c}\right),$
$\displaystyle\nabla\overline{p_{n}}^{n}$ $\displaystyle=\nabla
p+\nabla\left(\alpha_{w}p_{c}\right).$ (2.54)

We can then obtain an analytical form for $\boldsymbol{U}_{r}$ by using the
multi-phase Darcy equations presented above (Eqn. 2.51):

$\displaystyle\boldsymbol{U}_{r}$
$\displaystyle=\left(\overline{\boldsymbol{U}_{w}^{w}}-\overline{\boldsymbol{U}_{n}^{n}}\right),$
$\displaystyle=-\frac{M_{w}}{\phi_{f}\alpha_{w}}\left(\nabla\overline{p_{w}^{w}}-\rho_{w}\boldsymbol{g}\right)+\frac{M_{n}}{\phi\alpha_{n}}\left(\nabla\overline{p_{n}^{n}}-\rho_{n}\boldsymbol{g}\right),$
$\displaystyle=\phi^{-1}_{f}\left[-\frac{M_{w}}{\alpha_{w}}\nabla\overline{p_{w}^{w}}+\frac{M_{n}}{\alpha_{n}}\nabla\overline{p_{n}^{n}}+\left(\rho_{w}\frac{M_{w}}{\alpha_{w}}-\rho_{n}\frac{M_{n}}{\alpha_{n}}\right)\boldsymbol{g}\right],$
$\displaystyle=\phi^{-1}_{f}\left[\begin{array}[]{c}-\left(\frac{M_{w}}{\alpha_{w}}-\frac{M_{n}}{\alpha_{n}}\right)\nabla
p+\left(\rho_{w}\frac{M_{w}}{\alpha_{w}}-\rho_{n}\frac{M_{n}}{\alpha_{n}}\right)\boldsymbol{g}+\\\
\frac{M_{w}}{\alpha_{w}}\nabla\left(\alpha_{n}p_{c}\right)+\frac{M_{n}}{\alpha_{n}}\nabla\left(\alpha_{w}p_{c}\right)\end{array}\right],$
(2.57)
$\displaystyle={\phi_{f}}^{-1}\left[\begin{array}[]{c}-\left(M_{w}{\alpha}^{-1}_{w}-M_{n}{\alpha}^{-1}_{n}\right)\nabla
p+\left({\rho}_{w}M_{w}{\alpha}^{-1}_{w}-{\rho}_{w}M_{n}{\alpha}^{-1}_{n}\right)\boldsymbol{g}+\\\
\left(M_{w}{\alpha}_{n}{\alpha}^{-1}_{w}+M_{n}{\alpha}_{w}{\alpha}^{-1}_{n}\right)\nabla
p_{c}-\left(M_{w}{\alpha}^{-1}_{w}-M_{n}{\alpha}^{-1}_{n}\right)p_{c}\nabla{\alpha}_{w}\end{array}\right].$
(2.60)

A two-phase Darcy model for the single-field velocity $\boldsymbol{U}_{f}$ is
then formed to derive the continuum scale formulation of the drag force $\mu
k^{-1}\boldsymbol{U}_{f}$ and capillary force $\boldsymbol{F}_{c}$. This is
achieved by summing both phase velocities, Eq. 2.51, and using the pressure
gradient relationship, Eq. 2.54. We obtain:

$\displaystyle\boldsymbol{U}_{f}$
$\displaystyle=\overline{\boldsymbol{U}_{w}}+\overline{\boldsymbol{U}_{n}},$
$\displaystyle=-M_{n}\nabla\overline{p_{n}^{n}}-M_{w}\nabla\overline{p_{w}^{w}}+\left(\rho_{n}M_{n}+\rho_{w}M_{w}\right)\boldsymbol{g},$
(2.61) $\displaystyle=-\left(M_{n}+M_{w}\right)\nabla
p+\left(\rho_{n}M_{n}+\rho_{w}M_{w}\right)\boldsymbol{g}+\left[M_{w}\nabla\left(\alpha_{n}p_{c}\right)-M_{n}\nabla\left(\alpha_{w}p_{c}\right)\right],$

The previous equation can be recast into:

$0=-\nabla
p+\rho_{M}\boldsymbol{g}-M^{-1}\boldsymbol{U}_{f}+M^{-1}\left[M_{w}\nabla\left(\alpha_{n}p_{c}\right)-M_{n}\nabla\left(\alpha_{w}p_{c}\right)\right],$
(2.62)

where $M=M_{w}+M_{n}$ is the total mobility and
$\rho_{M}=\left(\rho_{w}M_{w}+\rho_{n}M_{n}\right)/\left(M_{w}+M_{n}\right)$
is a mobility-weighted average fluid density. This single-field two-phase
Darcy equation matches the two-phase micro-continuum momentum equation, Eq.
2.37, if this equation’s drag coefficient and the capillary force equal

$\mu
k^{-1}=M^{-1}=k_{0}^{-1}\left(\frac{\mu_{w}}{k_{r,w}}+\frac{\mu_{n}}{k_{r,n}}\right)^{-1},$
(2.63)

and

$\displaystyle\boldsymbol{F}_{c}$
$\displaystyle=M^{-1}\left[M_{w}\nabla\left(\alpha_{n}p_{c}\right)-M_{n}\nabla\left(\alpha_{w}p_{c}\right)\right],$
$\displaystyle=\left[M^{-1}\left(M_{w}\alpha_{n}-M_{n}\alpha_{w}\right)\left(\frac{\partial
p_{c}}{\partial\alpha_{w}}\right)-p_{c}\right]\nabla\alpha_{w},$ (2.64)

respectively. Here, the single-field relative permeability, Eq. 2.63, is a
harmonic average of the two-phase mobilities, in agreement with the proposal
of Wang & Beckermann, (1993) and Soulaine et al., (2019).

Finally, we note that in Eq. 2.62, the single-field fluid density $\rho_{M}$
in the buoyant term is a weighted average based on the fluid mobilities, or
more exactly, the fractional flow functions, $M_{i}M^{-1}$. This is a classic
concept in multiphase flow in porous media. As shown in Carrillo et al.,
(2020) a strictly equivalent solution can be derived where $\rho_{M}$ is
replaced by $\rho_{f}$ in Eq. 2.62 and the capillary force expression is
replaced by:

$\boldsymbol{F}_{c}=\boldsymbol{F}_{c,1}+\boldsymbol{F}_{c,2}=M^{-1}\left(M_{w}\alpha_{n}-M_{n}\alpha_{w}\right)[(\rho_{w}-\rho_{n})\boldsymbol{g}+\nabla
p_{c}]-p_{c}\nabla\alpha_{w}$ (2.65)

this is the expression shown in Sections 2.5 and 2.13.

###  2.7 Momentum Conservation Equation for a Linear Elastic Solid

We start with the differential form of the force balance for a solid body in
terms of the displacement vector ${\boldsymbol{d}}_{s}$. Fundamentally, it
consists of an elastic stress tensor ($\boldsymbol{\sigma})$ balanced by a
Terzaghi effective stress tensor ($\boldsymbol{\tau})$ (Carrillo & Bourg,,
2019; Jasak & Weller,, 2000).

$\frac{{\partial}^{2}\left({\rho}_{s}{\boldsymbol{d}}_{s}\right)}{\partial
t^{2}}-\nabla\cdot\boldsymbol{\sigma}=\nabla\cdot\boldsymbol{\tau}+{\rho}_{s}\boldsymbol{g}$
(2.66)

where $\boldsymbol{\tau}$ is a function of the confining pressure
$(P_{conf})$, the fluid pressure $(p)$, and the swelling pressure
$(p_{swell})$ such that: $\boldsymbol{\tau}=\
{\boldsymbol{P}}_{conf}-\boldsymbol{I}p-\boldsymbol{I}p_{swell}$. Applying the
volume averaging operators for a solid ($s$) in contact with a wetting ($w$)
and non-wetting fluid ($n$) we obtain:

$\overline{\frac{{\partial}^{2}\left({\rho}_{s}{\boldsymbol{d}}_{s}\right)}{\partial
t^{2}}}-\overline{\nabla\cdot\boldsymbol{\sigma}}=\overline{\nabla\cdot\boldsymbol{\tau}}\
+\overline{{\rho}_{s}\boldsymbol{g}}$ (2.67)

Expanding the first term from the left we get:

$\begin{split}=\frac{{\partial}^{2}\left(\overline{{\rho}_{s}{\boldsymbol{d}}_{s}}\right)}{\partial
t^{2}}-\frac{\partial}{\partial
t}\frac{1}{V}\int_{A_{s,w}}{{{\rho}_{s}{\boldsymbol{d}}_{s}\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,w}}{\frac{\partial\left({\rho}_{s}{\boldsymbol{d}}_{s}\right)}{\partial
t}{\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,w}dA}\\\
-\frac{\partial}{\partial
t}\frac{1}{V}\int_{A_{s,n}}{{{\rho}_{s}{\boldsymbol{d}}_{s}\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,n}dA}-\frac{1}{V}\int_{A_{s,n}}{\frac{\partial\left({\rho}_{s}{\boldsymbol{d}}_{s}\right)}{\partial
t}{\boldsymbol{U}}_{s}\cdot{\boldsymbol{n}}_{s,n}dA}\end{split}$ (2.68)

Given that linear elasticity theory only deals with infinitesimal
deformations, we can safely assume that the velocity of the solid-fluid
interface $({\boldsymbol{U}}_{s})$ is very close to zero. This approximation
yields:

$\overline{\frac{{\partial}^{2}\left({\rho}_{s}{\boldsymbol{d}}_{s}\right)}{\partial
t^{2}}}=\frac{{\partial}^{2}\left(\overline{{\rho}_{s}{\boldsymbol{d}}_{s}}\right)}{\partial
t^{2}}=\frac{{\partial}^{2}\left({\phi}_{s}{\rho}_{s}\overline{{\boldsymbol{d}}^{s}_{s}}\right)}{\partial
t^{2}}$ (2.69)

We now continue by expanding the stress terms,

$\begin{split}\frac{{\partial}^{2}\left({\phi}_{s}{\rho}_{s}\overline{{\boldsymbol{d}}^{s}_{s}}\right)}{\partial
t^{2}}-\nabla\cdot\overline{\boldsymbol{\sigma}}-\frac{1}{V}\int_{A_{s,n}}{\boldsymbol{\sigma}{\cdot\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,n}}{\boldsymbol{\sigma}{\cdot\boldsymbol{n}}_{s,n}dA}\\\
-\nabla\cdot\overline{\boldsymbol{\tau}}-\frac{1}{V}\int_{A_{s,w}}{\boldsymbol{\tau}{\cdot\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,n}}{\boldsymbol{\tau}{\cdot\boldsymbol{n}}_{s,n}dA}=0\end{split}$
(2.70)

Just as before, we can write the Terzaghi stress tensor (i.e. the fluid
pressure) as the sum of its intrinsic phase average and its sub-volume
deviations. Furthermore, we also expand the averaged Terzaghi term into its
phase averaged form,

$\begin{split}\frac{{\partial}^{2}\left({\phi}_{s}{\rho}_{s}\overline{{\boldsymbol{d}}^{s}_{s}}\right)}{\partial
t^{2}}-\nabla\cdot\overline{\boldsymbol{\sigma}}-\frac{1}{V}\int_{A_{s,w}}{{\boldsymbol{\sigma}\cdot\boldsymbol{n}}_{s,i}dA}-\frac{1}{V}\int_{A_{s,n}}{\boldsymbol{\sigma}{\cdot\boldsymbol{n}}_{s,n}dA}-\nabla\cdot{(\phi}_{s}{\overline{\boldsymbol{\tau}}}^{s})\\\
-\frac{1}{V}\int_{A_{s,w}}{\left(\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}+\widetilde{\boldsymbol{\tau}}\right){\cdot\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,n}}{(\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}+\widetilde{\boldsymbol{\tau}}){\cdot\boldsymbol{n}}_{s,n}dA}=0\end{split}$
(2.71)

Assuming that the intrinsic phase averages of the stress tensors are constant
along the integration surfaces and using the geometric relation between said
surfaces and the porosity gradient (Equation 2.10), we obtain:

$\frac{{\partial}^{2}\left({\phi}_{s}{\rho}_{s}\overline{{\boldsymbol{d}}^{s}_{s}}\right)}{\partial
t^{2}}-\nabla\cdot\overline{\boldsymbol{\sigma}}-\nabla\cdot{(\phi}_{s}\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}})+{\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}\cdot\nabla\phi}_{s}-\\\
\frac{1}{V}\int_{A_{s,w}}{{\left(\boldsymbol{\sigma}+\widetilde{\boldsymbol{\tau}}\right)\cdot\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,n}}{(\boldsymbol{\sigma}+\widetilde{\boldsymbol{\tau}}){\cdot\boldsymbol{n}}_{s,n}dA}=0$
(2.72)

This expression can be simplified into the following by the use of the chain
rule,

$\frac{{\partial}^{2}\left({\phi}_{s}{\rho}_{s}\overline{{\boldsymbol{d}}^{s}_{s}}\right)}{\partial
t^{2}}-\nabla\cdot\overline{\boldsymbol{\sigma}}-{\phi}_{s}\nabla\cdot\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}-\\\
\frac{1}{V}\int_{A_{s,w}}{{\left(\boldsymbol{\sigma}+\widetilde{\boldsymbol{\tau}}\right)\cdot\boldsymbol{n}}_{s,w}dA}-\frac{1}{V}\int_{A_{s,n}}{(\boldsymbol{\sigma}+\widetilde{\boldsymbol{\tau}}){\cdot\boldsymbol{n}}_{s,n}dA}=0$
(2.73)

Finally we represent the integral terms as two separate interaction terms
defined as
${\boldsymbol{B}}_{s,i}=\frac{1}{V}\int_{A_{s,i}}{{\left(\boldsymbol{\sigma}+\widetilde{\boldsymbol{\tau}}\right)\cdot\boldsymbol{n}}_{s,i}dA}$:

$\frac{{\partial}^{2}\left({\phi}_{s}{\rho}_{s}\overline{{\boldsymbol{d}}^{s}_{s}}\right)}{\partial
t^{2}}-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}+{{\phi}_{s}\rho}_{s}\boldsymbol{g}+{\boldsymbol{B}}_{s,w}+{\boldsymbol{B}}_{s,n}$
(2.74)

In the case of a linear elastic solid the stress tensor is a function of the
two lame coefficients $({\mu}_{s},\ {\lambda}_{s})$ and takes the following
form:

$\overline{\boldsymbol{\sigma}}={\phi}_{s}\mu_{s}\nabla\overline{{\boldsymbol{d}}^{s}_{s}}+{{\phi}_{s}{\mu}_{s}\left(\nabla\overline{{\boldsymbol{d}}^{s}_{s}}\right)}^{T}+{\phi}_{s}\lambda_{s}tr(\nabla\overline{{\boldsymbol{d}}^{s}_{s}})\boldsymbol{I}$
(2.75)

###  2.8 Momentum Conservation Equation for a Viscoplastic Solid

The elementary momentum conservation equation for a plastic solid can be
written as:

$\frac{\partial{\rho}_{s}{\boldsymbol{U}}_{s}}{\partial
t}+\nabla\cdot\left({\rho}_{s}{\boldsymbol{U}}_{s}{\boldsymbol{U}}_{s}\right)-\nabla\cdot\boldsymbol{\sigma}=\nabla\cdot\boldsymbol{\tau}+{\rho}_{s}\boldsymbol{g}$
(2.76)

where $\boldsymbol{\sigma}$ represents the plastic viscous stress tensor and
$\boldsymbol{\tau}$ is the Terzaghi effective stress tensor. Note how in this
case, the solid is essentially a viscous fluid, albeit a non-Newtonian one.
Fortunately, the averaging procedure for each term within this equation has
already been shown in Sections 2.5 and 2.7. Putting it all together, the
averaged momentum equation for a plastic solid is as follows:

$\frac{\partial{{\phi}_{s}\rho}_{s}{\overline{\boldsymbol{U}}}^{s}_{s}}{\partial
t}+\nabla\cdot\left({\phi}_{s}{\rho}_{s}{\overline{\boldsymbol{U}}}^{s}_{s}{\overline{\boldsymbol{U}}}^{s}_{s}\right)-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot{\overline{\boldsymbol{\tau}}}^{s}+{{\phi}_{s}\rho}_{s}\boldsymbol{g}+{\boldsymbol{B}}_{s,w}+{\boldsymbol{B}}_{s,n}$
(2.77)

For a viscoplastic solid the stress tensor is a function of the effective
plastic viscosity $({\mu}^{eff}_{s})$ and takes the following form:

$\overline{\boldsymbol{\sigma}}={\phi}_{s}\ {\mu}^{eff}_{s}\left(\
\nabla{\overline{\boldsymbol{U}}}^{s}_{s}+{\
\left(\nabla{\overline{\boldsymbol{U}}}^{s}_{s}\right)}^{T}-\frac{2}{3}\nabla\cdot\left({\overline{\boldsymbol{U}}}^{s}_{s}\boldsymbol{I}\right)\right)$
(2.78)

###  2.9 General Solid Momentum Conservation Equation

If we assume low solid velocities ($\mathrm{Re}_{s}\lesssim 1)$ we can
reasonably neglect the inertial terms and time derivatives in both the elastic
and plastic solid momentum equations (Eqns. 2.77 and 2.74 in order to arrive
to a general equation for solid mechanics:

$-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot{\overline{\boldsymbol{\tau}}}^{s}+{{\phi}_{s}\rho}_{s}\boldsymbol{g}+{\boldsymbol{B}}_{s,w}+{\boldsymbol{B}}_{s,n}$
(2.79)

where ${\overline{\boldsymbol{\sigma}}}^{s}_{s}$ is either the solid’s elastic
or viscoplastic stress tensor. The only thing left to do is to close the
${\boldsymbol{B}}_{s,i}$ interaction terms. This was done in Carrillo et al.,
(2020) and shown in the next section (Section 2.10), where:

${\boldsymbol{B}}_{s,w}+{\boldsymbol{B}}_{s,n}={\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)-{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{s}{\boldsymbol{F}}_{c,2}$
(2.80)

here, $\mu k^{-1}$ is the single-field mobility of the fluids (i.e. a
representation for drag) and ${\boldsymbol{F}}_{c,i}$ represents the capillary
forces arising from fluid-fluid interfacial dynamics within the porous medium.
Putting everything together we now obtain the averaged momentum equation for a
solid containing two immiscible phases:

$-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot{\overline{\boldsymbol{\tau}}}^{s}+{{\phi}_{s}\rho}_{s}\boldsymbol{g}+{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)-{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{s}{\boldsymbol{F}}_{c,2}$
(2.81)

###  2.10 Closure of the Solid’s Multiscale Interaction Terms

Just as we did for the fluid equations in Section 2.6, we will assume that the
sum of the averaged stresses at the solid-fluid interface can be expressed as
the sum of two independent terms: a drag force that captures shear-induced
momentum exchange $\left({\boldsymbol{B}}_{drag}\right)$ and a capillary force
originating from capillary pressure jumps across the integrated solid surfaces
within the porous media $\left({\boldsymbol{B}}_{cap}\right)$.

$\displaystyle{\boldsymbol{B}}_{drag}+{\boldsymbol{B}}_{cap}={\boldsymbol{B}}_{s,w}+{\boldsymbol{B}}_{s,n}$
(2.82)

We now seek closure of these two coupling terms. By conservation of momentum,
we know that any drag-induced momentum lost by the fluid must be gained by the
solid. Thus, by Eqn. 2.40,

$\displaystyle{\boldsymbol{B}}_{drag}={\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)$
(2.83)

Closure of the capillarity-induced interaction term ${\boldsymbol{B}}_{cap}$
is obtained by adding the solid and fluid momentum equations (Eqns. 2.40 and
2.79) within the porous medium at low Reynold numbers and low permeability,
which yields:

$0=-{\phi}_{f}\nabla p+{\phi}_{f}{\rho}_{f}\boldsymbol{g}-{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)+{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{f}{\boldsymbol{F}}_{c,2}\
$ $+$
$-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot{\overline{\boldsymbol{\tau}}}^{s}+{{\phi}_{s}\rho}_{s}\boldsymbol{g}+{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)+{\boldsymbol{B}}_{cap}$
$=$
$-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}-{\phi}_{f}\nabla
p+\left({{\phi}_{s}\rho}_{s}+{\phi}_{f}{\rho}_{f}\right)\boldsymbol{g}+{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{f}{\boldsymbol{F}}_{c,2}+{\boldsymbol{B}}_{cap}\boldsymbol{\
}$ (2.84)

In multiphase porous systems with incompressible grains and no swelling
pressure (i.e.
$\nabla\cdot\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}=-\nabla p)$, Biot
Theory states that $\nabla\cdot\overline{\boldsymbol{\sigma}}=\nabla
p-{\rho}^{*}\boldsymbol{g}+p_{c}\nabla{\alpha}_{w}$, where
${\rho}^{*}=({{\phi}_{s}\rho}_{s}+{\phi}_{f}{\rho}_{f})$ (Jha & Juanes,, 2014;
Kim et al.,, 2013). This expression is satisfied by the previous equation in
the absence of capillary forces (Carrillo & Bourg,, 2019). In multiphase
systems, however, it imposes the following equality,

$\displaystyle{\boldsymbol{B}}_{cap}=-({\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{f}{\boldsymbol{F}}_{c,2}+p_{c}\nabla{\alpha}_{w})$
(2.85)

Given that
${\boldsymbol{F}}_{c,1}=M^{-1}\left(M_{w}{\alpha}_{n}-M_{n}{\alpha}_{w}\right)\left(\nabla
p_{c}+\left({\rho}_{w}-{\rho}_{n}\right)\boldsymbol{g}\right)$ and that
${\boldsymbol{F}}_{c,2}=-p_{c}\nabla{\alpha}_{w}$ (see Eqns. 2.94 and 2.95),
the previous equation can be rearranged to obtain:

$\displaystyle{\boldsymbol{B}}_{cap}=-{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{s}{\boldsymbol{F}}_{c,2}$
(2.86)

Equation 2.86 gives closure to the last coupling parameter and marks the end
of this derivation. The result is the solid momentum conservation equation
shown in Sections 2.9 and 2.13.

###  2.11 Interfacial Conditions between Solid-Free Regions and Porous Regions

One of the most important features allowed by the equations presented above is
the existence of an interface between solid-free and microporous domains.
Although the creation of a rigorous un-averaged description of this interface
is still an open question, we approximate a solution to it by guaranteeing its
necessary components within our fluid and solid averaged equations.

An accurate description of fluid behavior at the interface requires three
components: 1) mass conservation across the interface, 2) continuity of
stresses across the interface, and 3) an interfacial wettability condition.
Components 1 and 2 are intrinsically fulfilled by our solver due to its
single-field formulation for velocity and pressure within the fluid
conservation equations (Eqns. 2.88 and 2.90). As shown in Neale & Nader,
(1974) and Carrillo & Bourg, (2019) these two components are necessary and
sufficient to model single-phase flow within a multiscale system. Furthermore,
these conditions have also been used for closure when modelling multiphase
flow in moving porous media (Lacis et al.,, 2017; Zampogna et al.,, 2019;
Carrillo et al.,, 2020). The required wettability condition at the porous
interface (Component 3) is included in our model through the implementation of
a penalized contact angle condition (Eqn. 2.96) following the steps outlined
in Horgue et al., (2014) and Carrillo et al., (2020).

The complementary solid conditions at the porous interface are very similar:
1) solid mass conservation across the interface, 2) continuity of fluid-
induced stresses across the interface, and 3) a discontinuity of solid
stresses at the interface. Just as before, the first two conditions are
intrinsically fulfilled through the use of a single set of mass and momentum
conservation equations across both domains and have also been used as closure
conditions in previous studies (Lacis et al.,, 2017; Zampogna et al.,, 2019).
The third condition is enforced by the use of volume-averaged solid rheology
models that tend towards infinitely deformable materials in solid-free
regions, as shown in Carrillo & Bourg, (2019). When volume-averaged, the
behavior of the solid’s stress tensor is domain dependent (i.e. solid fraction
dependent). Thus, in solid regions, the elasticity and viscosity of the porous
medium is determined by standard averaged rheological properties (the elastic
and viscoplastic moduli). Contrastingly, in solid-free regions, the solid
fraction tends to zero and, as such, said properties do as well. The result is
a stress-free “ghost” solid that does not apply resistance to the porous
region, creating the required stress discontinuity at the porous interface.

Although necessary, these conditions represent but an approximation to the
complete description of fluid and solid mechanics at the porous interface.
However, to the best of our knowledge, there does not exist an alternative set
of boundary conditions that can or have been used to model multiphase flow in
multiscale porous media.

###  2.12 Advection-Diffusion Equation for Tracer Particles

The derivation of the following advection-diffusion equation for tracer
particles flowing within a single fluid follows the same averaging procedure
as all the previous mass conservation equations.

$\frac{\partial\left({\phi}_{f}{\overline{C}}_{i,x}^{i}\right)}{\partial
t}+\nabla\left({\overline{C}}_{i,x}^{i}{\boldsymbol{U}}_{f}\right)-\nabla\cdot\left({\phi}_{f\
}D_{\mathrm{eff,}x}\nabla{\overline{C}}_{i,x}^{i}\right)=0$ (2.87)

where ${\overline{C}}_{i,x}^{i}$ and $D_{\mathrm{eff,}x}$ are the molar
concentration and effective diffusion coefficient of a given species x,
respectively. A reaction term can readily be added to Eqn. 2.87 to describe
reactive solutes (Soulaine et al.,, 2017). The expansion of Eqn. 2.87 into
multiphase flow is beyond the scope of our investigation, but can be found in
Maes & Soulaine, (2019).

###  2.13 Summary of Equations and Multiscale Parameters

The final set of conservation equations used in this framework now follows.
The combination of these solid and fluid conservation equations leads to a
model that tends towards multiphase Navier-Stokes in solid-free regions and
towards Biot Theory in porous regions (see Appendix D.4).

$\frac{\partial{\phi}_{f}}{\partial t}+\nabla\cdot{\boldsymbol{U}}_{f}=0$
(2.88)

$\frac{\partial{\phi}_{f}{\alpha}_{w}}{\partial
t}+\nabla\cdot\left({\alpha}_{w}{\boldsymbol{U}}_{f}\right)+\nabla\cdot\left({\phi}_{f}{\alpha}_{w}{\alpha}_{n}{\boldsymbol{U}}_{r}\right)=0$
(2.89)

$\frac{\partial{\rho}_{f}{\boldsymbol{U}}_{f}}{\partial
t}+\nabla\cdot\left(\frac{{\rho}_{f}}{{\phi}_{f}}{\boldsymbol{U}}_{f}{\boldsymbol{U}}_{f}\right)=-{\phi}_{f}\nabla
p+{\phi}_{f}{\rho}_{f}\boldsymbol{g}+\nabla\cdot\overline{\boldsymbol{S}}-\\\
{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)+{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{f}{\boldsymbol{F}}_{c,2}$
(2.90)

$\frac{\partial{\phi}_{s}}{\partial
t}+\nabla\cdot\left({\phi}_{s}{\overline{\boldsymbol{U}}}^{s}_{s}\right)=0$
(2.91)

$-\nabla\cdot\overline{\boldsymbol{\sigma}}={\phi}_{s}\nabla\cdot\overline{{\boldsymbol{\tau}}^{\boldsymbol{s}}}+{{\phi}_{s}\rho}_{s}\boldsymbol{g}+{\phi}_{f}\mu
k^{-1}\left({\boldsymbol{U}}_{f}-{\overline{\boldsymbol{U}}}_{s}\right)-{\phi}_{f}{\boldsymbol{F}}_{c,1}+{\phi}_{s}{\boldsymbol{F}}_{c,2}$
(2.92)

The closed-form expressions of the multi-scale parameters $\mu
k^{-1}\mathrm{,\ }\boldsymbol{U}_{r}\mathrm{\ ,and}\ {\boldsymbol{F}}_{c,i}$,
which are defined differently in each region now follow. A full derivation and
discussion of these parameters can be found in Sections 2.6 and 2.10 (Carrillo
et al.,, 2020).

$\mu k^{-1}=\begin{cases}0&\textnormal{in solid-free regions}\\\
k^{-1}_{0}{\left(\frac{k_{r,w}}{{\mu}_{w}}+\frac{k_{r,n}}{{\mu}_{n}}\right)}^{-1}&\textnormal{in
porous regions}\end{cases}$ (2.93)

${\boldsymbol{F}}_{c,1}=\begin{cases}-\frac{\gamma}{{\phi}_{f}}\nabla\cdot\left({\boldsymbol{n}}_{w,n}\right)\nabla{\alpha}_{w}&\textnormal{in
solid-free regions}\\\
M^{-1}\left(M_{w}{\alpha}_{n}-M_{n}{\alpha}_{w}\right)\left(\nabla
p_{c}+\left({\rho}_{w}-{\rho}_{n}\right)\boldsymbol{g}\right)&\textnormal{in
porous regions}\end{cases}$ (2.94)

${\boldsymbol{F}}_{c,2}=\begin{cases}0&\textnormal{in solid-free regions}\\\
-p_{c}\nabla{\alpha}_{w}&\textnormal{in porous regions}\end{cases}$ (2.95)

${\boldsymbol{n}}_{w,n}=\begin{cases}\frac{\nabla{\alpha}_{w}}{\left|\nabla{\alpha}_{w}\right|}&\textnormal{in
solid-free regions}\\\ {cos\left(\theta\right)\
}{\boldsymbol{n}}_{wall}+{sin\left(\theta\right)\
}{\boldsymbol{t}}_{wall}&\textnormal{at the interface between solid-free
porous regions}\end{cases}$ (2.96)

${\boldsymbol{U}}_{r}=\begin{cases}C_{\alpha\
}max\left(\left|{\boldsymbol{U}}_{f}\right|\right)\frac{\nabla{\alpha}_{w}}{\left|\nabla{\alpha}_{w}\right|}&\textnormal{in
solid-free regions}\\\
{\phi}^{-1}\left[\begin{array}[]{c}-\left(M_{w}{\alpha}^{-1}_{w}-M_{n}{\alpha}^{-1}_{n}\right)\nabla
p+\\\
\left({\rho}_{w}M_{w}{\alpha}^{-1}_{w}-{\rho}_{w}M_{n}{\alpha}^{-1}_{n}\right)\boldsymbol{g}+\\\
\left(M_{w}{\alpha}_{n}{\alpha}^{-1}_{w}+M_{n}{\alpha}_{w}{\alpha}^{-1}_{n}\right)\nabla
p_{c}-\\\
\left(M_{w}{\alpha}^{-1}_{w}-M_{n}{\alpha}^{-1}_{n}\right)p_{c}\nabla{\alpha}_{w}\end{array}\right]&\textnormal{in
porous regions}\end{cases}$ (2.97)

where $\gamma$ is the surface tension, $C_{\alpha}$ is an interface
compression parameter (traditionally set to values between 1-4 in the Volume-
of-Fluid method), $k_{0}$ is the absolute permeability, $k_{r,i}$ are the
fluids’ relative permeabilities, $M_{i}=k_{i,r}/{\mu}_{i}\ $are the fluids’
mobilities, $M=M_{w}+M_{n}$ is the single-field mobility, and $p_{c}$ is the
average capillary pressure between the two fluids in a given control volume.
The phase-specific parameters can be readily calculated from closed-form
relative permeability and capillary pressure models such as the Brooks-Corey
(Brooks & Corey,, 1964) and Van Genutchen (van Genuchten,, 1980) models.
Lastly, $\theta$ is the imposed contact angle at the porous wall, and
$\boldsymbol{n}_{wall}$ and $\boldsymbol{t}_{wall}$ are the normal and
tangential directions relative to said wall, respectively. For the reader’s
convenience, a full implementation of this model, complete with rheological,
relative permeability, and capillary pressure models, is included within the
hybridPorousInterFoam (for static porous media) and hybridBiotInterFoam (for
deformable porous media) open-source solvers.

###  2.14 Conclusion

A complete set of partial differential equations required to model multiphase
flow in multiscale deformable porous media is presented above. In the
following chapters, we focus on the numerical implementation of these
equations and their application to a wide range of different natural and
engineered systems. As such, for the rest of this dissertation, whenever any
partial differential equation is mentioned, we will be referring to the ones
described in this chapter. The main model limitations stemming from the
assumptions presented in this chapter are summarized and discussed in Chapter
9.

[75mm] Computers are useless. They can only give you answers. Pablo Picasso

##  Chapter 3 Numerical Implementation

Having rigorously derived the governing equations for multiphase flow in
deformable porous media in Chapter 2, we now turn to the task of implementing
them into a computational framework apt for performing numerical simulations.
This chapter will provide the computational tools required to perform all the
simulations shown in subsequent chapters. The work presented here is adapted
from Carrillo & Bourg, (2019), Carrillo et al., (2020), and Carrillo & Bourg,
2021b .

###  3.1 Numerical Platform

The implementation of the multiphase DBB model was done in OpenFOAM®, a free,
open-source, parallelizable, and widely used computational fluid mechanics
platform. This C++ code uses the Finite Volume Method to discretize and solve
partial differential equations in complex 3-D structured and unstructured
grids. Its object-oriented structure 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. The implementation described below
stems directly from OpenFOAM’s® standard two-phase incompressible flow solver
“interFoam”.

###  3.2 Solution Algorithm

The solution of the governing equations is done in a sequential manner,
starting with the fluid mechanics equations and following with the solid
mechanics equations for every time step. Of particular importance is the
handling and modification of the velocity-pressure coupling required for
modeling incompressible fluids in conjunction with a moving solid matrix. For
this step, we based our solution algorithm on the Pressure Implicit Splitting-
Operator (PISO) (Issa,, 1986). First, we explicitly solve the fluid saturation
equation (Eqn. 2.89) for ${\alpha}^{t+1}_{w}$ through the Multidimensional
Universal Limiter of Explicit Solution (MULES) algorithm (Márquez & Fich,,
2013). This allows for stable numerical advection of the saturation field by
the application of Flux Corrected Transport Theory (Rudman,, 1997). Then, we
update the boundary values of ${\boldsymbol{U}}_{f}$ and
${\boldsymbol{U}}_{r}$ in addition to the cell-centered values of the
permeability $k^{t+1}$, density ${\rho}_{f}^{t+1}$, and viscosity
${\mu_{f}}^{t+1}$ based on the newly calculated saturation field
${\alpha}^{t+1}_{w}$. The capillary forces ${\boldsymbol{F}}^{t+1}_{c,i}$ are
also updated accordingly. After that, a preliminary value of the fluid
velocity ${\boldsymbol{U}}^{*}_{f}$ is calculated by implicitly solving the
algebraically discretized form of the fluid momentum equation used in the
Finite Volume Method.

$\displaystyle
a_{p}{\boldsymbol{U}}^{*}_{f}=\boldsymbol{H}\left({\boldsymbol{U}}^{*}_{f}\right)+{\rho}^{t+1}_{f}\boldsymbol{g}+{\boldsymbol{F}}^{t+1}_{c,i}-\nabla
p^{t}$ (3.1)

where $\boldsymbol{H}\left({\boldsymbol{U}}^{*}_{f}\right)$ contains inertial,
convective, viscous, and drag source terms originating from neighboring cells
and $a_{p}$ represents these same terms but at the volume of interest. Note
that the ${\boldsymbol{U}}^{*}_{f}$ field does not follow mass conservation.
To account for this, we use the fluid continuity equation (Eqn. 2.88) in
conjunction with the previous equation (Eqn. 3.1) to update the velocity field
${\boldsymbol{U}}^{**}_{f}$ and calculate a preliminary mass-conservative
pressure field $p^{*}$. In other words, these fields must satisfy,

$\displaystyle{\boldsymbol{U}}^{**}_{f}=\frac{1}{a_{p}}\left(\boldsymbol{H}\left({\boldsymbol{U}}^{*}_{f}\right)+{\rho}^{t+1}_{f}\boldsymbol{g}+{\boldsymbol{F}}^{t+1}_{c,i}-\nabla
p^{*}\right)$ (3.2)

$\displaystyle\nabla\cdot{\boldsymbol{U}}^{**}_{f}=-\
\frac{\partial{\phi}_{f}}{\partial t}$ (3.3)

These equations can be recast into a single coupled equation which is then
used to implicitly solve for pressure. This step can be done through several
generalized matrix solvers that are standard in OpenFOAM®.

$\displaystyle\nabla\cdot\left(\frac{1}{a_{p}}\left(\boldsymbol{H}\left({\boldsymbol{U}}^{*}_{f}\right)+{\rho}^{t+1}_{f}\boldsymbol{g}+{\boldsymbol{F}}^{t+1}_{c,i}\right)\right)-\nabla\cdot\left(\frac{1}{a_{p}}\nabla
p^{*}\right)=-\ \frac{\partial{\phi}_{f}}{\partial t}$ (3.4)

After solving for pressure $p^{*}$, velocity can be re-calculated from
Equation 3.2. This semi-implicit pressure-velocity correction step is repeated
until the desired convergence is reached. It has been shown that at least two
pressure-velocity correction loops are required to ensure mass conservation
(Issa,, 1986).

At this point ${\boldsymbol{U}}^{t+1}_{f}$ and $p^{t+1}$ are set and used as
input values for updating the drag and pressure source terms present in the
solid mechanics momentum equation (Eqn. 2.92). Then, in the case of visco-
poro-plasticity, said equation is discretized in a similar way as the fluid
momentum equation (Eqn. 2.90) and used to implicitly solve for
${\boldsymbol{U}}^{t+1}_{s}$. In the case of poroelasticity, the solid
mechanics equation is solved through the algorithm presented in Jasak &
Weller, (2000). Here, the solid’s elastic equation (Eqn. 2.74) is discretized
and segregated into implicit and explicit components, after which it is
iteratively solved until convergence is reached. This segregated method not
only guarantees fast convergence but also memory efficiency. Finally, the
updated solid velocity is used to “advect” the solid fraction field
${\phi}_{s}$ by solving the mass conservation equation (Eqn. 2.91). At this
point the algorithm advances in time according to the imposed Courant-
Friedrichs-Lewy (CFL) number. Further discussion regarding the discretization
techniques and matrix-solution procedures can be found in Carrillo & Bourg,
(2019), Carrillo et al., (2020), Carrillo & Bourg, 2021b , Jasak, (1996), and
Jasak & Weller, (2000).

###  3.3 Open-Source Implementation

The complete set of governing equations and solution algorithms, along with
the necessary rheology, relative permeability, and capillary pressure models
(Appendix A and B) were implemented into the solvers hybridPorousInterFoam
(for static porous media) and  thybridBiotInterFoam (for deformable porous
media). These solvers, along with representative tutorial cases, automated
compilation and running procedures, and all the simulated cases presented in
this dissertation were incorporated into open-source CFD packages of the same
names. OpenFOAM® and our code are free to use under the GNU general public
license and can be found at https://openfoam.org/ and
https://github.com/Franjcf (Carrillo, 2020b, ; Carrillo, 2020a, ),
respectively.

[75mm] In God we trust, all others bring data. William Edwards Deming

##  Chapter 4 Applications to Single-Phase Flow in Multiscale Deformable
Porous Media

In this chapter we apply the equations derived in Chapter 2 and implemented in
Chapter 3 to model single phase flow through deformable porous media. Special
emphasis is given to the simulation of of swelling clays and their influence
on the permeability of sedimentary rocks. The work presented in this chapter
represents the first step towards verifying and showcasing the Multiphase DBB
model. This chapter is adapted from Carrillo & Bourg, (2019).

###  4.1 Introduction

How does the permeability of sedimentary media depend on porosity, mineralogy,
fluid chemistry, and stress history? This question has been a recurrent theme
in subsurface hydrology for over half a century (Berg,, 1970; Bourg & Ajo-
Franklin,, 2017; Brace,, 1980). It impacts a range of endeavors that shape
humanity’s energy landscape including hydrocarbon migration and recovery
(Alvarado & Manrique,, 2010), geothermal energy production (Barbier,, 1997),
geologic carbon sequestration (Klaus,, 2003), and radioactive waste storage
(Sellin & Leupin,, 2014). The last two technologies have the potential to
contribute up to half of the mitigation effort required to stabilize global
CO2 emissions (Metz et al.,, 2005; Socolow & Pacala,, 2004) but require the
ability to accurately predict the permeability evolution of ductile fine-
grained sedimentary rocks (shale, mudstone) over millennial time-scales in the
presence of geochemical and geomechanical disturbances (Bourg & Ajo-Franklin,,
2017; Neuzil,, 2013).

A major challenge associated with developing a predictive understanding of
flow in fine-grained sedimentary media is that these structures have two
characteristic length scales: a macroscale defined by the assemblages of
coarse grains of quartz, feldspar, or carbonate over distances of micrometers
and a microscale defined by the assemblages of clay minerals (primarily
smectite and Illite) over distances of nanometers (Bourg & Ajo-Franklin,,
2017). Figure 4.1a shows a conceptual model of the macroscale structure of
sedimentary rocks as a mixture of rigid coarse grains, a deformable
microporous clay matrix, and macropores (Crawford et al.,, 2008; Marion et
al.,, 1992; Revil & Cathles,, 1999). This model is consistent with electron
microscopy observations (Nadeev et al.,, 2013; Nole et al.,, 2016; Tuller &
Or,, 2003) as exemplified in Fig. 4.1b-d (Fiès & Bruand,, 1998; Peters,,
2009). Experimental data indicate that permeability in fine-grained soils,
sediments, and sedimentary rocks can be highly sensitive to the spatial
distribution of the clay matrix (Abichou et al.,, 2002; Nadeau,, 1998).

A second major challenge is that the microporous clay matrix is non-rigid: it
can swell or shrink in response to changes in salinity and deform in response
to fluid flow or external stresses in a manner that reflects the nanoscale
colloidal interactions between negatively-charged clay particles (Liu,, 2013;
Madsen & Muller-Vonmoos,, 1989; Suzuki et al.,, 2005; Teich-McGoldrick et
al.,, 2015). The resulting dynamics of the clay matrix give rise to
significant couplings between the hydrologic, chemical, and mechanical (HCM)
properties of clayey media (Carey et al.,, 2014; Erol,, 1979; Murad &
Cushman,, 1997). These couplings are particularly strong if the clay fraction
contains significant amounts of smectite (i.e., swelling clay minerals) and if
the pore water contains predominantly sodic salts. In these cases, swelling by
more than 1400 % and swelling pressures up to $\mathrm{\sim}$50 MPa have been
reported (Karnland et al.,, 2007; Norrish,, 1954). In short, fundamental
predictions of the permeability of fine-grained sediments and sedimentary
rocks require a model capable of describing coupled fluid flow in pores with
very different sizes (intergranular macropores and clay micropores with pore
widths on the order of micrometers and nanometers, respectively) while also
predicting the deformations of the microporous clay matrix induced by flow,
external stresses, and salinity changes (Bourg & Ajo-Franklin,, 2017).

Existing approaches to predicting the hydrology of fine-grained sedimentary
media have focused on addressing either one of the two challenges outlined
above, but not both simultaneously. Approaches focused on the existence of two
length scales have generally used an “ideal packing model” representation of
clay-rich sedimentary media as a microporous (clay) medium embedded within a
network of coarse grains (Crawford et al.,, 2008; Marion et al.,, 1992; Revil
& Cathles,, 1999). On this model, a threshold naturally emerges at a clay
mineral mass fraction of $\mathrm{\sim}$1/3 where the microporous clay matrix
becomes the load-bearing phase, in agreement with experimental observations on
the core-scale hydrologic and mechanical properties of sedimentary rocks
(Bourg & Tournassat,, 2015; Crawford et al.,, 2008). However, existing models
based on this framework invariably neglect the mechanics of the clay matrix by
assuming either that the clay has a fixed porosity or that it uniformly
occupies the space available within the network of coarse grains. One
consequence of this approximation is that these models do not capture the
influence of salinity on the permeability of clayey media (Kwon et al.,, 2004;
Quirk,, 1986).

Conversely, approaches focused on the HCM couplings such as Terzaghi’s
consolidation theory, Biot’s theory of poroelasticity, and Mixture Theory
(theories initially developed to describe clayey media and now widely applied
to other deformable porous media including hydrogels and biological tissues)
simplify the governing equations for fluid and solid dynamics into a single
macroscopic equation by assuming that both phases are superimposed continua
with negligible inertial forces (Auton & MacMinn,, 2017; Barry & Aldis,, 1992;
Jain & Juanes,, 2009; Santillán et al.,, 2018; Terzaghi,, 1964). This approach
results in an implicitly coupled momentum equation for deformable porous media
that can be paired to different body forces or solid deformation constitutive
relations to create system-specific models. A drawback of this approach in the
context of sedimentary rocks, however, is that it does not reflect the
existence of the two characteristic length scales illustrated in Fig. 4.1 or
consequences such as the permeability threshold at a finite clay content noted
above.

Figure 4.1: (a) Conceptual representation of the two key challenges posed by
fluid flow in clay-rich sedimentary rocks. The structure contains three
regions: rigid coarse grains, macropores, and a deformable microporous clay
matrix. Fluid flow (red arrows) occurs primarily in the bulk-fluid domain, but
the boundaries of this domain are influenced by the deformation of the
microporous clay domain (black arrows). Figures b-d are electron microscopy
images of clayey media, specifically (b) a Canadian shale and (c,d) a mixture
of clay and sand (denoted by the letters c and s, respectively) associated
with reductive and expansive clay swelling states (Fiès & Bruand,, 1998;
Peters,, 2009). Coarse grains, microporous clay, and macropores are shown in
gray, orange, and blue in Fig. 4.1a and in light gray, dark gray, and black in
Figs. 1b-d.

The advent of Digital Rock Physics (DRP) over the last decade—i.e., the
combined use of X-ray computed tomography (XCT) and computational methods to
generate numerical models of fluid flow in real rocks (Adler et al.,, 1992;
Fredrich,, 1999; Mirabolghasemi et al.,, 2015; Raeini et al.,, 2017)—provides
a potential route to simultaneously addressing both challenges outlined above.
In the last few years, a handful of DRP studies have examined rocks with two
characteristic length scales including fractured porous rocks, vuggy media,
and clay-rich sedimentary rocks (Keller et al.,, 2013; Li et al.,, 2016; Saif
et al.,, 2017). In these studies, the macroscale features (coarse grains,
macropores, and microporous regions) are fully resolved, while the microscale
features (i.e., particles and pores within the microporous regions) are below
the resolution of the XCT measurements. Many of these studies used a
computational framework based on a pore network model (PNM) with two
characteristic length scales (Mehmani & Prodanović,, 2014; Mehmani et al.,,
2013; Prodanović et al.,, 2015). Alternatively, at least two used an approach
based on computational fluid dynamics (CFD) where flow in the macropores and
microporous regions were coupled through slip-flow boundary conditions and
volume-averaging based on the Darcy-Brinkman formalism (Bijeljic et al.,,
2018; Guo et al.,, 2018). The resulting models demonstrate the possibility of
characterizing porous media with two characteristic length scales as mixtures
of three regions (as in Fig. 4.1a) and in coupling fluid flow in the
macropores and microporous regions (Golfier et al.,, 2015; Soulaine et al.,,
2017). However, these studies all assumed a static porous medium
configuration, a significant limitation in the case of clay-rich media, as
noted above.

Figure 4.2: Illustration of the three main theoretical frameworks combined in
the present study to model fluid flow in porous media with two characteristic
length scales and a deformable microporous matrix: (a) models based on the
Darcy-Brinkman formulation that describe coupled fluid flow in macropores
(clear region) and in a microporous medium (shaded region) (Ochoa-Tapia &
Whitaker,, 1995); (b) boundary methods that simulate flow around moving
impermeable solids (Breugmem & Verzicco,, 2013); and (c) poromechanics model
of a cylindrical conduit where the fluid and porous solid are in a single
domain (grey) and forces are imposed through boundary conditions (Auton &
MacMinn,, 2017). Here $Q$ represents the fluid flow rate, and $a(t)$ and
$b_{0}$ represent the inner and outer radii of the cylinder, respectively.

The main contribution of this chapter is the application of the system of
differential equations derived in Chapter 2 to single phase flow in deformable
porous media. The model is parameterized with a focus on viscoplastic clay-
rich sedimentary media (Section 4.2) and validated against experimental
results on the properties of these media (Section 4.3). With proper
parameterization, our framework should be applicable to other systems that
involve coupled fluid flow in macropores and in a deformable microporous
matrix such as soils (Murad & Cushman,, 1996; Vo & Heys,, 2011), hydrogels
(Datta et al.,, 2016), biological tissues (Dai et al.,, 2014), and fractures
(Noiriel et al.,, 2007). To illustrate this versatility, we also formulate and
apply the model for cases where the microporous matrix undergoes elastic
rather than viscoplastic deformation (Section 4.4).

###  4.2 Parametrization to Porous Media: Na-Smectite Clay

As noted in the introduction, a key challenge in subsurface hydrology is the
strong influence of clay minerals on fluid flow in sedimentary rocks (Bourg &
Ajo-Franklin,, 2017). This challenge is particularly profound in systems that
contain smectite clay, a microporous material whose permeability, swelling
pressure, and plastic rheology are highly sensitive to porosity and aqueous
chemistry (Aksu et al.,, 2015; Mondol et al.,, 2008; Spearman,, 2017).
Smectite is the predominant clay mineral in many unconsolidated sediments, in
bentonite (a mixture of sand and clay used in geotechnical applications), and
in soils formed in temperate weathering environments (Abichou et al.,, 2002;
Sposito,, 2008). Smectite and Illite clay minerals constitute roughly half of
the world’s sedimentary rock mass, with Illite being generally more abundant
on a mass basis but both clay minerals being roughly equally important on a
surface area basis. The present section parameterizes the set of equations
derived above for systems where the microporous matrix consists entirely of
smectite by choosing the appropriate permeability, salt diffusion, swelling
pressure, and rheology models. We focus on conditions where aqueous chemistry
is dominated by sodic salts such as NaCl, where smectite swelling and rheology
have been most extensively characterized. The result is a closed, fully
coupled system of equations.

####  4.2.1 Parametrization of the Permeability Function

First, we define the permeability of the clay matrix ($k$), the key variable
that determines the sub-REV-scale momentum interaction between water and clay.
For simplicity, we use the well-known Kozeny-Carman relation:

$k\left({\phi}_{f}\right)=\frac{1}{b\ t^{2}a^{2\ }_{s}}\
\frac{{\phi}^{3}_{f}}{{\left(1-{\phi}_{f}\right)}^{2}}$ (4.1)

Equation 4.1 describes the permeability of a bundle of capillary tubes of
uniform radius within an impermeable solid, where $b$ is a pore shape factor,
$t$ is the tortuosity of the pore network, and $a_{s}$ is the specific surface
area expressed as the ratio of surface area to solid volume. The parameter
values used in applying Eqn. 4.1 to Na-smectite and Illite are summarized in
Table 4.1. The shape factor was set to b = 3 $\mathrm{\pm}$ 1 based on
empirical observations in a range of porous media (Bourg & Ajo-Franklin,,
2017). Tortuosity was set to $t^{2}$ = 4.0 $\mathrm{\pm}$ 1.6 based on
measurements of water tracer diffusion in compacted Na-smectite (Bourg et
al.,, 2006). The specific surface area was set to 1.7 $\mathrm{\pm}$ 0.2
nm${}^{-1\ }$for smectite and 0.21 $\mathrm{\pm}$ 0.03 nm-1 for Illite based
on experimental values reported in previous studies ($a_{s}$ = 703
$\mathrm{\pm}$ 60 m2 g-1 for smectite and as = 78 $\mathrm{\pm}$ 10 m2 g-1 for
Illite, converted to units of area per volume using grain densities of 2400
$\mathrm{\pm}$ 140 kg m${}^{-3\ }$and 2700 $\mathrm{\pm}$ 70 kg m-3,
respectively) (Brooks,, 1955; Diamond & Kinter,, 1956; Mesri & Olson,, 1971;
Orchiston,, 1954; Quirk,, 1955). The values of specific surface area used here
are based on methods that probe the entire water-accessible surface area of
clay particles, including glycerol and ethylene glycol monoethyl ether (EGME)
retention techniques but not standard N2 adsorption techniques (Diamond &
Kinter,, 1956; Tournassat et al.,, 2015).

A well-known limitation of Eqn. 4.1 is that it does not capture the influence
of pore-size heterogeneity (Bourg & Ajo-Franklin,, 2017; Dixon et al.,, 1999)
and, therefore, does not accurately predict the macro-scale permeability of
clay-rich soils or sedimentary rocks (Mondol et al.,, 2008; Ren et al.,,
2016). The comparison with experimental data in Fig. 4.3A, however, indicates
that Eqn. 4.1 provides a reasonable permeability model for pure Na-smectite
and Illite without the need for fitting parameters.

Two important features of the permeability of pure clay are not described by
Eqn. 4.1. First, the permeability data in Fig. 4.3A relate exclusively to
sodium-exchanged clays. Other data indicate that smectite (and, to a smaller
extent, Illite) has a higher permeability at the same porosity values if
equilibrated with calcium or potassium electrolytes (Mesri & Olson,, 1971).
This permeability difference likely reflects the less uniform pore-size
distributions of K- and Ca-smectite (relative to Na-smectite) that results
from the stronger charge-screening capacity of K and Ca ions at the clay
surface. Second, clay-water mixtures exposed to one-dimensional consolidation
are known to develop a textural anisotropy that makes their permeability
tensor anisotropic (Hicher et al.,, 2000). Overall, Eqn. 4.1 is used here as a
simple, yet reasonably accurate permeability model in the range of conditions
of interest, i.e., fully-saturated Na-smectite or Illite in the absence of
excessive anisotropy. Modifications to account for clay anisotropy, other clay
counterions (such as calcium or potassium), multiphase flow effects (e.g., the
Klinkenberg Effect), or to implement other permeability models (Chapuis &
Aubertin,, 2003; Samarasinghe et al.,, 1982) can be readily carried out within
the present framework.

####  4.2.2 Parametrization of the Swelling Pressure Term

Second, we define the swelling pressure of the microporous clay matrix. For
this, we use the well-known Derjaguin-Landau-Verwey-Overbeek (DLVO) theory of
colloidal interactions. More precisely, we use a semi-empirical formulation
proposed by Liu, (2013) for the disjoining pressure in a water film of
thickness “h” between clay tactoids (i.e., stacks of clay particles separated
by $\mathrm{\leq}$ 3 water layers):

$\begin{split}p_{swell}=2{\hat{C}}_{\mathrm{NaCl}}&RT\mathrm{cosh}\left(y_{m}\
-1\right)\\\
&-\frac{A_{h}}{6\pi}\left(\frac{1}{h^{3}}-\frac{2}{{\left(h+D_{p}\right)}^{3}}+\frac{1}{{\left(h+2D_{p}\right)}^{3}}\right)+S_{0}{\mathrm{exp}\left(-\frac{h}{l}\right)\
}\ \end{split}$ (4.2)

In Eqn. 4.2, the first term on the right side is Langmuir’s model (Kemper &
Quirk,, 1972; Langmuir,, 1938) of the osmotic swelling pressure caused by
overlapping electrostatic double layers in a slit-shaped nanopore, where
${\hat{C}}_{\mathrm{NaCl}}$ is the salt concentration in mol m-3, R is the
ideal gas constant (8.31446 J K-1 mol-1), T is the temperature in Kelvin, and
$y_{m}$ is the scaled electrostatic potential at the mid-plane of the
nanopore, calculated here using the so-called “compression approach” for
solving the Poisson-Boltzmann equation between two charged parallel plates in
a symmetrical electrolyte solution (Liu & Neretnieks,, 2008). The salinity
dependence of $y_{m}$ gives rise to the well-known salinity dependence of clay
swelling. The second term is the contribution of London dispersion forces to
the swelling pressure modeled through Hamaker’s approach (Hamaker,, 1937),
where$\ D_{p}$ is the thickness of clay tactoids and $A_{h}$ is Hamaker’s
constant for clay tactoids separated by a water film. The third term is an
empirical description of the short-range “hydration repulsion” between clay
particles, where S0 and l are empirical coefficients. Parameter values and the
relation between ${\phi}_{f}$ and h were taken from Liu, (2013) and are listed
in Table 4.1. Equation 4.21 as parametrized by Liu, (2013) is used here as a
convenient parametric fit to experimental data on smectite swelling as a
function of ${\phi}_{f}$ and C.

We note that Eq. 4.2 makes significant simplifications and approximations. For
example, the hydration repulsion term is purely phenomenological, while the
optimal formulation and parameterization of the London dispersion term is
still unsettled in the case of smectite clay (Gilbert et al.,, 2015; Tester et
al.,, 2016). In addition, the first two terms in Eqn. 4.2 are based on mean-
field theories, i.e., they neglect short-range interactions between water,
ions, and clay surfaces (McBride,, 1997; Missana & Adell,, 2000). In the case
of Na-smectite, this last approximation is valid only at interparticle
distances greater than 3 nm (i.e. ${\phi}_{f}>0.75$) (Adair et al.,, 2001).
One consequence of this is that DLVO theory does not predict the existence of
stable “crystalline” swelling states with interparticle distances
$\mathrm{\leq}$ 1 nm that predominate at high salinity, high compaction, or in
aqueous chemistries dominated by divalent ions (Bourg & Ajo-Franklin,, 2017;
Pashley & Israelachvili,, 1984; Shen & Bourg,, 2021). Furthermore, this
formulation assumes that clay swelling pressure is controlled by a single
microstructural variable (interparticle distance h or, equivalently, porosity
${\phi}_{f}$); more complex formulations have been proposed in the case of
other microporous media with a deformable solid skeleton, such as activated
carbon or zeolites (Pijaudier-Cabot et al.,, 2011). Finally, the numerical
method used by Liu, (2013) to evaluate $y_{m}$ cannot be applied when C = 0;
for simplicity we use C = 0.001 M as an approximation when calculating
$p_{swell}$ in pure water (resulting in a $\mathrm{\sim}$0.05% error in said
value).

Figure 4.3: (a) Permeability of compacted Na-smectite and Illite as a function
of porosity; (b) swelling pressure of compacted bentonite as a function of
porosity at different salinities (error bars are smaller than the symbols);
and (c) porosity of Na-smectite as a function of salinity at a small constant
confining pressure. Symbols are experimental results (Mesri & Olson,, 1971;
Norrish,, 1954; SKB,, 2011). Lines are model predictions obtained using Eqns.
4.1 and 4.2 with no fitting parameters in Fig. 4.3A,B and with a single
fitting parameter (the unspecified small confining pressure used in the
experiments, which was set to 0.03 MPa; dashed lines show the sensitivity of
the model predictions to $\pm$ 0.01 MPa differences in this value) in Fig.
4.3C.

To test the accuracy of Eqn. 4.2 as parameterized by Liu, (2013), we used
macroscopic data on the swelling pressure of confined water-saturated Na-
smectite as a function of compaction and salinity (Fig. 4.3B) as well as X-ray
diffraction (XRD) measurements of h vs. salinity (Fig. 4.3C). Figure 4.3 shows
that Eqn. 4.2 overestimates porosity at salinities $\geq$ 0.5 M NaCl, where
crystalline swelling predominates, but accurately predicts experimental data
at lower salinities.

####  4.2.3 Parametrization of the Diffusion Coefficient of Dissolved Salts

Third, we define the effective diffusion coefficient of dissolved salt (x) in
the microporous clay using a form of Archie’s law (Blum et al.,, 2007; Boving
& Grathwohl,, 2001):

$D_{\mathrm{e,}x}={\phi}^{n}_{f}D_{x}$ (4.3)

Although there does not exist a universal Archie’s law exponent n for all
porous media, several studies have concluded that a value of 2.2 to 2.5
reasonably reproduces the diffusion of ions in pure compacted clay (Revil et
al.,, 2011; Shen & Chen,, 2007; Van Loon et al.,, 2003; Van Loon,, 2015).

Equation 4.3 is known to oversimplify salt diffusion in clayey media,
particularly by neglecting the salinity-dependence of $D_{\mathrm{e,}x}$ that
arises from anion repulsion in the electrical double layer on charged clay
surfaces (Malusis et al.,, 2003; Sherwood & Craster,, 2000; Van Loon et al.,,
2007; Underwood & Bourg,, 2020). A notable alternative to Eqn. 4.3 that
accounts for this effect is Kemper’s model (Kemper & van Schaik,, 1966; Kemper
& Rollins,, 1966).

$D_{\mathrm{e,}x}=\frac{1-e^{-y_{m}}}{{\ t}^{2}}D_{x}$ (4.4)

where $t$ is tortuosity as in Eqn. 4.1, ym is the same as in Eqn. 4.2, and
“$1-e^{-y_{m}}$” is the equilibrium ratio between anion concentration in the
clay nanopores and in bulk liquid water. Equation 4.4 is less empirical than
Eqn. 4.3 and more consistent with the permeability and swelling pressure
models used above. However, its use would require a more complex treatment of
the coupled advection and diffusion of ions in Eqn. 2.87, where the impact of
anion exclusion is neglected. We note that multiphase flow effects such as
Knudsen diffusion are neglected here due to the fact that we focus exclusively
in fully water-saturated media (Malek & Coppens,, 2003). Equation 4.3 is
therefore used here for simplicity.

####  4.2.4 Parametrization of Clay Rheology

Lastly, we define the effective viscosity of the microporous clay in a manner
that accounts for the impact of clay swelling and shearing on its plastic
viscosity and critical shear stress (Güven,, 1993; Maciel et al.,, 2009).
Spearman, (2017) developed a model based on floc fractal theory that predicts
the rheological properties of a wide variety of clays when sheared at
different solid-water ratios:

${\mu}^{eff}_{s}={\left[{\mathrm{\tau}}^{\frac{1}{4-D}}_{0}+{\left(1+\frac{1}{\beta{\gamma}^{m}}{\left(\frac{{\tau}_{0}}{{\mu}_{*}}\right)}^{m}\right)}^{\frac{1}{4-D}}{\left({\mu}_{*}\gamma\right)}^{\frac{1}{4-D}}\right]}^{4-D}{\gamma}^{-1}$
(4.5)

${\tau}_{0}={\tau}^{*}{\left[\frac{{\phi}_{s}/{\phi}^{max\
}_{s}}{(1-{\phi}_{s}/{\phi}^{max}_{s})}\right]}^{4-D}\ \ \ \ \ \vdots\ \ \ \ \
\
{\mu}_{*}=\frac{{\mu}_{f}}{{\left(1-{\phi}_{s}/{\phi}^{max}_{s}\right)}^{2}}\
\ $ (4.6)

In Eqns. 4.5-4.6, ${\phi}^{max}_{s}$ is the maximum solid fraction of clay,
${\tau}^{*}$ the critical shear stress at ${\phi}_{s}=0.5\ {\phi}^{max}_{s}$,
D is the fractal dimension, $\beta$ is a structural parameter, m is a
structural break-up parameter, and $\gamma$ is the strain rate (calculated at
each time step as the sum of symmetric components of $\nabla{\hat{U}}_{s}$).
With proper tuning, this model was shown to successfully describe the
rheologies of smectite, Illite, and kaolinite clays. Parameter values were
taken from Spearman, (2017) based on fitting rheological data pertaining to a
90% smectite-water mixture obtained by Coussot et al., (1993) and are listed
in Table 4.1.

We note that Eqns. 4.5 and 4.6 do not explicitly consider potential impacts of
salinity or dynamic changes in solid fraction (as in the case of shrinking and
swelling) on rheology. Instead, the Spearman model assumes that the clay is at
an equilibrated solid fraction when calculating the effective viscosity. For
simplicity, we assume that the resulting effective viscosity applies not only
to shear, but also to expansion and contraction.

The use of Eqns. 4.1-4.6 in the modeling framework derived in Chapter 2
provides closure for the coupled system of equations in the case of pure
smectite in sodium electrolytes such as NaCl at salinities up to 0.5 M. As
noted above, this closure relies on a number of assumptions and
approximations. Notably, it neglects clay fabric anisotropy (Hicher et al.,,
2000; Tessier,, 1990) and clay dispersion into the bulk water phase.
Furthermore, it relies on the assumption that DLVO theory, the Kozeny-Carman
equation, Spearman’s model, and the Fickian diffusion equation described above
are accurate in conditions beyond their validation in Fig. 4.3 (for example,
during dynamics shearing, shrinking, and swelling). The assumptions and
approximations listed above are not intrinsic to the modeling framework
derived in Chapter 2. They can be readily addressed in future studies as
additional information become available on the microscale permeability,
swelling pressure, rheology, and dispersion of microporous clay.

###  4.3 Model Validation and Application to Plastically Deformable Porous
Media (Clay)

####  4.3.1 Model Verification by Comparison to the Fluid-Driven Deformation
of a Clay Filter Cake

A quantitative validation of our model was realized by testing its ability to
predict the compaction of a clay plug caused by water flow through said plug.
To do this, we compared our numerical results to predictions of an analytical
model derived by Hewitt et al., (2016) at equivalent experimental conditions.
Said analytical model is based on Biot Theory and is able predict the 1-D
solid deformation of the plug as a function its initial porosity, the
deformation modulus of the solid, and the fluid pressure gradient across the
plug. Hewitt et al., (2016) validated their model against experimental results
obtained with an elastically deformable microporous medium placed within a 12
cm by 25 cm container based on the assumptions that permeability follows the
Kozeny-Carman equation, stresses follow the Terzaghi stress principle, and
fluid flow (and thus the resulting deformation-inducing pressure gradient) is
governed by Darcy’s law.

To replicate these conditions, our simulations where carried out on a 240 by
300 grid representing a 2-D container (12 cm by 15 cm) partially filled with
non-swelling clay. We induced fluid flow through the clay by a constant
pressure boundary condition at the top of the container, where fluid (and only
fluid) was allowed to leave through the lower boundary. In Fig. 4.4, we
compare our steady-state, non-dimentionalized numerical model predictions for
Na-smectite against Hewitt et al., (2016) analytical model at two initial
porosity values, showing good agreement between both models. Since Hewitt et
al., (2016) did not include a swelling pressure into their analytical
framework, we set $p_{swell}=0$ in our numerical simulations to ensure proper
comparison and verification of the flow-deformation mechanics. Please refer to
Table 4.1 for a detailed listing of the parameter used in this and all other
subsequent simulations.

Figure 4.4: Model prediction of the compaction associated with water flow
through a clay plug as a function of pressure drop (from 0 to 200 kPa) across
the clay plug: (a) schematic view of the simulated system; (b) equilibrated
system as predicted by our simulation framework (colors show clay fraction
${\phi}_{s}$ in the microporous region and fluid velocity ${\overline{U}}_{f}$
in the free fluid region); (c) comparison between Hewitt’s analytical model
(solid lines) and the predictions of our numerical model (symbols) for initial
clay porosities of 0.55 and 0.65; and (d) final configurations of the system
as a function of increasing pressure drop. The simulations were carried out in
a 1 by 300 grid (c) and a 240 by 300 grid (b, d). Results are reported as the
non-dimensional volume change of the clay plug vs. non-dimensional pressure
drop across the clay plug. Pressure in the simulations was non-dimentionalized
using a deformation modulus calculated from the clay’s critical stress at its
initial porosity [p* = p/$\sigma$, where $\sigma={\tau}^{initial}_{0\
}/({\phi}^{initial}_{f}/{\phi}^{final}_{f}-1)$]. For both simulations, mean
average error (MAE) was about 0.04 times the final measured volume change.

The resulting steady-state deformation profiles are a result of the balance
between the clay’s structural forces (critical stress) and the forces imposed
by the fluid (viscous drag): as compaction increases, the permeability of the
porous medium decreases, which in turn increases the viscous stresses imposed
on the medium. At the same time, given a constant pressure gradient, flow
through the plug decreases as permeability decreases, which then reduces the
magnitude of the viscous stresses. Steady-state is achieved once the viscous
drag is balanced by the porous medium’s structural forces, which explains why
systems with a higher initial compaction deform less than their counterparts
(Fig. 4.4). The good agreement shown in Fig. 4.4C is expected (both models
rely on similar assumptions such as the validity of the Kozeny-Carman equation
and Terzaghi and Biot’s principles) but also provides some confidence in the
ability of our model framework to represent the feedbacks between fluid flow
and solid deformation in microporous media.

####  4.3.2 Model Verification by Comparison to Oedometric Clay Swelling
Experiments

As a further quantitative validation of our model (and to verify the swelling
pressure effects disregarded in the previous section) we tested model
predictions against standard oedometric measurements of clay swelling driven
by a salinity change. Specifically, we used measurements by Di Maio, (1996) of
the volumetric swelling of water-saturated Na-smectite samples exposed to a
salinity shock (from 1 to 0 M NaCl) at different confining pressures (Fig.
4.5). The experimental conditions were straightforward: individual water-
saturated clay samples were first compressed within a oedometer chamber to
confining pressures of 40 kPa, 160 kPa, and 320 kPa. Subsequently, the samples
were placed in contact with a 1 M NaCl solution through a porous boundary
until equilibrated. Swelling was then induced by replacing the saltwater
solution with distilled water while maintaining a constant confining pressure.

The swelling portion of the experimentas was modeled by defining a 500 by 300
grid representing a 2-D container (5 cm by 3 cm) filled with 2 cm of 1M NaCl-
equilibrated compacted smectite clay, where the clay was only allowed to swell
in the positive y-direction as a result of salt diffusion out of the container
at the lower boundary (where salt concentration was set to zero). All other
boundary conditions on the container walls were set to replicate impermable
surfaces with no-slip boundary conditions. For simplicity, the confining
pressure was assumed uniform throughout the clay sample and was applied as a
constant in the Terzaghi stress tensor (see Section 2.7).

Figure 4.5: Macroscopic swelling of NaCl-saturated smectite clay (at confining
pressures 40, 160, and 320 kPa) initially equilibrated with a 1M NaCl solution
and then placed in contact with a reservoir of distilled water. Experimental
results (solid lines) and approximated error (shaded regions) were obtained
using data from Di Maio, (1996). Simulation predictions (symbols) were
obtained using two-dimensional 500 by 300 grids. For the three simulations,
MAE was 0.06 to 0.09 times the measured final volume change.

As shown in Figure 4.5, our model accurately captures the swelling kinetics of
Na-smectite clay driven by osmotic uptake of liquid water. The agreement
between our model and experimental data provides further confidence into the
ability of our model framework to capture feedbacks between hydrology,
mechanics, and salt transport in deformable microporous media. Overall,
Figures 4.4 and 4.5 show that our model, with no parameter fitting, yields
reasonably accurate predictions (with normalized MAEs below 0.09) of the
overall extent of clay swelling or compaction and of the associated kinetics.
Agreement is obtained regardless of whether the clay volume change is driven
by fluid flow (Fig. 4.4) or salinity changes (Fig. 4.5). This good agreement
suggests that our model captures key features of the coupled HCM behavior of
Na-smectite despite the assumptions and simplifications noted above.

####  4.3.3 Model Prediction of the Permeability of a Bead Pack Containing Na-
Smectite as a Function of Salinity and Clay Content

Having parametrized and tested our simulation framework, we used it to predict
the permeability of an idealized model of fine-grained soils and sedimentary
media: a bead pack containing coarse non-porous beads uniformly coated with
Na-smectite clay. The idealized system simulated here has previously been used
both as a conceptual model of the hydrology of sedimentary media and as an
idealized experimental model of the properties of engineered clay barriers
(Abichou et al.,, 2002; Revil & Cathles,, 1999; Tuller & Or,, 2003). To the
best of our knowledge, it has not been previously implemented in a numerical
framework that accounts for the HCM couplings in microporous clay. As
implemented here, this model system captures both key features of fine-grained
soils and sedimentary media identified in the introduction: the co-existence
of two characteristic length scales and the HCM coupled properties of the
microporous clay. As noted above, Na-smectite is used here because reasonably
accurate constitutive models exist to describe its swelling pressure,
permeability, and rheology at the microscopic scale, even though the
microporous regions in most soil and sedimentary media generally consist of
more complex mixtures of Ca/Na-smectite, other clay minerals, and organic
matter.

Our simulations are based on a representative, yet simplified, rock geometry
built upon the following assumptions: first, the medium’s coarse-grained,
load-bearing structure can be represented as a 2-D cross-section of a 3-D
randomly distributed spherical packed bed with a porosity of 0.64 obtained
from Finney, (1970); second, the initial distribution of clay is modeled as a
uniform film coating the coarse grains (a reasonable approximation given
previous imaging studies) (Abichou et al.,, 2002; Aksu et al.,, 2015; Peters,,
2009), and third, a two-dimensional geometry is sufficient to capture major
features of the simulated system (as in Quispe et al., (2005)). These are
significant approximations; in particular, feedbacks between fluid flow and
clay dynamics may create some degree of non-uniformity in the clay
distribution within sedimentary rocks (Song & Kovscek,, 2015) and percolation
thresholds associated with pore clogging may occur more readily in two-
dimensional than in three-dimensional systems.

Figure 4.6: Effect of clay mass fraction and NaCl concentration on the spatial
distribution of smectite clay and available macropore space in a 2-D cross
section of an ideal 3-D spherical packed bed (Finney,, 1970). The simulated
system is 3 by 3 mm. The grid resolution is 1020 by 1020 with a voxel size of
3 µm.

The system was simulated as a function of clay mass fraction (from 0 to 0.3)
and salinity (from 1 to 1000 mM NaCl). We chose this parameter space because,
as noted in the introduction, it represents the range over which micro and
macro-pores coexist and over which the coarse grains are load bearing.
Briefly, we initialized each set of simulations by populating the 3 by 3 mm
coarse grain structure with a uniform clay coating equilibrated at 1000 mM
NaCl until we reached the desired clay mass fraction. Clay parametrization was
consistent with the verification cases described in the previous sections. We
then applied a constant hydrostatic pressure gradient in the y direction and a
salinity step change across the sample and allowed the clay to swell until it
reached its equilibrium volume fraction based on the new salinity value. A
representative sample of the resulting steady-state configurations as a
function of clay mass fraction and salinity is shown in Fig. 4.6. Each sample
was put through several swelling-shrinking cycles to ensure consistent results
and to introduce heterogeneity approaching that present in natural systems.
Although not evident in the equations, heterogeneous clay distributions can
arise due to the combination of frictional forces imparted by the solid grains
and the clay’s non-Newtonian rheology. After we reached the steady state
configuration, the permeability of each sample was evaluated by integration of
the fluid velocities. Figures illustrating the predicted fluid velocity
distribution in a region containing both macropores and microporous clay are
shown in Fig. 4.7.

Figure 4.7: Steady-state fluid velocity profiles with a 10 % clay mass
fraction at different salinity values. Thin black lines surrounding the solid
grains represent the boundary between the clay matrix (within the lines) and
the free fluid (outside the lines). At high salt concentrations (A), flow
within the clay (greenish shades) does not control the overall permeability of
the medium and can be considered negligible when compared to flow around the
clay. Conversely, at lower salt concentrations and higher clay volumes (C),
flow is controlled by the internal permeability of the clay matrix due to the
absence of percolating flow paths through the macropores. The impact on the
overall permeability of the system is shown in Fig. 4.8.

Figure 4.8A shows that permeability varies by ten orders of magnitude within
the studied parameter space, with the largest permeability reduction (six
orders of magnitude) coming from changing the salinity at clay mass fractions
near 15%. This large permeability range reflects the very low permeability of
pure Na-smectite and the ability of the clay to block the main preferential
fluid flow paths, effectively dividing the medium’s macropores into isolated
regions. Because of the clay’s ability to swell, the location of the
percolation threshold depends on both salinity and clay mass fraction.

Figure 4.8: (a) Permeability of the simplified sedimentary rock model in Fig.
4.6 as a function NaCl concentration for a range of smectite clay mass
fractions. The unprobed parameter space in the lower left corner represents
conditions where clay more than fills the available space between the coarse
grains (i.e., clay becomes load bearing). (b) Data from Fig. 4.8A plotted as a
function of clay matrix volume fraction. The permeability data collapse into a
single curve that can be approximated with an error function (Eqn. 4.7 with A
= 9.6, S = 1.38 $\times$ A, and $V^{*}_{clay}$ =0.25). Note that the figure
also includes data obtained using Eqns. 4.1 and 4.2 parameterized for Illite
instead of smectite. The smectite and Illite permeabilities coincide below the
percolation threshold (where permeability is controlled by flow through
macropores) and differ by up to two orders of magnitude above the threshold
(where flow through the clay matrix predominates). (c) Comparison between Eqn.
4.7 (solid lines) and experimental datasets on the permeability of three
different types of siliciclastic sedimentary rocks (Gräfe et al.,, 2017; Heap
et al.,, 2017; Revil & Cathles,, 1999) and a series of glass bead-smectite