# Multiphase Flow Modelling in Multiscale Porous Media: An Open-Sourced Micro-
Continuum Approach

Francisco J. Carrillo  Ian C. Bourg  Cyprien Soulaine  Department of Chemical
and Biological Engineering, Princeton University, Princeton, NJ, USA
Department of Civil and Environmental Engineering, Princeton University,
Princeton, NJ, USA  Princeton Environmental Institute, Princeton University,
Princeton, NJ, USA  Earth Sciences Institute of Orléans, Université d’Orléans,
CNRS, BRGM, Orléans, France  French Geological Survey, BRGM, Orléans, France

###### Abstract

An open-sourced multiphase Darcy-Brinkman approach is proposed to simulate
two-phase flow in hybrid systems containing both solid-free regions and porous
matrices. This micro-continuum model is rooted in elementary physics and
volume averaging principles, where a unique set of partial differential
equations is used to represent flow in both regions and scales. The crux of
the proposed model is that it tends asymptotically towards the Navier-Stokes
volume-of-fluid approach in solid-free regions and towards the multiphase
Darcy equations in porous regions. Unlike existing multiscale multiphase
solvers, it can match analytical predictions of capillary, relative
permeability, and gravitational effects at both the pore and Darcy scales.
Through its open-source implementation, hybridPorousInterFoam, the proposed
approach marks the extension of computational fluid dynamics (CFD) simulation
packages into porous multiscale, multiphase systems. The versatility of the
solver is illustrated using applications to two-phase flow in a fractured
porous matrix and wave interaction with a porous coastal barrier.

###### keywords:

porous media , multi-scale , multiphase , micro-continuum, fracture , coastal
barrier

††journal: Journal of Computational Physics

##  1 Introduction

Virtually all aspects of subsurface engineering for energy and environmental
applications require in-depth understanding of multiphase flow within
heterogeneous porous media. Examples include enhanced hydrocarbon recovery,
geologic carbon sequestration, nuclear waste storage, geothermal energy
production, seasonal storage of natural gas in geologic formations, and gas
hydrate formation in sediments (Lake et al., 2014; Li and Benson, 2015; Rocco
et al., 2017; Yin et al., 2018). In addition, multiphase fluid dynamics in
heterogeneous porous media play key roles in the natural fluxes of water and
carbon in soils and sediments (Hassanizadeh et al., 2002; Or et al., 2013;
Maxwell et al., 2014; Scandella et al., 2017) as well as in a variety of
engineering processes (Baber et al., 2016; Jabbari et al., 2016). One largely
unresolved challenge in the field is the inability to predict and characterize
multiphase flow physics within inherently multiscale structures, particularly
in systems that contain both porous and solid-free domains (Helmig et al.,
2013). Although this challenge is widely recognized, there is increased
urgency in addressing it because of the need to sequester billions of tons of
CO2 and to efficiently extract hydrocarbons without causing extensive
environmental damage.

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 and Tartakovsky, 2009;
Datta et al., 2014). This coupling is intrinsically multiscale, as inertial
and viscous forces dominate at large pores or fractures while capillary forces
and interfacial energetics dominate within smaller porous micro-structures.
The complex linkage between microscopic geometric heterogeneities and
macroscopic processes makes it necessary to consider scale-dependent processes
across porous media in order to create truly predictive models: from the scale
of microscopic interfaces ($\sim\mu$m), to the scale of pore-networks and lab
columns ($\sim$cm), all the way up to the field scale ($\sim$km).

There are almost as many definitions of “multiscale” as authors that invoke
this concept. Nevertheless, multiscale modelling can be sorted in three main
categories (Helmig et al., 2010): _(i)_ the multiscale homogenization
strategy, _(ii)_ multiscale algorithm approach, and _(iii)_ the multiphysics
approach. The first of these, the multiscale homogenization strategy, aims at
deriving large scale models rooted in elementary physical principles by using
homogenization techniques including volume averaging, mixture theory, and
asymptotic expansions (Whitaker, 1999; Standnes et al., 2017; Battiato et al.,
2019; Starnoni and Pokrajac, 2020). A prime example is the seminal work of
Whitaker (1986a), which demonstrates that Darcy’s law arises from the
integration of Stokes equation over a porous Representative Elementary Volume
(REV). These upscaling techniques usually uncouple each scale’s relevant
physics through the scale-separation hypothesis. This way, effective
coefficients in large-scale models can be used to describe fine-scale
phenomena and geometric features. Such parameters are commonly estimated by
using complementary fine-scale simulations on REVs or sub-grid models. The
second strategy, the multiscale algorithm approach, solves flow physics on
interconnected grids with different degrees of refinement. This way, each
grid’s refinement level can be tuned to fit its respective scale of interest.
A portion of these algorithms primarily focus on fine-scale solutions, and
thus, use multi-scale finite volume/element solvers to speed up convergence in
fine grids (Jenny et al., 2003, 2006; Efendiev and Hou, 2007). Conversely,
alternative algorithms focus on large-scale behaviors and only solve for
small-scale behavior when needed (Tomin and Lunati, 2013, 2016). The third
strategy, the multiphysics approach, uses domain decomposition to solve
different physics within each scale’s sub-domain. In this method, sub-domains
have their own independent set of governing equations and only interact with
each other through the implementation of appropriate boundary conditions (Sun
et al., 2012; Baber et al., 2016). A popular implementation of this approach
uses the Beavers and Joseph (1967) conditions to couple a porous domain
governed by Darcy’s Law with a domain governed by the Navier-Stokes equations.

Here, we will implement concepts from all three strategies to propose an
alternative solution to the multiscale challenge. To do so, we will rely on
the micro-continuum approach (Soulaine and Tchelepi, 2016), whereby a single
equation is used to handle flow and transport in systems where a large scale
solid-free domain coexists with a small-scale porous domain (Figure 1). In the
case of single-phase flow and transport, this approach generally relies on the
well-known Darcy-Brinkman (DB) equation –also referred to as Darcy-Brinkman-
Stokes (DBS) equation– (Brinkman, 1947) that arises from volume averaging the
Stokes (or Navier-Stokes) equations in a control volume that contains both
fluids and solids (Vafai and Tien, 1981; Hsu and 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 1A) where the drag force term vanishes and
the DB equation turns into the Stokes (or Navier-Stokes) equation. In porous
regions (see Figure 1C), 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 and Nader, 1974), paving the path to hybrid
scale modeling (see Figure 1B). 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 and Joseph, 1967;
Neale and Nader, 1974).

Figure 1: 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 ability of the DB equation to handle two scales simultaneously has been
used to solve fluid flow in three-dimensional images of rock samples that
contain unresolved sub-voxel porosity (Knackstedt et al., 2006; Apourvari and
Arns, 2014; Scheibe et al., 2015; Soulaine et al., 2016; Guo et al., 2018;
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 and Tchelepi, 2016; Tomin and 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 and Tchelepi, 2016). This
approach tends to a full Navier-Stokes representation of the flow physics at
the pore scale and, hence, 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., 2020). Another example, presented in
Carrillo and Bourg (2019), leveraged this framework to create a Darcy-
Brinkman-Biot approach capable of predicting the coupled hydrology and
mechanics of soft porous media such as clays and elastic membranes.

The micro-continuum framework outlined above was limited, until recently, to
single-phase flow (Soulaine and Tchelepi, 2016). Horgue et al. (2014) and
later Soulaine et al. (2018) proposed the first two-phase micro-continuum
model by combining a two-phase variant of the DB equation with the volume of
fluid (VOF) approach (Hirt and Nichols, 1981) to two-phase flow in solid-free
regions. This formulation enabled multiphase flow in solid-free regions with
imposed wettability conditions at the solid surface while describing
microporous regions as impervious, fully-saturated porous domains (Soulaine et
al., 2018). Soulaine et al. (2019) later refined their formulation to enable
two-phase flow at both the pore and continuum scales, but with simplified flow
physics in the porous domain; in particular, the model could not describe the
interplay of gravity and capillarity effects within the microporous matrix.

In this paper, we expand upon Soulaine et al. (2019) to propose a fully
realized multiscale solver for two-phase flow in porous media rooted in
elementary physical principles and rigorously derived using the method of the
volume averaging (Whitaker, 1999). We show that there exists a single set of
partial differential equations that can be applied in pore, continuum, or
hybrid scale representation of multiphase flow in porous media. Particular
attention is paid to the rigorous derivation of gravity and capillary effects
in the porous domain. The resulting two-phase micro-continuum framework is
verified using a series of test cases where reference solutions exist. We show
that the multiscale solver converges to the standard Darcy scale solutions
(Buckley-Leverett, capillary-gravity equilibrium, drainage in a heterogeneous
reservoir) when used at the continuum scale 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 open-source
solver (hybridPorousInterFoam) accompanying the present article.

The paper is organized as follows. In Section 2, the multi-scale governing
equations are rigorously derived using the method of volume averaging. Multi-
scale parameters are then defined by asymptotic matching to the two-phase
Navier-Stokes and Darcy equations. In Section 3, we describe the numerical
algorithm used to solve the problem (governing equations, constitutive
relations, and boundary conditions) and present its numerical implementation
as an open-source simulation platform. In Section 4, we present the model
verification at the pore and continuum scales. In Section 5, we illustrate the
versatility of the proposed framework by describing two hybrid scale
applications: wave propagation in a coastal barrier and two-phase flow in a
fractured porous matrix. We close with a summary and conclusions.

##  2 Mathematical model

In this section, we derive the micro-continuum equations for two-phase flow.
First, we consider the conservation laws for multi-phase systems in the
continuous physical space. Then, the micro-continuum equations are formed by
volume averaging the continuous equations over each volume of an Eulerian
grid. Finally, information below the size of the grid cell (fluid-fluid
interface location and micro-structure geometry) is modelled with closure of
the multiscale parameters.

###  2.1 Governing equation in the continuous physical space

This section presents the basic hydrodynamic laws that govern multiphase flow
at the pore scale. The domain is decomposed into three disjoint subsets: a
solid phase $V_{s}$, a wetting liquid phase $V_{l}$, and a non-wetting gas
phase $V_{g}$ which is separated from $V_{l}$ by the interface $A_{lg}$ (see
Figure 2A). Although the fluids are referred to as liquid and gas (or wetting
and non-wetting), the derivation and resulting model are valid for any
incompressible, immiscible fluid pair including liquid-liquid and liquid-gas
systems.

Each fluid phase is assumed to be Newtonian and incompressible. Therefore,
mass conservation in each phase dictates

$\displaystyle\nabla\cdot\boldsymbol{v}_{i}=0$ $\displaystyle\text{in
$V_{i}$,}\quad i=l,g,$ (1)

where $\boldsymbol{v}_{i}$ is the velocity of phase $i$. Mass conservation at
the fluid-fluid interface yields

$\rho_{l}\left(\boldsymbol{v}_{l}-\boldsymbol{w}\right)\cdot\boldsymbol{n}_{lg}=\rho_{g}\left(\boldsymbol{v}_{g}-\boldsymbol{w}\right)\cdot\boldsymbol{n}_{lg}\hskip
14.22636pt\text{at $A_{ij}$},$ (2)

where $\rho_{i}$ is the density of phase $i$, $\boldsymbol{w}$ is the velocity
of the interface, and $\boldsymbol{n}_{lg}$ is the normal vector to the fluid-
fluid interface pointing from the wetting to the non-wetting phase. In the
absence of phase change,
$\boldsymbol{v}_{l}=\boldsymbol{v}_{g}=\boldsymbol{w}$ at the fluid/fluid
interface.

