Transport in Porous Media (2024) 151:1729–1754

https://doi.org/10.1007/s11242-024-02092-x

Computing Relative Permeability and Capillary Pressure

of Heterogeneous Rocks Using Realistic Boundary Conditions

AbdAllah A. Youssef1,2 · Qi Shao1 · S. K. Matthäi1

Received: 26 January 2024 / Accepted: 17 May 2024 / Published online: 19 June 2024
© The Author(s) 2024

Abstract
Relative permeability and capillary pressure are key parameters in multiphase flow mod-
elling. In heterogeneous porous media, flow direction- and flow-rate dependence result
from non-uniform saturation distributions that vary with the balance between viscous,
gravitational, and capillary forces. Typically, relative permeability is measured using con-
stant inlet fractional-flow—constant outlet fluid pressure conditions on samples mounted
between permeable porous plates to avoid capillary end-effects. This setup is replicated in
numeric experiments but ignores the extended geologic context beyond the sample size,
impacting the saturation distribution and, consequently, the upscaled parameters. Here, we
introduce a new workflow for measuring effective relative permeability and capillary pres-
sure at the bedform scale while considering heterogeneities at the lamina scale. We harness
the flexibility of numeric modelling to simulate continuum-REV-scale saturation distribu-
tions in heterogeneous rocks eliminating boundary artefacts. Periodic fluid flux boundary
conditions are applied in combination with arbitrarily oriented, variable-strength pressure
gradient fields. The approach is illustrated on a periodic model of cross-bedded sandstone.
Stepping saturation while applying variable-strength pressure-gradient fields with differ-
ent orientations, we cover the capillary-viscous force balance spectrum of interest. The
obtained relative permeability and capillary pressure curves differ from ones obtained with
traditional approaches highlighting that the definition of force balances needs considera-
tion of flow direction as an additional degree of freedom. In addition, we discuss when the
common viscous and the capillary limits are applicable and how they vary with flow direc-
tion in the presence of capillary interfaces.

Keywords Bedform · Capillary jump · Lamina scale · Cross flow · Effective properties ·
Capillary end effect

* AbdAllah A. Youssef
abdallah_youssef@ymail.com
1
Peter Cook Centre for CCS Research &, Department of Infrastructure Engineering, The University
of Melbourne, Parkville, VIC 3010, Australia
2
Centre for Integrative Petroleum Research, King Fahd University of Petroleum & Minerals
(KFUPM), 31261 Dhahran, Saudi Arabia

13
Vol.:(0123456789)
1730 A. A. Youssef et al.

Abbreviations

Subscript
w Water phase
α Phase ­(CO2 or water)
m Time level

Superscript
- Upscaled

Symbols
a Coefficients of discretized pressure equation
G Pressure gradient vector (ML−2T−2)
g Pressure gradient magnitude (ML−2T−2)
i Unit vector
k Permeability tensor (L2)
kr Relative permeability tensor (L2 L−2)
p, pc Pressure and capillary pressure (ML−1T−2)
R Right-hand side in discretized pressure equation
s Saturation (L3 L−3)
v Darcy velocity vector (L T−1)
vol Fluid volume
θ Polar angle (rad)
𝜑 Azimuth angle (rad)
∅ Porosity(L3 L−3)
𝜇 Fluid viscosity (ML−1T−1)

1 Introduction

Supporting global decarbonization by carbon geosequestration, the efficient utilization of

the pore space in the storage formations is essential. Storage design relies on performance
assessment conducted on digital twins of candidate sites and involves simulation of mul-
tiphase flow through highly heterogeneous rock masses (e.g. Ringrose 2020). Detailed field
characterization of outcrop analogues of storage sites (e.g. Willis and White 2000) reveals
a nested structure of geobodies with internal bedforms ranging from mm-scale lamina to
hm-scale foresets (Ringrose and Bentley 2015). The sedimentary sequences of interest typ-
ically have geometrically complex facies architectures (cf., Ahmed et al. 2014; Sech et al.
2009) with large aspect ratio horizons, lateral facies variations, and pinchouts (Fabuel-
Perez et al. 2010; Gross et al. 2023). High-resolution digital outcrop models reveal that,
with increasing scale, more and more geologic features become important (e.g. Willis and
White 2000). In sufficiently large outcrops, a nested arrangement of rock types is visible,
each with its unique multiphase flow properties and often separated from one another by
sharp interfaces. Such discontinuities act as capillary barriers, filtering non-wetting fluid
that collects upstream, eventually blocking any cross flow (Matthai and Burney 2018; Revil
et al. 1998; Shosa and Cathles 2001; Tran et al. 2020).
The spatial continuity of high permeability streaks within nested rock sequences
gives rise to a wide variation of flow velocity. Due to the nesting, and common

13
Computing Relative Permeability and Capillary Pressure of… 1731

power-law size—frequency distributions, which imply that high permeability streaks

exist on multiple length scales, these features cannot be neatly separated by scale. The
ensuing spatial correlation structure of permeability also implies a poor correlation, if
any, between permeability and flow velocity. This explains why fast, highly focused flow
can occur kilometres away from injector wells (Flemings et al. 2002; Shao et al. 2022).
While it was known already at the beginning of the twentieth century that nested
geo-heterogeneity gives rise to non-uniform macroscopic saturation distributions (e.g.
Krampert 1934), their importance for multiphase fluid flow has only recently received
more attention by studies like Perrin and Benson (2010), Zhang et al. (2013) and Jack-
son et al. (2020). Many open questions remain. Reasons for this lack of understanding
include that relevant heterogeneities are too large and saturation patterns within these
evolve too slowly to be amenable to laboratory studies. The relevant heterogeneities are
also too geometrically complex for conventional reservoir simulation, and they cannot
be separated into neat categories of scale due to their nested nature (cf., Oden et al.
2006).
Explicit representation of the nested heterogeneity of fluvio-deltaic sequences has been
attempted in two-dimensional high-resolution digital outcrop models (Willis et al. 1999;
Willis and White 2000). Advanced numeric simulation (Matthai and Burney 2018) reveals
the strong impact of this heterogeneity on sweep efficiency corroborating field observations
from Sleipner (Ringrose 2020). Simultaneously, coarse–fine intercalations in composite
rock types (e.g. Mishra et al. 2020) control micro-displacement efficiency (Jackson et al.
2020; Pickup and Hern 2002; Trevisan et al. 2017). Boon and Benson (2021) and Boon
et al. (2022) established that the laminated structure of composite quartz sandstones gives
rise to rate- and directional dependence of relative permeability. In field-data based simula-
tions, Shao et al. (2022) demonstrated that conditions triggering this flow-rate dependent
behaviour can occur far from injector wells. They also showed that rate-dependent relative
permeability has a marked effect on plume spreading and saturation distributions within
plumes, confirming earlier predictions by Pickup and Sorbie (1994) and Durlofsky (1997).
Yet, today, multiphase flow in heterogeneous media is still modelled without considering
rate dependence, attempting to account for its effects via permeability anisotropy or other
model adjustments (e.g. Li et al. 2023).
Among others, Virnovsky et al. (2004) and most recently Matthai and Tran (2023) dem-
onstrated that this flow-rate dependence is not due to capillary end effects but arises from
non-uniform saturation distributions caused by local variations in the balance of viscous-
gravitational and capillary forces. Since this behaviour is complex, analysis is still in its
infancy (Darman et al. 2001). However, indirectly acknowledging its practical importance,
oil and gas companies carry out special core analysis (SCAL, e.g. Huang and Honarpour
1998) at flow velocities deemed representative for the reservoir of interest (e.g. Breston and
Hughes 1949). Clearly, this is no satisfactory solution when flow velocity varies widely
across a reservoir. What can be done about this?
Notably, three distinct approaches exist for estimating ensemble saturation functions of
heterogeneous rocks: special core analysis (SCAL), mathematical analysis, and numeri-
cal simulation. SCAL is most widely used and porous plates serve to address capillary
end effects (e.g. Osoba et al. 1951; Huang and Honarpour 1998). Despite recent advances
in inverse modelling with regard to the handling of capillary end effects (An et al. 2023;
Moreno et al. 2021), there still are boundaries artefacts (cf., Jackson et al. 2018) and sam-
ple dimensions are limited (e.g. Berg and Ott 2012). Mathematical analysis can be used to
infer saturation patterns without boundary artefacts (e.g. Ekrann et al. 1996; Rabinovich
et al. 2016) but is limited to capillary limit (CL) or viscous limit (VL) assumptions.

13
1732 A. A. Youssef et al.

Numerical, flow-simulation-based analysis is becoming more and more attractive

