# The impact of sub-resolution porosity on numerical simulations of multiphase flow

Francisco Carrillo, Cyprien Soulaine, Ian Bourg

► To cite this version:

Francisco Carrillo, Cyprien Soulaine, Ian Bourg. The impact of sub-resolution porosity on numerical simulations of multiphase flow. Advances in Water Resources, 2022, 161, pp.104094. 10.1016/j.advwatres.2021.104094. insu-03576254

HAL Id: insu-03576254

https://insu.hal.science/insu-03576254

Submitted on 15 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# The impact of sub-resolution porosity on numerical simulations of multiphase flow

Francisco J. Carrillo<sup>a,*</sup>, Cyprien Soulaine<sup>b</sup>, Ian C. Bourg<sup>a,c</sup>

## ARTICLE INFO

**Keywords:**

Porous media

Sub-resolution porosity

Multiphase flow

Relative permeability

Multi-scale

Micro-continuum

## ABSTRACT

Sub-resolution porosity (SRP) is a ubiquitous, yet often ignored, feature in Digital Rock Physics. It embodies the trade-off between image resolution and field-of-view, and it is a direct result of choosing an imaging resolution that is larger than the smallest pores in a heterogeneous rock sample. In this study, we investigate the impacts of SRP on multiphase flow in porous rocks. To do so, we use our newly developed Multiphase Micro-Continuum model to perform first-of-a-kind direct numerical simulations of two-phase flow in porous samples containing SRP. We show that SRP properties (porosity, permeability, wettability) can impact predicted absolute permeabilities, fluid breakthrough times, residual saturations, and relative permeabilities by factors of 2, 1.5, 3, and 20, respectively. In particular, our results reveal that SRP can function as a persistent connector preventing the formation of isolated wetting fluid domains during drainage, thus dramatically increasing relative permeabilities to both fluids at low saturations. Overall, our study confirms previous evidence that flow within the SRP cannot be disregarded without incurring significant errors in numerical predictions or experimental analyses of multiphase flow in heterogeneous porous media.

## 1. Introduction

The emergence of Digital Rock Physics (DRP) has revolutionized the way we study porous media. It is now possible to directly characterize the pore structure of subsurface systems and perform three-dimensional direct numerical simulations of fluid flow in digital models of rock samples that approach the size of a Representative Elementary Volume (REV). As such, DRP has transformed our capacity to characterize and predict fluid flow in soils, sedimentary rocks, hydrocarbon reservoirs, and engineered porous systems (Mehmani et al., 2020; Han et al., 2020). The computation of rock transport parameters including absolute permeability (Spanne et al., 1994), dispersion coefficients (Bijeljic et al., 2013; Soulaine et al., 2021a), relative permeabilities, and capillary pressures (Raeini et al., 2014; Prodanović et al., 2014) has had direct impacts in the fields of reservoir engineering, hydrology, and CO₂ sequestration (Blunt et al., 2013; Soulaine et al., 2021b).

### 1.1. Rock imaging techniques and sub-resolution porosity

DRP is made possible by advances in high resolution imaging techniques, notably X-ray Microtomography (XCT) (Baker et al., 2012; Singh et al., 2018; Kohanpur and Valocchi, 2020) and focused ion

beam scanning electron microscopy (FIB-SEM) (Cnudde and Boone, 2013; Kelly et al., 2016; Welch et al., 2017; Ruspini et al., 2021). The first method, XCT, involves recording hundreds or thousands of two-dimensional (2-D) X-ray projections through a sample that are then computationally reconstructed to create a 3-D image. This method enables detailed volumetric representations of rock core samples spanning several cubic millimeters with a resolution of about 1 x 1 x 1 micrometers (Wildenschild and Sheppard, 2013; Blunt et al., 2013). The second method, FIB-SEM, involves repeated etching and imaging of a sample through alternating application of focused ion beams and scanning electron microscopy at considerably smaller scales. It yields images spanning ~5 cubic micrometers with an associated resolution of ~5 x 5 x 5 nanometers (Dewers et al., 2012). However, the associated repeated etching procedure means that FIB-SEM is a destructive imaging technique. These two techniques highlight an important limitation of current imaging techniques: the existence of an unavoidable trade-off between image resolution and field-of-view.

The inherent complexity of most natural rocks further complicates the imaging and characterization process. More often than not, rocks such as sandstones, carbonates, and shales exhibit heterogeneities that

<footnote><sup>a</sup> Department of Civil and Environmental Engineering, Princeton University, Princeton, NJ, USA</footnote>

<footnote><sup>b</sup> Institut des Sciences de la Terre d'Orléans, Université d'Orléans, CNRS, BRGM, Orléans, France</footnote>

<footnote><sup>c</sup> High Meadows Environmental Institute, Princeton University, Princeton, NJ, USA</footnote>

<footnote>* Corresponding author.
E-mail address: franjcf@outlook.com (F.J. Carrillo).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** (A) Conceptual representation of the Multiphase Micro-Continuum approach. The image shows an advancing fluid-fluid interface in a system that contains impermeable rock (top), free fluid (center), and a permeable porous medium (bottom). The two immiscible fluids are shown in different shades of blue (left and right) and the inset shows a sample REV over which the model's equations are averaged. (B) SEM image of a shaly sandstone obtained from Peters (2009) showing the distribution of its components: porous clay (dark gray), non-porous sand (light gray), and open pore space (black).</figcaption></figure>

span several length scales (Bear, 1988; Mousavi et al., 2013; Ak-barabadi et al., 2017; Beckingham et al., 2017), some of which cannot be properly resolved by the aforementioned imaging techniques. A common way to simplify the imaging process while accounting for these heterogeneities is to designate a “cutoff” voxel size that re-solves the largest pores within a given rock sample (or other features of interest such as cracks or fractures) while simultaneously acting as a “filter” for any pores smaller than that particular size. These small pores, which are not individually resolved, are then designated “sub-resolution porosity” (SRP) and labeled as a third phase in the rock-pore-SRP system during the eventual segmentation of sample images. The final result is a reconstructed image with an acceptable trade-off between resolution and field of view (Scheibe et al., 2015).

Until recently, and despite its abundance in reconstructed natural rock samples, SRP was generally assumed to have little influence on rock hydraulic properties predicted from flow simulations. Most computational models were based on the simplifying assumption that transport within the SRP is dominated by diffusion and thus contributes negligibly to fluid flow (Haggerty and Gorelick, 1995; Carrera et al., 1998; Gouze et al., 2008; Shabro et al., 2011; Gjetvaj et al., 2015). However, recent studies have shown that this assumption breaks down whenever SRP contributes significantly to the rock's percolating path by forming bridges between resolved pore spaces (Churcher et al., 1991; Tanino and Blunt, 2012; Wu et al., 2019). In these cases, SRP can impact permeability by a factor larger than 2 even when contributing only ~2% of the total porosity (Soulaine et al., 2016). In addition, recent evidence suggests that SRP also can have important impacts on multiphase flow, as shown by observations of dramatic changes in relative permeability curves and overall flow behavior associated with differences in SRP wetting properties in mixed-wet or strongly-wetting rocks (Zou et al., 2018; Rücker et al., 2019; Fan et al., 2020; Garfi et al., 2020).

## 1.2. Multiscale models

A potential path towards resolving the influence of SRP on hydro-logic processes is provided by sustained efforts to develop numerical techniques designed to account for this feature combined with steady advances in high-performance parallel computing (Arbogast, 1993a,b; Moctezuma-Berthier et al., 2004; Javadpour, 2009; Bauer et al., 2011; Jiang et al., 2013b). In particular, the development of multiscale/dual-porosity pore network models (D-PNM) has allowed for relatively fast and accurate assessment of the permeability of rocks containing multiscale heterogeneity (Békri et al., 1995). Classical PNM s rely on approximating the 3-dimensional resolved pore space through a series

of ideally-shaped pore “nodes” and “throats” (Fatt, 1956). The result is a system where the relevant fluid dynamics can be readily solved through idealized equations for flow (Dong and Blunt, 2009; Joekar-Niasar and Hassanizadeh, 2012; Jiang et al., 2012; Blunt et al., 2013; Huang et al., 2016; Suo et al., 2020). In D-PNMs, the presence of SRP is accounted for through the implementation of an additional fine-scale pore network (Ioannidis and Chatzis, 2000; Jiang et al., 2013; Prodanović et al., 2014; Sadeghnejad and Gostick, 2020; Moslemipour and Sadeghnejad, 2020) or through the creation of “micro-links” forming percolation paths between large pores (Bultreys et al., 2015; Xu et al., 2021). Accurate definition of SRP connectivity within these networks remains a challenge (Zhao et al., 2017; Petrovskyy et al., 2020) as, by definition, there are no discernible features from which to inform assignments of pore network topology within the SRP.