Momentum conservation in each fluid yields

$\displaystyle 0=-\nabla
p_{i}+\rho_{i}\boldsymbol{g}+\nabla\cdot\mathsf{S}_{i}$ $\displaystyle\text{in
$V_{i}$,}\quad i=l,g,$ (3)

where $\boldsymbol{g}$ is the gravity vector,
$\mathsf{S}_{i}=\mu_{i}\left(\nabla\boldsymbol{v}_{i}+\nabla\boldsymbol{v}_{i}^{T}\right)$
is the viscous stress tensor, and $p_{i}$ and $\mu_{i}$ are the pressure and
viscosity of phase $i$, respectively. In Eq. (3), the inertia terms have been
neglected and the momentum balance is described using the Stokes equation.
This simplification is common in models of subsurface fluid flow, where flow
rates are usually very low (Bear, 1972). The Stokes equation is adopted for
simplicity in the derivation of the micro-continuum momentum equation. For
simplicity, inertial effects will be integrated at the end of the derivation
based on the full Navier-Stokes equation.

Finally, momentum conservation at the fluid-fluid interface yields

$\left[p_{l}\mathsf{I}-\mathsf{S}_{l}\right]\cdot\boldsymbol{n}_{lg}=\left[p_{g}\mathsf{I}-\mathsf{S}_{g}\right]\cdot\boldsymbol{n}_{lg}+\sigma\kappa\boldsymbol{n}_{lg}\hskip
14.22636pt\text{at $A_{lg}$},$ (4)

where $\mathsf{I}$ is the unity tensor, $\sigma$ is the fluid-fluid
interfacial tension, and $\kappa=\nabla\cdot\boldsymbol{n}_{lg}$ is the
interface curvature.

###  2.2 Volume averaging: derivation of a single-field formulation

The mathematical model introduced in the former section is defined on a
continuous physical domain. Common computational procedures solve this system
of equations by discretizing the continuous domain into an ensemble of subset
volumes by using the Finite Volume Method (FVM) (Patankar, 1980). In the FVM
framework, all the physical variables are averaged over each discrete volume.
The averaging process and the discretization refinement level dictate that the
control volume can contain the following: one fluid, two fluids, one fluid and
a solid phase, or two fluids and a solid phase. Features with characteristic
length scales below that of the averaging volume (e.g., the geometry of solid-
fluid and fluid-fluid interfaces and the forces exerted onto them) must be
described using sub-grid scale representations. In this section, we use volume
averaging theorems to identify the form of the multiphase micro-continuum
equations.

Figure 2: Distribution of the fluid phases in (a) the continuous physical
domain, (b) the discrete Eulerian grid.

##### Volume averaging and single-field variables

In the FVM, the partial differential equations that describe conservation
laws, Eqs. (1) and (3), are transformed into discrete algebraic equations by
integrating them over each discrete volume $V$. This operation is carried out
using the volume averaging operator,

$\overline{\beta_{i}}=\frac{1}{V}\int_{V_{i}}\beta_{i}dV,$ (5)

where $\beta_{i}$ is a function defined in $V_{i}$ ($i=l,g$). As in standard
volume averaging theory, we also define a phase averaging operator,

$\overline{\beta_{i}}^{i}=\frac{1}{V_{i}}\int_{V_{i}}\beta_{i}dV.$ (6)

The averages defined by Eqs. (5) and (6) are related through the porosity
field $\phi$ and saturation field $\alpha_{l}$. The porosity field $\phi$ is
defined as $(V_{l}+V_{g})/V$, i.e., the volume occupied by both fluids divided
by the control volume $V$, such that:

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

The porosity field is the cornerstone of micro-continuum methods because it
delineates porous ($0<\phi<1$) and solid-free regions ($\phi=1$). It is
intrinsically related to the resolution of the simulation as illustrated in
Fig. 1. For example, in image-based flow simulations, the control volume size
corresponds to the imaging instrument resolution and the porosity field
obtained from the gray-scale is used to model sub-voxel micro-structures
(Apourvari and Arns, 2014; Soulaine et al., 2016; Scheibe et al., 2015; Guo et
al., 2018; Singh, 2019; Abu-Al-Saud et al., 2020). By construction of micro-
continuum models, all cells must have non-zero porosity (Soulaine and
Tchelepi, 2016). Hence, a pure solid phase ($\phi=0$) in the micro-continuum
framework is described instead as a very low-porosity, very low-permeability
domain ($\phi\approx 0$).

The saturation field $\alpha_{l}$ is defined as $V_{l}/(V_{l}+V_{g})$, i.e.,
the volume of liquid divided by the volume occupied by both fluids within the
control volume, such that

