Multiphase flow modeling in multiscale porous media: An open-source micro-continuum approach

A multiphase Darcy-Brinkman approach is proposed to simulate two-phase ﬂow 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 ﬂow in both regions and scales. The crux of the proposed model is that it tends asymptotically towards the Navier-Stokes volume-of-ﬂuid 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 ﬂuid dynamics (CFD) simulation packages into porous multiscale, multiphase systems. The versatility of the solver is illustrated using applications to two-phase ﬂow in a fractured porous matrix and wave interaction with a porous coastal barrier. © 2020 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (


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 CO 2 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 (∼ µm), to the scale of pore-networks and lab columns (∼cm), all the way up to the field scale (∼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(Jenny et al., , 2006;;Efendiev and Hou, 2007).Conversely, alternative algorithms focus on large-scale behaviors and only solve for small-scale behavior when needed (Tomin andLunati, 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).
The ability of the DB equation to handle two scales simultaneously has been used to solve fluid flow in threedimensional 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 microcontinuum 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 twophase 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, capillarygravity 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.

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.

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 where v i is the velocity of phase i. Mass conservation at the fluid-fluid interface yields where ρ i is the density of phase i, w is the velocity of the interface, and 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, v l = v g = w at the fluid/fluid interface.
Momentum conservation in each fluid yields where g is the gravity vector, S i = µ i ∇v i + ∇v T i is the viscous stress tensor, and p i and µ 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 where I is the unity tensor, σ is the fluid-fluid interfacial tension, and κ = ∇ • n lg is the interface curvature.

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.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, where β i is a function defined in V i (i = l, g).As in standard volume averaging theory, we also define a phase averaging operator, The averages defined by Eqs. ( 5) and ( 6) are related through the porosity field φ and saturation field α l .The porosity field φ is defined as (V l + V g )/V, i.e., the volume occupied by both fluids divided by the control volume V, such that: 1, in solid-free regions, ]0; 1[ , in porous regions. (7) The porosity field is the cornerstone of micro-continuum methods because it delineates porous (0 < φ < 1) and solid-free regions (φ = 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 (φ = 0) in the micro-continuum framework is described instead as a very low-porosity, very low-permeability domain (φ ≈ 0).
The saturation field α 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 0, in regions saturated with gas, ]0; 1[ , in unsaturated regions, 1, in regions saturated with liquid. (8) A saturation field α 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 α l +α g = 1 is always valid and α g is deduced from the knowledge of α l .The averaging operators defined by Eqs. ( 5) and (6) are related by 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 p and velocity v are defined as weighted sums of the pressure and Darcy velocity in each fluid phase: ) 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, v = v l + v g .
The governing equations solving for p and 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 α i , v i i , and 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), where A i j is the surface area between the two fluids, A is is the surface area between fluid i and the solid phase, n i j is the normal vector at the fluid-fluid interface pointing from i to j, and 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): 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 v, p, and α l .Because the volume-averaged continuity equations, Eq. ( 12), involve averaged phase velocities v i they must be transformed into equations in terms of the microcontinuum 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: 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: From the definitions of single-field and relative velocities, we can show that v l the saturation equation can be expressed as: In equation ( 15), the wetting phase saturation α l is advected by the single-field velocity v.The third term on the left-hand side is an additional convection term involving the relative velocity 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 v r and this term must be closed.In the forthcoming discussion, we will see that different descriptions of v r are derived for solid-free (φ = 1) and porous regions (0 ≤ φ < 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, 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 singlefield 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): where the two last terms on the right-hand side, are the drag forces exerted by phase k on phase i.In short, D is reflects the friction of fluid i on the solid surface and D i j 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 ≤ φ < 1) or fluids only (φ = 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: where S is the single-field shear stress S = µ ∇v + ∇v T and µ is the average fluid viscosity µ = α l µ l + α g µ g .
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 µk −1 v and a surface tension force F c : Eventually, if the porous structure is immobile, the porosity φ can be removed from the derivatives and the multiscale single-field momentum equation becomes: 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 p, the single-field velocity v, and the wetting fluid saturation α l : 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 v r , the drag force µk −1 v, and the surface tension force 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.

Closure and multi-scale parameters
In the following, we show how the multiscale parameters v r , µk −1 , and 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 α l defined in the two-phase multiscale microcontinuum model.In cells saturated by the wetting phase, α l = 1.In cells that contain the non-wetting phase only, α l = 0. Finally, 0 < α 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 (φ = 1).The mass balance and saturation equations, Eqs. ( 13) and ( 15), remain, therefore, unchanged.The saturation equation with φ = 1 is the equation used in algebraic VOF solvers such as interFoam, the VOF solver of the open-source CFD code OpenFOAM R .
There, the relative velocity 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 where C α is a model parameter used to control the compression of the interface and n lg is mean normal vector.
For low values of C α , 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 α is often chosen between 0 and 4. The mean normal vector n lg is computed by using the gradient of the phase indicator function α l .The relation between these two vectors can be obtained by applying Eq. ( 11) to the liquid phase indicator function 1 l (a function equal to 1 in V l and 0 elsewhere) in solid-free regions such that (Quintard and Whitaker, 1994), Therefore, 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 D ls = D gs = 0. Therefore, the Darcy term in the momentum equation vanishes: 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: 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 κ is approximated by a mean interface curvature κ.Brackbill et al. (1992) assumes that the mean curvature of the interface can be approximated by calculating the divergence of the mean normal vector, κ = ∇• n lg .
Because κ and σ 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): 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 (v i with i = g, l) with separate pressure fields (p i with i = g, l).
The incompressible, immiscible two-phase Darcy model consists of a saturation equation for the wetting phase, a mass balance equation, and two momentum balance equations, one for each phase, these can also be written as, 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 µ 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, ∇ • α i S i i , 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 : 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 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 p, Eq. ( 9), along with the definition of capillary pressure, Eq. ( 32), results in: Based on the equations presented above, the multi-phase Darcy model implies the following expression for v r : 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 v is then formed to derive the continuum scale formulation of the drag force µk −1 v and capillary force F c .This is achieved by summing both phase velocities, Eq. ( 30), and using the pressure gradient relationship, Eq. ( 33).We obtain: The previous equation can be recast into: where M = M l + M g is the total mobility and ρ * = ρ l M l + ρ g M g / M l + M g 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 and 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 ρ * 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 ρ * is replaced by ρ in Eq. ( 36) and the capillary force expression is replaced by: 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 θ 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., 2018Soulaine et al., , 2019)).This is achieved by replacing the mean normal vector n lg at the interface between the clear fluid and the porous regions by a locally modified normal, nlg , that satisfies the condition, nlg = cos θn wall + sin θt wall , where n wall and 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 vr , µk −1 , and 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 /µ i ) and capillary pressure p c (α l ).The relative velocity follows the relation: The single-field relative permeability is given by: , in porous regions. (42) The body force F c describes the capillary forces within a computational cell using: where the modified normal at the fluid-fluid interface is: in solid-free regions, cos θn wall + sin θt wall , at the interface between solid-free and porous regions. (44) Finally, the single-field fluid density is expressed as: in solid-free regions, ρ g M g + ρ l M l M −1 , in porous regions. (45)