The expansion of multiscale models into multiphase flow further complicates matters, as the effects of capillarity and wettability need to be modeled through representative relative permeability and capillary pressure models in order to obtain accurate flow representations within the SRP (Carrillo et al., 2020). For this reason, the few studies that implemented multiphase D-PNMs have relied on the assumption of quasi-static fluid displacement, an assumption valid for simulating flow at low capillary numbers (Mehmani et al., 2013; Bultreys et al., 2015; Xu et al., 2021) and where both phases are effectively set at a given saturation. These studies have leveraged D-PNMs to study how the amount and distribution of SRP affects the relative permeability behavior of artificial rock samples (Mehmani and Prodanović, 2014) and how SRP characterization and connectivity affect the wetting properties of natural rocks (Bultreys et al., 2016; Song et al., 2021; Isah et al., 2020). Unfortunately, due to the simplifying assumptions of D-PNMs outlined above, extension of these studies to dynamic systems with mixed-wet SRP or systems with viscously-dominated flow remains impossible.

The Micro-Continuum approach presents an alternative route to simulating dynamic flow processes in systems with SRP. This approach relies on locally-averaged Navier-Stokes equations that asymptotically approach Darcy's law in regions with SRP and the Navier-Stokes equations in fully resolved pores. This model has proven fairly flexible and has been used to evaluate the effects of static (Knackstedt et al., 2006; Apourvari and Arns, 2014; Scheibe et al., 2015; Soulaine and Tchelepi, 2016; Kang et al., 2019; Singh, 2019), reactive (Soulaine et al., 2017; Noiriel and Soulaine, 2021; Trinchero et al., 2021), and deformable (Carrillo and Bourg, 2019) SRP on the permeability of heterogeneous porous media. Furthermore, through careful consideration of capillary and viscous effects within the SRP (i.e., fluid mobility, relative permeabilities, and capillary pressures), recent investigations have successfully expanded and validated the Micro-Continuum Approach for situations involving the flow of multiple fluids in multiscale
porous media (Soulaine et al., 2018; Carrillo et al., 2020; Carrillo and Bourg, 2021b,a). In this approach, the impact of simplifying model assumptions is greatly reduced relative to the D-PNM approach at the expense of relatively high computational costs. As such, this approach allows for the simulation of dynamic multiscale systems in domain sizes that approach that of an REV.

### 1.3. Objective of this paper

In this study, we leverage the capabilities of the Multiphase Micro-Continuum Approach to systematically examine the influence of SRP properties (permeability, porosity, wettability) on Direct Numerical Simulation predictions of multiphase flow in a digital model of a carbonate rock. In particular, we characterize the rock's absolute permeability, relative permeability curves, residual permeabilities, and fluid breakthrough times on the ~30 mm³ scale of an OCT image. We hypothesize that the SRP properties outlined above have even greater impacts on multiphase flow than on single phase flow, such that their neglect or misrepresentation leads to inaccurate predictions of rock hydraulic properties. To the best of our knowledge, this is the first application of Direct Numerical Simulations to multiphase flow in rock samples containing unresolved porosity and the first computational effort to systematically examine the impacts of SRP wetting properties on the aforementioned rock flow properties.

## 2. Materials and methods

### 2.1. Mathematical model

The Multiphase Micro-Continuum framework for incompressible immiscible flow in rigid porous media consists of three volume-averaged partial differential equations. They describe the conservation and transport of fluid mass (Eq. (1)), fluid saturation (Eq. (2)), and fluid momentum (Eq. (3)). Once implemented in a suitable numerical solver, these equations are used to solve for the single-field pressure ($\rho$), the single-field fluid velocity ($U$), and the wetting-fluid saturation ($\alpha_w$). A full description of the model can be found in Carrillo et al. (2020). Here, we have:

$$ \nabla \cdot U = 0, \qquad (1) $$

$$ \frac{\partial \phi \alpha_w}{\partial t} + \nabla \cdot (\alpha_w U) + \nabla \cdot (\phi \alpha_w \alpha_n U_r) = 0, \qquad (2) $$

$$ \frac{1}{\phi} \left( \frac{\partial \rho U}{\partial t} + \nabla \cdot \left( \frac{\rho}{\phi} U U \right) \right) = - \nabla p + \nabla \cdot S - \mu k^{-1} U + F_c, \qquad (3) $$

where the subscripts w and n refer to the wetting and non-wetting fluids, $\phi$ is the cell porosity, $\rho$ is the single-field density, $\mu k^{-1}$ is the drag coefficient of the unresolved porous media (a function of the cell permeability, saturation, and fluid viscosities), $F_c$ are the capillary forces, and $S = \mu(\nabla U + (\nabla U)^T)$ is the averaged single-field shear stress tensor. Lastly, $U_r$ is the fluid relative velocity, i.e. the cell-averaged difference in velocity between the wetting and non-wetting fluids. Here, gravity is neglected and the phrase “single-field” refers to averaged variables that depend on the local saturation, density, and viscosity of both fluids within a given cell (Maes and Soulaine, 2019).

A key feature of Eqs. (1)–(3) is that they are valid in control volumes that contain any combination of the three relevant phases (porous solid, wetting fluid, non-wetting fluid), meaning that they can be applied to systems that contain both solid-free ($\phi = 1$) and porous regions ($\phi < 1$). Due to the scale separation hypothesis (Whitaker, 1986), this unique set of equations tends towards distinct solutions in solid-free and porous regions. Notably, the single-field momentum equation tends to a solution that can be asymptotically matched to the two-phase

Navier–Stokes equations in solid-free regions and two-phase Darcy’s law in porous regions (Carrillo et al., 2020):

$$ \left\{ \begin{array}{l} \frac{\partial \rho U}{\partial t} + \nabla \cdot (\rho U U) = -\nabla p + \nabla \cdot S + F_c, \quad \text{if } \phi = 1, \\ U = -\frac{k}{\mu} (\nabla p - F_c), \quad \text{if } \phi < 1. \end{array} \right. \qquad (4) $$

As such, the Multiphase Micro-Continuum model is ideally suited for simulating multiphase flow in OCT images that contain SRP, as illustrated schematically in Fig. 1.

The asymptotic matching noted above requires appropriate definitions of the relative velocity $U_r$, drag force $\mu k^{-1}U$, and capillary forces $F_c$. These variables reflect the influence of sub-grid-scale structure and dynamics, including the fluid distribution and the impact of porous micro-structure on flow within the SRP. For this reason, these parameters are defined differently in the solid-free region ($\phi = 1$) and porous regions ($\phi < 1$). In particular, the single-field drag force is negligible in solid-free regions and, in porous regions, depends on absolute ($k_0$) and relative permeabilities ($k_{r,i}$) within the SRP:

$$ \mu k^{-1} = \begin{cases} 0, & \text{if } \phi = 1, \\ k_0^{-1} \left( \frac{k_{r,w}}{\mu_w} + \frac{k_{r,n}}{\mu_n} \right)^{-1}, & \text{if } \phi < 1. \end{cases} \qquad (5) $$

The capillary forces within the solid-free region are proportional to the surface tension $\gamma$ and the curvature of the fluid-fluid interface as described by the Continuum Surface Force formulation (Brackbill et al., 1992). In the porous region, capillary forces are a function of the fluid mobilities ($M_i = k_0 k_{i,r}/\mu_i$; $M = M_w + M_n$) and the average capillary pressure $p_c$:

$$ F_c = \begin{cases} -\gamma \nabla \cdot (\hat{n}_{w,n}) \nabla \alpha_w, & \text{if } \phi = 1, \\ \left[ M^{-1} (M_w \alpha_n - M_n \alpha_w) \left( \frac{\partial p_c}{\partial \alpha_w} \right) - p_c \right] \nabla \alpha_w, & \text{if } \phi < 1, \end{cases} \qquad (6) $$

where the normal at the fluid-fluid interface, $\hat{n}_{w,n}$, is given by

$$ \hat{n}_{w,n} = \begin{cases} -\frac{\nabla \alpha_w}{|\nabla \alpha_w|}, & \text{if } \phi = 1, \\ \cos \theta_p \mathbf{n}_{wall} + \sin \theta_p \mathbf{t}_{wall}, & \text{at the SRP surface.} \end{cases} \qquad (7) $$

Eq. (7) imposes a contact angle $\theta_p$ at the SRP surface following the approach developed by Horgue et al. (2014), where $\mathbf{t}_{wall}$ and $\mathbf{n}_{wall}$ are the tangential and normal directions relative to the SRP surface. The specification of the contact angle at non-porous rock surfaces, $\theta_r$, follows a similar implementation.

The relative fluid velocity is given by:

$$ U_r = \begin{cases} C_a \max(|U|) \frac{\nabla \alpha_w}{|\nabla \alpha_w|}, & \text{if } \phi = 1, \\ \phi^{-1} \begin{bmatrix} -\left(\frac{M_w}{\alpha_w} - \frac{M_n}{\alpha_n}\right) \nabla p \\ + \left(\frac{M_w}{\alpha_w} + \frac{M_n}{\alpha_n}\right) \nabla p_c \\ -\left(\frac{M_w}{\alpha_w} - \frac{M_n}{\alpha_n}\right) p_c \nabla \alpha_w \end{bmatrix}, & \text{if } \phi < 1, \end{cases} \qquad (8) $$

where $C_a$ is an interface compression parameter used in the Volume-of-Fluid method (typically set to values between 1 and 4), and the expression within the SRP is imposed by asymptotic matching to two-phase Darcy’s law (Carrillo et al., 2020).

Closure of the system of equations presented above requires appropriate constitutive models to solve for $p_c$ and $k_{r,i}$ within the SRP. For simplicity, we use the well-known Van Genuchten model (van Genuchten, 1980):

$$ k_{r,n} = (1 - \alpha_w)^{1/2} (1 - \alpha_w^{1/m})^{2m}, \qquad (9) $$

$$ k_{r,w} = \alpha_w^{1/2} (1 - (1 - \alpha_w^{1/m})^m)^2, \qquad (10) $$

$$ p_c = p_{c,0} \left( (\alpha_w)^{-\frac{1}{m}} - 1 \right)^{-m}, \qquad (11) $$
where *m* is a wetting parameter that controls the internal wettability of the SRP and *p*<sub>c,0</sub> is the entry capillary pressure of the SRP. The SRP is internally water-wet if *m* > 1, intermediate wet at *m* = 1, and oil-wet if *m* < 1. Note that the sign of the entry capillary pressure was changed for values of *m* > 1 to prevent unphysical parameterizations where the SRP is both water-wet and oil-wet at the same time.

Finally, we note that the formulation of an accurate interfacial
condition between a solid-free region and an unsaturated porous region
(see Fig. 1) is still an active area of research. Our formulation proposes
a suitable approximation of said condition by: 1) Ensuring mass and
stress continuity across the porous interface, and 2) upscaling the
relevant capillary forces and surface tension force discontinuity into
an apparent contact angle θp. The effects of pore surface roughness,
aperture, interfacial energies, thin-films and adsorption dynamics can
be potentially included and upscaled into θp by using the theory de-
veloped in Wenzel (1936), Whyman et al. (2008), Cassie and Baxter
(1944), Zou et al. (2015) and/or (AlRatrout et al., 2018). However,
the integration of such models is outside the scope of this investigation.
Further discussion on the multiphase interfacial condition at the porous
interface can be found in Carrillo et al. (2020), Soulaine et al. (2018),
Carrillo and Bourg (2021b).