because image segmentation and model building techniques are making it easier to build
digital twins of rock masses with a complex internal structure. Simultaneously, the physi-
cal realism of numeric methods is improving. The geologic context of the sample volume
can also be accounted for by non-local boundary conditions (e.g. Durlofsky 2005).
Numerical approaches for generating non-uniform saturation distributions fall into two
main categories: steady-state and unsteady-state methods. Steady-state ensemble param-
eters are evaluated once the saturation distribution has reached a time-invariant state. Sim-
ilar to dynamic core flooding experiments, unsteady analysis quantifies parameters from
evolving saturation patterns and fronts. These parameters often bear the prefix pseudo
(Barker and Dupouy 1996; Barker and Thibeau 1997).
The aim of the present study is to establish a more physically realistic simulation frame-
work and setup for the computation of non-uniform saturation distributions in nested het-
erogeneous porous media including material interfaces. Relative permeability and capillary
pressure curves ought to be computed for realistic in situ conditions and force balances
which arise locally from injection measures and fluid properties. For this purpose, the
promising use of periodic boundary conditions for the saturation equation (Virnovsky et al.
2004; Matthai and Tran 2023) is extended to the pressure equation, choosing a solution
approach that admits arbitrarily oriented far-field fluid pressure gradients. Rather than try-
ing to replicate physical experiments, this approach builds on the versatility of numeric
simulation methods and the possibilities they offer to control the simulation setup.
After a brief review of numerical relative permeability upscaling, this article pro-
ceeds with a rationale for the proposed periodic boundary conditions. Then, test cases are
described, and simulation results presented focussing on cross-bedded sandstone. A final
discussion focusses on less intuitive features of the results obtained. The simulation models
and results are included as digital supplementary material.

2 Methodology

The numeric steady-state upscaling technology employed here builds on the embedded dis-
continuity method for multiphase flow modelling presented by Nick and Matthäi (2011a,
2011b), Tran et al., (2020) and Matthai and Tran (2023). Thus, we only detail the math-
ematical modelling and implementation of the periodic boundary conditions in this col-
located finite element – finite volume discretisation with embedded discontinuities (Mat-
thai and Tran 2023) here as it underpins this new upscaling approach. We start by a brief
review of relevant aspects of numerical upscaling.

2.1 Background

Steady-state numeric upscaling methods can be distinguished by the applied boundary con-
ditions (Odsæter et al. 2015). Constant inlet fractional flow methods go back to physical
experimentation done at Penn-State University in the 50’s: Co-injection of fluid phases is per-
formed at a constant fractional flow while outlet pressure is fixed, all remaining boundaries are
sealed (Geffen et al. 1951; Osoba et al. 1951). While simplifying measurement, this method’s
attractiveness is limited by a potential transition zone downstream of the inlet depending how
inflow is implemented (Ekrann and Aasen 2000). The thickness of this zone affects the results
and might even cause a rate dependency of ensemble parameters for finite size samples (cf.,

13
Computing Relative Permeability and Capillary Pressure of… 1733

Benham et al. 2021; Boon and Benson 2021). Impervious sample jackets acting as no-flow
boundaries will channel the flow, complicating the determination of anisotropy, or directional
variation of the flow behaviour (e.g. Durlofsky 2005).
In their efforts to circumvent such boundary artefacts, Dale et al. (1997) proposed an upscal-
ing method using periodic boundary conditions. This natural extension of the periodic flux
boundary conditions from single-phase upscaling (Durlofsky 1991) to multiphase flow was dem-
onstrated by Virnovsky et al. (2004). Matthai and Tran (2023) further explored the effects of
material interfaces in such numeric experiments as modelled with the method of Van Duijn et al.
(1995) and Van Duijn & De Neef (1998). The periodic boundaries are relevant to modelling
flow in flow cells as reservoir representation of geological beds (Mikes et al. 2006). However, a
remaining limitation of this method lies in the uniform boundary pressure it produces. In a realis-
tic setting, fractional flow and fluid pressure in the region of interest vary with the broader hetero-
geneity structure (S. Jonoud and Jackson 2008a, b), i.e. interdigitating geological features that are
cut short by the model boundaries. Furthermore, the angle-weighted averaging of flow functions
obtained in orthogonal flow directions (e.g. Eq. 18, Boon et al. 2022) introduces inconsistencies
by averaging upscaled relative permeability curves derived for flow parallel and perpendicular to
layers although these do not necessarily emerge from the same fine-scale saturation distribution.
Fully periodic pressure boundary conditions permit a more comprehensive representation of the
flow behaviour while eliminating the aforementioned boundary artefacts.
Similar to the flux-periodic method of Matthai and Tran (2023), the new upscaling proce-
dure with fully periodic boundary conditions presented here involves two steps: 1) determina-
tion of stationary non-uniform saturation distributions, and 2) computation of upscaled flow
functions via homogenisation.

2.2 Determination of Saturation Distributions

Considering the heterogenous domain Ω (see Fig. 1), we determine the spatial steady-state
saturation distribution, s(x ∈ Ω), herein referred as saturation state, by solving the incom-
pressible, gravity-free two-phase flow problem,
𝜕s𝛼
−∇.𝐯𝛼 = 𝜙 in (t > 0) ∩ Ω (1)
𝜕t
where v is the Darcy velocity, 𝜙 and t are porosity and time, respectively. The subscript α
refers to the fluid phase either ­CO2 or brine. The domain boundaries are subjected to Dir-
ichlet pressure conditions (Durlofsky 1991; Wen et al. 2003),

Fig. 1  Continuous rotation of a pressure field (Eq. 4) around a symmetric unit cell of a heterogene-
ous domain (Ω) representing a cross-bedded sandstone as is discussed in more detail further below. Col-
oured arrows show periodicity of the variable flux across the domain boundaries 𝜕 Ω (Eq. 3) which is not
necessarily uniform across boundaries

13
1734 A. A. Youssef et al.

p = p𝐱̂ 𝐨 + 𝐆.(𝐱̂ − 𝐱̂ 𝐨 ) (2)

where p is pressure at point ̂

𝐱 ∈ 𝜕Ω, while p̂𝐱0 is the pressure at a reference point ̂
𝐱0 ∈ 𝜕Ω
and G is the pressure gradient vector. The boundary conditions are completed by assigning
full flux periodicity,

(𝐧.𝐯𝛼 )||𝜕Ω𝜂,1 = −(𝐧.𝐯𝛼 )||𝜕Ω𝜂,2 (3)

where n is a unit normal vector perpendicular to 𝜕Ω and pointing outward, and 𝜕Ω𝜂,1 and
𝜕Ω𝜂,2 are two opposite boundaries (Fig. 1). By contrast to single-phase flow, the selection
of the boundary condition 2 controls s such that s = s(x,G). In other words, the saturation
state relies on both, direction and the magnitude of the pressure field. Now, we define G for
a two-dimensional flow as
( )
𝐆 = g cos(𝜃)𝐢x + sin(𝜃)𝐢y ∀ 𝜃 ∈ [0, 2𝜋) (4)

and for a three-dimensional flow,

( )
𝐆 = g sin(𝜑) cos(𝜃)𝐢x + sin(𝜑) sin(𝜃)𝐢y + cos(𝜑)𝐢z ∀ 𝜃 ∈ [0, 2𝜋), 𝜑 ∈ [0, 𝜋) (5)

where g the pressure gradient magnitude, i is a unit vector, and θ and φ are polar and azi-
muthal angles, respectively, defined in spherical coordinates. Figure 1 illustrates how these
boundary conditions provoke fluxes for different flow angles in 2D. Crucially, this defini-
tion enables the flux to vary freely across boundaries.
Dale et al. (1997) discusses the non-uniqueness of the steady-state solution that arises
when capillary jumps exist at the layer interfaces. Since their focus was on the boundary
value problem, they did not elaborate on the initial saturation state of the system. Here
we prescribe a unique initial saturation distribution as additional initial-value constraint
for each calculation. For example, the capillary-limit saturation distribution as computed
with the method of Ekrann et al. (1996). Importantly, it will be shown that this does not
uniquely determine the final saturation distribution (Fig. 2).

2.3 Estimation of Ensemble Properties

Integrating the saturation distribution over the flow domain, we recover the average satu-
ration, s , verifying that it matches the prescribed value. This takes into account the pore
volume

∫ s d(volp )
Ω
s= (6)
∫ d(volp )
Ω

The saturation dependent ensemble capillary pressure, pc , and relative permeability, 𝐤r ,

are no longer single-argument functions, but depend on flow direction and local force bal-
ances (Boon et al. 2022). Both pc and 𝐤r are field variables, with pc being a scalar and
𝐤r being a vector. In two dimensions, they are functions of g and θ. In three dimensions,
they also depend on φ. These functions are defined over a continuous spectrum of satura-
tion states between the CL and the VL. To map out this range, the workflow in Fig. 2 is
followed.

13
Computing Relative Permeability and Capillary Pressure of… 1735

Fig. 2  Proposed ensemble upscaling workflow. Δg is the incremental increase of the pressure gradient mag-
nitude while incremental changes in polar and azimuth angles are Δ𝜃 = 2π
N
and Δ𝜑 = M 𝜋
, respectively

13
1736 A. A. Youssef et al.

For each saturation state, 𝐤r is estimated from the tensor operation