Numerical implementation
The two-phase multi-scale micro-continuum model proposed in this paper is implemented in the open-source CFD platform OpenFOAM R 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.

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: In equation ( 46), V and δt stand for the cell volume and time step, respectively.The subscript P indicates values at the cell center.The coefficients a NP account for the influence of neighboring control volumes and primarily include convective and diffusive fluxes across cell faces.K f s 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, S = Vρ k v k P δt .Eq. ( 46) can then be rearranged as: This equation forms a matrix system that results from the momentum equation discretization.The term a P = Vρ k+1 δt + a P + K f s represents the diagonal term of this matrix.Following OpenFOAM R internal notations (Jasak, 1996), the operator H (X) = NP a NP X NP + S is introduced and Eq. ( 47) becomes: 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 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: Further details regarding the discretization procedure in OpenFOAM R 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.

Solution algorithm
The discretized equations are solved using OpenFOAM R in a segregated way.In particular, the pressurevelocity 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 R VOF solver interFoam.The numerical scheme uses the following sequence of steps.First, the saturation equation, Eq. ( 15), is solved explicitly using the OpenFOAM R 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 v and v r are updated according to Eqs. ( 10) and ( 14).Third, the single-field relative permeability k k+1 , density ρ k+1 , and viscosity µ k+1 are updated using the new value of the saturation field, α k+1 l .The surface tension force, F k+1 c , is computed using Eq. ( 43).Fourth, the velocity field v * is calculated by solving implicitly the momentum equation, 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 v * (which does not satisfy the continuity equation, Eq. ( 13)) is corrected.This is achieved by finding (v * * , p * ) that obeys, Based on these two equations, the pressure equation is formulated as and solved implicitly with a generalized method of Geometric-Algebraic Multi-Grid (GAMG) embedded in OpenFOAM R .The corrected velocity 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 (v, p) coupling satisfies mass conservation.
The resulting values are set to (v k+1 , p k+1 ) and, then, the algorithm marches in time as dictated by the imposed Courant-Friedrichs-Lewy (CFL) number.

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).

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.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.

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 1 and 2.