$\alpha_{l}=\begin{cases}0,&\textnormal{in regions saturated with gas,}\\\
\left]0;1\right[,&\textnormal{in unsaturated regions,}\\\ 1,&\textnormal{in
regions saturated with liquid.}\\\ \end{cases}$ (8)

A saturation field $\alpha_{l}$ such as that described by Eq. (8) is used in
continuum scale simulations of multiphase flow in porous media (where it
represents actual saturation) and in pore scale simulations of multiphase flow
in solid-free regions that rely on the VOF representation (where it is used to
track the evolution of the immiscible fluid-fluid interface). The relationship
$\alpha_{l}+\alpha_{g}=1$ is always valid and $\alpha_{g}$ is deduced from the
knowledge of $\alpha_{l}$. The averaging operators defined by Eqs. (5) and (6)
are related by $\overline{\beta_{i}}=\phi\alpha_{i}\overline{\beta_{i}}^{i}$
($i=l,g$).

The two-phase micro-continuum approach relies on single-field variables, i.e.,
unique fluid pressure and velocity fields that are defined throughout the
entire grid regardless of the nature of the phases that occupy the cells. The
single-field pressure $\overline{p}$ and velocity $\overline{\boldsymbol{v}}$
are defined as weighted sums of the pressure and Darcy velocity in each fluid
phase:

$\overline{p}=\alpha_{l}\overline{p}_{l}^{l}+\alpha_{g}\overline{p}_{g}^{g},$
(9)

and

$\overline{\boldsymbol{v}}=\phi\left[\alpha_{l}\overline{\boldsymbol{v}_{l}}^{l}+\alpha_{g}\overline{\boldsymbol{v}_{g}}^{g}\right],$
(10)

respectively. We note that the use of porosity-weighted values in Eq. (10)
yields a single-field velocity equal to the sum of the filtration (Darcy)
velocities in each phase,
$\overline{\boldsymbol{v}}=\overline{\boldsymbol{v}_{l}}+\overline{\boldsymbol{v}_{g}}$.

The governing equations solving for $\overline{p}$ and
$\overline{\boldsymbol{v}}$ are obtained using a two-step strategy. First, the
volume averaging operator, Eq. (5), is applied to the continuity equations,
Eq. (1), and to the momentum equations, Eq. (3), leading to two pairs of
partial differential equations solving for $\alpha_{i}$,
$\overline{\boldsymbol{v}_{i}}^{i}$, and $\overline{p}_{i}^{i}$ ($i=l,g$).
Second, pairs of phase-averaged equations are added to each other to form the
governing equations for the single-field variables. In the averaging process,
the volume averaging operator is applied to spatial differential operators
(gradient and divergence). This operation is not straightforward because
integrals and derivatives can not be interchanged in volumes that contain
interfaces including fluid-fluid and fluid-solid interfaces. This is achieved
using the spatial volume averaging theorems (Howes and Whitaker, 1985;
Whitaker, 1999),

$\displaystyle\overline{\nabla\beta_{i}}$
$\displaystyle=\nabla\overline{\beta_{i}}+\frac{1}{V}\int_{A_{ij}}\beta_{i}\boldsymbol{n}_{ij}dA+\frac{1}{V}\int_{A_{is}}\beta_{i}\boldsymbol{n}_{is}dA,$
(11) $\displaystyle\overline{\nabla\cdot\boldsymbol{\beta}_{i}}$
$\displaystyle=\nabla\cdot\overline{\boldsymbol{\beta}_{i}}+\frac{1}{V}\int_{A_{ij}}\boldsymbol{\beta}_{i}\cdot\boldsymbol{n}_{ij}dA+\frac{1}{V}\int_{A_{is}}\boldsymbol{\beta}_{i}\cdot\boldsymbol{n}_{is}dA,$

where $A_{ij}$ is the surface area between the two fluids, $A_{is}$ is the
surface area between fluid $i$ and the solid phase, $\boldsymbol{n}_{ij}$ is
the normal vector at the fluid-fluid interface pointing from $i$ to $j$, and
$\boldsymbol{n}_{is}$ is the normal vector at the solid surface pointing from
the fluid to the solid. The surface integral terms in these equations
transform the boundary conditions at the discontinuity between the fluid
phases and at the solid surface into body forces. In others words, the
interfacial conditions are included directly in the partial differential
equations that describe the conservation laws in the Eulerian grid.

##### Mass balance and saturation equations

The application of the volume averaging theorem, Eq. (11), along with the
continuity equations, Eq. (1), yields (Whitaker, 1986b):

$\frac{\partial\phi\alpha_{i}}{\partial
t}+\nabla.\overline{\boldsymbol{v}_{i}}=0,\quad i=g,l.$ (12)

The two-phase micro-continuum framework developed in this paper consists of a
set of partial differential equations that solve for the single-field
variables $\overline{\boldsymbol{v}}$, $\overline{p}$, and $\alpha_{l}$.
Because the volume-averaged continuity equations, Eq. (12), involve averaged
phase velocities $\overline{\boldsymbol{v}_{i}}$ they must be transformed into
equations in terms of the micro-continuum single-field variables, namely a
total fluid conservation equation and a saturation equation.

The total fluid conservation equation is obtained by summing the two
continuity equations and assuming that the porous structure is immobile with
time, such that:

$\nabla\cdot\overline{\boldsymbol{v}}=0.$ (13)

Equation (13) is a divergence-free velocity that is commonly used together
with the momentum equation to derive the pressure equation.

The saturation equation is obtained by first introducing the concept of
relative velocity:

$\overline{\boldsymbol{v}_{r}}=\left(\overline{\boldsymbol{v}_{l}}^{l}-\overline{\boldsymbol{v}_{g}}^{g}\right).$
(14)

From the definitions of single-field and relative velocities, we can show that
$\overline{\boldsymbol{v}_{l}}^{l}=\phi^{-1}\overline{\boldsymbol{v}}+\alpha_{g}\overline{\boldsymbol{v}_{r}}$.
Because
$\overline{\boldsymbol{v}_{l}}=\phi\alpha_{l}\overline{\boldsymbol{v}_{l}}^{l}$,
the saturation equation can be expressed as:

$\frac{\partial\phi\alpha_{l}}{\partial
t}+\nabla\cdot\left(\alpha_{l}\overline{\boldsymbol{v}}\right)+\nabla\cdot\left(\phi\alpha_{l}\alpha_{g}\overline{\boldsymbol{v}_{r}}\right)=0.$
(15)

In equation (15), the wetting phase saturation $\alpha_{l}$ is advected by the
single-field velocity $\overline{\boldsymbol{v}}$. The third term on the left-
hand side is an additional convection term involving the relative velocity
$\overline{\boldsymbol{v}_{r}}$. The saturation equation, Eq. (15), is exact,
i.e., it is derived from elementary physical principles without any
assumptions. However, there is no conservation law to solve for
$\overline{\boldsymbol{v}_{r}}$ and this term must be closed. In the
forthcoming discussion, we will see that different descriptions of
$\overline{\boldsymbol{v}_{r}}$ are derived for solid-free ($\phi=1$) and
porous regions ($0\leq\phi<1$). In the first case, the convection term
involving the relative velocity serves to compress the fluid-fluid interface
and ensures a sharp transition between the immiscible phases. In the second
case, $\overline{\boldsymbol{v}_{r}}$ is closed by matching Eq. (15) to the
standard saturation equation used in multi-phase Darcy flow solvers.

##### Momentum equation

A similar procedure is used to form the multiscale momentum equation. First,
the volume averaged equations are derived for each fluid. Then, the two
resulting equations are combined to form the single-field conservation law.

The application of the volume averaging theorem, Eq. (11), to the Stokes
momentum conservation equation for fluid $i$, Eq. (3), yields (Whitaker,
1986b; Lasseux et al., 1996; Ishii and Hibiki, 2011):

$0=-\nabla\left(\phi\alpha_{i}\overline{p_{i}}^{i}\right)+\phi\alpha_{i}\rho_{i}\boldsymbol{g}+\nabla\cdot\left(\phi\alpha_{i}\overline{\mathsf{S}_{i}}^{i}\right)+\boldsymbol{D}_{is}+\boldsymbol{D}_{ij},\quad
i=g,l,$ (16)

where the two last terms on the right-hand side,

$\boldsymbol{D}_{ik}=\frac{1}{V}\int_{A_{ik}}\boldsymbol{n}_{ik}\cdot\left(-p_{i}\mathsf{I}+\mathsf{S}_{i}\right)dA,$
(17)

are the drag forces exerted by phase $k$ on phase $i$. In short,
$\boldsymbol{D}_{is}$ reflects the friction of fluid $i$ on the solid surface
and $\boldsymbol{D}_{ij}$ reflects interfacial shear between the two fluids.
These terms accounts for shear that occurs at scales below that of a control
volume; therefore, the description of these terms must differ depending on
whether the computational cells contain a porous solid structure
($0\leq\phi<1$) or fluids only ($\phi=1$). These drag forces will be derived
later on in Section 2.3; for the time being, they are kept in their integral
forms.

The sum of the two phase-averaged momentum conservation equations yields:

$0=-\nabla\left(\phi\overline{p}\right)+\phi\rho\boldsymbol{g}+\nabla\cdot\left(\phi\overline{\mathsf{S}}\right)+\boldsymbol{D}_{ls}+\boldsymbol{D}_{gs}+\boldsymbol{D}_{lg}+\boldsymbol{D}_{gl},$
(18)

where $\overline{\mathsf{S}}$ is the single-field shear stress
$\left[\overline{\mathsf{S}}=\mu\left(\nabla\overline{\boldsymbol{v}}+\nabla\overline{\boldsymbol{v}}^{T}\right)\right]$
and $\mu$ is the average fluid viscosity
$\left[\mu=\alpha_{l}\mu_{l}+\alpha_{g}\mu_{g}\right]$. 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: a
drag force $\mu k^{-1}\overline{\boldsymbol{v}}$ and a surface tension force
$\boldsymbol{F}_{c}$:

$-\mu
k^{-1}\overline{\boldsymbol{v}}+\boldsymbol{F}_{c}=\phi^{-1}\left[\boldsymbol{D}_{ls}+\boldsymbol{D}_{gs}+\boldsymbol{D}_{lg}+\boldsymbol{D}_{gl}\right].$
(19)

Eventually, if the porous structure is immobile, the porosity $\phi$ can be
removed from the derivatives and the multiscale single-field momentum equation
becomes:

$0=-\nabla\overline{p}+\rho\boldsymbol{g}+\nabla\cdot\overline{\mathsf{S}}-\mu
k^{-1}\overline{\boldsymbol{v}}+\boldsymbol{F}_{c}.$ (20)

##### Summary of the derivation

The single-field micro-continuum model for incompressible, immiscible two-
phase flow in a rigid porous medium, derived above using volume averaging
theory, consists of a set of three partial differential equations, namely a
total mass balance equation, Eq. (13), a saturation equation, Eq. (15), and a
momentum equation, Eq. (20), that can be solved for the single-field pressure
$\overline{p}$, the single-field velocity $\overline{\boldsymbol{v}}$, and the
wetting fluid saturation $\alpha_{l}$:

$\displaystyle\nabla\cdot\overline{\boldsymbol{v}}=0,$
$\displaystyle\frac{\partial\phi\alpha_{l}}{\partial
t}+\nabla\cdot\left(\alpha_{l}\overline{\boldsymbol{v}}\right)+\nabla\cdot\left(\phi\alpha_{l}\alpha_{g}\overline{\boldsymbol{v}_{r}}\right)=0,$
(21)
$\displaystyle\frac{1}{\phi}\left(\frac{\partial\rho\overline{\boldsymbol{v}}}{\partial
t}+\nabla\cdot\left(\frac{\rho}{\phi}\bar{\boldsymbol{v}}\bar{\boldsymbol{v}}\right)\right)=-\nabla\overline{p}+\rho\boldsymbol{g}+\nabla\cdot\overline{\mathsf{S}}-\mu
k^{-1}\overline{\boldsymbol{v}}+\boldsymbol{F}_{c}.$

In Eq. (21), the momentum equation has been modified from Eq. (20) to include
the inertial effects that were neglected above for clarity. The derivation
with these inertial effects follows the same averaging procedure as described
above, starting from the Navier-Stokes (rather than Stokes) equation (Vafai
and Tien, 1981; Hsu and Cheng, 1990; Bousquet-Melou et al., 2002; Goyeau et
al., 2003). The numerical implementation described in Section 3 for solving
Eq. (21) accounts for the inertial effects.

The set of equations presented above is valid throughout the computational
domain regardless of the content of a cell. This characteristic is a
fundamental aspect of our multiscale solver. It means that the same equations
for multiphase flow and transport can be used in both solid-free and porous
regions, unlike in the case of multi-physics solvers that involve mortars (Sun
et al., 2012; Baber et al., 2016). This feature allows the proposed solver to
be applied in media where the pore space is fully resolved and flow is
described using the Navier-Stokes equations (pore scale modeling), in media
where pores are not resolved and flow is described using Darcy’s law
(continuum scale modeling), and in intermediate situations that include both
fully resolved solid-free regions and porous regions (hybrid scale modelling)
as illustrated in Figure 1.

A critical feature of the multiscale solver developed in this paper is that it
tends asymptotically to the solution of the two-phase Navier-Stokes equations
when used as a pore scale model and to the solution of the two-phase Darcy
equations when used as a continuum scale model. This is achieved by defining
the relative velocity $\overline{\boldsymbol{v}_{r}}$, the drag force $\mu
k^{-1}\overline{\boldsymbol{v}}$, and the surface tension force
$\boldsymbol{F}_{c}$. These terms are referred to as multiscale parameters
because they describe sub-grid scale information such as the location of the
fluid-fluid interface and the hydrodynamic impact of the porous micro-
structure. They have a different meaning and a different formulation depending
on whether the computational grid blocks contain solid material or not.

###  2.3 Closure and multi-scale parameters

In the following, we show how the multiscale parameters
$\overline{\boldsymbol{v}_{r}}$, $\mu k^{-1}$, and $\boldsymbol{F}_{c}$ can be
derived by matching Eq. (21) to its two desired asymptotic models: in the pore
scale limit, the algebraic Volume-of-Fluid method; in the continuum scale
limit, the multiphase form of Darcy’s law.

##### Algebraic Volume-of-Fluid model in the pore scale limit

In CFD, the Volume of Fluid (VOF) method (Hirt and Nichols, 1981) is a
standard approach to track the interface movement of two immiscible fluids in
a fixed Eulerian grid. This approach is known to approximate the solution of
the physical problem, Eqs. (1)-(4), using a Finite-Volume grid. In the VOF
approach, a phase indicator representing the volume of fluid in each grid
block is used to track the distribution of the fluid phases in the
computational domain as illustrated in the upper part of Figure 2B. This phase
indicator has the same form as the saturation field $\alpha_{l}$ defined in
the two-phase multiscale micro-continuum model. In cells saturated by the
wetting phase, $\alpha_{l}=1$. In cells that contain the non-wetting phase
only, $\alpha_{l}=0$. Finally, $0<\alpha_{l}<1$ in cells containing the
immiscible interface between both fluids. The VOF approach relies on a single-
field formulation of the Navier-Stokes equations to compute the two-phase
flow. If a cell of the Finite-Volume grid is considered as a control volume,
then all the derivation introduced in the previous section can be used to
derive the VOF momentum, mass balance, and saturation equations (Maes and
Soulaine, 2019).

In the standard VOF approaches, the cells do not contain solid ($\phi=1$). The
mass balance and saturation equations, Eqs. (13) and (15), remain, therefore,
unchanged. The saturation equation with $\phi=1$ is the equation used in
algebraic VOF solvers such as _interFoam_ , the VOF solver of the open-source
CFD code OpenFOAM®. There, the relative velocity
$\overline{\boldsymbol{v}_{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 $\overline{\boldsymbol{v}}$:

$\overline{\boldsymbol{v}_{r}}=C_{\alpha}\max\left(\left|\overline{\boldsymbol{v}}\right|\right)\overline{\boldsymbol{n}_{lg}},$
(22)

where $C_{\alpha}$ is a model parameter used to control the compression of the
interface and $\overline{\boldsymbol{n}_{lg}}$ is 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}_{lg}$ is
computed by using the gradient of the phase indicator function $\alpha_{l}$.
The relation between these two vectors can be obtained by applying Eq. (11) to
the liquid phase indicator function $\mathds{1}_{l}$ (a function equal to 1 in
$V_{l}$ and 0 elsewhere) in solid-free regions such that (Quintard and
Whitaker, 1994),

$\nabla\alpha_{l}=-\frac{1}{V}\int_{A_{lg}}\boldsymbol{n}_{lg}dA.$ (23)

Therefore,

$\overline{\boldsymbol{n}_{lg}}=-\frac{\nabla\alpha_{l}}{\left|\nabla\alpha_{l}\right|},$
(24)

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}_{ls}=\boldsymbol{D}_{gs}=0$. Therefore, the Darcy
term in the momentum equation vanishes:

$\mu k^{-1}\overline{\boldsymbol{v}}=0.$ (25)

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

$\boldsymbol{D}_{lg}+\boldsymbol{D}_{gl}=\phi\boldsymbol{F}_{c}=\frac{1}{V}\int_{A_{lg}}\boldsymbol{n}_{lg}\cdot\sigma\kappa
dA.$ (26)

This equation cannot be used directly, because the terms under the volume
integral require 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}_{lg}}$. Because
$\overline{\kappa}$ and $\sigma$ are constant within a control volume, they
can be pulled out of the integral in Eq. (26) to obtain (after applying Eq.
(23)) the so-called Continuum Surface Force (CSF) formulation (Brackbill et
al., 1992):

$\boldsymbol{F}_{c}=\phi^{-1}\sigma\nabla\cdot\left(\frac{\nabla\alpha_{l}}{\left|\nabla\alpha_{l}\right|}\right)\nabla\alpha_{l}.$
(27)

##### Standard two-phase Darcy model in the continuum scale limit

In this section, 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 (Muskat, 1949; Miller et al., 1998; Pinder and Gray,
2008). The model can be derived by applying the volume averaging operators on
a Representative Elementary Volume of the porous structure (Whitaker, 1986b;
Lasseux et al., 1996), along the same lines of the derivation in Section 2.2.
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{v}_{i}}$ with $i=g,l$) with
separate pressure fields ($\overline{p_{i}}$ with $i=g,l$).

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