## 2.2. Numerical implementation

The mathematical model presented in Section 2.1 was numerically implemented in OpenFOAM®, a free, parallel, C++ simulation platform that uses the Finite Volume Method to discretize and solve partial differential equations in three-dimensional grids. Mass conservation and incompressibility (Eqs. (1) and (3)) were ensured through the Pressure Implicit Splitting-Operator (PISO) algorithm (Issa, 1986). The evolution of the fluid–fluid interface (Eq. (2)) was solved using the Multidimensional Universal Limiter of Explicit Solution (MULES) algorithm (Márquez and Fich, 2013) and a Piecewise-Linear Interface Calculation (PLIC) compression scheme, both which help sharpen the interface and keep the simulation numerically stable. Numerical stability was also ensured by limiting the global Courant number to a maximum of 0.15 in all our simulations and by using a second order linear-upwind scheme to discretize the remaining equations. Grid convergence analysis results are shown in the Supplemental Information. Extensive validation of the modeling framework is presented in Carrillo et al. (2020) and the open-source implementation is available from the author's GitHub repository (Carrillo, 2020).

## 2.3. *Studied rock sample*

Simulations were performed on a reconstructed 3-D OCT sample of an Estaillades Carbonate rock obtained from Bultreys (2016) through the DigitalRockPortal. This sample was 7 mm in diameter and was scanned with a UGCT's HECTOR micro-CT scanner. This set of images has been used in several previous D-PNM studies (Bultreys et al., 2015, 2016). The sample (1000 by 1000 by 1000 voxels, 3.1 µm per voxel) is ideally suited for our purposes, as it is a mono-mineralic calcite rock containing both intergranular macropores and unresolved intragranular micropores (i.e., SRP). Voxels containing solid rock, resolved pores, and unresolved pores were identified through a 3-phase segmentation procedure following the steps outlined in Bultreys et al. (2015). This yielded a sample with 56.2% solid rock voxels, 11.8% porous voxels, and 32% microporous voxels (Fig. 2).

Due to the computational cost associated with performing direct numerical simulations on such a large physical space, we extracted a 200 by 200 by 200 voxel sub-sample from the original sample in order to perform our simulations. This representative sub-sample was extracted from the Cartesian center of the sample, starting in projection number 750 and ending in projection 950. The computational cost was further reduced by removing all grid cells corresponding to solid rock voxels in the resulting computational mesh, yielding a sample of

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Representative cross-section of the OCT projection of Estaillades carbonate rock used in this study (Bultreys, 2016). (A) Full 2-D view of the sample, which is 7 mm in diameter with a resolution of 3.1 µm per voxel. (B) 500 by 500 pixels cropped sample. (C) Corresponding segmented image. For all figures, black is open pore space, dark gray corresponds to domains that contain SRP, and the lightest color is solid calcite.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Spatial distribution of the SRP (red), pore space (blue), and solid rock (transparent) within the extracted 3-D rock representation (200 by 200 by 200 cells). A) The complete computational mesh. B) The corresponding SRP, which accounts for 21% of the voxels. C) The associated open pore space, which accounts for 40% of the voxels.</figcaption></figure>

about 3.2 million cells (see Fig. 3). In order to maintain adequate mesh resolution while properly representing the mobile fluid–fluid interface within the open pore space, we implemented a dynamic mesh refinement algorithm that allowed the mesh to become up to 16 times finer at said interface. No mesh refinement was carried out within the SRP. Lastly, as is customary for these types of simulations and to properly control the flow rate into the sample, we added two “buffer” regions at the inlet (top) and outlet (bottom) boundaries of our samples. The other 4 sample sides were defined with no-flow boundary conditions.

3. Base simulation setup and upscaling

3.1. Base simulation parameterization

The main purpose of our simulations is to perform a sensitivity analysis of the impact of SRP properties on single and multiphase flow at the scale of the full digital rock image, hereafter referred to as the macroscopic scale. For this, we first introduce a “base” simulation that will be parameterized using experimental values and then used as a template for the systematic variation of SRP properties. This workflow is conceptually similar to the one performed in Hashemi et al. (2021).

Our base simulation involves the injection of oil into a fully-water-
saturated rock sample at a constant rate of 0.1 µL s⁻¹ until the simula-
tion reaches a steady-state. The choice of the labels “oil” and “water”
is arbitrary: our main goal is to examine the flow of two immiscible
and incompressible fluids. The advancing fluid is non-wetting (θp and
θr > 90°) in our base simulation, but the wettability of the solid
by the two fluids is reversed in some of our simulations. The rock
and fluid properties are summarized in Tables 1–2. The subsequent
<table><caption>**Table 1**
Simulated fluid properties. These were kept constant for all simulations.</caption><thead><tr><th>Property</th><th>Value</th></tr></thead><tbody><tr><td>ρ<sub>w</sub></td><td>1000 kg m<sup>-3</sup></td></tr><tr><td>μ<sub>w</sub></td><td>0.001 Pa s</td></tr><tr><td>ρ<sub>n</sub></td><td>800 kg m<sup>-3</sup></td></tr><tr><td>μ<sub>n</sub></td><td>0.1 Pa s</td></tr><tr><td>γ</td><td>0.03 kg s<sup>-2</sup></td></tr></tbody></table>

<table><caption>**Table 2**
Simulated rock and SRP parameters. The second column represents each parameter value used in the base simulation and the third column shows the range over which each parameter was varied. These ranges where chosen to create representative samples of the associated parameter space: from strongly hydrophobic to strongly hydrophilic systems, and from systems with permeable (or impermeable) SRP to systems with no SRP.</caption><thead><tr><th>Property</th><th>Base Value</th><th>Range</th></tr></thead><tbody><tr><td>φ</td><td>0.5</td><td>0 – 1</td></tr><tr><td>k<sub>0</sub></td><td>10<sup>-13</sup> m<sup>2</sup></td><td>10<sup>-12</sup> – 10<sup>-17</sup> m<sup>2</sup></td></tr><tr><td>θ<sub>r</sub></td><td>30°</td><td>30 – 150°</td></tr><tr><td>θ<sub>p</sub></td><td>30°</td><td>30 – 150°</td></tr><tr><td>m</td><td>1</td><td>0.2 – 1.5</td></tr><tr><td>p<sub>c,0</sub></td><td>± 1.35 × 10<sup>4</sup> Pa</td><td>n/a</td></tr></tbody></table>