𝐤r𝛼 = 𝐤−1 𝐤𝛼 (7)

where both 𝐤𝛼 and 𝐤 are estimated from the solution of a single-phase upscaling prob-
lem with periodic boundary conditions. This involves conducting at least two separate flow
simulations in 2D and three in 3D. Typically, but not necessarily, these simulations are per-
formed in orthogonal directions. For each run, the volumetrically averaged pressure gradi-
ent, ⟨𝜕p⟩, and the volumetrically averaged velocity, ⟨𝜕𝐕⟩ are determined on Ω. The tensors
appearing in (7) are obtained by solving the overdetermined system (Durlofsky 1991; Lang
et al. 2014),

⎡ kxx ⎤
⎢ kxy ⎥
⎢ ⎥
kxz
⎡ ⟨𝜕𝐩⟩1 ⟨𝜕𝐩⟩2 ⟨𝜕𝐩⟩3 ⎤⎢ k ⎥ ⎡ ⟨𝐯⟩ ⎤
1 1 1 1

⎢ ⟨𝜕𝐩⟩1 ⟨𝜕𝐩⟩2 ⟨𝜕𝐩⟩3 ⎥⎢

2 2 2 ⎢ yx ⎥ ⎢ ⟨𝐯⟩2 ⎥
⎢ ⎥ (8)
3 ⎥ kyy = −⎢ 3⎥
⎢ ⟨𝜕𝐩⟩1 ⟨𝜕𝐩⟩2 ⟨𝜕𝐩⟩3 ⎥⎢⎢ kyz ⎥⎥ ⎢ ⟨𝐯⟩ ⎥
3 3

⎣ 𝐈1 𝐈2 𝐈2 ⎦ ⎣ 𝟎 ⎦
⎢ kzx ⎥
⎢k ⎥
⎢ zy ⎥
⎣ kzz ⎦

⎡ 𝜕x pi 𝜕y pi 𝜕z pi ⎤ ⎡ 0 0 0 ⎤ ⎡ 0 0 0 ⎤
where ⟨𝜕𝐩⟩i1 =⎢ 0 0 0 ⎥, ⟨𝜕𝐩⟩i2 = ⎢ 𝜕x pi 𝜕y pi 𝜕z pi ⎥, ⟨𝜕𝐩⟩i3 = ⎢ 0 0 0 ⎥,
⎢ ⎥ ⎢ ⎥ ⎢ i i i⎥
⎣ 0 0 0 ⎦ ⎣ 0 0 0 ⎦ ⎣ 𝜕x p 𝜕y p 𝜕z p ⎦
⎡ ⟨vx ⟩i ⎤
⎢ � �i ⎥ [ ]T
⟨𝐯⟩i = ⎢ vy ⎥, 𝟎 = 0 0 0 . The superscript reflects the simulation number, and the sub-
⎢ � v �i ⎥
⎣ z ⎦
script is the index of the row containing non-zero elements in the pressure gradient matri-
ces. This overdetermined system is solved using least square fitting. pc is evaluated as the
pore volume weighted capillary pressure which is calculated from

∫ pc d(volp )
Ω
pc = (9)
∫ d(volp )
Ω

We thus estimate pc as the pore-volume-weighted value rather than using the differ-
ence between the average of the non-wetting phase pressure and the average wetting phase
pressure as proposed by Kueper and McWhorter (1992). This aligns with the microscopic
view of the porous medium where pc represents pressure differences across phase inter-
faces (e.g. Adler and Brenner 1988).

2.4 Implementation of Periodic Pressure Boundary Conditions in Collocated Finite

Element‑finite Volume Framework

We implement the periodic pressure boundary conditions in the collocated finite element-
finite volume scheme with embedded discontinuities (Nick & Matthai, 2011a, b; Tran et al.

13
Computing Relative Permeability and Capillary Pressure of… 1737

2020) so that we can capture pressure and saturation discontinuities arising from capillary
heterogeneity within the rock volume of interest. The space–time solution method consid-
ered here is operator splitting, solving the pressure and saturation equations sequentially
using an IMPES (implicit for pressure and explicit for saturation) approach. The pressure
equation that is treated as elliptic at the small sample scale, is solved with the finite element
(FE), phase velocities are computed, and then the hyperbolic saturation transport equation
is solved with the finite volume (FV) method, see Geiger et al. (2004) and Paluszny et al.
(2007). Various iterative solution schemes are available to deal with the non-linearity of
this system (Aziz and Settari 1979).
The resulting system of linear algebraic equations for pressure needs to be modified to
incorporate the periodic boundary conditions. Since the pressure equation is solved with
the FEM, no-flow boundary conditions are “natural”. In such case, the pressure system
takes the form
∑
ai,j pj = Rj ∀i ∈ {1, 2, ..., N}
(10)
j

and no modification is required. By contrast, the concept of periodicity implies domain

repetition. Consequently, adjustments are required for the equations containing shared
nodes located on 𝜕Ω. For instance, for i = 1, the nodal contributions of finite elements
associated with the mirrored part of the domain must be reintroduced (Fig. 3). These coef-
ficients appear in the equation of the corresponding periodic node i = 10. Specifically,
a10,10 is added to a1,1, a10,12 is added to a1,2 and a10,15 is added to a1,5. A similar treatment is
applied to the right-hand side of the algebraic system. However, this superposition forms
only one part of the boundary treatment. Thus, Fig. 3 shows that the mirrored boundary
part contains a node corresponding to node 16, which shares certain parent elements with
node 1. To recover this element’s contribution, we need to add the term involving p16 in
i = 10 to i = 1, adjusting pressure as indicated by boundary condition (2). Finally, for i = 1,
Eq. (10) transforms into

Fig. 3  Triangular hybrid finite element (FE)-finite volume (FV) mesh with material properties that are
piecewise constant on the finite elements. The triangular FE are delimited by solid black lines connecting
their nodes (black filled circles). The polygonal FV are delimited by dashed lines. They span multiple FE’s
dividing them into FV sectors. Their part missing from the left periodic inflow boundary is shown in faint
colour. A unique numbering relates the node IDs on the inflow with those on the outflow boundary

13
1738 A. A. Youssef et al.

(a1,1 + a10,10 )p1 + (a1,2 + a10,12 )p2 + a1,3 p3 + a1,4 p4 + (a1,5 + a10,15 )p5
( ) (11)
+ a10,16 p16 + 𝐆.(𝐱̂ 1 − 𝐱̂ 10 ) = R1 + R10

This approach is readily generalized to three dimensions. It is crucial to emphasize that

the coefficients for interior nodes 3 and 4 remain unchanged throughout this modification.
Now returning to the saturation equation discretized with the FVM and explicitly in
time
Δt ∑
sm+1 = sm
𝛼 + (𝐯𝛼 .𝐧)j ∀i ∈ {1, 2, ..., N} (12)
𝛼i i Ai 𝜙i j

where the superscripts denote the time level, and the summation on the right-hand side
accounts for the flux accumulation from all FV sectors (Fig. 3), to accommodate periodic-
ity, pairs of corresponding node-centred finite volumes situated on opposite sides of the
periodic boundary 𝜕Ω are treated as one. The net flux through such recombined FVs is
fully defined permitting the computation of saturation at the upstream node as if it were
on the inside of the model domain. Corner-node finite volumes require special considera-
tion as they are recombined from three or more incomplete finite volumes. Generally, the
number of corresponding unified periodic boundary nodes, denoted as NF, shared between
F boundaries is determined using the formula

∑
i=F−1
NF = F + Ni (13)
i=1

The first term in this series is 1.

All computations presented in this paper were performed with OpenCSMP, an open-
source C +  + 17 library that implements the hybrid finite element – finite volume frame-
work discussed, including the use of split boundaries to model the dynamic pressure—sat-
uration (embedded) discontinuities with the extra degrees of freedom required to model the
potential jump discontinuities in pressure and saturation that will arise in the simulations
for certain parameter combinations. The results presented here adhere the split boundaries
criteria except at the model boundaries due to the difficulties brought by periodicity. We
hope to solve this in the future versions.

3 Demonstration Models

We present two demonstration models that represent REVs of cross-bedded sandstone. The
first is a three-dimensional sedimentological (bedform) model (Fig. 4). To demonstrate the
effects of its layered heterogeneity on the saturation pattern, it is sliced parallel to-, and
perpendicular to the sediment transport direction. The second model (Fig. 5) is a much-
simplified layered idealization of a slice parallel to the sediment transport direction, ame-
nable to analytical solution. Presenting simulations done with these models, we showcase
our approach and workflow, also aiming to improve the general understanding of direc-
tional and rate-dependent multiphase flow in heterogeneous porous media.

13
Computing Relative Permeability and Capillary Pressure of… 1739

Fig. 4  Cross-bedded sandstone model. a 3D model; b Section model PERPENDICULAR to the sediment
transport direction; c) Section model PARALLEL to the sediment transport direction. Reduced permeabil-
ity beds are in black and permeable beds in grey

Fig. 5  OBLIQUE layer model