$\frac{\partial\phi\alpha_{l}}{\partial
t}+\nabla.\overline{\boldsymbol{v}_{l}}=0,$ (28)

a mass balance equation,

$\nabla.\overline{\boldsymbol{v}}=0,$ (29)

and two momentum balance equations, one for each phase,

$\displaystyle\overline{\boldsymbol{v}_{i}}=\phi\alpha_{i}\overline{\boldsymbol{v}_{i}}^{i}$
$\displaystyle=-\frac{k_{0}k_{r,i}}{\mu_{i}}\left(\nabla\overline{p_{i}}^{i}-\rho_{i}\boldsymbol{g}\right),\quad
i=g,l,$
$\displaystyle=-M_{i}\left(\nabla\overline{p_{i}}^{i}-\rho_{i}\boldsymbol{g}\right),\quad
i=g,l,$ (30)

these can also be written as,

$0=-\nabla\overline{p_{i}}^{i}+\rho_{i}\boldsymbol{g}-M_{i}^{-1}\overline{\boldsymbol{v}_{i}},\quad
i=g,l,$ (31)

where $k_{0}$ is the absolute permeability of the porous structure, $k_{r,l}$
and $k_{r,g}$ are the relative permeabilities with respect to each fluid
(classically represented here as functions of water saturation; more complex
formulations exist that account for viscous coupling between the two fluids or
for the Klinkenberg effect in the gas phase (Standnes et al., 2017; Picchi and
Battiato, 2018; Guo et al., 2018)), and $M_{i}=\frac{k_{0}k_{r,i}}{\mu_{i}}$
are the fluid mobilities. These momentum equations arise from further
simplification of the volume averaged Stokes equations, Eq. (16), where the
drag forces are combined and described as a Darcy term. Moreover, by relying
on scale separation arguments, Whitaker (1986a) showed that the viscous
dissipative term,
$\nabla\cdot\left(\alpha_{i}\overline{\mathsf{S}}_{i}^{i}\right)$, is
negligible in comparison to the drag forces. This feature is a fundamental
aspect of the multiscale micro-continuum framework 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 and Nader, 1974).

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

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