sensitivity analysis was performed by independently modifying the SRP's porosity ($\phi = 0$ to 1), absolute permeability ($k_0 = 10^{-12}$ to $10^{-17}$ m$^2$), internal wetting properties ($m = 0.2$ to 1.5 in Eqs. (9) and (11)), and the contact angles formed by fluid-fluid interfaces on the external surface of SRP and impermeable rock domains ($\theta_p$ and $\theta_r = 30°$ to $150°$). The decoupling of the contact angle at the SRP and impermeable rock surfaces allows the investigation of mixed-wet systems (Song and Kovscek, 2015; Huang et al., 2016; Akbarabadi et al., 2017) and establishes the possibility of defining a roughness- or saturation-dependent contact angle in future studies (Wenzel, 1936; Whyman et al., 2008). The decoupling of internal ($m$) and surface ($\theta$) wetting properties allows us to differentiate between macroscopic and microscopic wetting effects, where $\theta$ impacts multiphase flow in the solid-free domain and $m$ impacts multiphase flow within the SRP.

Lastly, we carried out additional single-phase flow simulations for each case where we varied the SRP porosity and permeability. This was necessary to calculate each case's absolute permeability and relative permeability curves (see Section 3.2). On average, each multiphase simulation ran for approximately 120 h on ten 28-core Broadwell Xeon nodes.

## 3.2. Calculation of absolute permeability and relative permeability curves

Relative permeabilities were calculated through modification of the upscaling approach presented in Raeini et al. (2014), where macroscopic relative permeability $K_{r,i}$ is defined as the ratio between the apparent permeability $K_i$ calculated from transient, multi-phase flow experiments and the upscaled absolute permeability $K_0$ calculated from steady-state, single-phase flow experiments:

$$K_{r,i} = \frac{K_i}{K_0} = \frac{Q_i / \Delta P_i}{Q_{i,s} / \Delta P_{i,s}} \qquad (12)$$

In Eq. (12), the subscript $i$ identifies properties pertaining to either fluid and the subscript $s$ refers to quantities obtained from single-phase experiments. Furthermore, $Q_i = \int U \cdot n \alpha_i dA$ is the volumetric fluid flow rate of phase $i$ passing through an area $A$ into the porous medium, and $\Delta P_i$ is the pressure drop in phase $i$ across said medium. The latter is defined as follows:

$$\Delta P_i = - \frac{1}{Q_i} \int_{V_f} (-\nabla p + F_c) \cdot U dV_{f,i}, \qquad (13)$$

## SRP Effects on Absolute Permeability

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Sample absolute permeability as a function of SRP properties. Each label shows the parameter that was modified relative to the base simulation. Values indicated to the right of each bar show the percent change in absolute permeability relative to the base simulation described in Section 3.1.</figcaption></figure>

$$= -\frac{1}{Q_i} \int_{V_f} \left( \frac{D}{Dt}(\rho U) - \nabla \cdot S + \mu k^{-1} U \right) \cdot U dV_{f,i},$$

where $V_f$ is the fluid volume of the sample excluding the buffer zones. A drag term ($\mu k^{-1} U$) is included Eq. (13) to account for the momentum dissipation (i.e. pressure drop) induced by the presence of SRP in the sample. The calculation of $\Delta P_{i,s}$ follows Eq. (13) sans the capillary force term.

Relative permeability curves were constructed by matching each $K_{r,i}$ value to the corresponding saturation in the porous medium at a specific point in time. This so-called “unsteady” approach, where $K_{r,i}$ values are not calculated at steady state (Amaefule and Handy, 1982; Johnson et al., 1959), enables calculating relative permeability curves without needing to carry out a distinct steady-state multiphase simulation for each data point, a current necessity across numerical frameworks for rock models with realistically complex pore structures given current computational capabilities. However, this comes at the expense of accuracy or, more precisely, at a risk that the resulting relative permeability curves may be sensitive to fluid flow rate (Diamantopoulos and Durner, 2012). To minimize the impact of this approximation, we focus on characterizing the sensitivity of $K_{r,i}$ to different SRP properties, as opposed to absolute values of $K_{r,i}$.

## 4. Impact of SRP on absolute permeability

In the following three sections, we quantify the effects of SRP properties on the rock's overall absolute permeability (this section), relative permeability curves (Section 5), and time-dependent saturation profiles (Section 7). For the remainder of this study, each simulation case is identified by the variable that is changed with respect to the base simulation established in Section 3.1 and parameterized according to Tables 1 and 2. We now start by evaluating the effect of SRP on the rock's absolute permeability.

Fig. 4 shows that the sample's absolute permeability is overestimated by 57% if the SRP is neglected and assumed to be open pore space ($\phi = 1$) and underestimated by 34% if it is ignored and assumed to be impermeable ($\phi = 0$), where the former's permeability more than doubles the latter. The overall trend in Fig. 4 is fairly intuitive: as the SRP's porosity and/or permeability increase, so does the rock's absolute permeability. This is in line with the findings of Mehmani and Prodanović (2014), and Soulaine et al. (2016). However, whereas some previous studies have observed that SRP can have a disproportionately large impact on permeability, implying that it forms key percolation pathways for single-phase flow (Soulaine et al., 2016), the factor of ~2 impact of SRP on absolute permeability observed here is consistent with the predictions of the well-known Kozeny-Carman (KC) equation,
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Sensitivity of drainage and imbibition relative permeability curves to different SRP properties. (A) Sensitivity to SRP absolute permeability, from $k_0 = 10^{-17}$ to $10^{-12}$ m$^2$. (B) Sensitivity to SRP internal wettability, from oil-wetting ($m < 1$) to water-wetting ($m > 1$). (C) Sensitivity to the external wettability of rock and SRP domains, from water-wetting ($\theta_r = 30^\circ, \theta_p = 30^\circ$), to mixed-wetting ($\theta_r = 150^\circ, \theta_p = 30^\circ$ and $\theta_r = 30^\circ, \theta_p = 150^\circ$), to oil-wetting ($\theta_r = 150^\circ, \theta_p = 150^\circ$). (D) Sensitivity to SRP porosity, from $\phi = 0$ to $\phi = 1$. Unless specified, all parameterized values not indicated in the legend are held constant and equal to the values described in Tables 1 and 2. Each color pair represents the oil (top) and water (bottom) relative permeability curves for a given simulated case. The base simulation is shown in black for all cases.</figcaption></figure>

where $k_{KC}/k_{base} = 0.58, 1$, and $1.68$ for $\phi = 0, 0.5$, and $1$, respectively. This implies that the SRP is relatively uniformly distributed in the studied rock sample (in close agreement with Bultreys et al. (2016)). As examined in the following sections, greater impacts of SRP are observed in systems with multiple fluid phases, where SRP wettability and relative permeability become key factors controlling fluid flow.

# 5. Impact of SRP on relative permeability curves

Changes in sample relative permeability as a function of SRP porosity, wetting properties, and absolute permeability are not particularly intuitive. These often involve non-linear behaviors brought about by the combination of capillary forces and the sample's geometry. Throughout the following discussion we will see that the SRP has two primarily competing effects: it can enhance flow by connecting otherwise-isolated macroscopic flow paths, but it also can reduce flow by being less permeable than the open pore space. We will show that the balance between these two roles is strongly dictated by SRP properties.

The four sets of relative permeability curves present in Fig. 5 exhibit two distinct behaviors reflecting different responses to changes in SRP properties. In one observed behavior, the curves for both fluids shift up (or down) in the same direction with respect to the y-axis. This implies that the sample becomes more (or less) permeable to both phases simultaneously. In the other observed behavior, the water and oil relative permeability curves shift up or down in opposite directions, indicating that an increased permeability to one fluid is associated with a decreased permeability to the other fluid.

## 5.1. Sensitivity to SRP absolute permeability

Fig. 5A demonstrates that an increase in SRP absolute permeability enhances the relative permeability curves of both oil and water. This enhancement occurs in addition to the enhancement in absolute permeability presented in Fig. 4. The enhancement of water relative permeability is entirely expected as the SRP is water-wet (and mostly water-saturated) in this scenario, such that greater SRP permeability naturally facilitates the flow of water. The enhancement of oil permeability is less intuitive. Since oil minimally accesses the SRP in this scenario, this enhancement is likely indirect, i.e., greater SRP permeability facilitates water drainage from the open pore space, which in turn facilitates the flow of oil. In other words, displacement of water from large pores can now occur through the SRP. This process, though relatively inefficient, is facilitated by a higher SRP permeability.

We note, that both effects essentially disappear at SRP permeabilities below $\sim 10^{-17}$ m$^2$ as shown in Fig. SI1 in the Supporting Information. In short, SRP permeability is only important if it is sufficiently large that flow can actually occur within the SRP.

## 5.2. Sensitivity to SRP's internal wettability