consisting of two layers inclined
at 45°. The model is meshed
with a double periodic triangular
mesh. The mesh consists of 2721
nodes and 5224 elements

13
1740 A. A. Youssef et al.

Table 1  Properties of the Property Bed 1 Bed 2

constitutive rocktypes of the
cross-bedded sandstone
Porosity 0.28 0.286
Permeability (mD) 100 1000
Brooks-Corey parameter 1.5 2
Residual water saturation 0.27 0.105
Entry pressure (kPa) 3.162 1.0

Table 2  Mesh parameters of Model PERPENDICULAR PARALLEL

cross-bedded sandstone sections
Number of nodes 46,217 149,900
Number of elements 67,718 230,080

3.1 Cross‑bedded Sandstone

Sandstones are coarse-grained clastic sediments deposited in high-energy environments

such as the foreshore, estuaries, or river channels. These different environments of dep-
osition manifest in characteristic types of cross stratification. Cross-bedding is often
marked by rhythmic alternations in grain size inversely correlated with permeability and
capillary entry pressure (e.g. Boggs 2012; Trevisan et al. 2017). Our model of cross-
bedded sandstone (Fig. 5) was created using the algorithm of Rubin (1987), expressing
sedimentary structures using the superposition of waves. Three bedforms are propagat-
ing in the z-direction (sediment transport direction), controlling the interference pat-
terns created (cf., Rubin and Carter 2005). The bedforms are truncated in the y-direction
(deposition direction). Here the bedform with the highest level, y, at each (x,z) point is
selected. Perpendicular to the bedform propagation direction (x-direction), the shapes of
the bedforms are due to the interference of two propagating waves. Wave characteristics
and sedimentation rate are set such that a fully periodic REV with the dimensions of
30 × 6.4 × 80 cm results with reference to the x, y and z axes, respectively. The volumes
between the surfaces are populated with two alternating rock types with distinct prop-
erties, see Table 1. For flow simulation, fluid properties are chosen to model viscosity
and density contrasts representative of water and ­CO2 at a subsurface depth of 1.5 km
for a normal geothermal gradient. It follows that the supercritical C ­ O2 is 22 times less
viscous than water.
Since the three-dimensional structure of these bedforms is complex (Youssef et al.
2022), we simulate and illustrate the flow along cross sections. Two perpendicular cross
sections are considered (Fig. 4): one perpendicular to the sediment transport direction
(Model PERPENDICULAR), and the other parallel to the sediment transport direction
(Model PARALLEL). Both sections are meshed with an unstructured triangular mesh.
Table 2 summarizes the mesh parameters of each section.

13
Computing Relative Permeability and Capillary Pressure of… 1741

Table 3  Rock properties of Property Fine sand Coarse sand

OBLIQUE layer model
Porosity 0.28 0.286
Permeability (mD) 365 2534
Brooks–Corey parameter 1.8 2.2
Residual water saturation 0.159 0.104
Entry pressure (kPa) 1.0 0.75

3.2 Oblique Layer Model

Square model OBLIQUE has a side length of 10 cm and consists of a repetitive sequence
of regularized fine and coarse-grained sandstone layers (Fig. 5, Table 3). The layers are
inclined by 45°; each is ~ 1.77 cm thick. The finite-element mesh is fully periodic, repre-
senting an REV with layers matching at opposing boundaries. The model has an analytical
solution described in the Results section. Its purpose is to illustrate force balance variation
with flow direction. The model idealizes the cross-beds in Model PARALLEL and is simu-
lated with the same fluid properties.

4 Results

Model OBLIQUE is ideally suited to illustrate the benefits of the periodic pressure bound-
ary conditions and the difference they make to the saturation distribution as compared with
the constant fractional flow setup.

4.1 Effective Properties of Model OBLIQUE

Due to its symmetry, it is sufficient to investigate flow directions for a single quadrant with
𝜃 ∈ [ −𝜋 , 𝜋 ]. Remaining angles can be extrapolated from this range. The effective intrinsic
4 4
permeability of the model is anisotropic with a principal axis aligned with its layers, and a
value equivalent to their arithmetic average of 1450 mD. The minimum permeability axis
is oriented perpendicular to the layers with a value equivalent to their harmonic mean of
639 mD. Instead of conducting multiple numerical simulations for different flow directions
and pressure gradients, we can construct the ensemble saturation functions analytically, see
Appendix. The effective 𝐤r and pc (Fig. 6), illustrate the flow direction dependence. They
are also rate dependent and this feature is most pronounced for flow perpendicular to the
layers (𝜃 = −𝜋 4
). For the case of 𝜃 = 𝜋4 , capillary forces dominate over the viscous ones,
and the CL persists across the entire range of pressure gradients applied. Importantly, this
implies that the rate dependence estimated for the case of flow along the layers (e.g. Boon
and Benson 2021) merely is an artifact of the boundary conditions (cf., Ekrann and Aasen
2000). On the other hand, for −𝜋 4
< 𝜃 < 𝜋4 , the balance between viscous and capillary forces
varies, with a wider range observed when θ is close to −𝜋 4
and a narrower range when θ is
close to 𝜋4 .
Another interesting feature is a variation of the viscous-capillary force balance with sat-
uration. While the rate-dependent velocity (green curve, Fig. 6d) maintains the same ratio

13
1742 A. A. Youssef et al.