Buckley-Leverett
We first consider the well-established Buckley-Leverett semi-analytical solution for two-phase flow in a horizontal one-dimensional with no capillary effects (4 m long, 2000 cells, φ = 0.5, k −1 0 = 1 × 10 11 m −2 ).In this case, water is injected into an air-saturated reservoir at a constant flow rate with the following boundary conditions: v water = 1 × 10 −5 m s −1 , ∂p inlet ∂x = 0 Pa m −1 , and p outlet = 0 Pa.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.

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): where the symbols are consistent with the ones presented in previous sections.Given the parameters presented in Tables 1 and 2, Eq. ( 54) is solved iteratively to obtain α f ront l = 0.467 and α f ront l = 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.

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, φ = 0.5, k −1 0 = 1 × 10 11 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: p bottom = 0 Pa.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): which can be rearranged to yield: 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 ∂p c ∂α 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.C show the steady state saturation profiles and the resulting equilibrium saturation gradients for both implemented capillary pressure models, respectively.

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 × 10 −4 m s −1 , p outlet = 0 Pa ).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 × 10 −13 to 4 × 10 −13 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 hybridPorousInter-Foam 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.

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 R .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).

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 −1 0 = 1 × 10 20 m −2 ) and allowing them to reach equilibrium for different prescribed contact angles (θ water = 60 • , 90 • , 150 • ).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 θ 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 5 • 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.

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, θ water = 45 • ) 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): where R is the tube's radius.We then numerically simulate the system with hybridPorousInterFoam and in-terFoam, using impermeable porous boundaries with the former (k −1 0 = 1 × 10 20 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 52 • , a small yet significant difference to the imposed contact angle of 45 • .We are not the first to note that interFoam (the standard pore scale multiphase flow solver in OpenFOAM R 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.

Taylor film
We now model the drainage of ethanol (µ eth.=1.2 × 10 −3 Pa s, ρ eth.= 789 kg m −3 ) by air in a 2-D microchannel (800 by 100 µm, 280 by 116 cells, θ eth.= 20 • , injection velocity "U" = 0.4 m/s, p outlet = 0 Pa).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), where Ca is the capillary number defined as Ca = µ eth.U σ .We can solve Eq. ( 58) with the given simulation parameters to obtain a film thickness of 4.35 µm.Simulations of this system performed using hybridPorousInterFoam with impermeable porous boundaries (k −1 0 = 1 × 10 20 m −2 ) and interFoam with conventional boundaries yield a value of 4.50 µm, representing a relative error of about 3% or 0.15 µm.These tests and their results are consistent with numerical simulations reported by Graveleau et al. (2017) and Maes and Soulaine (2019) using interFoam.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, θ oil =45 • , p outlet = 0 Pa) 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 −1 0 = 1 × 10 20 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) .

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.

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).
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 −1 0 = 5 × 10 7 m −2 , φ 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 = ∆ρ(Length Scale) 2 g y σ >> 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).

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% watersaturated microfracture in the former and we inject water into a 90 % air-saturated microfracture in the latter 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.
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.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(Soulaine et al., , 2018)), and solid mechanics (Carrillo and Bourg, 2019) into this framework will be the focus of subsequent papers.

Figure 1 :
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.

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

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

Figure 4 :
Figure 4: Comparison of the time-dependent saturation profiles calculated from our numerical framework and Buckley-Leverett's semianalytical 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 andCshow the semi-analytical (lines) and numerical (symbols) solutions of the system when using the Brooks-Corey and Van Genuchten relative permeability models, respectively.

Figure 5 :
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.

Figure 6 :
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.

Figure 7 :
Figure 7: Simulation setup for oil drainage within a heterogeneous reservoir.The different colored blocks represent the spatially variable permeability field.

Figure 8 :
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 × 10 −13 to 4 × 10 −13 m 2 as shown in the previous figure.

Figure 9 :
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.

Figure 10 :
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 φ =0.001 and k −1 0 = 1 × 10 20 m −2 in hybridPorousInterFoam) and the blue and red colors represent oil and water, respectively.

Figure 11 :
Figure 11: Coastal barrier simulation at different times.The thin white line represents the boundary of the porous solid: φ and k −1 0 are set to 0.5 and 5 × 10 7 m −2 below said line and to 1 and 0 above it, respectively.

Figure 12 :
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, (φ = 1, k −1 0 = 0) from the porous fracture walls (φ = 0.5, k −1 0 = 4 × 10 12 m −2 ) located above and below it.
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, ourmodel 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.
expression (instead of Eq. 38)F c = (ρ * − ρ) g + M −1 M l ∇ α g p c − M g ∇ (α l p c ) = M −1 M l α g − M g α l (ρ l − ρ g ) g + M −1 M l α g − M g α l ∇p c − p c ∇α l = M −1 M l α g − M g α l [(ρ l − ρ g ) g + ∇p c ] − p c ∇α l (B.2) At equilibrium (v l = v g = 0),the multiphase Darcy's law yields (ρ l − ρ g ) g + ∇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 = −∇p + ρg − p c ∇α l .The preceding derivation suggests an alternative formulation where density is defined identically in the clear fluid and porous regions (ρ = ρ l α l + ρ g α g ) while the expression for F c becomes ∇α l , in the clear fluid regions, M −1 M l α g − M g α l [(ρ l − ρ g ) g + ∇p c ] − p c ∇α l , in the porous regions.(B.3)

Table 1 :
Table of Fluid Properties

Table 2 :
Table of Model Parameters