Fig. 5B shows that an increase in SRP internal wettability, from oil-wetting ($m < 1$) to water-wetting ($m > 1$), also enhances the flow potential of both fluids. This effect is likely analogous to that observed for SRP absolute permeability: a more hydrophilic SRP should remain more fully water-saturated, and hence more permeable to water (because of the impact of saturation on relative permeability within the SRP, represented using Eqs. (9)–(10)). As in Fig. 5A, this greater ability of water to flow through the SRP indirectly facilitates oil flow, likely by aiding water drainage from the open pore space. The reasoning of this last point is analogous to the one above.
## 5.3. Sensitivity to SRP and rock surface contact angles

Fig. 5C shows that the relative permeability curves shift in opposite directions in response to changes in the external wettability of the rock or SRP surfaces. Specifically, as the pore walls become more hydrophobic, permeability to water decreases, while permeability to oil increases. The impact on oil flow is relatively small, likely because of the partial cancellation of two competing effects: more hydrophobic surfaces should inhibit oil flow by causing this flow to occur preferentially in smaller pores or closer to the pore walls; simultaneously, more hydrophobic surfaces should enhance oil flow by minimizing the tendency towards trapping of oil droplets through capillary effects. Therefore, we posit that a decrease in capillary number (Ca) and/or a decrease in sample homogeneity would likely enhance the trapping effect and may reverse the order of the oil relative permeability curves.

The impact on water flow is larger, a counter-intuitive observation. If water flows predominantly within the SRP, the impact of surface contact angles on water flow should be minimal. Alternatively, if water flows predominantly in the open pore space, surface contact angles should have relatively minor impact on relative permeability to water because of the competing effects noted above in the case of oil. In fact, an increase in water relative permeability with $\theta$ (opposite to that observed here) was reported by Fan et al. (2020), Bakhshian and Hosseini (2019), Zhao et al. (2018). A possible explanation of our results is that residual water flow in our simulated system relies on the combination of SRP and residual macropore water flow, whereas the studies noted above consider systems with no SRP. In such systems, water retained in the open pore space through capillary forces, such as in capillary film coatings on rough pore walls (Tokunaga and Wan, 1997; Khishvand et al., 2016), may represent a key residual flow path that would be eliminated in systems with hydrophobic pore walls.

## 5.4. Sensitivity to SRP porosity

The effects of modulating SRP internal porosity between 0 to 1 are shown in Fig. 5D. The overall magnitude of the relative permeability changes is in close agreement with Mehmani and Prodanović (2014), where the authors found that the addition of pore-clogging SRP can modify the relative permeability of the wetting and non-wetting phases by about a factor of 2. We note, again, that this effect occurs in addition to the significant impact of SRP porosity on absolute permeability presented in Fig. 4.

In addition to this significant influence of SRP porosity on relative permeability, our results also show unexpected complexity. In particular, the impact of SRP porosity on water flow is non-monotonous, with minimum water relative permeabilities observed at either $\phi = 0$ or 1 and larger water relative permeabilities observed at intermediate $\phi$ values. This observation is consistent with the expected trend if residual water flow relies on a combination of both SRP and residual macropore water as suggested above: values $\phi = 0$ or 1 would inhibit water flow by eliminating the SRP water component of these residual flow paths.

# 6. Impact of SRP on residual relative permeability

As noted above, our results strongly suggest that the SRP can function as an efficient and persistent connector between otherwise-disconnected water bodies, particularly at low water saturations. We call this increase in permeability the ‘SRP-enhanced relative permeability’. A key manifestation of this is the persistence of significant relative permeability in the water phase at water saturations below 0.5, in agreement with experimental observations for rocks with significant microporosity (Bennion and Bachu, 2010). In contrast, pore network model simulations of multiphase flow generally predict that relative permeability to water is nearly zero at water saturations below ~0.2 to 0.5 (Prodanović et al., 2014; Huang et al., 2016).

## SRP Effects on Residual Relative Permeability

<figure><img src="image_9.png" /><figcaption>Fig. 6. Steady-state water relative permeability for select cases. Each label shows the only varied parameter with respect to the base simulation. The percentages to the right of the bars show the percent change in residual permeability with respect to the aforementioned base simulation specified in Section 3.1.</figcaption></figure>

A convenient way to characterize this effect is by ranking the relative permeabilities of water once each system has achieved a steady state, as seen in Fig. 6. The overarching trend is clear: increasing the SRP’s permeability to water (increasing absolute permeability or SRP wettability with higher values of $m$) also increases the steady-state relative permeability of said fluid (up to 20 times). The reason for this is not obvious, as higher SRP permeability should lead to higher displacement of the defending fluid, lower residual saturations, and thus, lower (not higher!) steady-state permeabilities. This leads us to believe that increasing the flow capability of the SRP also leads to the creation of enhanced percolation pathways that are persistent and remain connected throughout the sample, even at low water saturations. This enhanced permeability effect is heavily influenced by the SRP pore size, wettability and the resulting entry capillary pressure, where water-favorable entry pressures ($m > 1, p_{c,0} < 0$) increase this effect and non-favorable capillary pressures ($m < 1, p_{c,0} > 0$) diminish it. This phenomenon is consistent with experiments in mixed-wet porous media (AlRatrout et al., 2018) and somewhat analogous to thin-film flow in soils, where small amounts of water facilitate transport above the soil’s water table (Tokunaga and Wan, 1997; Lebeau and Konrad, 2010) through capillary effects. However, our current simulation setup does not allow us to comment on what would happen to this effect as a result of dynamic changes in sample wettability over time. This persistence of significant residual relative permeability to water has potentially important implications in the physics of soil drying (Or et al., 2013) and in hydrocarbon recovery from tight sandstone formations (Tian et al., 2019).

# 7. Impact of SRP on dynamic saturation evolution

The presence of SRP has the following competing effects on the evolution of oil saturation within the sample during oil-flooding: 1) It increases the residual saturation of its wetting phase (be it oil or water) by acting as a fluid reservoir that “defends” itself against the non-wetting phase. 2) It decreases the residual saturation of the defending fluid phase by adding additional inter-pore connectivity and outflow routes (Mehmani and Prodanović, 2014). The balance between these two effects is dictated by the flow properties of the SRP.

Fig. 7 shows that fluid injection into the sample follows two characteristic behaviors: 1) An initial linear increase in saturation, where the slope is primarily dictated by the injection rate. 2) A non-linear plateauing slope dictated by the slow drainage of the defending fluid through the SRP and flow of the injected fluid into the SRP, which are influenced by the SRP’s flow properties. The transition point between these two primary flow mechanisms is dictated by the “breakthrough
<figure><img src="image_10.png" /><figcaption>Fig. 7. The evolution of oil saturation vs time for all studied cases. Each label shows the parameter that was varied with respect to the base simulation. Note that oil breakthrough occurs when about 1/2 of the macropore space still contains water, suggesting that much of the later (and slower) increase in oil content corresponds to water drainage from open pore space. The insets at the top left and bottom right show the final configurations for the cases with the highest and lowest final oil saturations, respectively. In said insets, the blue phase represents oil within open pore space and red represents the SRP that has been invaded by oil.</figcaption></figure>

time”, the point at which the injected fluid first reaches the sample’s
outlet boundary. The next two sections will leverage the information
within the oil-flooding saturation curves in Fig. 7 to study the effects
of the SRP on the dynamic and static properties of these experiments.

7.1. Impact on breakthrough time

We now present a general ranking of the breakthrough times for oil flooding as a function of SRP properties obtained from the results in Fig. 7 (and Figure SI2). The samples are well distributed around the standard base case and obey the following trends: The slowest breakthrough times correspond to cases with oil-wetting surface contact angles, where the oil explores more of the porous medium before reaching the outlet, in agreement with experimental observations of multiphase flow in bead-packs and micromodels with no SRP (Zhao et al., 2016; Hu et al., 2017). These are followed by the sample case with no SRP, where the reasoning is the same as above. Sample cases with a less water-wet SRP (decreasing *m*) or with lower SRP permeability or porosity further decrease the breakthrough times by limiting the ability of water to drain through the SRP, such that the oil explores less of the sample before reaching the outlet. The effects of dynamic contact angle hysteresis is not evaluated here and will be investigated in future work. Overall, our results show that oil breakthrough times are sensitive to SRP parameters (±30%) even though drainage occurs predominantly in the larger pores, a result that has potentially important implications in enhanced oil recovery and geologic CO₂ sequestration.

7.2. Impact on residual saturations

Finally, we observe in Fig. 8 (and Fig. SI3) that residual water saturations are highly correlated with oil breakthrough times: samples with faster breakthrough times generally have higher residual water saturations at steady state. The reasoning behind this behavior is very similar to the one developed to explain the difference in oil breakthrough times: if more residual defending fluid is present, the invading fluid explores less of the available space and hence travels through the sample more rapidly. Overall, this analysis indicates that SRP has a considerable impact on a sample's residual saturations (±400%), strongly implying that it should not be neglected during the design of subsurface fluid extraction and sequestration processes.

<figure><img src="image_11.png" /><figcaption>Fig. 8. Residual water saturations vs breakthrough times for all studied cases. Note that the longer it takes for oil to break through the sample, the lower the final water saturation. This trend holds for all cases with an SRP porosity of $\phi = 0.5$ (blue circles). Changes in SRP porosity have the additional effect of changing the water storage capacity of the sample (red stars).</figcaption></figure>