For simplicity, we follow the classical approximation that $p_{c}$ depends
only on saturation (Leverett, 1940; Brooks and Corey, 1964; Van Genutchen,
1980). Alternative formulations have been proposed to account for observed
disequilibrium and hysteretic effects in the macroscopic capillary pressure
(Hassanizadeh et al., 2002; Gray et al., 2015; Li et al., 2018; Miller et al.,
2019; Starnoni and Pokrajac, 2020).

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
$\overline{\boldsymbol{v}}_{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 $\overline{p}$, Eq. (9), along with
the definition of capillary pressure, Eq. (32), results in:

$\displaystyle\nabla\overline{p_{l}}^{l}$
$\displaystyle=\nabla\overline{p}-\nabla\left(\alpha_{g}p_{c}\right),$
$\displaystyle\nabla\overline{p_{g}}^{g}$
$\displaystyle=\nabla\overline{p}+\nabla\left(\alpha_{l}p_{c}\right).$ (33)

Based on the equations presented above, the multi-phase Darcy model implies
the following expression for $\overline{\boldsymbol{v}}_{r}$:

$\displaystyle\overline{\boldsymbol{v}}_{r}$
$\displaystyle=\left(\overline{\boldsymbol{v}_{l}}^{l}-\overline{\boldsymbol{v}_{g}}^{g}\right),$
$\displaystyle=-\frac{M_{l}}{\phi\alpha_{l}}\left(\nabla\overline{p_{l}}^{l}-\rho_{l}\boldsymbol{g}\right)+\frac{M_{g}}{\phi\alpha_{g}}\left(\nabla\overline{p_{g}}^{g}-\rho_{g}\boldsymbol{g}\right),$
$\displaystyle=\phi^{-1}\left[-\frac{M_{l}}{\alpha_{l}}\nabla\overline{p_{l}}^{l}+\frac{M_{g}}{\alpha_{g}}\nabla\overline{p_{g}}^{g}+\left(\rho_{l}\frac{M_{l}}{\alpha_{l}}-\rho_{g}\frac{M_{g}}{\alpha_{g}}\right)\boldsymbol{g}\right],$
$\displaystyle=\phi^{-1}\left[-\left(\frac{M_{l}}{\alpha_{l}}-\frac{M_{g}}{\alpha_{g}}\right)\nabla\overline{p}+\left(\rho_{l}\frac{M_{l}}{\alpha_{l}}-\rho_{g}\frac{M_{g}}{\alpha_{g}}\right)\boldsymbol{g}+\frac{M_{l}}{\alpha_{l}}\nabla\left(\alpha_{g}pc\right)+\frac{M_{g}}{\alpha_{g}}\nabla\left(\alpha_{l}pc\right)\right],$
$\displaystyle=\phi^{-1}\left[-\left(\frac{M_{l}}{\alpha_{l}}-\frac{M_{g}}{\alpha_{g}}\right)\nabla\overline{p}+\left(\rho_{l}\frac{M_{l}}{\alpha_{l}}-\rho_{g}\frac{M_{g}}{\alpha_{g}}\right)\boldsymbol{g}+\left(M_{l}\frac{\alpha_{g}}{\alpha_{l}}+M_{g}\frac{\alpha_{l}}{\alpha_{g}}\right)\nabla
p_{c}-\left(\frac{M_{l}}{\alpha_{l}}-\frac{M_{g}}{\alpha_{g}}\right)p_{c}\nabla\alpha_{l}\right].$
(34)

In Soulaine et al. (2019) only the first term of the right-hand side was
considered, such that the model could not account for gravity or capillary
effects in the porous domain. The comprehensive formulation presented in Eq.
(34) overcomes these limitations.

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

$\displaystyle\overline{\boldsymbol{v}}$
$\displaystyle=\overline{\boldsymbol{v}_{l}}+\overline{\boldsymbol{v}_{g}},$
$\displaystyle=-M_{g}\nabla\overline{p_{g}}^{g}-M_{l}\nabla\overline{p_{l}}^{l}+\left(\rho_{g}M_{g}+\rho_{l}M_{l}\right)\boldsymbol{g},$
(35)
$\displaystyle=-\left(M_{g}+M_{l}\right)\nabla\overline{p}+\left(\rho_{g}M_{g}+\rho_{l}M_{l}\right)\boldsymbol{g}+\left[M_{l}\nabla\left(\alpha_{g}p_{c}\right)-M_{g}\nabla\left(\alpha_{l}p_{c}\right)\right],$

The previous equation can be recast into:

$0=-\nabla\overline{p}+\rho^{*}\boldsymbol{g}-M^{-1}\overline{\boldsymbol{v}}+M^{-1}\left[M_{l}\nabla\left(\alpha_{g}p_{c}\right)-M_{g}\nabla\left(\alpha_{l}p_{c}\right)\right],$
(36)

where $M=M_{l}+M_{g}$ is the total mobility and
$\rho^{*}=\left(\rho_{l}M_{l}+\rho_{g}M_{g}\right)/\left(M_{l}+M_{g}\right)$
is a mobility-weighted average fluid density. This single-field two-phase
Darcy equation matches the two-phase micro-continuum momentum, Eq. (20), if
the drag coefficient and the capillary force equal

$\mu
k^{-1}=M^{-1}=k_{0}^{-1}\left(\frac{\mu_{l}}{k_{rl}}+\frac{\mu_{g}}{k_{rg}}\right)^{-1},$
(37)

and

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

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

Finally, we note that in Eq. (36), the single-field fluid density $\rho^{*}$
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 Appendix B, a strictly
equivalent solution can be derived where $\rho^{*}$ is replaced by $\rho$ in
Eq. (36) and the capillary force expression is replaced by:

$\boldsymbol{F}_{c}=M^{-1}\left(M_{l}\alpha_{g}-M_{g}\alpha_{l}\right)[(\rho_{l}-\rho_{g})\boldsymbol{g}+\nabla
p_{c}]-p_{c}\nabla\alpha_{l}$ (39)

##### Condition at the interface between a clear fluid region and a porous
domain

The multiscale parameters are derived above for solid-free and porous regions.
In hybrid scale simulations, however, both regions can exist concomitantly in
the computational grid (see Fig. 1B). Here, a condition at the interface
between porous and solid-free domains is proposed.

First, we note that for single-phase flow the DB equation captures the slip
length induced by the continuity of stresses between the two regions (Neale
and Nader, 1974). If the porous matrix has sufficiently low permeability,
fluid velocities in the porous domain are near zero and a no-slip condition is
recovered at the interface between solid-free and porous regions (Angot et
al., 1999; Khadra et al., 2000; Soulaine and Tchelepi, 2016). This enables the
use of micro-continuum simulations at the pore-scale using a penalized
approach, i.e., the solid phase is described as a low-permeability porous
medium.

For two-phase flow, the discontinuity in porosity leads to a change in the
form of the surface tension force. Here, we treat this discontinuity by
assuming that the fluid-fluid interface of a droplet on a porous substrate
forms a contact angle $\theta$ with the solid surface (see Fig. 3). The
contact angle is an upscaled parameter that depends on various sub-grid scale
properties including interfacial energies, surface roughness, and the presence
of thin precursor films (Wenzel, 1936; Cassie and Baxter, 1944; Meakin and
Tartakovsky, 2009). In the present model, the contact angle is imposed by
locally modifying the orientation of the fluid-fluid interface relative to the
solid surface (Horgue et al., 2014; Soulaine et al., 2018, 2019). This is
achieved by replacing the mean normal vector $\overline{\boldsymbol{n}_{lg}}$
at the interface between the clear fluid and the porous regions by a locally
modified normal, $\hat{\boldsymbol{n}}_{lg}$, that satisfies the condition,

$\hat{\boldsymbol{n}}_{lg}=\cos\theta\boldsymbol{n}_{wall}+\sin\theta\boldsymbol{t}_{wall},$
(40)

where $\boldsymbol{n}_{wall}$ and $\boldsymbol{t}_{wall}$ are the normal and
tangent vectors to the porous surface, respectively. The numerical strategy to
implement Eq. (40) is described in details in Horgue et al. (2014) and
Soulaine et al. (2018). The effectiveness of this interfacial condition is
demonstrated in Section 4.2.

##### Summary of the multiscale parameters

The multiscale parameters $\bar{\boldsymbol{v}}_{r}$, $\mu k^{-1}$, and
$\boldsymbol{F}_{c}$ were derived by asymptotic matching to the VOF method in
solid-free regions and to the multiphase Darcy model in porous regions. The
resulting parameters, therefore, have different forms in different regions. In
the porous domains, the multiscale parameters depend on concept such as
relative permeability $k_{r,i}$ (also described in terms of fluid mobility,
$M_{i}=k_{r,i}/\mu_{i}$) and capillary pressure
$p_{c}\left(\alpha_{l}\right)$.

The relative velocity follows the relation:

$\bar{\boldsymbol{v}}_{r}=\begin{cases}C_{\alpha}\max\left(\left|\overline{\boldsymbol{v}}\right|\right)\frac{\nabla\alpha_{l}}{\left|\nabla\alpha_{l}\right|},&\textnormal{in
clear fluid regions,}\\\
\phi^{-1}\left[-\left(\frac{M_{l}}{\alpha_{l}}-\frac{M_{g}}{\alpha_{g}}\right)\nabla\overline{p}+\left(\frac{\rho_{l}M_{l}}{\alpha_{l}}-\frac{\rho_{g}M_{g}}{\alpha_{g}}\right)\boldsymbol{g}+\left(\frac{M_{l}\alpha_{g}}{\alpha_{l}}+\frac{M_{g}\alpha_{l}}{\alpha_{g}}\right)\nabla
p_{c}-\left(\frac{M_{l}}{\alpha_{l}}-\frac{M_{g}}{\alpha_{g}}\right)p_{c}\nabla\alpha_{l}\right],&\textnormal{in
porous regions.}\end{cases}$ (41)

The single-field relative permeability is given by:

$\mu k^{-1}=\begin{cases}0,&\textnormal{in solid-free regions,}\\\
k_{0}^{-1}\left(\frac{k_{r,l}}{\mu_{l}}+\frac{k_{r,g}}{\mu_{g}}\right)^{-1},&\textnormal{in
porous regions.}\end{cases}$ (42)

The body force $\boldsymbol{F}_{c}$ describes the capillary forces within a
computational cell using:

$\boldsymbol{F}_{c}=\begin{cases}-\phi^{-1}\sigma\nabla.\left(\hat{\boldsymbol{n}}_{lg}\right)\nabla\alpha_{l},&\textnormal{in
solid-free regions,}\\\
\left[M^{-1}\left(M_{l}\alpha_{g}-M_{g}\alpha_{l}\right)\left(\frac{\partial
p_{c}}{\partial S}\right)-p_{c}\right]\nabla\alpha_{l},&\textnormal{in porous
regions,}\end{cases}$ (43)

where the modified normal at the fluid-fluid interface is:

$\hat{\boldsymbol{n}}_{lg}=\begin{cases}-\frac{\nabla\alpha_{l}}{\left|\nabla\alpha_{l}\right|},&\textnormal{in
solid-free regions,}\\\
\cos\theta\boldsymbol{n}_{wall}+\sin\theta\boldsymbol{t}_{wall},&\textnormal{at
the interface between solid-free and porous regions.}\end{cases}$ (44)

Finally, the single-field fluid density is expressed as:

$\rho=\begin{cases}\rho_{l}\alpha_{l}+\rho_{g}\alpha_{g},&\textnormal{in
solid-free regions,}\\\
\left(\rho_{g}M_{g}+\rho_{l}M_{l}\right)M^{-1},&\textnormal{in porous
regions.}\end{cases}$ (45)

Figure 3: Conceptual Representation of the multiphase DB micro-continuum
approach. Here $\theta$ represents the contact angle and REV is a
Representative Elementary Volume.

##  3 Numerical implementation

The two-phase multi-scale micro-continuum model proposed in this paper is
implemented in the open-source CFD platform OpenFOAM® version 7.0 from
https://www.openfoam.org. This code is a C++ library that solves partial
differential equations with the finite-volume method. It handles complex
structured and unstructured three dimensional grids by default and has
demonstrated a good scalability for parallel computing of flow in porous media
(Orgogozo et al., 2014; Horgue et al., 2015; Guibert et al., 2015). One of its
features is that it solves the coupled equations using sequential approaches.
The present section details the solution algorithm developed in this paper.
Particularly close attention is paid to the description of the velocity-
pressure coupling.

###  3.1 Discretization of the equations

The momentum equation, Eq. (21), is transformed into a set of algebraic
equations after application of the finite-volume discretization procedure. The
nonlinearity introduced by the advection term is dealt with by linearizing
around the latest velocity field. The momentum equation is expressed in semi-
discrete form (with successive time levels denoted by $k$ and $k+1$) using a
Euler implicit difference scheme:

$\mathcal{V}\left(\frac{\rho^{k+1}\mathbf{v}_{P}^{k+1}-\rho^{k}\mathbf{v}_{P}^{k}}{\delta
t}\right)=-a_{P}^{{}^{\prime}}\textbf{v}_{P}^{k+1}+\sum_{NP}a_{NP}^{{}^{\prime}}\textbf{v}_{NP}^{k+1}-\nabla
p+\rho\boldsymbol{g}+\boldsymbol{F}_{c}-K_{fs}\textbf{v}_{P}^{k+1}.$ (46)

In equation (46), $\mathcal{V}$ and $\delta t$ stand for the cell volume and
time step, respectively. The subscript $P$ indicates values at the cell
center. The coefficients $a_{NP}^{{}^{\prime}}$ account for the influence of
neighboring control volumes and primarily include convective and diffusive
fluxes across cell faces. $K_{fs}$ corresponds to the exchange of momentum of
the fluids with regard to the solid structure, i.e., the Darcy term in Eq.
(21). The pressure gradient, buoyancy term, and capillary force are not
discretized at this stage.

All explicit source terms other than the pressure gradient, buoyancy term, and
capillary force are combined into a single vector,
$\textbf{S}=\frac{\mathcal{V}\rho^{k}\mathbf{v}_{P}^{k}}{\delta t}$. Eq. (46)
can then be rearranged as:

$\left(\frac{\mathcal{V}\rho^{k+1}}{\delta
t}+a_{P}^{{}^{\prime}}+K_{fs}\right)\textbf{v}_{P}^{k+1}=\sum_{NP}a_{NP}^{{}^{\prime}}\textbf{v}_{NP}^{k+1}+\mathbf{S}-\nabla
p+\rho\boldsymbol{g}+\boldsymbol{F}_{c}.$ (47)

This equation forms a matrix system that results from the momentum equation
discretization. The term $a_{P}=\left(\frac{\mathcal{V}\rho^{k+1}}{\delta
t}+a_{P}^{{}^{\prime}}+K_{fs}\right)$ represents the diagonal term of this
matrix. Following OpenFOAM® internal notations (Jasak, 1996), the operator
$\textbf{H}\left(\textbf{X}\right)=\sum_{NP}a_{NP}^{{}^{\prime}}\textbf{X}_{NP}+\textbf{S}$
is introduced and Eq. (47) becomes:

$a_{P}\textbf{v}_{P}^{k+1}=\textbf{H}\left(\textbf{v}^{k+1}\right)-\nabla
p+\rho\boldsymbol{g}+\boldsymbol{F}_{c}.$ (48)

This semi-discretized form of the momentum balance is used to form the
pressure equation. This is usually achieved by dividing Eq. (48) by the
diagonal coefficient, $a_{P}$, and substituting the semi-discretized form of
$\textbf{v}^{k+1}$ into the overall mass balance, Eq. (13), which is a
divergence free velocity in the absence of phase change. Finally, the pressure
equation can be written as:

$\nabla.\left(\frac{\textbf{H}\left(\textbf{v}^{k+1}\right)+\rho\boldsymbol{g}+\boldsymbol{F}_{c}}{a_{P}}\right)-\nabla.\left(\frac{1}{a_{P}}\nabla
p^{k+1}\right)=0.$ (49)

Further details regarding the discretization procedure in OpenFOAM® can be
found in Jasak (1996) and Weller et al. (1998). The saturation equation, Eq.
(15), is discretized with a Van Leer limiter function for the convection term
and a forward Euler scheme for time discretization.

###  3.2 Solution algorithm

The discretized equations are solved using OpenFOAM® in a segregated way. In
particular, the pressure-velocity coupling formed by Eqs. (48) and (49) is
handled by a predictor-corrector algorithm along the same lines as the
Pressure Implicit Splitting-Operator (PISO) algorithm originally designed by
Issa (1985) to solve the transient Navier-Stokes equations. It is built on the
top of the OpenFOAM® VOF solver _interFoam_. The numerical scheme uses the
following sequence of steps. First, the saturation equation, Eq. (15), is
solved explicitly using the OpenFOAM® implementation of the Flux Corrected
Transport (FCT) theory (Rudman, 1997) called Multidimensional Universal
Limiter with Explicit Solution (MULES). Details regarding the MULES algorithm
can be found in the Chapter 5 of Damian (2013). Second, the boundary values of
$\boldsymbol{v}$ and $\boldsymbol{v_{r}}$ are updated according to Eqs. (10)
and (14). Third, the single-field relative permeability $k^{k+1}$, density
$\rho^{k+1}$, and viscosity $\mu^{k+1}$ are updated using the new value of the
saturation field, $\alpha_{l}^{k+1}$. The surface tension force,
$\boldsymbol{F}_{c}^{k+1}$, is computed using Eq. (43). Fourth, the velocity
field $\mathbf{v}^{*}$ is calculated by solving implicitly the momentum
equation,

$a_{P}\mathbf{v}_{P}^{*}=\textbf{H}\left(\mathbf{v}^{*}\right)+\rho^{k+1}\boldsymbol{g}+\boldsymbol{F}_{c}^{k+1}-\nabla
p^{k},$ (50)

where the gradient of the pressure field is evaluated from the values computed
at the previous time step. This stage is called the momentum predictor. Fifth,
the predicted velocity $\textbf{v}^{*}$ (which does not satisfy the continuity
equation, Eq. (13)) is corrected. This is achieved by finding
$\left(\mathbf{v}^{**},p^{*}\right)$ that obeys,

$\mathbf{v}_{P}^{**}=\frac{1}{a_{P}}\left[\textbf{H}\left(\mathbf{v}^{*}\right)+\rho^{k+1}\boldsymbol{g}+\boldsymbol{F}_{c}^{k+1}-\nabla
p^{*}\right],$ (51) $\nabla.\mathbf{v}^{**}=0.$ (52)

Based on these two equations, the pressure equation is formulated as

$\nabla.\left(\frac{\textbf{H}\left(\textbf{v}^{*}\right)+\rho^{k+1}\boldsymbol{g}+\boldsymbol{F}_{c}^{k+1}}{a_{P}}\right)-\nabla.\left(\frac{1}{a_{P}}\nabla
p^{*}\right)=0,$ (53)

and solved implicitly with a generalized method of Geometric-Algebraic Multi-
Grid (GAMG) embedded in OpenFOAM®. The corrected velocity $\textbf{v}^{**}$ is
then computed point-wise from Eq. (51). This step (the PISO loop) may be
repeated several times to ensure convergence. Issa (1985) demonstrated that at
least two iterations are required to ensure that the solution of the pressure-
velocity $\left(\textbf{v},p\right)$ coupling satisfies mass conservation. The
resulting values are set to $(\textbf{v}^{k+1},p^{k+1})$ and, then, the
algorithm marches in time as dictated by the imposed Courant-Friedrichs-Lewy
(CFL) number.

###  3.3 Open-Source Toolbox: hybridPorousInterFoam

The accompanying open-source toolbox follows the implementation described
above and consists of four distinct parts: a main directory that includes the
licence files, instructional files, release notes, and automated compilation
procedures along with three main toolbox sub-directories. The three sub-
directories consist of a Solver sub-directory that includes the code for the
hybridPorousInterFoam solver; a Tutorials sub-directory that includes all the
verification and example cases presented in this paper and; and a Libraries
sub-directory that includes both the dynamically linked libraries used in the
implementation of the penalized contact angle and, also, the Brooks and Corey
(1964) and Van Genutchen (1980) porous media models used to calculate the
required sub-voxel description of the fluid-fluid interface in terms of
relative permeability and capillary pressures (see Appendix A). This last
library was obtained from the open-sourced toolbox published in Horgue et al.
(2015). The hybridPorousInterFoam toolbox can be accessed from the author’s
repository (https://github.com/Franjcf).

##  4 Verification

In this section, the two-phase micro-continuum model is used in various
situations for which reference solutions exist. The objective is to verify
that the multiscale solver converges effectively towards its two asymptotic
limits, namely the two-phase Darcy model at the continuum scale and the VOF
formulation at the pore scale.

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

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

###  4.1 Darcy scale validation

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

####  4.1.1 Buckley-Leverett

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

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

####  4.1.2 Gravity dominated Buckley-Leverett

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

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

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

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

####  4.1.3 Gravity-capillarity equilibrium

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

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

which can be rearranged to yield:

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

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

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

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

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

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

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

###  4.2 Pore scale validation

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

####  4.2.1 Contact angle on a flat plate

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

####  4.2.2 Capillary rise

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

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

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

####  4.2.3 Taylor film

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

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

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

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

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

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

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

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

##  5 Hybrid Scale Applications

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

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

###  5.1 Wave propagation in a coastal barrier

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

Figure 11: Coastal barrier simulation at different times. The thin white line
represents the boundary of the porous solid: $\phi$ and $k_{0}^{-1}$ are set
to 0.5 and $5\text{\times}{10}^{7}\text{\,}\mathrm{m}^{-2}$ below said line
and to 1 and 0 above it, respectively.

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

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

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

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

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

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

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

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

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

##  6 Conclusion

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

The core derivation of our micro-continuum framework relies solely on
fundamental principles and uses the the methods of volume averaging and
asymptotic matching to modify and expand the Navier-Stokes equations. Through
this study, we showed that our model can successfully simulate multiphase
Darcy and Navier-Stokes flow to the same standard as conventional single-scale
solvers. The coupling between the two scales at porous interfaces is handled
by ensuring mass conservation and continuity of stresses at said boundary, as
well as by implementing a constant contact angle wettability condition. We
then leveraged all these features to show that our model can be used to model
hybrid scale systems such as wave interaction with a porous coastal barrier
and drainage and imbibition in a fractured porous matrix.

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

## Acknowledgments

This work was supported by the National Science Foundation, Division of Earth
Sciences, Early Career program through Award EAR-1752982. C.S was partially
sponsored by BRGM through the TRIPHASIQUE project and by the French National
Agency for Research that funded the FraMatI project under contract
ANR-19-CE05-0002.

## References

## References

  * Abu-Al-Saud et al. (2020) Abu-Al-Saud, M., Gmira, A., Al-Enezi, S., Yousef, A., et al., 2020. Pore-scale simulation of fluid flow in carbonates using micro-ct scan images. In: International Petroleum Technology Conference. International Petroleum Technology Conference.
  * Alvarado and Manrique (2010) Alvarado, V., Manrique, E., 8 2010. Enhanced oil recovery: An update review. Energies 3 (9), 1529–1575.
URL http://www.mdpi.com/1996-1073/3/9/1529

  * Angot et al. (1999) Angot, P., Bruneau, C.-H., Fabrie, P., 1999. A penalization method to take into account obstacles in incompressible viscous flows. Numerische Mathematik 81 (4), 497–520.
  * Apourvari and Arns (2014) Apourvari, S. N., Arns, C. H., 2014. An assessment of the influence of micro-porosity for effective permeability using local flux analysis on tomographic images. International Petroleum Technology Conference, 19-22 January, Doha, Qatar.
  * Auriault (2009) Auriault, J.-L., 2009. On the domain of validity of brinkman’s equation. Transport in porous media 79 (2), 215–223.
  * Aussillous et al. (2000) Aussillous, P., Quéré, D., Qué, D., 2000. Quick deposition of a fluid on the wall of a tube. Physics of Fluids 12 (10), 2367–32107.
URL
https://doi.org/10.1063/1.1289396http://aip.scitation.org/toc/phf/12/10https://doi.org/10.1063/1.1289396%0Ahttp://aip.scitation.org/toc/phf/12/10

  * Baber et al. (2016) Baber, K., Flemisch, B., Helmig, R., 2016. Modeling drop dynamics at the interface between free and porous-medium flow using the mortar method. International Journal of Heat and Mass Transfer 99, 660–671.
  * Battiato et al. (2019) Battiato, I., O’Malley, D., Miller, C. T., Takhar, P. S., Valdés-Parada, F. J., Wood, B. D., et al., 2019. Theory and applications of macroscale models in porous media. Transport in Porous Media, 1–72.
  * Bear (1972) Bear, J., 1972. Dynamics of fluids in porous media. Elsevier, New York.
URL Dynamicsoffluidsinporousmedia

  * Beavers and Joseph (1967) Beavers, G. S., Joseph, D. D., 1967. Boundary conditions at a naturally permeable wall. Journal of Fluid Mechanics 30, 197–207.
  * Bousquet-Melou et al. (2002) Bousquet-Melou, P., Goyeau, B., Quintard, M., Fichot, F., Gobin, D., 2002. Average momentum equation for interdendritic flow in a solidifying columnar mushy zone. International Journal of Heat and Mass Transfer 45 (17), 3651 – 3665\.
URL http://www.sciencedirect.com/science/article/pii/S0017931002000777

  * Brackbill et al. (1992) Brackbill, J., Kothe, D., Zemach, C., 1992. A continuum method for modeling surface tension. Journal of Computational Physics 100 (2), 335 – 354.
URL http://www.sciencedirect.com/science/article/pii/002199919290240Y

  * Brinkman (1947) Brinkman, H. C., 1947. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl. Sci. Res. A1, 27–34.
  * Brooks and Corey (1964) Brooks, R., Corey, A., 1964. Hydraulic properties of porous media. Colorado State University, Hydro Paper 3, 27.
  * Carrillo and Bourg (2019) Carrillo, F. J., Bourg, I. C., 2019. A darcy-brinkman-biot approach to modeling the hydrology and mechanics of porous media containing macropores and deformable microporous regions. Water Resources Research 55, 8096–8121.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019WR024712

  * Cassie and Baxter (1944) Cassie, A., Baxter, S., 1944. Wettability of porous surfaces. Transactions of the Faraday Society 40, 546–551.
  * Chen and Meiburg (1998) Chen, C. Y., Meiburg, E., 1998. Miscible porous media displacements in the quarter five-spot configuration. Part 1. the homogeneous case. Journal of Fluid Mechanics 371, 233–268.
URL https://doi.org/10.1017/S0022112098002195

  * Chen and Wilkinson (1985) Chen, J. D., Wilkinson, D., 1985. Pore-scale viscous fingering in porous media. Physical Review Letters 55 (18), 1892–1895.
  * Cox (1986) Cox, R., 1986. The dynamics of the spreading of liquids on a solid surface. part 1. viscous flow. Journal of Fluid Mechanics 168, 169–194.
  * Damian (2013) Damian, S. M., 2013. An extended mixture model for the simultaneous treatment of short and long scale interfaces. Ph.D. thesis, FACULTAD DE INGENIERÍA Y CIENCIAS HÍDRICAS UNIVERSIDAD NACIONAL DEL LITORAL.
  * Datta et al. (2014) Datta, S. S., Dupin, J.-B., Weitz, D. A., 2014. Fluid breakup during simultaneous two-phase flow through a three-dimensional porous medium. Physics of Fluids 26, 062004.
  * Di Donato et al. (2003) Di Donato, G., Huang, W., Blunt, M., 4 2003. Streamline-Based Dual Porosity Simulation of Fractured Reservoirs. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
URL http://www.onepetro.org/doi/10.2118/84036-MS

  * Douglas et al. (1991) Douglas, J., Hensley, J. L., Arbogast, T., 1991. A dual-porosity model for waterflooding in naturally fractured reservoirs. Computer Methods in Applied Mechanics and Engineering 87 (2-3), 157–174.
  * Efendiev and Hou (2007) Efendiev, Y., Hou, T., 2007. Multiscale finite element methods for porous media flows and their applications. Applied Numerical Mathematics 57 (5-7), 577–596.
  * Ferrari and Lunati (2013) Ferrari, A., Lunati, I., 7 2013. Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy. Advances in Water Resources 57, 19–31.
  * Golfier et al. (2002) Golfier, F., Zarcone, C., Bazin, B., Lenormand, R., Lasseux, D., Quintard, M., 2002\. On the ability of a darcy-scale model to capture wormhole formation during the dissolution of a porous medium. Journal of fluid Mechanics 457, 213–254.
  * Goyeau et al. (2003) Goyeau, B., Lhuillier, D., Gobin, D., Velarde, M., 2003. Momentum transport at a fluid–porous interface. International Journal of Heat and Mass Transfer 46 (21), 4071–4081.
  * Graveleau et al. (2017) Graveleau, M., Soulaine, C., Tchelepi, H. A., 2017. Pore-scale simulation of interphase multicomponent mass transfer for subsurface flow. Transport in Porous Media 120 (2), 287–308.
  * Gray et al. (2015) Gray, W. G., Dye, A. L., McClure, J. E., Pyrak-Nolte, L. J., Miller, C. T., 2015\. On the dynamics and kinematics of two-fluid-phase flow in porous media. Water Resources Research 51, 5365–5381.
  * Gründing et al. (2019) Gründing, D., Smuda, M., Antritter, T., Fricke, M., Rettenmaier, D., Kummer, F., Stephan, P., Marschall, H., Bothe, D., 2019. Capillary Rise – A Computational Benchmark for Wetting Processes. arXiv, 1–38.
URL http://arxiv.org/abs/1907.05054

  * Guibert et al. (2015) Guibert, R., Nazarova, M., Horgue, P., Hamon, G., Creux, P., Debenest, G., 2015\. Computational permeability determination from pore-scale imaging: Sample size, mesh and method sensitivities. Transport in Porous Media 107 (3), 641–656.
  * Guo et al. (2018) Guo, B., Ma, L., Tchelepi, H. A., 2018. Image-based micro-continuum model for gas flow in organic-rich shale rock. Advances in Water Resources 122, 70–84.
  * Hassanizadeh et al. (2002) Hassanizadeh, S. M., Celia, M. A., Dahle, H. K., 2002. Dynamic effects in the capillary pressure-saturation relationship and its impacts on unsaturated flow. Vadose Zone Journal 1, 38–57.
  * Helmig et al. (2013) Helmig, R., Flemisch, B., Wolff, M., Ebigbo, A., Class, H., 2013. Model coupling for multiphase flow in porous media. Advances in Water Resources 51, 52–66.
  * Helmig et al. (2010) Helmig, R., Niessner, J., Flemisch, B., Wolff, M., Fritz, J., 2010. Handbook of Geomathematics. Springer-Verlag Berlin Heidelberg, Ch. Efficient Modeling of Flow and Transport in Porous Media Using Multiphysics and Multiscale Approaches, pp. 418 – 457.
  * Higuera et al. (2013) Higuera, P., Lara, J. L., Losada, I. J., 1 2013. Realistic wave generation and active wave absorption for Navier-Stokes models. Application to OpenFOAM®. Coastal Engineering 71, 102–118.
  * Hirt and Nichols (1981) Hirt, C., Nichols, B., 1981. Volume of fluid (vof) method for the dynamics of free boundaries. Journal of Computational Physics 39 (1), 201 – 225.
URL http://www.sciencedirect.com/science/article/pii/0021999181901455

  * Holzbecher (2009) Holzbecher, E., 2009. Modeling of Viscous Fingering. The COMSOL Conference 2009 Milan Modeling, 6.
  * Horgue et al. (2014) Horgue, P., Prat, M., Quintard, M., 2014. A penalization technique applied to the volume-of-fluid method: Wettability condition on immersed boundaries. Computers & Fluids 100, 255 – 266.
URL http://www.sciencedirect.com/science/article/pii/S0045793014002254

  * Horgue et al. (2015) Horgue, P., Soulaine, C., Franc, J., Guibert, R., Debenest, G., 2015. An open-source toolbox for multiphase flow in porous media. Computer Physics Communications 187 (0), 217– 226.
URL http://www.sciencedirect.com/science/article/pii/S0010465514003403

  * Howes and Whitaker (1985) Howes, F. A., Whitaker, S., 1985. The spatial averaging theorem revisited. Chemical Engineering Science 40 (8), 1387 – 1392.
URL http://www.sciencedirect.com/science/article/pii/0009250985800786

  * Hsu and Cheng (1990) Hsu, C., Cheng, P., 1990. Thermal dispersion in a porous medium. International Journal of Heat and Mass Transfer 33 (8), 1587–1597.
  * Ishii and Hibiki (2011) Ishii, M., Hibiki, T., 2011. Thermo-Fluid Dynamic of Two-Phase Flow, 2nd Edition. Smart Energy Systems: Nanowatts to Terawatts. Springer.
  * Issa (1985) Issa, R. I., 1985. Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics 62, 40–65.
  * Jabbari et al. (2016) Jabbari, M., Jambhekar, V. A., Hattel, J. H., Helmig, R., 2016. Drying of a tape-cast layer: Numerical modelling of the evaporation process in a graded/layered material. International Journal of Heat and Mass Transfer 103, 1144–1154.
  * Jasak (1996) Jasak, H., 1996. Error analysis and estimation for the finite volume method with applications to fluid flows. Ph.D. thesis, Department of Mechanical Engineering Imperial College of Science, Technology and Medicine.
  * Jenny et al. (2003) Jenny, P., Lee, S., Tchelepi, H. A., 2003. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. Journal of Computational Physics 187 (1), 47–67.
  * Jenny et al. (2006) Jenny, P., Lee, S. H., Tchelepi, H. A., 2006. Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media. Journal of Computational Physics 217 (2), 627–641.
  * Jurin (1719) Jurin, J., 1 1719. II. An account of some experiments shown before the Royal Society; with an enquiry into the cause of the ascent and suspension of water in capillary tubes. Philosophical Transactions of the Royal Society of London 30 (355), 739–747.
  * Kang et al. (2019) Kang, D. H., Yang, E., Yun, T. S., 2019. Stokes-brinkman flow simulation based on 3-d $\mu$-ct images of porous rock using grayscale pore voxel permeability. Water Resources Research 55 (5), 4448–4464.
  * Karimi-Fard and Durlofsky (2016) Karimi-Fard, M., Durlofsky, L. J., 2016. A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features. Advances in water resources 96, 354–372.
  * Khadra et al. (2000) Khadra, K., Angot, P., Parneix, S., Caltagirone, J.-P., 2000. Fictitious domain approach for numerical modelling of navier-stokes equations. International Journal for Numerical Methods in Fluids 34 (8), 651–684.
URL http://dx.doi.org/10.1002/1097-0363(20001230)34:8<651::AID-
FLD61>3.0.CO;2-D

  * Knackstedt et al. (2006) Knackstedt, M., Sheppard, A., Nguyen, V., Arns, C., Sok, R., et al., 2006. 3d imaging and characterization of the pore space of carbonate core; implications to single and two phase flow properties. In: SPWLA 47th Annual Logging Symposium. Society of Petrophysicists and Well-Log Analysts.
  * Lacis et al. (2017) Lacis, U., Zampogna, G. A., Bagheri, S., 3 2017. A computational continuum model of poroelastic beds. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473 (2199).
  * Lake et al. (2014) Lake, L., Johns, R. T., Rossen, W. R., Pope, G. A., 2014. Fundamentals of enhanced oil recovery. Society of Petroleum Engineers, Richardson, TX.
  * Larsen et al. (2019) Larsen, B. E., Fuhrman, D. R., Roenby, J., 7 2019. Performance of interFoam on the simulation of progressive waves. Coastal Engineering Journal 61 (3), 380–400.
URL https://www.tandfonline.com/doi/full/10.1080/21664250.2019.1609713

  * Lasseux et al. (1996) Lasseux, D., Quintard, M., Whitaker, S., 1996. Determination of permeability tensors for two-phase flow in homogeneous porous media: Theory. Transport in Porous Media 24 (2), 107–137.
  * Lenormand et al. (1988) Lenormand, R., Touboul, E., Zarcone, C., 1988. Numerical models and experiments on immiscible displacements in porous media. Journal of Fluid Mechanics 189, 165–187.
  * Leverett (1940) Leverett, M. C., december 1940. Capillary behavior in porous solids. Transactions of the AIME 142 (1), 152–169.
  * Li and Benson (2015) Li, B., Benson, S. M., 2015. Influence of small-scale heterogeneity on upward co2 plume migration in storage aquifers. Advances in Water Resources 83, 389–404.
  * Li et al. (2018) Li, T., Schlüter, S., Dragila, M. I., Wildenschild, D., 2018. An improved method for estimating capillary pressure from 3d microtomography images and its application to the study of disconnected nonwetting phase. Advances in Water Resources 114, 249–260.
  * Liu et al. (2019) Liu, J., Ju, Y., Zhang, Y., Gong, W., 12 2019. Preferential Paths of Air-water Two-phase Flow in Porous Structures with Special Consideration of Channel Thickness Effects. Nature Scientific Reports 9 (1), 16204.
URL http://www.nature.com/articles/s41598-019-52569-9

  * Liu et al. (1997) Liu, X., Ormond, A., Bartko, K., Ying, L., Ortoleva, P., 1997. A geochemical reaction-transport simulator for matrix acidizing analysis and design. Journal of Petroleum Science and Engineering 17 (1), 181–196.
  * Maes and Soulaine (2019) Maes, J., Soulaine, C., 2019. A unified single-field volume-of-fluid-based formulation for multi-component interfacial transfer with local volume changes. Journal of Computational Physics.
  * Maxwell et al. (2014) Maxwell, R. M., Putti, M., Meyerhoff, S., Delfs, J.-O., Ferguson, I. M., Ivanov, V., Kim, J., Kolditz, O., Kollet, S. J., Kumar, M., Lopez, S., Niu, J., Paniconi, C., Park, Y.-J., Phanikumar, M. S., Shen, C., Sudicky, E. A., Sulis, M., 2014. Surface-subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feebacks. Water Resources Research 50, 1531–1549.
  * Meakin and Tartakovsky (2009) Meakin, P., Tartakovsky, A. M., 2009. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Reviews of Geophysics 47, RG3002.
  * Miller et al. (2019) Miller, C., Bruning, K., Talbot, C., McClure, J., Gray, W., 2019. Nonhysteretic capillary pressure in two-fluid porous medium systems: Definition, evaluation, validation, and dynamics. Water Resources Research 55, 6825–6849.
  * Miller et al. (1998) Miller, C. T., Christakos, G., Imhoff, P. T., McBride, J. F., Pedit, J. A., Trangenstein, J. A., 1998. Multiphase flow and transport modeling in heterogeneous porous media: challenges and approaches. Advances in Water Resources 21 (2), 77 – 120.
URL https://www.sciencedirect.com/science/article/pii/S030917089600036X

  * Molins et al. (2020) Molins, S., Soulaine, C., Prasianakis, N., Abbasi, A., Poncet, P., Ladd, A., Starchenko, V., Roman, S., Trebotich, D., Tchelepi, H., Steefel, C., 2020. Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: Review of approaches and benchmark problem set. Computational Geosciences, 1–34.
  * Morton (2002) Morton, R. A., 2002. Factors controlling storm impacts on coastal barriers and beaches - A preliminary basis for near real-time forecasting. Journal of Coastal Research 18 (3), 486–501.
URL https://www.jstor.org/stable/4299096

  * Muskat (1949) Muskat, M., 1949. Physical principles of oil production. McGraw-Hill, New York.
URL Physicalprinciplesofoilproduction

  * Nandlal and Weijermars (2019) Nandlal, K., Weijermars, R., 6 2019. Drained rock volume around hydraulic fractures in porous media: planar fractures versus fractal networks. Petroleum Science.
  * Neale and Nader (1974) Neale, G., Nader, W., 1974. Practical significance of brinkman’s extension of darcy’s law: coupled parallel flows within a channel and a bounding porous medium. The Canadian Journal of Chemical Engineering 52 (4), 475–478.
  * Or et al. (2013) Or, D., Lehmann, P., Shahraeeni, E., Shokri, N., 2013. Advances in soil evaporation physics–a review. Vadose Zone Journal 12 (4).
  * Orgogozo et al. (2014) Orgogozo, L., Renon, N., Soulaine, C., Hénon, F., Tomer, S., Labat, D., Pokrovsky, O., Sekhar, M., Ababou, R., Quintard, M., 2014. An open source massively parallel solver for richards equation: Mechanistic modelling of water fluxes at the watershed scale. Computer Physics Communications 185 (12), 3358 – 3371.
URL http://www.sciencedirect.com/science/article/pii/S0010465514002719

  * Patankar (1980) Patankar, S. V., 1980. Numerical Heat Transfer And Fluid Flow. Taylor & Francis, Washington, DC.
  * Picchi and Battiato (2018) Picchi, D., Battiato, I., 2018. The impact of pore-scale flow regimes on upscaling of immiscible two-phase flow in porous media. Water Resources Research 54, 6683–6707.
  * Pinder and Gray (2008) Pinder, G. F., Gray, W. G., 2008. Essentials of Multiphase Flow and Transport in Porous Media. John Wiley & Sons.
  * Quintard and Whitaker (1994) Quintard, M., Whitaker, S., 1994. Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media. Advances in Water Resources 17 (4), 221 – 239.
  * Riaz and Tchelepi (2006) Riaz, A., Tchelepi, H. A., 2006. Numerical simulation of immiscible two-phase flow in porous media. Physics of Fluids 18 (1).
  * Rocco et al. (2017) Rocco, S., Woods, A. W., Harrington, J., Norris, S., 2017. An experimental model of episodic gas release through fracture of fluid confined within a pressurized elastic reservoir. Geophysical Research Letters 44, 751–759.
  * Roman et al. (2016) Roman, S., Soulaine, C., AlSaud, M. A., Kovscek, A., Tchelepi, H., 2016. Particle velocimetry analysis of immiscible two-phase flow in micromodels. Advances in Water Resources 95, 199–211.
URL http://ac.els-
cdn.com/S0309170815002018/1-s2.0-S0309170815002018-main.pdf?_tid=028cf0f2-39bb-11e7-bfcf-00000aacb35e&acdnat=1494886282_d68eea9c1f173d770eb2063395074657

  * Rudman (1997) Rudman, M., 1997. Volume-tracking methods for interfacial flow calculations. International journal for numerical methods in fluids 24 (7), 671–691.
  * Rusche (2002) Rusche, H., 2002. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. Ph.D. thesis, Imperial College of Science, Technology & Medicine Department of Mechanical Engineering Exhibition Road, London SW7 2BX.
  * Saffman and Taylor (1988) Saffman, P., Taylor, G., 1988. The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Dynamics of Curved Fronts 245 (1242), 155–174.
URL
http://rspa.royalsocietypublishing.org/content/royprsa/245/1242/312.full.pdf

  * Scandella et al. (2017) Scandella, B. P., Delwiche, K., Hemond, H. F., Juanes, R., 2017. Persistence of bubble outlets in soft, methane-generating sediments. Journal of Geophysical Research: Biogeosciences 122, 1298–1320.
  * Scheibe et al. (2015) Scheibe, T. D., Perkins, W. A., Richmond, M. C., McKinley, M. I., Romero-Gomez, P. D. J., Oostrom, M., Wietsma, T. W., Serkowski, J. A., Zachara, J. M., 2015\. Pore-scale and multiscale numerical simulation of flow and transport in a laboratory-scale column. Water Resources Research 51 (2), 1023–1035.
URL http://dx.doi.org/10.1002/2014WR015959

  * Singh (2019) Singh, K., 2019. How hydraulic properties of organic matter control effective liquid permeability of mudrocks. Transport in Porous Media, 1–17.
  * Soulaine et al. (2019) Soulaine, C., Creux, P., Tchelepi, H. A., 2019. Micro-continuum framework for pore-scale multiphase fluid transport in shale formations. Transport in Porous Media.
  * Soulaine et al. (2016) Soulaine, C., Gjetvaj, F., Garing, C., Roman, S., Russian, A., Gouze, P., Tchelepi, H., May 2016. The impact of sub-resolution porosity of x-ray microtomography images on the permeability. Transport in Porous Media 113 (1), 227–243.
  * Soulaine et al. (2017) Soulaine, C., Roman, S., Kovscek, A., Tchelepi, H. A., 2017. Mineral dissolution and wormholing from a pore-scale perspective. Journal of Fluid Mechanics 827, 457–483.
URL
https://www.cambridge.org/core/product/identifier/S0022112017004992/type/journal_article

  * Soulaine et al. (2018) Soulaine, C., Roman, S., Kovscek, A., Tchelepi, H. A., 2018. Pore-scale modelling of multiphase reactive flow. Application to mineral dissolution with production of CO2. Journal of Fluid Mechanics 855, 616–645.
  * Soulaine and Tchelepi (2016) Soulaine, C., Tchelepi, H. A., 2016. Micro-continuum approach for pore-scale simulation of subsurface processes. Transport In Porous Media 113, 431–456.
  * Standnes et al. (2017) Standnes, D. C., Evje, S., Andersen, P. O., 2017. A novel relative permeability model based on mixture theory approach accounting for solid-fluid and fluid-fluid interactions. Transport in Porous Media 119, 707–738.
  * Starnoni and Pokrajac (2020) Starnoni, M., Pokrajac, D., 2020. On the concept of macroscopic capillary pressure in two-phase porous media flow. Advances in Water Resources 135, 103487\.
  * Sun et al. (2012) Sun, T., Mehmani, Y., Balhoff, M. T., 2012. Hybrid multiscale modeling through direct substitution of pore-scale models into near-well reservoir simulators. Energy & Fuels 26 (9), 5828–5836.
  * Tam (1969) Tam, C. K., 1969. The drag on a cloud of spherical particles in low reynolds number flow. Journal of Fluid Mechanics 38 (03), 537–546.
  * Tomin and Lunati (2013) Tomin, P., Lunati, I., 2013. Hybrid multiscale finite volume method for two-phase flow in porous media. Journal of Computational Physics 250 (0), 293 – 307.
URL http://www.sciencedirect.com/science/article/pii/S0021999113003513

  * Tomin and Lunati (2016) Tomin, P., Lunati, I., 2016. Investigating darcy-scale assumptions by means of a multiphysics algorithm. Advances in water resources 95, 80–91.
  * Tomin and Voskov (2018) Tomin, P., Voskov, D., 2018. Robust and accurate formulation for modeling of acid stimulation. In: ECMOR XVI-16th European Conference on the Mathematics of Oil Recovery.
  * Vafai and Tien (1981) Vafai, K., Tien, C., 1981. Boundary and inertia effects on flow and heat transfer in porous media. International Journal of Heat and Mass Transfer 24 (2), 195–203.
  * Van Genutchen (1980) Van Genutchen, M. T., 1980. A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal (02), 1–8.
  * Voinov (1976) Voinov, O., 1976. Hydrodynamics of wetting. Fluid dynamics 11 (5), 714–721.
  * Wang and Beckermann (1993) Wang, C.-Y., Beckermann, C., 1993. A two-phase mixture model of liquid-gas flow and heat transfer in capillary porous media-i. formulation. International journal of heat and mass transfer 36, 2747–2747.
  * Weller et al. (1998) Weller, H. G., Tabor, G., Jasak, H., Fureby, C., NOV/DEC 1998. A tensorial approach to computational continuum mechanics using object-oriented techniques. COMPUTERS IN PHYSICS 12 (6), 620–631.
  * Wenzel (1936) Wenzel, R. N., 8 1936. Resistance of solid surfaces to wetting by water. Industrial and Engineering Chemistry 28 (8), 988–994.
  * Whitaker (1986a) Whitaker, S., 1986a. Flow in porous media i: A theoretical derivation of darcy’s law. Transport in Porous Media 1, 3–25, 10.1007/BF01036523.
URL http://dx.doi.org/10.1007/BF01036523

  * Whitaker (1986b) Whitaker, S., 1986b. Flow in porous media ii: The governing equations for immiscible, two-phase flow. Transport in porous media 1 (2), 105–125.
  * Whitaker (1999) Whitaker, S., 1999. The method of volume averaging, theory and applications of transport in porous media. Dorderecht: Kluwer Academic.
  * Whyman et al. (2008) Whyman, G., Bormashenko, E., Stein, T., 1 2008. The rigorous derivation of Young, Cassie-Baxter and Wenzel equations and the analysis of the contact angle hysteresis phenomenon. Chemical Physics Letters 450 (4-6), 355–359.
  * Xu et al. (2017) Xu, Z., Liu, H., Valocchi, A. J., 5 2017. Lattice Boltzmann simulation of immiscible two-phase flow with capillary valve effect in porous media. Water Resources Research 53 (5), 3770–3790.
URL http://doi.wiley.com/10.1002/2017WR020373

  * Yan et al. (2016) Yan, B., Wang, Y., Killough, J. E., 2 2016. Beyond dual-porosity modeling for the simulation of complex flow mechanisms in shale reservoirs. Computational Geosciences 20 (1), 69–91.
  * Yin et al. (2018) Yin, Z., Khurana, M., Tan, H. K., Linga, P., 2018. A review of gas hydrate growth kinetic models. Chemical Engineering Journal 342, 9–29.
  * Zacharoudiou et al. (2018) Zacharoudiou, I., Boek, E. S., Crawshaw, J., 12 2018. The impact of drainage displacement patterns and Haines jumps on CO2 storage efficiency. Scientific Reports 8 (1).
  * Zampogna et al. (2019) Zampogna, G. A., Lācis, U., Bagheri, S., Bottaro, A., 1 2019. Modeling waves in fluids flowing over and through poroelastic media. International Journal of Multiphase Flow 110, 148–164.

##  Appendix A Relative Permeability and Capillary Pressure Models

###  A.1 Relative Permeability Models

The two implemented relative permeability models rely on a definition of the
effective saturation of the wetting fluid, $\alpha_{l,eff}$, as a function of
each fluid’s irreducible saturation $\alpha_{i,irr}$:

$\alpha_{l,eff}=\frac{\alpha_{l}-\alpha_{l,irr}}{1-\alpha_{g,irr}-\alpha_{l,irr}}$
(59)

The model proposed by Brooks and Corey (1964) relates the relative
permeability of each phase to the effective saturation through the following
relation, where $m$ is a non-dimensional coefficient dictated by porous media
properties:

$k_{rg}=(1-\alpha_{l,eff})^{m}$ (60) $k_{rl}=\alpha_{l,eff}^{m}$ (61)

Alternatively, the Van Genutchen (1980) model relates the relative
permeabilities to the wetting fluid’s effective saturation in the following
way:

$k_{r,g}=(1-\alpha_{l,eff})^{1/2}(1-\alpha_{l,eff}^{1/m})^{2m}$ (62)
$k_{r,l}=\alpha_{l,eff}^{1/2}(1-(1-\alpha_{l,eff}^{1/m})^{m})^{2}$ (63)

###  A.2 Capillary Pressure Models

The implemented capillary pressure models rely on slightly different
formulations for the effective wetting fluid saturation, $\alpha_{l,pc}$:

$\alpha_{l,pc}=\frac{\alpha_{l}-\alpha_{pc,irr}}{\alpha_{pc,max}-\alpha_{pc,irr}}$
(64)

where $\alpha_{pc,max}$ is the maximum saturation of the wetting fluid and
$\alpha_{pc,irr}$ is its irreducible saturation. The Brooks and Corey (1964)
model uses the following expression to calculate capillary pressure within a
porous medium:

$p_{c}=p_{c,0}\alpha_{l,pc}^{-\beta}$ (65)

where $p_{c,0}$ is the entry capillary pressure and $1/\beta$ is a parameter
that can be calculated from the pore size distribution. Alternatively, the Van
Genutchen (1980) model uses the relation:

$p_{c}=p_{c,0}(\alpha_{l,pc}^{-1/m}-1)^{1-m}$ (66)

##  Appendix B Alternative formulation of multi-scale parameters

For consistency with Eq. 20, we can recast Eq. 36 as follows

$0=-\nabla\overline{p}+\rho\boldsymbol{g}-M^{-1}\overline{\boldsymbol{v}}+(\rho^{*}-\rho)\boldsymbol{g}+M^{-1}\left[M_{l}\nabla\left(\alpha_{g}p_{c}\right)-M_{g}\nabla\left(\alpha_{l}p_{c}\right)\right],$
(67)

The capillary force $\boldsymbol{F}_{c}$, then, includes a term in
$(\rho^{*}-\rho)\boldsymbol{g}$. After some manipulation, we obtain the
following expression (instead of Eq. 38)

$\displaystyle\boldsymbol{F}_{c}$
$\displaystyle=(\rho^{*}-\rho)\boldsymbol{g}+M^{-1}\left[M_{l}\nabla\left(\alpha_{g}p_{c}\right)-M_{g}\nabla\left(\alpha_{l}p_{c}\right)\right]$
$\displaystyle=M^{-1}\left(M_{l}\alpha_{g}-M_{g}\alpha_{l}\right)(\rho_{l}-\rho_{g})\boldsymbol{g}+M^{-1}\left(M_{l}\alpha_{g}-M_{g}\alpha_{l}\right)\nabla
p_{c}-p_{c}\nabla\alpha_{l}$
$\displaystyle=M^{-1}\left(M_{l}\alpha_{g}-M_{g}\alpha_{l}\right)[(\rho_{l}-\rho_{g})\boldsymbol{g}+\nabla
p_{c}]-p_{c}\nabla\alpha_{l}$ (68)

At equilibrium
$(\overline{\boldsymbol{v}_{l}}=\overline{\boldsymbol{v}_{g}}=0)$, the
multiphase Darcy’s law yields $(\rho_{l}-\rho_{g})\boldsymbol{g}+\nabla
p_{c}=0$. Therefore, the equation presented above is consistent with the
expectation that at equilibrium the capillary term should be independent of
the fluid mobilities and the overall momentum equation should reduce to
$0=-\nabla\overline{p}+\rho\boldsymbol{g}-p_{c}\nabla\alpha_{l}$. The
preceding derivation suggests an alternative formulation where density is
defined identically in the clear fluid and porous regions
($\rho=\rho_{l}\alpha_{l}+\rho_{g}\alpha_{g}$) while the expression for
$\boldsymbol{F}_{c}$ becomes

$\boldsymbol{F}_{c}=\begin{cases}\sigma\nabla.\left(\hat{\boldsymbol{n}}_{lg}\right)\nabla\alpha_{l},&\textnormal{in
the clear fluid regions,}\\\
M^{-1}\left(M_{l}\alpha_{g}-M_{g}\alpha_{l}\right)[(\rho_{l}-\rho_{g})\boldsymbol{g}+\nabla
p_{c}]-p_{c}\nabla\alpha_{l},&\textnormal{in the porous regions.}\end{cases}$
(69)