Fig. 6  Model OBLIQUE: rate dependence of 𝐤r across layers (a), parallel to layers (b), and pc (c). Relative
permeabilities of the constituent rocktypes are shown by dashed black lines. Curve colours in a,b,c cor-
respond to one the force balance represented in graph (d indicating variation of viscous-capillary force bal-
ance between CL and VL as represented by flow velocity variation with flow direction.. Red curves reflect
VL behaviour and blue ones the CL. The dashed line represents a fixed velocity of (3.5 × ­10−6m/s, which is
a typical front spreading velocity of 1 ft/day as observed during water flooding of oil fields

between the CL and VL limits, Fig. 6 reveals that 𝐤r is close to the CL for average satura-
tions near the residual water saturation before shifting towards VL at high water satura-
tions. This phenomenon can be explained considering capillary spreading near the layer
interfaces: at a low water saturation, viscous forces encounter a significant resistance trying
to overcome capillary spreading upstream of the interface. This resistance increases with
| dp |
the capillary gradient (∝ | ds c |) when saturation builds up ahead of the interface (Fig. 7). At
| w|
high water saturations, the strength of the capillary gradient diminishes, and viscous forces
dominate (Fig. 7). This behaviour matches with Boon and Benson’s (2021) observations.
In this scenario, multiple saturation states exhibit a similar pc as shown in Fig. 6c.
This is due to the model characteristics: each layer has an equivalent pore volume,
resulting in an equal contribution to the pc curve, see Eq. (9). In addition, the high vis-
cosity contrast between C ­ O2 and water promotes similar fractional flow curves. Thus,
both layers have similar saturations at the VL, thereby narrowing the window of rate-
dependent behaviour of the pc curve. Crucially, this behaviour is due to the unique

13
Computing Relative Permeability and Capillary Pressure of… 1743

Fig. 7  Model OBLIQUE: steady-

state saturation distributions
along the diagonal perpendicular
to the layers at the viscous-cap-
illary force balance represented
by the green line in Fig. 6. The
interface between the two layers
is represented by a vertical black
dashed line in the centre

fluid and material properties of this model and should therefore not be interpreted as
rate-independent.
The ratio between viscous and capillary forces can be estimated in terms of capil-
lary number Nc = 1Dkp CO2 , where l is period length, µ is viscosity, k is the effective per-
v l𝜇
c
meability and pc is the capillary pressure evaluated at certain average saturation. This
definition considers the variation from CL to VL by only accounting for the velocity
perpendicular to interfaces. Based on the model parameters, the rate-dependent win-
dow lies in the range of Nc = 10−7 ∶ 10−3 where the lower limit represents CL, and the
upper limit corresponds to VL. The application of this range shall be restricted to
media possessing the same heterogeneity structure, permeability contrast, viscosity
contrast, relative correlation length and aspect ratio.

4.2 Effective Properties of Three‑dimensional Flow in Model OBLIQUE

Extruding Model OBLIQUE along the z-axis perpendicular to the xy-plane, the princi-
pal axes of the absolute permeability tensor remain the same, but there is the additional
third component aligned with the z-axis. This implies a semi-isotropic absolute perme-
ability tensor, where the permeability along the z-axis is equivalent to the arithmetic
average of the layers. Similarly, kr is semi-anisotropic. However, |vT |3D that would be
required to achieve the same saturation state as in two dimensions is larger. For the
same saturation distribution, |vT |3D can be estimated by scaling |vT |2D by a factor of
csc(φ). When flow is exactly aligned with the z-axis (𝜑 = 0 ), the flux across the layers
diminishes, and we arrive at the CL. By contrast, when 𝜑 = 𝜋2 , the 3D case collapses to
the 2D one and the results in Fig. 6. This conclusion matches the observations reported
by Ringrose et al. (1993), that flow through domains extruded to 3D differs from flow
through corresponding cross sections.

13
1744 A. A. Youssef et al.

Fig. 8  Steady-state water saturation obtained with the Penn-State method (left), water saturation profile
along the diagonal perpendicular to the interfaces (white line), plotted from top-left to bottom-right corner
using Penn-State method (centre), and periodic boundaries (right)

4.3 Influence of Boundary Conditions

To investigate the influence of boundary conditions, the saturation distribution in model

OBLIQUE is re-evaluated with the boundary conditions of the Penn-State method. In
this setup, co-injection occurs from a highly permeable region located adjacent to the
left boundary, while the right boundary is maintained at a fixed pressure. The remaining
boundaries are sealed. Figure 8 shows the saturation distribution observed after injecting
more than 45 pore volumes (PV) and a water cut of 0.8 at the same injection velocity of
2.13 × 10−8 m/s. These numeric results reveal a completely different asymmetric saturation
pattern as compared with the repetitive symmetric saturation cycles resulting from periodic
boundaries. This comparison corroborates earlier studies emphasizing the strong influence
of boundary conditions on saturation distributions (cf., Ekrann and Aasen 2000).

4.4 Upscaled Properties of Cross‑Bedded Sandstone Models PERPENDICULAR

and PARALLEL

Due to the geometric complexity of these models, there is no analytical solution for their
ensemble saturation functions. Instead, methodology (Fig. 2) is applied to extract extended
saturation functions from this model. For their analysis, it is practical also to define vis-
cous-capillary force balance in terms of G rather than 𝐯T . The periodic boundary setup
allows total flux to vary freely along boundaries as determined by the local rock properties,
structure, and the applied fluid-pressure gradient. In the multiphase flow analysis, simula-
tions for two fixed far-field pressure gradients were analysed, applying gradients of 500
and 14,000 Pa/m in three different flow directions: 𝜃 = 0, 𝜃 = 𝜋4 and 𝜃 = 𝜋2 . The analyti-
cal CL saturation distributions were used to initialize the models. Then, the fluid mixture

Table 4  Eigen values and corresponding eigen vectors of intrinsic permeability measured for models PER-
PENDICULAR and PARALLEL
Model PERPENDICULAR Model PARALLEL

kmax kmin kmax kmin

o o o
463 mD − 16.69 224.6 mD 73.31 536.7 mD 11.35 186.4 mD − 78.65o

13
Computing Relative Permeability and Capillary Pressure of… 1745

Fig. 9  Model PERPENDICULAR: steady-state saturation distribution in cross-bedded sandstone computed

for sw = 0.7. The white lines represent equally spaced pressure contours. Each row of images corresponds
to the flow direction indicated by the black arrow on the left. The left column represents the elevated pres-
sure gradient of 14 kPa/m and the right-column 0.5 kPa/m: Row (1) horizontal flow from left to right, row
(2) oblique upward flow with θ = 45o, row (3) vertical upward flow from bottom to top

Fig. 10  Model PARALLEL: steady-state saturation distribution computed for sw = 0.7 with equally spaced
pressure contours in white. Flow direction is indicated by black arrow on the left; row (1) horizontal flow
from left to right, row (2) oblique flow with θ = 45o, row (3) vertical flow from bottom to top. Top three
plots represent steep pressure gradient (14 kPa/m), bottom ones weak gradient (0.5 kPa/m)

13
1746 A. A. Youssef et al.

was cycled through these models until a pseudo steady-state saturation distribution was
obtained for the given force-balance.
For both models PARALLEL and PERPENDICULAR, single-phase flow analysis indi-
cates that the principal axes of their intrinsic permeability tensors are oriented oblique to
the coordinate axes (Table 4). Figures 9 and 10 show stationary saturation distributions
obtained for the resulting six scenarios computed at sw = 0.7. The obtained saturation dis-
tributions, highlight the importance of the macroscopic heterogeneity for residual trapping.
The saturation distribution is highly non-uniform as C­ O2 collects upstream of the lamina-
tions and gets trapped where fine-grained layers converge on unconformity surfaces. Heter-
ogeneity-related trapping is most pronounced at intermediate to high C ­ O2 saturations, but
leads to high local saturations already at small average saturations. As expected, trapping
is inversely correlated with viscous forces breaking down capillary barriers. Trapping also
is more pronounced when the far-field fluid pressure gradient is oriented perpendicular to
the sediment transport direction by comparison with when it is aligned with the sediment
transport direction. This observation will have to be re-evaluated once we will be able to
mesh and run the REV-scale model in three dimensions.
Changing G can be seen to have a significant influence on the saturation pattern
obtained for model PERPENDICULAR (Fig. 9). Its effect is less pronounced on model
PARALLEL (Fig. 10). By contrast with the Penn-State setup, the saturation distributions
do not vary systematically between inlet and outlet, indicating that this variation results
from the different boundary conditions. Note also, that the non-uniformity of saturation
within each rock-type indicates that the typical CL and VL assumptions cannot be justified
except for extreme values of 𝐆 not encountered in the real world.
Figure 11 shows sample of eigen value of 𝐤r . As indicated previously, the change in
g has a marked influence on 𝐤r of model PERPENDICULAR as compared with model
PARALLEL. For the sake of an effective visualization of 𝐤r , and including the effect of
the off-diagonal tensor components, we present the eigen vectors of 𝐤𝛼 , scaled by the eigen
values determined from 𝐤 (Table 4). Figure 12 shows 𝐤r evaluated for the horizontal flow
case in both slices. Similar to model OBLIQUE, the orientation of 𝐤𝛼 of model PARAL-
LEL hardly deviates from the orientation of 𝐤 . This indicates that the flow is controlled by
geometry and the non-wetting phase is filtrated homogenously by the coarse–fine interfaces

Fig. 11  Eigen value-based 𝐤r curves of model PERPENDICULAR (left) and model PARALLEL (right) for
the 2 horizontal pressure gradients

13
Computing Relative Permeability and Capillary Pressure of… 1747

Fig. 12  𝐤r of model PERPENDICULAR (left) and model PARALLEL (right) for the 2 horizontal pressure
gradients. The radial extent corresponds to the scaled maximum and minimum of phase permeability ten-
sors, while the angle records the orientation of the eigen vectors. The curves are parameterized by water
saturation. direction

( )
Fig. 13  Model PERPENDICULAR: ­CO2 saturation distribution for 𝐆 = 14 kPa∕m, π . The grayscale saturation
4
map is overlain by flow vectors scaled by pressure gradient magnitude

Fig. 14  Model PERPENDICULAR: C ­ O2 saturation distribution for g = 500 Pa/m along y-direction from
bottom to top. The saturation colour map is masked with flow field (vectors) scaled to pressure gradient
magnitude. The dashed yellow ellipses show selected locations where the pressure gradient is weak

13
1748 A. A. Youssef et al.

for both pressure gradients considered. In contrast, the orientation of 𝐤r of model PER-
PENDICULAR varies with saturation and 𝐆. Moreover, the orientation of ­CO2 and water
relative permeability tensors do not necessary align for the same saturation. This confirms
that the capillary interfaces in these models exert a strong influence on two-phase flow.
Figure 13 shows that these interfaces may be characterized by infinite capillary pressure
gradients (∇pCO2 → ∞) as pressure discontinuities arise when they act as membrane seals
where ­CO2 collects until the entry pressure of the interface is exceeded, there are jump
discontinuities in fluid pressure (the gradient becomes infinity) that are supported by the
special discretisation of the interfaces with split boundaries.
To gain a deeper understanding of the saturation patterns, ( we analyse
) the saturation pat-
tern obtained from model PERPENDICULAR using 𝐆 = 500 Pa ,
m 2
𝜋
and the parameters
in Table 3. Figure 14 shows C ­ O2 saturation distribution and pressure gradient vectors.
While the overall gradient is small, locally, pressure varies strongly (long vectors). Areas
with flat pressure gradients are marked by ellipses. They coincide with high ­CO2 satura-
­ O2 accu-
tions and capillary sealing of the layer interfaces downstream. In these settings, C
mulates in the high- permeability laminae until the flow stagnates.

5 Discussion

Using two-dimensional models of cross-bedded sandstone, this study demonstrates a novel

steady-state upscaling methodology including dynamic capillary sealing and pressure
compartmentalisation as well as doubly periodic boundary conditions for fluid pressure
and saturation, respectively. Full periodicity (Eq. 3) was already proposed as essential by
Virnovsky et al. (2004), extending the one-dimensional analysis of Dale et al. (1997). Here
we use a more elaborate simulation framework (e.g. Matthai and Tran 2023; Shao et al.
2022; Tran et al. 2020), and also consider continuously varying flow directions. Among
other features, this model correctly represents discontinuities in both saturation and pres-
sure, going beyond previous studies (e.g. Lohne et al. 2006; Pickup et al. 2000; Ringrose
et al. 1993). Our motivation for considering continuously varying directions is that flow
alignment with layers and coordinate axes (e.g. Boon et al. 2022; Jackson et al. 2018;
Jonoud and Jackson 2008a, b) is only one end-member scenario within a spectrum because
the local force balance depends on the complex interplay between the spatial arrangement
of layers and their internal property variations. Although computationally more costly, our
“full physics” approach allows for an analysis of flow-rate effects that is not possible with
invasion percolation-based models that cannot analyse the time dependence of the process
(cf., Cavanagh and Haszeldine 2014; Trevisan et al. 2017). While this analysis only con-
siders capillary and viscous forces (Eq. 1) which may be justifiable at the sub-meter scale
(e.g. Pickup and Sorbie 1994; Virnovsky et al. 2004), Matthai and Tran (2023) showed
that the underpinning simulation approach also handles buoyancy-driven flow and gravita-
tional phase segregation. In addition, the simple mathematical nature and robustness make
the steady-state upscaling more attractive than unsteady state especially at the small length
scale where numerical dispersion is insignificant (Pickup and Stephen 2000). However, at
larger scales, a combination of both steady state and unsteady state seems advisable as it
permits modeling of coarsening instabilities (Youssef et al. 2023).
In the numeric experiments shown, the global viscous-capillary force balance is defined
in terms of the pressure gradient as illustrated in Fig. 2. However, we can use flux or the

13
Computing Relative Permeability and Capillary Pressure of… 1749

capillary number by iterating the pressure gradient to match required flow velocity. One
concern related to capillary number is that it has many definitions (e.g. Berg and Ott 2012;
Boon and Benson 2021; Dale et al. 1997; Lohne et al. 2006). It has an intrinsic length
scale, and it is questionable whether it can be readily transferred to complex nested flow
geometries (Jonoud and Jackson 2008a, b). Thus, how should one define the drainage-cap-
illary pressure gradient in the case where there is a series of jump discontinuities defined
by entry pressure in actuality?
The saturation patterns observed in Cross-bedded sandstone or OBLIQUE are largely
due to capillary heterogeneity as opposed to permeability heterogeneity. Therefore, from a
hydrodynamic point of view, building realistic models involving different rocktypes criti-
cally depends on the physical realism of the model representation of capillary interfaces.
The input parameters for these building blocks can be estimated either from core flooding
on homogenous samples or pore scale simulation (e.g. Øren et al. 2019).
Our analysis indicates that the final saturation state is unique and independent of the ini-
tial saturation distribution. Is this conclusion valid only when input imbibition and drainage
saturation curves for the constituent rock types coincide?—This remains an open question
that needs to be addressed when considering the hysteresis of rate- and direction-dependent
saturation functions, to be investigated perhaps via physical experiments. In addition, the
approach proposed states that a form of steady state is achieved between capillary and vis-
cous forces. Achieving the steady state means that there are no further changes in both
pressure and saturation. It is expected that phase pressure reaches the steady-state faster
than saturation, hence, we propose a steady-state criterion as
Vp 1 𝜕s
max( ∑ )≤𝜀 (14)
Vp s 𝜕t

where Vp is the pore volume associated with each control volume in the model and 𝜀 is a
threshold value that can be determined by trial and error. This definition gives more weight
to larger FVs compared to smaller ones.
The concept described here can be easily modified to more general flow scenarios via
the border region approach (Wen et al. 2003). The border region yields equivalent parame-
ters for models that are lacking separation of scales and the determined parameters collapse
to effective if separation of scales exists. In th former case, the resultant equivalent tensor
system 8 is not necessarily symmetric nor positive definite and tuning may be required
(Durlofsky 2005). The border region approach simplifies the upscaling process as com-
pared to the problematic alignment of grids with the flow field (cf., Wen et al. 2000).

6 Conclusions

We present a steady-state upscaling study of multiphase fluid flow through heterogene-

ous composite rock types, forming geometrically complex facies associations, such as
cross-bedding or foresets as are typical for fluvial sandstones. To overcome the issue of
unrealistically uniform boundary fractional flows and saturation conditions, complicating
the analysis of such rocks in physical experiments, we employ doubly periodic boundary
conditions for both, saturation and pressure. This setup allows us to obtain non-uniform
saturation distributions free from boundary artefacts and suitable for the analysis of flow
at arbitrary far-field pressure gradients and force balances, including cases where the flow
is locally or globally blocked by dynamic capillary seals formed at coarse–fine interfaces

13
1750 A. A. Youssef et al.

due to saturation build-up upstream and choking of wetting phase flow as well resulting
in pressure and saturation jump discontinuities. The demonstrated ability to model these
blockages allows us to study macroscopic heterogeneity-induced capillary trapping as an
additional contribution to the amount of immobilized ­CO2 inside of the injection plume.
Our results further corroborate the finding that the directionality and rate dependence of
relative permeability is not an artifact of the boundary conditions but a significant and
reproducible feature of heterogeneous porous media. Thus, in the investigated cross-bed-
ded sandstone model, capillary number varies significantly with flow direction, ranging
over four orders of magnitude ranging from 10−7 to 10−3, dependent on whether the flow
is aligned with or perpendicular to the bedding planes, respectively. Our new simulation/
analysis method forms part of an already partially automated albeit computationally expen-
sive workflow for the analysis of scanned and segmented core samples (applicable also
to outcrop analogues). By contrast with invasion percolation approaches, this full-physics
method permits a rigorous analysis of rate- and pressure gradient-direction dependent rela-
tive permeability / capillary pressure functions for the modelling and simulation of plume
­ O2 storage sites.
spreading in C

Appendix

Inspecting the 1D analytical solution of Dale et al. (1997), we make two important
observations:

1. with all other parameters constant, the flux across the layer interfaces plays a crucial
role in controlling the transition between CL and VL conditions;
2. the phenomenon of rate dependency detected in steady-state upscaling performed using
periodic boundaries is a direct consequence of the capillary spreading around interfaces.

For the Model OBLIQUE, we deduce the following corollaries:

1. the effective parameters can be evaluated using the projection of the total velocity on
the interfaces. In other words, the total flow velocity for any prescribed flow direction is
|𝐯 |
|𝐯T | = | T |1D 𝜋 𝜋
∀𝜃 ∈ [− , ] (A1)
| |2D cos(𝜃 + 𝜋 ) 4 4
4

where |𝐯T |1D is the total velocity magnitude obtained from the 1D analytical solution, and
|𝐯T |2D is the total flow velocity magnitude when the flow is oriented at angle 𝜃;
2. due to the homogeneous capillary spreading around interfaces, 𝐤r is aligned with the
principal axes of absolute permeability.

Supplementary Information The online version contains supplementary material available at https://​doi.​
org/​10.​1007/​s11242-​024-​02092-x.

Acknowledgements This study is a contribution to GeoCquest, a BHP-supported collaborative project

of The University of Melbourne (Australia), the University of Cambridge (UK), and Stanford University

13
Computing Relative Permeability and Capillary Pressure of… 1751

(USA), aimed at developing a better understanding of small-scale heterogeneity and its influence on ­CO2
flow and trapping mechanisms. The authors thank Mohamed Batekh for editing Python scripts for Rhino.

Funding Open Access funding enabled and organized by CAUL and its Member Institutions.

Declarations
Conflict of interest The authors have no competing interests to declare that are relevant to the content of this
article.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long
as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Com-
mons licence, and indicate if changes were made. The images or other third party material in this article
are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article’s Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly
from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References
Adler, P.M., Brenner, H.: Multiphase flow in porous media. Annu. Rev. Fluid Mech. 20, 35–59 (1988).
https://​doi.​org/​10.​1146/​annur​ev.​fl.​20.​010188.​000343
Ahmed, S., Bhattacharya, J.P., Garza, D.E., Li, Y.: Facies architecture and stratigraphic evolution of a river-
dominated delta front, turonian ferron sandstone, Utah, U.S.A. J. Sediment. Res. 84, 97–121 (2014).
https://​doi.​org/​10.​2110/​jsr.​2014.6
An, S., Wenck, N., Manoorkar, S., Berg, S., Taberner, C., Pini, R., Krevor, S.: Inverse modeling of core
flood experiments for predictive models of sandstone and carbonate rocks. Water Resour. Res. (2023).
https://​doi.​org/​10.​1029/​2023W​R0355​26
Aziz, K., Settari, A.: Petroleum reservoir simulation. Appl. Sci. (1979). https://​doi.​org/​10.​2118/​97816​13999​
646
Barker, J.W., Dupouy, P., 1996. An Analysis of Dynamic Pseudo Relative Permeability Methods, in:
ECMOR V - 5th European Conference on the Mathematics of Oil Recovery. European Association of
Geoscientists & Engineers. https://​doi.​org/​10.​3997/​2214-​4609.​20140​6870
Barker, J.W., Thibeau, S.: A Critical Review of the Use of Pseudorelative Permeabilities for Upscaling. SPE
Reserv. Eng. 12, 138–143 (1997). https://​doi.​org/​10.​2118/​35491-​PA
Benham, G.P., Bickle, M.J., Neufeld, J.A.: Upscaling multiphase viscous-To-capillary transitions in hetero-
geneous porous media. J. Fluid Mech. (2021). https://​doi.​org/​10.​1017/​jfm.​2020.​1134
Berg, S., Ott, H.: Stability of CO2-brine immiscible displacement. Int. J. Greenhouse Gas Control 11, 188–
203 (2012). https://​doi.​org/​10.​1016/j.​ijggc.​2012.​07.​001
Boggs, S., 2012. Principles of sedimentology and stratigraphy.
Boon, M., Benson, S.M.: A physics-based model to predict the impact of horizontal lamination on CO2
plume migration. Adv. Water Resour. (2021). https://​doi.​org/​10.​1016/j.​advwa​tres.​2021.​103881
Boon, M., Matthäi, S.K., Shao, Q., Youssef, A.A.A., Mishra, A., Benson, S.M.: Anisotropic rate-dependent
saturation functions for compositional simulation of sandstone composites. J. Pet. Sci. Eng. (2022).
https://​doi.​org/​10.​1016/j.​petrol.​2021.​109934
Breston, J.N., Hughes, R.V.: Relation between pressure and recovery in long core water floods. J. Petrol.
Technol. 1, 100–110 (1949). https://​doi.​org/​10.​2118/​949100-G
Cavanagh, A.J., Haszeldine, R.S.: The sleipner storage site: capillary flow modeling of a layered CO2 plume
requires fractured shale barriers within the Utsira formation. Int. J. Greenhouse Gas Control 21, 101–
112 (2014). https://​doi.​org/​10.​1016/j.​ijggc.​2013.​11.​017
Dale, M., Ekrann, S., Mykkeltveit, J., Virnovsky, G.: Effective relative permeabilities and capillary pressure
for one-dimensional heterogeneous media. Transp. Porous Media 26, 229–260 (1997)
Darman, N.H., Durlofsky, L.J., Sorbie, K.S., Pickup, G.E.: Upscaling immiscible gas displacements: quan-
titative use of fine-grid flow data in grid-coarsening schemes. SPE J. 6, 47–56 (2001). https://​doi.​org/​
10.​2118/​69674-​PA
Durlofsky, L.J.: Numerical calculation of equivalent grid block permeability tensors for heterogeneous
porous media. Water Resour. Res. 27, 699–708 (1991). https://​doi.​org/​10.​1029/​91WR0​0107

13
1752 A. A. Youssef et al.

Durlofsky, L.J.: Use of higher moments for the description of upscaled, process independent relative perme-
abilities. SPE J. 2, 474–484 (1997). https://​doi.​org/​10.​2118/​37987-​PA
Durlofsky, L.J., 2005. Upscaling and Gridding of Fine Scale Geological Models for Flow Simulation, in: 8th
International Forum on Reservoir Simulation. Stresa, Italy, pp. 1–59.
Ekrann, S., Aasen, J.O.: Steady-state upscaling. Transp. Porous Media 41, 245–262 (2000)
Ekrann, S., Dale, M., Langaas, K., Mykkeltveit, J., 1996. Capillary Limit Effective Two-Phase Properties
for 3D Media, in: European 3-D Reservoir Modelling Conference, Stavanger, Norway. SPE. https://​doi.​
org/​10.​2118/​35493-​MS
Fabuel-Perez, I., Hodgetts, D., Redfern, J.: Integration of digital outcrop models (DOMs) and high reso-
lution sedimentology – workflow and implications for geological modelling: Oukaimeden Sand-
stone Formation, High Atlas (Morocco). Pet. Geosci. 16, 133–154 (2010). https://​doi.​org/​10.​1144/​
1354-​079309-​820
Flemings, P.B., Stump, B.B., Finkbeiner, T., Zoback, M.: Flow focusing in overpressured sandstones: the-
ory, observations, and applications. Am. J. Sci. 302, 827–855 (2002). https://​doi.​org/​10.​2475/​ajs.​302.​
10.​827
Geffen, T.M., Member, J., Me, A.I., Owens, W.W., Parrish, D.R., Morse, R.A.: Experimental investigation
of factors affecting laboratory relative permeability measurements. J. Petrol. Technol. 3, 99–110 (1951)
Geiger, S., Roberts, S., Matthäi, S.K., Zoppou, C., Burri, A.: Combining finite element and finite volume
methods for efficient multiphase flow simulations in highly heterogeneous and structurally complex
geologic media. Geofluids 4, 284–299 (2004). https://​doi.​org/​10.​1111/j.​1468-​8123.​2004.​00093.x
Gross, E.C., Carr, M., Jobe, Z.R.: Three-dimensional bounding surface architecture and lateral facies hetero-
geneity of a wet aeolian system: Entrada Sandstone, Utah. Sedimentology 70, 145–178 (2023). https://​
doi.​org/​10.​1111/​sed.​13035
Huang, D.D., Honarpour, M.M.: Capillary end effects in coreflood calculations. J. Petrol. Sci. Eng. 19(1–2),
103–117 (1998)
Jackson, S.J., Agada, S., Reynolds, C.A., Krevor, S.: Characterizing drainage multiphase flow in heteroge-
neous sandstones. Water Resour. Res. 54, 3139–3161 (2018). https://​doi.​org/​10.​1029/​2017W​R0222​82
Jackson, S.J., Lin, Q., Krevor, S.: Representative elementary volumes, hysteresis, and heterogeneity in mul-
tiphase flow from the pore to continuum scale. Water Resour. Res. (2020). https://​doi.​org/​10.​1029/​
2019W​R0263​96
Jonoud, S., Jackson, M.D.: Validity of steady-state upscaling techniques. SPE Reservoir Eval. Eng. 11, 405–
416 (2008a). https://​doi.​org/​10.​2118/​100293-​PA
Jonoud, S., Jackson, M.D.: New criteria for the validity of steady-state upscaling. Transp. Porous Media 71,
53–73 (2008b). https://​doi.​org/​10.​1007/​s11242-​007-​9111-x
Krampert, E.W., 1934. Geological Characteristics of Producing Oil and Gas Fields in Wyoming, Colorado,
and Northwestern New Mexico, in: Problems of Petroleum Geology. American Association of Petro-
leum Geologists. https://​doi.​org/​10.​1306/​SV633​4C33
Kueper, B.H., McWhorter, D.B.: The use of macroscopic percolation theory to construct large-scale capil-
lary pressure curves. Water Resour. Res. 28, 2425–2436 (1992). https://​doi.​org/​10.​1029/​92WR0​1176
Lang, P.S., Paluszny, A., Zimmerman, R.W.: Permeability tensor of three-dimensional fractured porous
rock and a comparison to trace map predictions. J. Geophys. Res. Solid Earth 119, 6288–6307 (2014).
https://​doi.​org/​10.​1002/​2014J​B0110​27
Li, B., Li, L., Wen, X.H., Sun, T., 2023. Bridging Computational Stratigraphy and Reservoir Simulation for
Geologically Realistic High-Resolution Reservoir Modeling, in: Society of Petroleum Engineers - SPE
Reservoir Simulation Conference, RSC 2023. Society of Petroleum Engineers. https://​doi.​org/​10.​2118/​
212244-​MS
Lohne, A., Virnovsky, G., Durlofsky, L.J.: Two-stage upscaling of two-phase flow: from core to simulation
scale. SPE J. 11, 304–316 (2006). https://​doi.​org/​10.​2118/​89422-​PA
Matthäi, S.K., Burney, C., 2018. Scale-transgressive simulation of the impact of heterogeneity on CO 2
migration using high-resolution outcrop-analog models, in: 14th Greenhouse Gas Control Technolo-
gies Conference Melbourne. pp. 21–26.
Matthai, S.K., Tran, L.K.: Numeric determination of relative permeability of heterogeneous porous media
with capillary discontinuities. Adv. Water Resour. (2023). https://​doi.​org/​10.​1016/j.​advwa​tres.​2023.​
104430
Mikes, D., Barzandji, O.H.M., Bruining, J., Geel, C.R.: Upscaling of small-scale heterogeneities to flow
units for reservoir modelling. Mar. Pet. Geol. 23, 931–942 (2006). https://​doi.​org/​10.​1016/j.​marpe​tgeo.​
2005.​06.​005
Mishra, A., Kurtev, K.D., Haese, R.R.: Composite rock types as part of a workflow for the integration of
mm-to cm-scale lithological heterogeneity in static reservoir models. Mar. Pet. Geol. 114, 104240
(2020). https://​doi.​org/​10.​1016/j.​marpe​tgeo.​2020.​104240

13
Computing Relative Permeability and Capillary Pressure of… 1753

Moreno, Z., Anto-Darkwah, E., Rabinovich, A.: Semi-analytical modeling of rate-dependent relative perme-
ability in heterogeneous formations. Water Resour. Res. (2021). https://​doi.​org/​10.​1029/​2021W​R0297​
10
Nick, H.M., Matthäi, S.K.: Comparison of three FE-FV numerical schemes for single- and two-phase flow
simulation of fractured porous media. Transp. Porous Media 90, 421–444 (2011a). https://​doi.​org/​10.​
1007/​s11242-​011-​9793-y
Nick, H.M., Matthäi, S.K.: A hybrid finite-element finite-volume method with embedded discontinuities for
solute transport in heterogeneous media. Vadose Zone J. 10, 299–312 (2011b). https://​doi.​org/​10.​2136/​
vzj20​10.​0015
Oden, J.T., Belytschko, T., Fish, J., Hughes, T.J.R., Johnson, C., Keyes, D., Laub, A., Petzold, L., Srolovitz,
D., Yip, S., 2006. Simulation-Based Engineering Science.
Odsæter, L.H., Berg, C.F., Rustad, A.B.: Rate dependency in steady-state upscaling. Transp. Porous Media
110, 565–589 (2015). https://​doi.​org/​10.​1007/​s11242-​015-​0573-y
Øren, P.E., Ruspini, L.C., Saadatfar, M., Sok, R.M., Knackstedt, M., Herring, A.: In-situ pore-scale imaging
and image-based modelling of capillary trapping for geological storage of CO2. Int. J. Greenhouse Gas
Control 87, 34–43 (2019). https://​doi.​org/​10.​1016/j.​ijggc.​2019.​04.​017
Osoba, J.S., Richardson, J.G., Kerver, J.K., Hafford, J.A., Blair, P.M.: Laboratory measurements of relative
permeability. J. Petrol. Technol. 1(4), 369–382 (1951)
Paluszny, A., Matthäi, S.K., Hohmeyer, M.: Hybrid finite element?finite volume discretization of complex
geologic structures and a new simulation workflow demonstrated on fractured rocks. Geofluids 7, 186–
208 (2007). https://​doi.​org/​10.​1111/j.​1468-​8123.​2007.​00180.x
Perrin, J.C., Benson, S.: An experimental study on the influence of sub-core scale heterogeneities on CO2
distribution in reservoir rocks. Transp. Porous Media 82, 93–109 (2010). https://​doi.​org/​10.​1007/​
s11242-​009-​9426-x
Pickup, G.E., Hern, C.Y.: The development of appropriate upscaling procedures. Transp. Porous Media 46,
119–138 (2002)
Pickup, G.E., Sorbie, K.S.: The scaleup of two-phase flow in porous media using phase permeability ten-
sors. SPE J. 1(04), 369–382 (1994)
Pickup, G.E., Stephen, K.D.: An assessment of steady-state scale-up for small-scale geological models.
Pet. Geosci. (2000). https://​doi.​org/​10.​1144/​petgeo.​6.3.​20
Pickup, G., Ringrose, P.S., Sharif, A.: Steady-state upscaling: from lamina-scale to full-field model. SPE
J. 5, 208–217 (2000). https://​doi.​org/​10.​2118/​62811-​PA
Rabinovich, A., Li, B., Durlofsky, L.J.: Analytical approximations for effective relative permeability
in the capillary limit. Water Resour. Res. 52, 7645–7667 (2016). https://​doi.​org/​10.​1002/​2016W​
R0192​34
Revil, A., Cathles, L.M., Shosa, J.D., Pezard, P.A., de Larouzière, F.D.: Capillary sealing in sedimen-
tary basins: a clear field example. Geophys. Res. Lett. 25, 389–392 (1998). https://​doi.​org/​10.​1029/​
97GL0​3775
Ringrose, P.: How to Store CO2 Underground: Insights from early-mover CCS Projects. Springer, Cham
(2020). https://​doi.​org/​10.​1007/​978-3-​030-​33113-9
Ringrose, P., Bentley, M.: Reservoir Model Design. Springer, Dordrecht (2015). https://​doi.​org/​10.​1007/​
978-​94-​007-​5497-3
Ringrose, P.S., Sorbie, K.S., Corbett, P.W.M., Jensen, J.L.: Immiscible flow behaviour in laminated
and cross-bedded sandstones. J. Pet. Sci. Eng. 9, 103–124 (1993). https://​doi.​org/​10.​1016/​0920-​
4105(93)​90071-L
Rubin, D.M., Carter, C.L., 2005. Bedforms 4.0 : MATLAB Code for Simulating Bedforms and
Cross-Bedding.
Rubin, D.M., 1987. Cross-Bedding, Bedforms, and Paleocurrents: Concepts in Sedimentology and Pale-
ontology, Society of economic paleontologists and mineralogists.
Sech, R.P., Jackson, M.D., Hampson, G.J.: Three-dimensional modeling of a shoreface-shelf parase-
quence reservoir analog: Part 1. Surface-based modeling to capture high-resolution facies archi-
tecture. Am. Assoc. Pet. Geol. Bull. 93, 1155–1181 (2009). https://​doi.​org/​10.​1306/​05110​908144
Shao, Q., Boon, M., Youssef, A.A., Kurtev, K., Benson, S.M., Matthai, S.K.: Modelling CO2 plume
spreading in highly heterogeneous rocks with anisotropic, rate-dependent saturation functions: a
field-data based numeric simulation study of Otway. Int. J. Greenhouse Gas Control (2022). https://​
doi.​org/​10.​1016/j.​ijggc.​2022.​103699
Shosa, J.D., Cathles, L.M., 2001. Experimental Investigation of Capillary Blockage of Two Phase Flow
in Layered Porous Media, in: Petroleum Systems of Deep-Water Basins: Global and Gulf of Mex-
ico Experience: 21st Annual. Society of economic paleontologists and mineralogists, pp. 725–740.
https://​doi.​org/​10.​5724/​gcs.​01.​21.​0725

13
1754 A. A. Youssef et al.

Tran, L.K., Kim, J.C., Matthäi, S.K.: Simulation of two-phase flow in porous media with sharp material
discontinuities. Adv. Water Resour. (2020). https://​doi.​org/​10.​1016/j.​advwa​tres.​2020.​103636
Trevisan, L., Krishnamurthy, P.G., Meckel, T.A.: Impact of 3D capillary heterogeneity and bedform
architecture at the sub-meter scale on CO2 saturation for buoyant flow in clastic aquifers. Int. J.
Greenhouse Gas Control 56, 237–249 (2017). https://​doi.​org/​10.​1016/j.​ijggc.​2016.​12.​001
Van Duijn, C.J., Molenaar, J., De Neef, M.J.: The Effect of Capillary Forces on Immiscible Two-Phase Flow
in Heterogeneous Porous Media, Transport in Porous Media. Kluwer Academic Publishers (1995)
Virnovsky, G.A., Friis, H.A., Lohne, A.: A steady-state upscaling approach for immiscible two-phase
flow. Transp. Porous Media 54, 167–192 (2004)
Wen, X.-H., Durlofsky, L.J., Edwards, M.G.: Use of border regions for improved permeability upscaling.
Math. Geol. (2003). https://​doi.​org/​10.​1023/A:​10262​30617​943
Wen, X.-H., Durlofsky, L.J., Lee, S.H., Edwards, M.G., 2000. Full Tensor Upscaling of Geologically
Complex Reservoir Descriptions, in: SPE Annual Technical Conference and Exhibition, Dallas,
Texas. SPE, pp. 129–139. https://​doi.​org/​10.​2118/​62928-​MS
Willis, B.J., White, C.D.: Quantitative outcrop data for flow simulation. J. Sediment. Res. 70, 788–802
(2000). https://​doi.​org/​10.​1306/​2DC40​938-​0E47-​11D7-​86430​00102​C1865D
Willis, B.J., Bhattacharya, J.P., Gabel, S.L., White, C.D.: Architecture of a tide-influenced river delta in
the Frontier Formation of central Wyoming, USA. Sedimentology 46, 667–688 (1999)
Youssef, A., Tertois, A., Shao, Q., Matthai, S., 2022. Simulating Unsteady CO2 Flow through Brine-
Saturated Cross-Bedded Sandstones: Towards Relative Permeability Curves for Unstable Displace-
ment, in: ECMOR 2022. European Association of Geoscientists & Engineers, pp. 1–19. https://​doi.​
org/​10.​3997/​2214-​4609.​20224​4019
Youssef, A.A.A., Shao, Q., Matthäi, S.K., 2023. Two-Step Upscaling of Sub-Seismic Geo-Heterogeneity
with Flow-Rate-And Direction Dependent Saturation Functions, in: Society of Petroleum Engineers
- SPE Reservoir Simulation Conference, RSC 2023. Society of Petroleum Engineers. https://​doi.​
org/​10.​2118/​212215-​MS
Zhang, Y., Kogure, T., Chiyonobu, S., Lei, X., Xue, Z.: Influence of heterogeneity on relative permeability
for CO 2/Brine: CT observations and numerical modeling. Energy Procedia (2013). https://​doi.​org/​10.​
1016/j.​egypro.​2013.​07.​241

Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and
institutional affiliations.

13