8. Conclusions

In this paper we studied the effects of XCT Sub-Resolution Porosity (SRP) on a rock's absolute permeability, relative permeability, residual saturations, and fluid breakthrough times. Our results quantify how these four properties react to changes in the porosity, permeability, and wettability of the SRP. One notable finding is that SRP can function as a persistent connector between otherwise-isolated fluid clusters during multiphase flow, even at low saturations. These results were obtained from numerical simulations performed with our newly-developed Multiphase Micro-Continuum framework. To the best of our knowledge, this is the first two-phase flow model and study to take into account SRP without having to rely on a quasi-static assumption or simplified pore-network models.

As such, this investigation establishes a framework for performing two-phase flow simulations in digital rock systems that have two characteristic length-scales. Potential improvements to our methodology include the simulation of larger and more diverse rock samples, a very attainable task due to the current continuous and rapid growth of high-performance computing resources. Finally, our results suggest
potentially fruitful opportunities for future work aimed at quantifying the effects of SRP on upscaled capillary pressure curves, and broadening the investigated parameter space to different types of rocks involving different geometries, different amounts of SRP, and different SRP-induced connectivity. These avenues will more extensively test the conclusions presented in this study and lead the way towards greater understanding of multiscale rock physics and the development of more accurate and predictive upscaled permeability models.

## Nomenclature

<table><tr><td>ρ<sub>i</sub></td><td>Density of phase i (kg/m<sup>3</sup>)</td></tr><tr><td>ρ</td><td>Single-field fluid density (kg/m<sup>3</sup>)</td></tr><tr><td>U</td><td>Single-field fluid velocity (m/s)</td></tr><tr><td>U<sub>r</sub></td><td>Relative fluid velocity (m/s)</td></tr><tr><td>m</td><td>Van-Genuchten wettability parameter</td></tr><tr><td>p</td><td>Single-field fluid pressure (Pa)</td></tr><tr><td>p<sub>c</sub></td><td>Average capillary pressure (Pa)</td></tr><tr><td>S</td><td>Single-field fluid viscous stress tensor (Pa)</td></tr><tr><td>Q</td><td>Volumetric fluid flow rate (m<sup>3</sup>/s)</td></tr><tr><td>ΔP</td><td>Macroscopic pressure difference (Pa)</td></tr><tr><td>γ</td><td>Fluid-fluid interfacial tension (Pa.m)</td></tr><tr><td>φ</td><td>Porosity field</td></tr><tr><td>α<sub>w</sub></td><td>Saturation of the wetting phase</td></tr><tr><td>α<sub>n</sub></td><td>Saturation of the non-wetting phase</td></tr><tr><td>μ<sub>i</sub></td><td>Viscosity of phase i (Pa.s)</td></tr><tr><td>k<sub>0</sub></td><td>SRP absolute permeability (m<sup>2</sup>)</td></tr><tr><td>K<sub>0</sub></td><td>Sample absolute permeability (m<sup>2</sup>)</td></tr><tr><td>k<sub>r,i</sub></td><td>SRP relative permeability for fluid i</td></tr><tr><td>K<sub>r,i</sub></td><td>Sample relative permeability for fluid i</td></tr><tr><td>F<sub>c</sub></td><td>Average capillary forces (Pa/m)</td></tr><tr><td>C<sub>α</sub></td><td>Parameter for the compression velocity model</td></tr><tr><td>M<sub>i</sub></td><td>Mobility of phase i (m<sup>3</sup>/kg.s)</td></tr><tr><td>M</td><td>Total mobility (m<sup>3</sup>/kg.s)</td></tr><tr><td>θ<sub>r</sub></td><td>Rock surface contact angle</td></tr><tr><td>θ<sub>p</sub></td><td>SRP surface contact angle</td></tr><tr><td>n<sub>wall</sub></td><td>Normal vector to the porous surface</td></tr><tr><td>t<sub>wall</sub></td><td>Tangent vector to the porous surface</td></tr><tr><td>p<sub>c,0</sub></td><td>Entry capillary pressure (Pa)</td></tr><tr><td>V<sub>f</sub></td><td>Total volume of fluid in the sample (m<sup>3</sup>)</td></tr><tr><td>V<sub>f,i</sub></td><td>Total volume of fluid i in the sample (m<sup>3</sup>)</td></tr></table>

## CRediT authorship contribution statement

**Francisco J. Carrillo:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. **Cyprien Soulaine:** Conceptualization, Funding acquisition, Project administration, Supervision, Writing – original draft, Writing – review & editing. **Ian C. Bourg:** Conceptualization, Funding acquisition, Project administration, Supervision, Writing – original draft, Writing – review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work was supported by the National Science Foundation, United States, Division of Earth Sciences, Early Career program through Award EAR-1752982. F.J.C. was additionally supported by a Mary and Randall Hack '69 Research Award from the High Meadows Environmental Institute at Princeton University. C.S was sponsored by the French Agency for Research (Agence Nationale de la Recherche, ANR)

through the labex VoltaireANR-10-LABX-100-01 and the grant FraMatI, France ANR-19-CE05-0002. We do not report any conflicts of interest. The code for the computational model used in this manuscript is archived at https://doi.org/10.5281/zenodo.4013969 (Carrillo, 2020) and can also be found at https://github.com/Franjcf. The Estaillades Carbonate rock sample was obtained from Bultreys (2016) through the DigitalRockPortal.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.advwatres.2021.104094.

## References

Akbarabadi, M., Saraji, S., Piri, M., Georgi, D., Delshad, M., 2017. Nano-scale experimental investigation of in-situ wettability and spontaneous imbibition in ultra-tight reservoir rocks. Adv. Water Resour. 107, 160–179.

AlRatrout, A., Blunt, M.J., Bijeljic, B., 2018. Wettability in complex porous materials, the mixed-wet state, and its relationship to surface roughness. Proc. Natl. Acad. Sci. USA 115 (36), 8901–8906.

Amaeule, J.O., Handy, L.L., 1982. Effect of interfacial tensions on relative oil/water permeabilities of consolidated porous media. Soc. Petrol. Eng. J. 22 (3), 371–381.

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. In: International Petroleum Technology Conference, 19-22 January, Doha, Qatar.

Arbogast, T., 1993a. Gravitational forces in dual-porosity systems: I. Model derivation by homogenization. Transp. Porous Media 13 (2), 179–203.

Arbogast, T., 1993b. Gravitational forces in dual-porosity systems: II. Computational validation of the homogenized model. Transp. Porous Media 13 (2), 205–220.

Baker, D.R., Mancini, L., Polacci, M., Higgins, M.D., Gualda, G.A., Hill, R.J., Rivers, M.L., 2012. An introduction to the application of X-ray microtomography to the three-dimensional study of igneous rocks. Lithos 148, 262–276.

Bakhshian, S., Hosseini, S.A., 2019. Pore-scale analysis of supercritical CO2-brine immiscible displacement under fractional-wettability conditions. Adv. Water Resour. 126, 96–107.

Bauer, D., Youssef, S., Han, M., Bekri, S., Rosenberg, E., Fleury, M., Vizika, O., 2011. From computed microtomography images to resistivity index calculations of heterogeneous carbonates using a dual-porosity pore-network approach: Influence of percolation on the electrical transport properties. Phys. Rev. E 84 (1), 1–12.

Bear, J., 1988. Dynamics of Fluids in Porous Media. In: Dover Civil and Mechanical Engineering Series, Dover.

Beckingham, L.E., Steefel, C.I., Swift, A.M., Voltolini, M., Yang, L., Anovitz, L.M., Sheets, J.M., Cole, D.R., Kneafsey, T.J., Mitnick, E.H., Zhang, S., Landrot, G., Ajo-Franklin, J.B., DePaolo, D.J., Mito, S., Xue, Z., 2017. Evaluation of accessible mineral surface areas for improved prediction of mineral reaction rates in porous media. Geochem. Cosmochim. Acta 205, 31–49.

Békri, S., Thovert, J.F., Adler, P.M., 1995. Dissolution of porous media. Chem. Eng. Sci. 50 (17), 2765–2791.

Bennion, D.B., Bachu, S., 2010. Drainage and imbibition CO2/brine relative permeability curves at reservoir conditions for carbonate formations. In: SPE Annual Technical Conference and Exhibition, Paper SPE 134028.

Bijeljic, B., Raeini, A., Mostaghimi, P., Blunt, M.J., 2013. Predictions of non-Fickian solute transport in different classes of porous media using direct simulation on pore-scale images. Phys. Rev. E 87 (1).

Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A., Pentland, C., 2013. Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216.

Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100 (2), 335–354.

Bultreys, T., 2016. Estaillades Carbonate, Vol. 2. Digital Rocks Portal, http://www.digitalrocksportal.org/projects/58.

Bultreys, T., Stappen, J.V., Kock, T.D., Boever, W.D., Boone, M.A., Hoorebeke, L.V., Cnudde, V., 2016. Investigating the relative permeability behavior of microporosity-rich carbonates and tight sandstones with multiscale pore network models. J. Geophys. Res. Solid Earth 121 (11), 7929–7945.

Bultreys, T., Van Hoorebeke, L., Cnudde, V., 2015. Multi-scale, micro-computed tomography-based pore network models to simulate drainage in heterogeneous rocks. Adv. Water Resour. 78, 36–49.

Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guinerà, J., 1998. On matrix diffusion: Formulations, solution methods and qualitative effects. Hydrogeol. J. 6 (1), 178–190.

Carrillo, F.J., 2020. HybridPorousInterFoam.

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 Resour. Res. 55 (10), 8096–8121.
Carrillo, F.J., Bourg, I.C., 2021a. Capillary and viscous fracturing during drainage in porous media. Phys. Rev. E 103, 063106.

Carrillo, F.J., Bourg, I.C., 2021b. Modeling multiphase flow within and around deformable porous materials: A Darcy-Brinkman-Biot approach. Water Resour. Res. 57 (2), e2020WR028734.

Carrillo, F.J., Bourg, I.C., Soulaine, C., 2020. Multiphase flow modeling in multiscale porous media: An open-source micro-continuum approach. J. Comput. Phys. X 8, 100073.

Cassie, A., Baxter, S., 1944. Wettability of porous surfaces. Trans. Faraday Soc. 40, 546–551.

Churcher, P.L., French, P.R., Shaw, J.C., Schramm, L.L., 1991. Rock properties of berea sandstone, baker dolomite, and Indiana limestone. In: SPE International Symposium on Oilfield Chemistry. Society of Petroleum Engineers, pp. 431–446.

Cnudde, V., Boone, M.N., 2013. High-resolution X-ray computed tomography in geosciences: A review of the current technology and applications. Earth-Sci. Rev. 123, 1–17.

Dewers, T.A., Heath, J., Ewy, R., Duranti, L., 2012. Three-dimensional pore networks and transport properties of a shale gas formation determined from focused ion beam serial imaging. Int. J. Oil Gas Coal Technol. 5, 229–248.

Diamantopoulos, E., Durner, W., 2012. Dynamic nonequilibrium of water flow in porous media: A review. Vadose Zone J. 11, v2j2011.0197.

Dong, H., Blunt, M.J., 2009. Pore-network extraction from micro-computerized-tomography images. Phys. Rev. E 80 (3), 1–11.

Fan, M., McClure, J.E., Armstrong, R.T., Shabaninejad, M., Dalton, L.E., Crandall, D., Chen, C., 2020. Influence of clay wettability alteration on relative permeability. Geophys. Res. Lett. 47 (18), 1–10.

Fatt, I., 1956. The network model of porous media. Petrol. Trans. AIME 207, 144–181.

Garfi, G., John, C.M., Lin, Q., Berg, S., Krevor, S., 2020. Fluid surface coverage showing the controls of rock mineralogy on the wetting state. Geophys. Res. Lett. 47 (8), 1–9.

Gjetvaj, F., Russian, A., Gouze, P., Dentz, M., 2015. Dual control of flow field heterogeneity and immobile porosity on non-fickian transport in berea sandstone. Water Resour. Res. 51 (10), 8273–8293.

Gouze, P., Melean, Y., Le Borgne, T., Dentz, M., Carrera, J., 2008. Non-fickian dispersion in porous media explained by heterogeneous microscale matrix diffusion. Water Resour. Res. 44 (11).

Haggerty, R., Gorelick, S.M., 1995. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour. Res. 31 (10), 2383–2400.

Han, J., Han, S., Kang, D.H., Kim, Y., Lee, J., Lee, Y., 2020. Application of digital rock physics using X-ray CT for study on alteration of macropore properties by CO₂ EOR in a carbonate oil reservoir. J. Pet. Sci. Eng. 189, 107009.

Hashemi, L., Blunt, M., Hajibeygi, H., 2021. Pore-scale modelling and sensitivity analyses of hydrogen-brine multiphase flow in geological porous media. Nat. News 11, 8348.

Horgue, P., Prat, M., Quintard, M., 2014. A penalization technique applied to the "volume-of-fluid" method: Wettability condition on immersed boundaries. Comput. & Fluids 100, 255–266.

Hu, R., Wan, J., Kim, Y., Tokunaga, T.K., 2017. Wettability effects on supercritical CO₂-brine immiscible displacement during drainage: Pore-scale observation and 3D simulation. Int. J. Greenh. Gas Control 60, 129–139.

Huang, X., Bandilla, K.W., Celia, M.A., 2016. Multi-physics pore-network modeling of two-phase shale matrix flows. Transp. Porous Media 111, 123–141.

Ioannidis, M.A., Chatzis, I., 2000. A dual-network model of pore structure for vuggy carbonates. In: Proceedings of the Annual Symposium of the Society of Core Analysts.

Isah, A., Adebayo, A.R., Mahmoud, M., Babalola, L.O., El-Husseiny, A., 2020. Drainage mechanisms in gas reservoirs with bimodal pores – a core and pore scale study. J. Nat. Gas Sci. Eng. 86, 103652.

Issa, R.I., 1986. Solution of the implicitly discretised fluid flow equations by operator-splitting. J. Comput. Phys. 62 (1), 40–65.

Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J. Can. Pet. Technol. 48 (8), 16–21.

Jiang, Z., van Dijke, M.I., Geiger, S., Kronbauer, D., Mantovani, I.F., Fernandes, C.P., 2013. Impact of the spatial correlation of microporosity on fluid flow in carbonate rocks. In: Society of Petroleum Engineers - SPE Reservoir Characterisation and Simulation Conference and Exhibition, RCSC 2013: New Approaches in Characterisation AndModelling of Complex Reservoirs, Vol. 2. pp. 637–646.

Jiang, Z., van Dijke, M.I., Sorbie, K.S., Couples, G.D., 2013b. Representation of multiscale heterogeneity via multiscale pore networks. Water Resour. Res. 49 (9), 5437–5449.

Jiang, Z., van Dijke, M.I., Wu, K., Couples, G.D., Sorbie, K.S., Ma, J., 2012. Stochastic pore network generation from 3D rock images. Transp. Porous Media 94 (2), 571–593.

Joekar-Niasar, V., Hassanizadeh, S.M., 2012. Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. Crit. Rev. Environ. Sci. Technol. 42 (18), 1895–1976.

Johnson, E.F., Bossler, D., Bossler, V.N., 1959. Calculation of relative permeability from displacement experiments. Trans. AIME 216 (01), 370–372.

Kang, D.H., Yang, E., Yun, T.S., 2019. Stokes-brinkman flow simulation based on 3-D μ-CT images of porous rock using grayscale pore voxel permeability. Water Resour. Res. 55 (5), 4448–4464.

Kelly, S., El-Sobky, H., Torres-Verdín, C., Balhoff, M.T., 2016. Assessing the utility of FIB-SEM images for shale digital rock physics. Adv. Water Resour. 95, 302–316.

Khishvand, M., Alizadeh, A.H., Piri, M., 2016. In-situ characterization of wettability and pore-scale displacements during two- and three-phase flow in natural porous media. Adv. Water Resour. 97, 279–298.

Knackstedt, M.A., Arns, C.H., Ghou, A., Sakellariou, A., Senden, T.J., Sheppard, A.P., Sok, R.M., Nguyen, V., Pinczewski, W.V., 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, pp. 1–15.

Kohanpur, A.H., Valocchi, A.J., 2020. Pore-network stitching method: A pore-to-core upscaling approach for multiphase flow. Transp. Porous Media 135, 659–685.

Lebeau, M., Konrad, J.M., 2010. A new capillary and thin film flow model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 46 (1), 1–15.

Maes, J., Soulaine, C., 2019. A unified single-field volume-of-fluid-based formulation for multi-component interfacial transfer with local volume changes. J. Comput. Phys. 402, 109024.

Márquez, S., Fich, D., 2013. An Extended Mixture Model for the Simultaneous Treatment of Short and Long Scale Interfaces (Ph.D. thesis). Universidad Nacional del Lateral.

Mehmani, A., Kelly, S., Torres-Verdín, C., 2020. Leveraging digital rock physics workflows in unconventional petrophysics: A review of opportunities, challenges, and benchmarking. J. Pet. Sci. Eng. 190, 107083.

Mehmani, A., Prodanović, M., 2014. The effect of microporosity on transport properties in porous media. Adv. Water Resour. 63, 104–119.

Mehmani, A., Prodanović, M., Javadpour, F., 2013. Multiscale, multiphysics network modeling of shale matrix gas flows. Transp. Porous Media 99 (2), 377–390.

Moctezuma-Berthier, A., Vizika, O., Thovert, J.F., Adler, P.M., 2004. One- and two-phase permeabilities of vugular porous media. Transp. Porous Media 56 (2), 225–244.

Moslemipour, A., Sadeghnejad, S., 2020. Dual-scale pore network reconstruction of vugular carbonates using multi-scale imaging techniques. Adv. Water Resour. 147, 103795.

Mousavi, M., Prodanovic, M., Jacobi, D., 2013. New classification of carbonate rocks for process-based pore-scale modeling. SPE J. 18 (2), 243–263.

Noiriel, C., Soulaine, C., 2021. Pore-scale imaging and modelling of reactive flow in evolving porous media: tracking the dynamics of the fluid-rock interface. Transp. Porous Media.

Or, D., Lehmann, P., Shahraeeni, E., Shokri, N., 2013. Advances in soil evaporation physics-a review. Vadose Zone J. 12, v2j2012–0163.

Peters, C.A., 2009. Accessibilities of reactive minerals in consolidated sedimentary rock: An imaging study of three sandstones. Chem. Geol. 265 (1–2), 198–208.

Petrovskyy, D., Dijke, M.v., Jiang, Z., Geiger, S., 2020. Phase connectivity in pore-network models for capillary-driven flow. Adv. Water Resour. 147, 103776.

Prodanović, M., Mehmani, A., Sheppard, A.P., 2014. Imaged-based multiscale network modelling of microporosity in carbonates. Geol. Soc. Lond. Special Publ. 406 (1), 95–113.

Raeini, A.Q., Blunt, M.J., Bijeljic, B., 2014. Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces. Adv. Water Resour. 74, 116–126.

Rücker, M., Bartels, W.B., Singh, K., Brussee, N., Coorn, A., van der Linde, H.A., Bonnin, A., Ott, H., Hassanizadeh, S.M., Blunt, M.J., Mahani, H., Georgiadis, A., Berg, S., 2019. The effect of mixed wettability on pore-scale flow regimes based on a flooding experiment in ketton limestone. Geophys. Res. Lett. 46 (6), 3225–3234.

Ruspini, L.C., Øren, P.E., Berg, S., Masalmeh, S., Bultreys, T., Taberner, C., Sorop, T., Marcelis, F., Appel, M., Freeman, J., et al., 2021. Multiscale digital rock analysis for complex rocks. Transp. Porous Media 2, 1573–1634.

Sadeghnejad, S., Gostick, J., 2020. Multiscale reconstruction of vuggy carbonates by pore-network modeling and image-based technique. SPE J. 25, 253–267.

Scheibe, T.D., Perkins, W.A., Richmond, M.C., McKinley, M.I., Romero-Gomez, P.D.J., Oostrom, M., Wietisma, 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 Resour. Res. 51 (2), 1023–1035.

Shabro, V., Torres-Verdín, C., Javadpour, F., 2011. Numerical simulation of shale-gas production: From pore-scale modeling of slip-flow, knudsen diffusion, and langmuir desorption to reservoir modeling of compressible fluid. In: Society of Petroleum Engineers - SPE Americas Unconventional Gas Conference 2011, UGC 2011. Society of Petroleum Engineers (SPE), pp. 706–716.

Singh, K., 2019. How hydraulic properties of organic matter control effective liquid permeability of mudrocks. Transp. Porous Media 19 (3), 761–777.

Singh, K., Menke, H., Andrew, M., Rau, C., Bijeljic, B., Blunt, M.J., 2018. Data descriptor: Time-resolved synchrotron X-ray microtomography datasets of drainage and imbibition in carbonate rocks. Sci. Data 5 (1), 1–8.

Song, W., Kovscek, A.R., 2015. Functionalization of micromodels with kaolinite for investigation of low salinity oil-recovery processes. Lab Chip 15 (16), 3314–3325.
Song, W., Yao, J., Zhang, K., Sun, H., Yang, Y., 2021. The impacts of pore structure and relative humidity on gas transport in shale: A numerical study by the image-based multi-scale pore network model. Transp. Porous Media 2, 1573-1634.

Soulaine, C., Girolami, L., Arbaret, L., Roman, S., 2021a. Digital rock physics: Computation of hydrodynamic dispersion. Oil Gas Sci. Technol. 76, 51.

Soulaine, C., Gjetvaj, F., Garing, C., Roman, S., Russian, A., Gouze, P., Tchelepi, H.A., 2016. The impact of sub-resolution porosity of X-ray microtomography images on the permeability. Transp. Porous Media 113 (1), 227-243.

Soulaine, C., Maes, J., Roman, S., 2021b. Computational microfluidics for geosciences. Front. Water 3, 1-11.

Soulaine, C., Roman, S., Kovscek, A., Tchelepi, H.A., 2017. Mineral dissolution and wormholing from a pore-scale perspective. J. Fluid Mech. 827, 457-483.

Soulaine, C., Roman, S., Kovscek, A., Tchelepi, H.A., 2018. Pore-scale modelling of multiphase reactive flow: application to mineral dissolution with production of J. Fluid Mech. 855, 616-645.

Soulaine, C., Tchelepi, H.A., 2016. Micro-continuum approach for pore-scale simulation of subsurface processes. Transp. Porous Media 113 (3), 431-456.

Spanne, P., Thovert, J.F., Jacquin, C.J., Lindquist, W.B., Jones, K.W., Adler, P.M., 1994. Synchrotron computed microtomography of porous media: Topology and transports. Phys. Rev. Lett. 73 (14), 2001-2004.

Suo, S., Liu, M., Gan, Y., 2020. An LBM-PNM framework for immiscible flow: With applications to droplet spreading on porous surfaces. Chem. Eng. Sci. 218, 115577.

Tanino, Y., Blunt, M.J., 2012. Capillary trapping in sandstones and carbonates: Dependence on pore structure. Water Resour. Res. 48 (8), 8525.

Tian, L., Feng, B., Zheng, S., Gu, D., Ren, X., Yang, D., 2019. Performance evaluation of gas production with consideration of dynamic capillary pressure in tight sandstone reservoirs. J. Energy Resour. Technol. 141, 022902.

Tokunaga, T.K., Wan, J., 1997. Water film flow along fracture surfaces of porous rock. Water Resour. Res. 33, 1287-1295.

Trinchero, P., Iraola, A., Bruines, P., Gylling, B., 2021. Water-mineral reactions in a translated single realistic fracture: Consequences for contaminant uptake by matrix diffusion. Water Resour. Res. 57 (10), e2021WR030442.

van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Am. J. 44 (5), 892-898.

Welch, N.J., Gray, F., Butcher, A.R., Boek, E.S., Crawshaw, J.P., 2017. High-resolution 3D FIB-SEM image analysis and validation of numerical simulations of nanometre-scale porous ceramic with comparisons to experimental results. Transp. Porous Media 118 (3), 373-392.

Wenzel, R.N., 1936. Resistance of solid surfaces to wetting by water. Ind. Eng. Chem. 28 (8), 988-994.

Whitaker, S., 1986. Flow in porous media I: A theoretical derivation of Darcy's law. Transp. Porous Media 1 (1), 3-25.

Whyman, G., Bormashenko, E., Stein, T., 2008. The rigorous derivation of Young, Cassie-Baxter and Wenzel equations and the analysis of the contact angle hysteresis phenomenon. Chem. Phys. Lett. 450 (4-6), 355-359.

Wildenschild, D., Sheppard, A.P., 2013. X-ray Imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems. Adv. Water Resour. 51, 217-246.

Wu, Y., Tahmasebi, P., Lin, C., Munawar, M.J., Cnudde, V., 2019. Effects of micropores on geometric, topological and transport properties of pore systems for low-permeability porous media. J. Hydrol. 575, 327-342.

Xu, Z., Lin, M., Jiang, W., Ji, L., Cao, G., 2021. Rapid multiscale pore network modeling for drainage in tight sandstone. J. Pet. Sci. Eng. 204, 108682.

Zhao, J., Kang, Q., Yao, J., Viswanathan, H., Pawar, R., Zhang, L., Sun, H., 2018. The effect of wettability heterogeneity on relative permeability of two-phase flow in porous media: A lattice Boltzmann study. Water Resour. Res. 54 (2), 1295-1311.

Zhao, B., MacMinn, C.W., Juanes, R., 2016. Wettability control on multiphase flow in patterned microfluidics. Proc. Natl. Acad. Sci. 113, 10251-10256.

Zhao, B., Shang, Y., Jin, L., Jia, B., 2017. Characterizing connectivity of multiscale pore structure in unconventional reservoirs by the complex network theory. OnePetro.

Zou, S., Armstrong, R.T., Arns, J.Y., Arns, C.H., Hussain, F., 2018. Experimental and theoretical evidence for increased ganglion dynamics during fractional flow in mixed-wet porous media. Water Resour. Res. 54 (5), 3277-3289.

Zou, L., Jing, L., Cvetkovic, V., 2015. Roughness decomposition and nonlinear fluid flow in a single rock fracture. Int. J. Rock Mech. Min. Sci. 75, 102-118.