# Probabilistic Forecast of Multiphase Transport under Viscous and Buoyancy
Forces in Heterogeneous Porous Media

Farzaneh Rajabi  Hamdi A. Tchelepi  Department of Energy Resources
Engineering, Stanford University, Stanford, CA 94305, USA

###### Abstract

In this study, we develop a probabilistic approach to map the parametric
uncertainty to the output state uncertainty in first-order hyperbolic
conservation laws. We analyze this problem for nonlinear immiscible two-phase
transport in heterogeneous porous media in the presence of a stochastic
velocity field. The uncertainty in the velocity field can arise from the
incomplete description of either porosity field, injection flux, or both. The
uncertainty in the total-velocity field leads to the spatiotemporal
uncertainty in the saturation field. Given information about the
spatial/temporal statistics of the correlated heterogeneity, we leverage
method of distributions to derive deterministic equations that govern the
evolution of single-point CDF of saturation. Unlike Buckley Leverett equation,
the equation for the raw CDF function is linear in space and time. Hereby, we
give routes to circumventing the computational cost of Monte Carlo scheme
while obtaining the full statistical description of saturation. We conduct a
set of numerical experiments and compare statistics of saturation computed
with the method of distributions, against those obtained using the statistical
moment equations approach and kernel density estimation post-processing of
high-resolution Monte Carlo simulations. This comparison demonstrates that the
CDF equations remain accurate over a wide range of statistical properties,
i.e. standard deviation and correlation length of the underlying random
fields, while the corresponding low-order statistical moment equations
significantly deviate from Monte Carlo results, unless for very small values
of standard deviation and correlation length.

###### keywords:

uncertainty quantification, stochastic partial differential equations, gravity
segregation, multiphase flow, Buckley Leverett equation, method of
distributions, nonlinear hyperbolic conservation laws

##  1 Introduction

Long-term resource allocation, perspective planning, and risk mitigation of
available energy resources rely on the accurate modeling of physical
phenomena, including fluid flow and transport systems in the complex large-
scale geological formations, e.g. oil reservoirs as well as carbon capture and
sequestration equipment. Deterministic modeling of flow and transport in such
systems has been widely studied for both single-phase and multiphase flow by
[9], [41], [42]. Rigorous modeling of these systems primarily depends on
sophisticated strategies for uncertainty quantification and stochastic
treatment of the system under study. Such an uncertainty is inherent to, and
critical for any physical modeling, essentially due to the incomplete
knowledge of state of the world, noisy observations, and limitations in
systematically recasting physical processes in a suitable mathematical
framework. These sources of uncertainty, render stochastic modeling of partial
differential equations particularly crucial. To this end, accurate predictions
of outputs (e.g. saturation fields) from reservoir simulations guarantee
precise oil recovery forecasts. These quantitative predictions rely on the
quality of the input measurements/data, such as the reservoir permeability and
porosity fields as well as forcings, such as initial and boundary conditions.
However, the available information about a particular geologic formation (e.g.
from well logs and seismic data of an outcrop) is usually sparse and
inaccurate compared to the size of the natural system and the complexity
arising from multi-scale heterogeneity of the underlying system. Eventually,
the uncertainty in the flow prediction can have a huge impact on the oil
recovery. Therefore, we can no longer rely on a naive deterministic model, and
we need to develop probabilistic flow predictions. Thereupon, because of
inherent uncertainties in the input data, originating from sparsity or
inaccuracy in the measurements, a stochastic framework is usually employed to
describe their variability. Such a probabilistic framework is the key building
block for making quantitative analysis of the uncertainty propagation from
system inputs to the output state variable of interest. Full probability
distribution of the random variables are then the essential metric to quantify
the uncertainty in the outputs. Within this stochastic framework, predicting
the saturation fields consists of determining the saturation probability
distributions, by leveraging the knowledge of probabilistic distribution of
the spatiotemporal uncertain random inputs.

Probabilistic analysis of single-phase flow and transport has been explored
and well-addressed by many researchers. Different approaches have been adopted
for this purpose, among which earlier works of [11], and [34] have presented a
comprehensive study for the stochastic hydrology. Such stochastic frameworks
quantify the associated uncertainty of the previously studied deterministic
models. Later works have focused on deriving deterministic equations for
leading statistical moments of the system states, called Statistical Moments
Equations (SME) method. This approach employed for single-phase contaminant
transport has been investigated by [17], [12], [39], [24], [8]. Even though
this scheme provides the first two statistical moments (ensemble mean and
variance) to forecast a system’s average response, they fail to predict the
probability of rare events, which is critical for risk assessment.

This deficiency in predictive capabilities of SME methods can be tackled by
using the probabilistic method of distributions, i.e. probability density
function (PDF) or cumulative density function (CDF) methods. The scheme was
first adopted in the turbulence flow literature by [33]. The PDF/CDF methods
provide a full probabilistic distribution by deriving deterministic equations
for PDF or CDF of the system states. PDF methods, sometimes require dealing
with closures that are problem dependent and sometimes not conveniently
attainable. This method has been utilized for the particle-base single-phase
transport of stochastic total-velocity field using a Markovian approximation
by [28], [27], and [29]. A More recent study by [40] has presented a
comparison between CDF method, moment equations (SME) approach, and Monte
Carlo scheme. There is a key advantage in using CDF method over their PDF
counterparts, as formulation of boundary conditions for PDF method is not
unique, and the closure problems might require burdensome manipulation. In
contrast, the boundary conditions for CDF scheme tend to be uniquely
formulated.

Stochastic treatment of immiscible multiphase transport i.e. Buckley Leverett
model in highly heterogeneous media, poses additional challenges compared to
their single-phase analogue. This originates from nonlinearity of the
transport equation and coupling of the total-velocity, pressure and saturation
fields. The nonlinearity stems from the fact that the uncertain velocity field
depends on the system state itself. Previous studies focused on using
streamline-base strategies for solving flow (pressure and velocity), and
transport (saturation) problems, i.e. Buckley Leverett equation, and provided
the statistics for the travel-time metric. Several studies in this context
were conducted by [13] for single-point distribution of one-dimensional
highly-heterogeneous porous media, in which the cumulative distribution
function obtained from the streamline-based simulation approach was presented
to be more efficient than Monte Carlo modeling of the nonlinear transport
equation. Moreover, [15] investigated single-point multi-dimensional
heterogeneous formations by using the distribution scheme called frozen (time
independent) streamline method (FROST). Another study by [14] extended the
single-point CDF approach to multi-point distribution of saturation in two-
phase transport, where they have computed the covariance function of
saturation in order to build more robust confidence intervals. Another study
by [10] extended the argument of [13] to estimate the statistics of fluid
saturations in the channelized porous systems by using the frozen streamline
assumption of FROST scheme. The SME methods were extended to be used in the
multiphase flow problems as well. Numerous studies were performed by [44] for
a Lagrangian treatment of moment equations subject to a stochastic
permeability field, and also for the Eulerian moment equations by [16], for
forward and inverse modeling by [26]. Method of distributions, discussed
previously for the uncertainty quantification of single-phase flow, remains as
one of the most robust and accurate approaches when extended to multiphase
flow models. Distribution methods provide full probabilistic description of
the saturation field. [37] proposed a framework to estimate one-point CDF of
saturation for one-dimensional Buckley Leverett model, in which the total
velocity field is random with a known given distribution, due to the uncertain
total flux as a function of time. They derive a general deterministic equation
for the single-point CDF of saturation for a horizontal reservoir, in which
gravity segregation plays no role. In this framework, a deterministic equation
is derived for a quantity called raw CDF, which conceptualizes the Heaviside
function over the saturation domain per se.

More theoretical approaches have used polynomial chaos expansion method for
uncertainty assessment of both single-phase and multiphase flow physical
frameworks. These methods have attracted interest because of their faster
predictive power compared to the more commonly used Mote Carlo scheme.
Stochastic spectral methods have been also utilized to estimate the statistics
of the state variable with a direct probabilistic collocation method by [23].
Besides the collocation methods, Pettersson and Tchelepi have recently applied
stochastic Galerkin method to derive statistics of the first moment for
saturation. The method is based on truncating the Karhunen-Loeve (KL)
expansion of random processes. A major advantage of these methods is the fact
that convergence analysis can be conducted. However, in practice, they suffer
from two disadvantages; one is the curse of dimensionality, the other is the
fact that they require further sophistication for highly nonlinear problems as
proposed by [32].

Monte Carlo method has been widely used to solve stochastic differential
equations describing two-phase immiscible flow in heterogeneous media. A large
number of equiprobable realizations of the reservoir description serve as
input for the corresponding deterministic governing equation. MC methods, even
though suffering from slow convergence rate (inversely proportional to the
square root of the number of realizations) of statistics that define dynamic
quantities of interest, are one of the simplest methods to implement. PDF and
CDF of the output variable from MC simulations is then computed by some
statistical post-processing. The robustness and easy implementation makes them
one of the first candidates to be chosen. Also their convergence is guaranteed
by the law of large numbers. Monte Carlo simulations can be accelerated using
Multilevel Monte Carlo (MLMC) approach. [30] used MLMC to accelerate the
uncertainty quantification of streamline simulation-based two-phase flow model
with the stochastic underlying permeability field.

Numerous studies including [35], [22], [20], [21], [18], [38], [43], [25], [3]
have investigated deterministic behavior of buoyancy forces interacting with
viscous and capillary forces in multiphase hyperbolic transport models.

Moreover, while uncertainty assessment of the hyperbolic conservation laws
(kinematic wave models) has been explored in several studies by [6], [43] and
[19], to the best of our knowledge, no study has been conducted on evaluating
the single-point distribution of saturation for stochastic nonlinear first-
order hyperbolic conservation laws in two-phase flow models, in which buoyancy
forces have been taken into account. Our study extends the idea of [37] to
propose an efficient distribution method for estimating the single-point
cumulative distribution function (CDF) of saturation for the Buckley Leverett
model subject to the gravitational forces. More broadly, we have investigated
the horizontal, updip and downdip flooding as well as pure gravity segregation
scenarios. This framework results in a full probabilistic description of the
saturation field for different physical set ups. While high heterogeneity is
traditionally regarded as an obstacle to designing stochastic methods, this
distribution method exploits the physics of highly heterogeneous porous
formation. The method identifies the underlying random fields that explain the
saturation distribution, together with the nonlinear mapping from these random
fields to the saturation domain. We show that the statistics and distributions
of these random fields can be estimated very efficiently and simply, leading
to a method that is computationally inexpensive compared to the exhausting
Monte Carlo (MC) simulations.

In section 2, we review the governing equations for the Buckley Leverett
model, and explain the differences in flux functions of horizontal vs.
gravitational multiphase case studies. We then explain the CDF method in
section 3, and derive a general deterministic equation for the (single-point)
CDF of saturation for the case studies with gravitational forces, subject to a
random porosity field and/or random injection flux. In section 4, we
illustrate the numerical aspects of the nonlinear hyperbolic laws with three
different approaches, Monte Carlo (MC), method of distributions (MD) and
statistical moment equation (SME) approach. We perform the Monte Carlo
simulation for the Riemann problem of nonlinear BL equation using the mass
conservative finite volume Godunov scheme, and then compare its solution for
several test cases against those of our MD method and the low-order SME
approximation, in order to validate our MD framework. A more comprehensive
analysis of low-order approximation for different random quantities is offered
in section 5. We conclude this work by analyzing the convergence and accuracy
of our method in sections 6,7, followed by the conclusion section in 8.

##  2 Problem formulation

We consider the nonlinear, incompressible, two-phase (Darcy) flow in a
heterogeneous reservoir with intrinsic heterogeneous permeability $K$ and
spatially-varying porosity $\phi(x)$ fields, where there are two immiscible
phases; a wetting phase (e.g. water) is injected to displace a non-wetting
phase (e.g. oil) through an imbibition process such as waterflooding. The key
underlying physical mechanisms are viscous and gravitational forces, whereas
capillarity, chemical reactions, diffusion, and changes of state are
neglected. We are primarily interested in estimating saturation fields while
quantifying the uncertainty associated with their description for horizontal
domains as well as vertical reservoirs in which gravity segregation plays a
major role throughout the displacement process.

The inevitable uncertainty in the response parameters stems from several
sources of uncertainty in the given input information of the reservoir. These
input specifications include initial saturations of the wetting and non-
wetting phases ($S_{wi}$ and $S_{oi}$), boundary conditions for the wetting
phase ($S_{w}(x=0,t)$ at the inlet of domain), parameters representing the
fluid properties, such as wetting and non-wetting phase viscosities ($\mu_{w}$
and $\mu_{o}$), parameters describing rock properties, such as porosity and
permeability ($\phi(x)$ and $K(x)$), parameters outlining the fluid-rock
interactions, such as relative permeabilities ($K_{rw}(S_{w})$,
$K_{ro}(S_{o})$), as well as the injection and production pressures of the
wetting phase ($P_{inj}$ and $P_{prod}$), where they are all assumed to be
known parameters of the problem under study.

An ideal prediction scenario calls for the complete knowledge of the
aforementioned parameters. However, data sparsity and heterogeneity of the
input definitions render such a perfect framework unattainable. The primary
justification to consider the inputs as stochastic variables is lack of
knowledge, i.e. having limited measurements which are only available at well
locations, for the static properties (porosity, permeability) as well as
dynamic properties (saturation, injection/production rates) of the reservoir.
The uncertainty can also arise from random initial/boundary conditions, as
well as fluid-rock interactions. Unlike rock properties that are inherent to
the specific subsurface of evaluation, and therefore may not be quantified
deterministically, the fluid properties are usually accessible through lab
experiments and are often taken as deterministic values. To this end, we focus
on the uncertainty in porosity field $\phi(x)$ and total Darcy flux $q(t)$
(the coefficients of the PDE are stochastic) in this study, whereas the rest
are assumed to be deterministically known parameters. More precisely, we
describe these random parameters using a prescribed PDF, i.e. by assigning a
specific probabilistic description $p_{\varphi}(\phi)$ and $p_{q}(Q)$ to them,
fully characterized by their known mean ($\mu$), standard deviation
$(\sigma^{2})$, and correlation structure $(Cr)$. Once such a probabilistic
setting for the input parameters is established, we can propagate the
uncertainty from inputs to output quantities of interest, essentially by
solving a coupled system of nonlinear transport equation for saturation which
satisfies mass conservation, along with the Darcy equation for pressure field.
This comprehensive framework enables us to identify a probabilistic spatio-
temporal description for the water saturation.

The nonlinear hyperbolic conservation law with a non-convex flux is
conventionally studied by employing the Buckley Leverett (BL) equation which
is the result of combining mass balance for each phase $\alpha=w,o$ with
Darcy’s law by [5]. The two-phase extension of the single-phase Darcy’s law
for wetting and non-wetting phases leads to,

$\displaystyle\textbf{u}_{w}=-\dfrac{kk_{rw}}{\mu_{w}}\nabla\varphi_{w}=-\dfrac{kk_{rw}}{\mu_{w}}\nabla(P_{w}+\rho_{w}gh)=-\dfrac{kk_{rw}}{\mu_{w}}\nabla(P_{o}-P_{c}+\rho_{w}gh)$
(1)
$\displaystyle\textbf{u}_{o}=-\dfrac{kk_{ro}}{\mu_{o}}\nabla\varphi_{o}=-\dfrac{kk_{ro}}{\mu_{o}}\nabla(P_{o}+\rho_{o}gh)=-\dfrac{kk_{ro}}{\mu_{o}}\nabla(P_{o}+\rho_{o}gh)$
(2)

Where capillary pressure relating the pressures of two phases is
$P_{c}(S_{w})=P_{o}-P_{w}$. In the above-mentioned equations,
$\textbf{u}_{\alpha}(\alpha=\\{w,o\\})$ are the phase Darcy velocities of the
wetting (water) and nonwetting (oil) phases. $k_{r\alpha}$, $\mu_{\alpha}$ and
$\rho_{\alpha}$ $(\alpha=\\{w,o\\})$ represent relative permeability,
viscosity and density of each phase, respectively. In the case of vertical
reservoirs, $h$ is the vertical distance travelled by each phase. Applying
mass conservation on a control volume results in the transport equation
(saturation equation),

$\displaystyle\phi(\textbf{x})\dfrac{\partial S_{\alpha}}{\partial
t}+\nabla\cdot\textbf{u}_{\alpha}=0\quad\text{for $\alpha=\\{w,o\\}$}$ (3)

Where the phase saturations, $S_{w}$ and $S_{o}$ are the volume fraction of
the pore space occupied by the corresponding fluid phase. Consequently,
$S_{w}+S_{o}=1$ simplifies the two mass conservation equations for the wetting
and nonwetting phases, to the incompressiblity condition of
$\nabla\cdot\textbf{u}_{T}=0$ for total Darcy flux. Subsequently, we obtain an
equation in terms of pressures,

$\displaystyle\textbf{u}_{T}=\textbf{u}_{w}+\textbf{u}_{o}=-\lambda_{T}\nabla
P_{o}+\lambda_{w}\nabla P_{c}-(\lambda_{o}P_{o}+\lambda_{w}P_{w})g\nabla h$
(4)

Where the relative mobilities are defined as
$\lambda_{\alpha}=\dfrac{kk_{r\alpha}}{\mu_{\alpha}}$ $(\alpha=\\{w,o\\})$ for
each phase and the total mobility reads $\lambda_{T}=\lambda_{o}+\lambda_{w}$.
The three terms on the right-hand side of Eq. (4) represent viscous, capillary
and buoyancy forces, respectively. The dimensionless flux function of the
wetting phase, referred to as the fractional flow of water is defined as the
wetting phase velocity divided by the total velocity,

$\displaystyle
f_{w}=\dfrac{\textbf{u}_{w}}{\textbf{u}_{T}}=\dfrac{\lambda_{w}}{\lambda_{T}}\bigg{(}1+\dfrac{\lambda_{o}}{\textbf{u}_{T}}\left(\nabla
P_{c}-\Delta\rho g\nabla h\right)\bigg{)}$ (5)

In the absence of capillarity, the two phase pressures are equal and identical
to the global pressure. For such a scenario, we consider one-dimensional flow
in an inclined reservoir with a constant dip angle $\theta$, where
$0\leq\theta\leq 90$, ranging from 0 for the horizontal reservoirs to 1 for
the vertical case studies, resulting in $\nabla h=sin\theta$. The
dimensionless gravity number $N_{g}$ is defined as the ratio of buoyancy to
viscous forces,

$\displaystyle N_{g}=\dfrac{kg(\rho_{w}-\rho_{o})}{\mu_{o}\textbf{u}_{T}}$ (6)

Subsequently, defining the viscosity ratio $m$ as
$m=\dfrac{\mu_{w}}{\mu_{o}}$, while employing the definitions for mobilities,
the fractional flow of water can be recast in the following form,

$\displaystyle
f_{w}(S_{w},\theta)=\dfrac{k_{rw}}{k_{rw}+mk_{ro}}\bigg{(}1-N_{g}k_{ro}sin\theta\bigg{)}$
(7)

Rewriting the mass conservation of water saturation (Eq.(3)) using the
fractional flow $f_{w}$, while applying the chain rule and incompressibility
condition leads to,

$\displaystyle\dfrac{\partial S_{w}}{\partial
t}+\textbf{v}_{T}(x,t)\cdot\nabla f_{w}(S_{w})=0$ (8)

Where $\textbf{v}_{T}(x,t)$ is defined as the interstitial velocity (total
seepage velocity field)
$\textbf{v}_{T}(\textbf{x},t)=\dfrac{q_{T}(t)}{\phi(\textbf{x})}$. Where
$q_{T}$ is the total volumetric flow rate, obtained by multiplying the total
Darcy velocity by the area A (A is assumed to be 1 in this work). we assume
the total flux is time-dependent. This is a feasible assumption because
incompressibility condition imposes $q_{T}(x,t)$ to be x-independent and hence
equal to the injection flow rate at the inlet boundary which could be time-
dependent. That is, $q_{T}(x,t)=q_{T}(x=0,t)=q_{T_{inlet}}(t)$. In this study,
we assume the total velocity field $\textbf{v}_{T}(x,t)$ is a random variable,
where its uncertainty arises from two different sources/physical phenomena;
the spatially-varying porosity field is uncertain in space, whereas the
temporally-varying total flux is a random variable in time. We will
investigate scenarios in which either or both of these phenomena are taken to
be stochastic.
It should be noted that fractional flow $f(S_{w})$ is a continuous smooth, and
hence differentiable function. To this end, $\nabla
f_{w}(S_{w})=f_{w}^{\prime}(S_{w})\nabla S_{w}$. Thereupon, supplemented with
the following initial/boundary conditions, Buckley Leverett equation defines
an IBVP for water saturation which reads,

$\displaystyle\dfrac{\partial S_{w}}{\partial
t}+\textbf{v}_{T}(x,t)f^{\prime}_{w}(S_{w})\cdot\nabla
S_{w}=0,\qquad\qquad\text{where}\quad\textbf{v}_{T}=\dfrac{q_{T}(t)}{\phi(x)}$
$\displaystyle
S_{w}(\textbf{x},t)=S_{B}=1-S_{oi},\quad\textbf{x}\in\Gamma_{i},t>0,$
$\displaystyle S_{w}(\textbf{x},t=0)=S_{wi},\quad\textbf{x}\in\Omega$ (9)

Where $S_{wi}$ is the initial water saturation in the heterogeneous reservoir
domain $\Omega$. In other words, the reservoir is assumed to be initially oil-
saturated with a small irreducible water saturation $S_{wi}$. $S_{B}$ is water
saturation at the injection boundary $\Gamma_{i}$, i.e. $x=0$, where the
wetting phase is being injected. $S_{oi}$ is the irreducible oil saturation.
$\Gamma_{o}$ is the boundary along which the flow outlet happens. Although
initial/boundary conditions may be in general random variables, we assume they
are both deterministic constants. Also, in order to guarantee self-similarity
in solutions of saturation, we impose $S_{wi}$ and $S_{B}$ as uniform in space
and time-independent quantities. For the test cases with gravitational force
included, the initial condition is non-uniform. Consequently, this equation
admits analytical solutions in terms of $x/t$. Moreover, the following
Dirichlet boundary conditions for the pressure field (pressure control) is
considered,

$\displaystyle P(\textbf{x})=p_{inj},\qquad\textbf{x}\in\Gamma_{i}$
$\displaystyle P(\textbf{x})=p_{prod},\qquad\textbf{x}\in\Gamma_{o}$ (10)

Therefore, predictions of two-phase flow entails solving the first-order
nonlinear hyperbolic Buckley Leverett equation with a nonconvex flux (the flux
function changes its convexity within the given saturation interval). Solving
the BL equation provides the distribution of saturation as a function of space
and time. The evolution of the saturation field is described in terms of
nonlinear kinematic waves, where different complex wave propagation scenarios
happen depending on the initial/boundary conditions as well as behavior of the
flux function (i.e. whether the flux function includes gravity or not, and if
it does, how the convexity of flux function changes throughout the domain).
These scenarios include formation of shocks and rarefaction zones as well as
reflection of the shocks from endpoints of the domain. The latter case which
is more involved, is the physical phenomena behind the so called pure gravity
segregation in a sealed gravity column, in which gravitational forces are the
only driving mechanism, as will be elaborated more in the upcoming sections.

Illustrating behavior of the flux function requires the dip angle as well as
relative permeabilities to be fully specified. According to Eq. (7), the
fractional flow is indeed a function of saturation and the dip angle, and
hence can result in three different scenarios depending on the inclination
angle $\theta$,

  * 1.

$N_{g}sin\theta=0$ turns Eq. (7) to that of a horizontal reservoir, meaning
gravity plays no role and viscous forces are the only driving mechanism.

  * 2.

$N_{g}sin\theta>0$ represents an updip flooding, in which gravity retards the
water flow, thereby reducing velocity of the front movement. It is noteworthy
that $sin\theta>0$ can result in a fractional flow smaller than 0. Assuming
that the wetting phase is heavier than the nonwetting one, $f_{w}<0$
represents counter-current flow of the heavier wetting phase moving downdip,
while the lighter nonwetting phase is moving updip. The $0<f_{w}<1$ region
indicates the co-current flow of water and oil moving updip.

  * 3.

$N_{g}sin\theta<0$ results in a downdip flooding scenario, where gravity
increases the water flow, and hence causes a faster movement of the saturation
front. This scenario corresponds to the counter-current flow of oil flowing
updip, while water is flowing downdip in a tilted reservoir. $sin\theta<0$
leads to a fractional flow larger than 1, thereby causing the $f_{w}>1$ region
its fractional flow curve. In contrast, in the $0<f_{w}<1$ region, the two
phases are moving downdip co-currently.

(a) Horizontal flooding (b) Upslope background drift

(c) Downslope background drift (d) Gravity column

Figure 1: Continuous (analytical) fractional flow curves for the four cases
studies of our problem (horizontal, upslope background drift, downslope
background drift and pure gravity segregation). The dashed curves are the flux
function $f_{w}(S_{w})$. In cases (b) and (c) $f(S_{w})$ accounts for a
combination of gravitational forces and viscous background drift. The blue
lines show the concave envelopes corresponding to imbibition at the left edge
of the domain, while the red lines are the convex envelopes corresponding to
drainage at the right edge. In case (d), the inner envelope manifest
reflection of waves from boundaries. All cases have adopted Brooks and Corey
relative permeabilities, and have a viscosity ratio of m=0.5.

The flux function corresponding to these scenarios are schematically
represented in Fig. (1). Horizontal case ($N_{g}=0$) with an S-shaped curved
is characterized by a shock and a trailing rarefaction zone where viscous
forces are the only driving mechanism, whereas the vertical updip injection
($N_{g}=10$) is characterized by a shock and a rarefaction zone, and vertical
downdip injection ($N_{g}=-10$) case is characterized by two shocks and a
rarefaction zone in between. The two former cases represent the combined
effect of viscous and buoyancy forces. The case in sub-figure (d) is
physically similar to the vertical case, except that it has closed boundaries
at the top and bottom. Therefore, unlike the vertical flooding in which we
inject at one boundary, and recover oil from the other boundary, we assume
that in this case, buoyancy leads the whole physics with no injection/viscous
forces helping it. To this end, we consider a gravity column in a
heterogeneous domain $[0,3]$ confined by a sealing medium at the top and
bottom. We define an initial condition with a heavy fluid on top of a lighter
fluid, described by setting $S_{L}=1$ in the top half of the domain and
$S_{R}=0$ in the bottom half of the domain. The fluids are separated by a
sharp interface, and initially the interface is at $Z=1$. Then, the two fluids
start segregating, where the heavier fluid would like to move downward
(positive $z$ direction) and reside below the lighter, and hence the lighter
fluid moves upward (negative $z$ direction). While the shock regions are
representative of a single-phase fluid, the rarefaction area is indeed two-
phase. Therefore, the outer envelope shows the two waves before hitting the
boundaries for the first time. These waves move in the opposite directions and
will be eventually reflected from each boundary, schematically shown by the
inner envelope. Note that $m\neq 1$, and the corresponding flux function is
asymmetric, and hence, the two waves move asynchronously, i.e. one may hit its
boundary much earlier than the other, depending on $m$. Once reflected from
both boundaries, the waves follow a very small rarefaction region. Afterwards,
at $S^{\star\star}_{R}$ and $S^{\star\star}_{L}$ two shocks form again, start
moving towards each other and meet at a point $z_{j}$, and then reach to the
equilibrium at $t_{eq}$, where the initial condition has been completely
flipped.

The current study only investigates the imbibition process (indicated by blue
lines in Fig. (1)). Nonetheless, extending the manipulations of this work to
the drainage process (depicted by red lines) is unambiguously straightforward;
the horizontal flooding case would still entail inspection of a rarefaction
followed by a shock. While the shock-rarefaction scenarios for updip and
downdip would be flipped in the drainage process. That is, contrary to the
imbibition scenario, the downdip case gives rise to one shock, while its updip
counterpart consists of two shocks with a rarefaction fan in between.

To shed more light on the special case of pure gravity segregation, we will
elaborate more on the distinctions between this scenario and the formerly
explained downdip flooding, where we assumed injection is happening at the
inlet. Strictly speaking, we assume a vertical column in a heterogeneous
domain $[0,3]$ with sealing top and bottom boundaries, so $q_{T}=0$ in pure
gravity segregation. Assuming the heavier fluid on top and the lighter fluid
beneath gives rise to an unstable situation and needs to reach a vertical
equilibrium. Therefore, even though $\textbf{u}_{T}=0$, we have
$\textbf{u}_{o}\neq 0,\textbf{u}_{w}\neq 0$, where the continuity condition
results in $\textbf{u}_{o}+\textbf{u}_{w}=\text{constant}$. Assuming boundary
conditions $\textbf{u}_{o}=\textbf{u}_{w}=0$ at the top and bottom of the
domain, we conclude the constant is zero, and hence
$\textbf{u}_{o}=-\textbf{u}_{w}$. Consequently, $\textbf{u}_{T}=0$ entails a
revisited definition for $f_{w}$, as we can no longer use the conventional
definition in Eq. (5). Note that we are using u for the vertical velocity with
a small abuse of notation, however, we keep in mind u is referring to a
horizontal or vertical velocity depending on the case study.
In order to find the corresponding flux function, we start from Eq. (1) for
the vertical velocity of water and oil. After rearranging and subtracting
them, and applying $\text{u}_{o}=-\text{u}_{w}$ in the absence of capillarity,
we end up with,

$\displaystyle\textbf{u}_{w}=\dfrac{g(\rho_{o}-\rho_{w})k}{\mu_{w}}\dfrac{k_{rw}}{\bigg{(}1+\dfrac{k_{rw}\mu_{o}}{k_{ro}\mu_{w}}\bigg{)}}\equiv\textbf{u}_{c}G(S_{w})$
(11)

Where the first fraction is assumed to be a characteristic velocity resulted
from gravity $\textbf{u}_{c}$, and the second fraction is the dimensionless
flux function $G(S_{w})$. Hence, fractional flow for pure gravity segregation
is defined as $G(S_{w})=\dfrac{\textbf{u}_{w}}{\textbf{u}_{c}}$, and is equal
to,

$\displaystyle
G(S_{w})=\dfrac{k_{rw}}{\bigg{(}1+\dfrac{k_{rw}\mu_{o}}{k_{ro}\mu_{w}}\bigg{)}}$
(12)

Therefore, the transport equation is recast as,

$\displaystyle\dfrac{\partial S_{w}}{\partial
t}+\dfrac{\textbf{u}_{c}}{\phi(z)}G^{\prime}_{w}(S_{w})\cdot\nabla S_{w}=0$
$\displaystyle
S_{w}(\textbf{z},t)=S_{B}=1-S_{oi},\quad\textbf{z}\in\Gamma_{i},t>0,$
$\displaystyle
S_{w}(\textbf{z},t=0)=1-S_{oi},\quad\textbf{z}\in\Omega^{-}\equiv[-1,0]$
$\displaystyle
S_{w}(\textbf{z},t=0)=S_{wi},\quad\textbf{z}\in\Omega^{+}\equiv[0,+1]$ (13)

Numerical treatment of first order hyperbolic conservation laws poses
additional challenges, mainly due to their nonlinearity lying solely in the
flux function. Such a nonlinear behavior in fractional flow is essentially
introduced by relative permeabilities of each phase. There are several
empirical models for defining the relative permeabilities as a function of
water saturation. Brooks and Corey [4], and Van Genuchten [36] are the most
commonly employed empirical models. We employ the Brooks and Corey model to
define relative permeabilities given by the following expressions,

$\displaystyle
k_{rw}(S_{w})=\left(\dfrac{S_{w}-S_{wi}}{S_{B}-S_{wi}}\right)^{2},\qquad
k_{ro}(S_{w})=\left(\dfrac{S_{B}-S_{w}}{S_{B}-S_{wi}}\right)^{2}$ (14)

A closer inspection of Eq. (2) immediately reveals the fact that hyperbolic
Buckley Leverett equation is indeed in the form of an advection equation with
velocity $\textbf{v}_{s}(x,t,s)\equiv\dfrac{q(t)}{\phi(x)}f^{\prime}(S_{w})$.
Combining this relation with Eq. (7) and Eq. (14) leads to an expression for
the advection speed (for simplicity, we refer to $S_{w}$ as $S$),

$\displaystyle\textbf{v}_{s}(x,t,S)=\dfrac{q(t)}{\phi(x)}\dfrac{\partial
f_{w}(S)}{\partial S}$
$\displaystyle=\dfrac{q(t)}{\phi(x)}\dfrac{2m(S-S_{wi})(1-S-S_{oi})(1-S_{oi}-S_{wi})}{\bigg{[}(S-S_{wi})^{2}+m(1-S-S_{oi})^{2}\bigg{]}^{2}}$
$\displaystyle+N_{g}\dfrac{q(t)}{\phi(x)}\dfrac{2(1-S-S_{oi})(S-S_{wi})\bigg{[}(S-S_{wi})^{3}+m(1-S-S_{oi})^{3}\bigg{]}}{(1-S_{wi}-S_{oi})^{2}\bigg{[}(S-S_{wi})^{2}+m(1-S-S_{oi})^{2}\bigg{]}^{2}}$
(15)

Where the second term drops for the horizontal cases with a gravity number of
$(N_{g}=0)$. Similarly, for the pure segregation scenario,

$\displaystyle\textbf{v}_{s}(z,S)=\dfrac{\textbf{u}_{c}}{\phi(z)}\dfrac{\partial
G(S)}{\partial
S}=\dfrac{\textbf{u}_{c}}{\phi(z)}\dfrac{2m(1-S-S_{oi})(S-S_{wi})\bigg{[}(S-S_{wi})^{3}+m(1-S-S_{oi})^{3}\bigg{]}}{(1-S_{wi}-S_{oi})^{2}\bigg{[}(S-S_{wi})^{2}+m(1-S-S_{oi})^{2}\bigg{]}^{2}}$
(16)

A hyperbolic conservation law with piece-wise constant initial condition
$S(z,t=0)=S_{0}(z)=\begin{cases}S_{L}\quad z\leq z_{d}\\\ S_{R}\quad
z>z_{d}\end{cases}$ forms a Riemann problem, where $S_{L}$ and $S_{R}$
represent the left and right values of the initial condition, separated at a
discontinuity location $z_{d}$. While in general, $S_{L}$ and $S_{R}$ can be
random, we assume they are deterministic constants. As the flux function
changes its convexity throughout the saturation domain, the solution of this
Riemann problem can be comprised of shocks and rarefaction zones, depending on
whether the flux function is convex or concave over that specific interval.
The exact analytical solution to Eq. (2) subject to the given boundary/initial
conditions, in one-dimensional domain of [a,b], and a final simulation time T,
for a given realization of the total velocity $\textbf{v}_{T}(x,t)$, can be
expresses as following ([1], [31]),

$\displaystyle S_{w}(x,t;v)=\begin{cases}S_{B},&x<X_{1}\\\
S_{\text{r}}=(f_{w}^{\prime})^{(-1)}\Bigg{[}x\bigg{(}\dfrac{\bigintss_{0}^{t}q(t^{\prime})dt^{\prime}}{\varphi(x)}\bigg{)}^{-1}\Bigg{]},&X_{1}\leq
x<X_{2}\\\ S_{wi}&x\geq X_{2}\end{cases}$ (17)

Where $X_{1}=a$ and $X_{2}=x_{f}(t;v)$ for the no-gravity case. Whereas the
solution for downdip flooding and gravity column before reflection of the
waves takes into account the presence of two shocks, resulting in the same
solution above, however, $X_{1}=x_{f_{2}}(t;v)$ and $X_{2}=x_{f_{1}}(t;v)$,
where $x_{f_{i}}(t;v)$ for $i=1,2$ describes the front location at time t for
the given realization of total velocity. This quantity can be found by
satisfying the entropy constraint through Rankine-Hugoniot condition at the
location of discontinuities (shocks),

$\displaystyle\dfrac{dx_{f_{i}}}{dt}=\dfrac{q(t)}{\phi(x)}\dfrac{f_{w}(S^{\star}_{i})-f_{w}(S_{wi})}{S^{\star}_{i}-S_{wi}}$
(18)

Where the continuity condition imposed on the speed of the shocks determines
$S^{\star}_{i}$,

$\displaystyle\dfrac{f_{w}(S^{\star}_{i})-f_{w}(S_{wi})}{S^{\star}_{i}-S_{wi}}=f^{\prime}_{w}(S^{\star}_{i})$
(19)

For the special case of horizontal reservoir domain (with fractional flow
$f_{w}=\dfrac{k_{rw}}{k_{rw}+mk_{ro}}$ supplemented by Brooks and Corey
relative permeabilities), $S^{\star}$ has an analytical expression (Eq. 4.7 in
[14]). Nevertheless, we resort to a general numerical root finding scheme
which handles $N_{g}\neq 0$ cases as well.

###  2.1 Nondimensionalization

We will investigate the impact of uncertainty in the inputs on the output
saturation field, first by assuming that $\phi(x)$ is random and $q$ is a
deterministic constant, second by assuming that $q_{T}(t)$ is random and
$\phi$ is a deterministic constant, and third by assuming both are stochastic
fields. For the case in which $q_{T}(t)$ is random with the given mean
$\mu_{q}$ and correlation length $\tau_{q}$, we can make the transport Eq. (2)
dimensionless by defining,

$\displaystyle x_{D}=\dfrac{x}{L},\qquad
t_{D}=\dfrac{1}{L}\bigintssss_{0}^{t}q(t^{\prime})dt^{\prime}$ (20)

This rescaling will recast Eq. (2) as the dimensionless form in the one-
dimensional domain $\Omega$,

$\displaystyle\phi(x)\dfrac{\partial S_{w}}{\partial t_{D}}+\dfrac{\partial
f(S_{w})}{\partial x_{D}}=0$ (21)

For the cases in which $q$ is a deterministic constant, we can simply use the
dimensionless quantities below to make Eq. (2) dimensionless,

$\displaystyle x_{D}=\dfrac{x}{L},\qquad t_{D}=\dfrac{qt}{L}$ (22)

This definitions lead to the dimensionless form of Eq. (2),

$\displaystyle\phi(x)\dfrac{\partial S_{w}}{\partial t_{D}}+\dfrac{\partial
f(S_{w})}{\partial x_{D}}=0$ (23)

In order to make transport Eq. (2) for pure gravity segregation dimensionless,
since there is no injection flux for this scenario, we will only investigate
the effect of uncertainty in porosity field on the water saturation.
Therefore, we introduce the dimensionless time and distance parameters,

$\displaystyle
t_{D}=t\dfrac{g(\rho_{o}-\rho_{w})k}{\mu_{w}H}=\dfrac{t}{\tau},\qquad\qquad
z_{D}=\dfrac{z}{H}$ (24)

Where $H$ is the total length of the vertical domain and
$\tau=\dfrac{\mu_{w}H}{g(\rho_{o}-\rho_{w})k}$ is the characteristic time
scale for gravity segregation. That is, if $\tau=1$, the particle has enough
time to travel along the whole $H$. Such definitions allow us to work with the
dimensionless form of Eq. (2) as,

$\displaystyle\phi(z)\dfrac{\partial S_{w}}{\partial t_{D}}+\dfrac{\partial
G(S_{w})}{\partial z_{D}}=0$ (25)

Figure 2: Characteristics plots for horizontal, updip and downdip injection
(from left to right).

Figure 3: (left) Characteristic plot in the $z-t$ plane, (right) the
corresponding finite volume solution of the saturation profiles $S_{w}(z,t)$
for a gravity column with a heavy fluid on top of a light fluid initially
separated by a sharp interface at $z_{d}$. After reflection from boundaries
the two shocks join each other at an equilibrium point $z_{j}$. Asymmetric
pattern of the right and left branches arises from $m\neq 1$ and the
asymmetric non-uniform initial condition. Solid black lines on the
characteristics plot correspond to the shocks. Regions I, II, III, IV
represent the right-moving shock before and after reflection, left-moving
shock before and after reflection, respectively.

##  3 Single-point CDF equation for Saturation

In order to provide full probabilistic distribution of water saturation, we
will extend the idea of [37] to solve a deterministic equation that governs
the space-time evolution of the CDF of water saturation
$F_{s}(S_{w};\textbf{x},t)$ for more general cases of gravitational forces
taken into consideration. This approach overcomes the limitations of solving
the nonlinear hyperbolic Eq. (2), by introducing a linear hyperbolic equation
for the so-called random “raw” CDF function of water saturation,

$\displaystyle\Pi(\theta,S_{w};\textbf{x},t)=\mathcal{H}(\theta-
S_{w}(\textbf{x},t))$ (26)

where $\mathcal{H}(.)$ is the Heaviside step function and $\theta$ is a
deterministic value (outcome) that the random water saturation $S_{w}$ takes
at a space-time point (x, t).

For continuous solutions of Eq. (2), the linear stochastic hyperbolic equation
for the raw CDF reads [37],

$\displaystyle\dfrac{\partial\Pi}{\partial
t}+\textbf{v}_{s}(\theta,x,t)\cdot\nabla\Pi=0\qquad x\in\Omega,t>0$
$\displaystyle\Pi(\theta,S_{w};x,t=0)=\mathcal{H}(\theta-S_{wi})$
$\displaystyle\Pi(\theta,S_{w};x=0,t)=\mathcal{H}(\theta-1+S_{oi})$ (27)

Where the boundary condition corresponds to the injection of wetting phase
fluid, i.e. imbibition process at the inlet $x=0$. The advection velocity is
defined using Eq. (2) for the first three cases, and it follows Eq. (16) for
pure gravity segregation.

As mentioned before, Eq. (3) holds only for the continuous solutions of the
raw CDF, in which we are indeed excluding the treatment of shocks
(discontinuities). A general form of Eq. (3) comprises of imposing the non-
smooth solutions to fulfill entropy conditions at the location of
discontinuities, allowing for a physically meaningful solution. This is
primarily done by adding a kinetic defect term $\mathcal{M}(\theta,x,t)$ to
Eq. (3), making it $\dfrac{\partial\Pi}{\partial
t}+\textbf{v}_{s}(\theta,x,t)\cdot\nabla\Pi=\mathcal{M}(\theta,x,t)$, which
essentially results in a solution valid for the the entire domain of
saturation. Although, the higher dimensional problems require solving the
later equation as an alternative to Eq. (3) , their one-dimensional
counterparts do not suffer from additional complications emerged from
interaction of waves at higher dimensions, and hence, the problem in one-
dimension could be studied by employing Eq. (3) for the separate smooth parts
of the domain. To this end, for horizontal displacement and updip flooding, we
divide the solution of Eq. (3), into two sub-regions similarly to [37],

$\displaystyle\Pi(\theta,x,t)=\begin{cases}\Pi_{a}=\mathcal{H}(\theta-
S_{wi}),\qquad\hfill&S_{wi}<\theta<S^{\star},\qquad\hfill x_{f}(t;v)<x<b\\\
\Pi_{b}=\mathcal{H}(\theta-
S_{r}),\qquad\hfill&S^{\star}<\theta<1-S_{oi},\qquad\hfill
a<x<x_{f}(t;v)\end{cases}$ (28)

Whereas for the downdip flooding case as well as the pure segregation scenario
before the fastest-moving wave reflects from its boundary, we deal with three
sub-regions. That is, for $0<t<t_{ref_{1}}$,

$\displaystyle\Pi(\theta,x,t)=\begin{cases}\Pi_{a}=\mathcal{H}(\theta-
S_{wi}),\qquad\hfill&S_{wi}<\theta<S_{R}^{\star},\qquad\hfill
x_{f_{R}}(t;v)<x<b\\\ \Pi_{b}=\mathcal{H}(\theta-
S_{r}),\qquad\hfill&S_{R}^{\star}<\theta<S_{L}^{\star},\qquad\hfill
x_{f_{L}}(t;v)<x<x_{f_{R}}(t;v)\\\
\Pi_{c}=\mathcal{H}(S_{B}-\theta),\qquad\hfill&S_{L}^{\star}<\theta<1-S_{oi},\qquad\hfill
a<x<x_{f_{L}}(t;v)\end{cases}$ (29)

It should be noted that $\Pi_{a}$ and $\Pi_{c}$ are independent of $x$ and
$t$. After the waves reflect from boundaries in the pure segregation case,
three separate rarefaction fans along with two shocks in between develop and
move towards each other. As is represented in Fig 3, the rarefaction region in
the middle is indeed the continuation of the characteristics in the
rarefaction fan of the region $\pi_{b}$ in Eq. 29, up to the time $t_{j}$.
Eventually, when the two shocks meet at time $t_{j}$, they form one shock
along with two rarefaction zones. Therefore, at $t=t_{j}$, the middle
rarefaction area has totally disappeared and we have
$x_{f_{L}}(t;v)=x_{f_{R}}(t;v)=x_{eq}$. That being said, after the slowest-
moving wave has reflected from its boundary, i.e. for $t_{ref_{2}}<t<t_{j}$,

$\displaystyle\Pi(\theta,x,t)=\begin{cases}\Pi_{a}=\mathcal{H}(\theta-
S_{r_{L}}),\quad S_{wi}<\theta<S_{L}^{\star\star},\quad a<x<x_{f_{L}}(t;v)\\\
\Pi_{b}=\mathcal{H}(\theta-
S_{r_{m}})\mathcal{H}(t_{j}-t)+0\cdot\mathcal{H}(t-t_{j}),\quad
S_{R}^{\star}<\theta<S_{L}^{\star},\quad\ x_{f_{L}}(t;v)<x<x_{f_{R}}(t;v)\\\
\Pi_{c}=\mathcal{H}(\theta-S_{r_{R}}),\qquad
S_{R}^{\star\star}<\theta<1-S_{oi},\qquad x_{f_{R}}(t;v)<x<b\end{cases}$ (30)

The ensemble average of $\Pi$ over all possible values of $S(x,t)$ yields
$F_{s}(S_{w};x,t)$ which is a deterministic quantity. Using the single-point
PDF of water saturation, $p_{s}(S_{w};x,t)$ we have,

$\displaystyle
F_{s}(\theta;x,t)=\big{\langle}\Pi(\theta,S_{w}(x,t))\big{\rangle}=\int_{-\infty}^{\infty}\mathcal{H}(\theta-
S_{w}^{\prime})p_{s}(S_{w}^{\prime};x,t)dS_{w}^{\prime}$ (31)

Therefore, by solving a linear equation for $\pi$, followed by finding only
the first ensemble moment of $\pi$, we find spatio-temporal evolution of the
full single-point CDF of saturation. This is a more straightforward approach
than first numerically solving the nonlinear Eq. (2) for an exhaustive number
of Monte Carlo realizations of the random variables, and subsequently
discovering the numerical CDF and PDF using methods like kernel density
estimation ([2]).

It is noteworthy that similarly to the past literature ([37], [14]), we prefer
employing a CDF framework rather than working with PDF equations, due to the
smoothness of CDF as well as the trivial formulation of boundary conditions
for CDF methodology, i.e. $F_{s}(\theta=S_{wi};x,t)=0$ and
$F_{s}(\theta=1-S_{oi};x,t)=1$. The boundary conditions can not be uniquely
defined for the PDF scheme. Therefore, using CDF obviates the need for dealing
with closure approximations inherent to the PDF scheme.

The method of characteristics can be utilized to form an analytical solution
for the one-dimensional linear hyperbolic Eq. (3), specifically in the
continuous rarefaction zone where there is no discontinuity present. While the
analytical solution for a horizontal reservoir has been studied by [37], the
more general case studies comprising of buoyancy forces demand manipulations
for handling two shocks introduced in this work. The raw CDF solution for the
horizontal reservoirs reads [37],

$\displaystyle\Pi_{b}=\mathcal{H}(\theta-
S_{wi})\mathcal{H}(x-C)+\mathcal{H}(\theta-1+S_{oi})\mathcal{H}(C-x)$ (32)

In order to define a family of characteristics $x=x(t;x_{0})$ along which the
original hyperbolic PDE becomes an ODE $\dfrac{d\pi}{dt}=0$, we start by using
equation $\dfrac{dx}{dt}=\textbf{v}(\theta,x(t),t),\qquad x(t=0)=x_{0}$, where
$x_{0}$ specifies where the characteristic line has originated from. Thus
$x=\bigintssss_{0}^{t}\textbf{v}(\theta,x(t^{\prime}),t^{\prime})dt^{\prime}+x_{0}$,
For the aforementioned setting of our problem, i.e. the initial and boundary
conditions of the gravity column, the characteristic solution is depicted in
Fig (3). The non-unity viscosity ratio causes an asymmetry and non-synchronism
in the movement of two waves. To this end, the solution will be discussed in
four distinct regions, i.e. for the right-moving branch before being reflected
from right boundary and after being reflected, and Similarly, for the left-
moving branch before and after being reflected from left boundary.

As Fig. (3) displays, before the faster wave (left-moving wave in
$(S_{D},S^{\star}_{L})$) hits into its boundary ($t<t_{{ref}_{1}}$), all the
characteristics are emanating from $z$ axis where $t=0$. Therefore, initial
condition is used to find the solution. Whereas, at $t=t_{{ref}_{1}}$, the
left-moving wave gets reflected from left boundary and the characteristics
start originating from t axis where $z=0$, while the direction of
characteristics has been reversed by the virtue of the reflection from
boundary. Thus, after reflection of the left-moving wave, and before the two
waves meet each other at time $t_{j}$ and stabilize at $t_{eq}$ somewhere in
the middle of domain identified by $z_{eq}$, we use the boundary condition at
$z=a$ to find the solution. The analysis of the right-moving wave
$(S^{\star}_{R},S_{D})$ follows a similar trend. Before hitting the right
boundary $t<t_{{ref}_{2}}$, characteristics originate from $z$ axis, and
initial condition is leveraged to find the solution. However, during
$t_{{ref}_{2}}<t<t_{j}$, characteristics develop from $t$ axis at the right
boundary $(z=L)$ and continues up to $t_{j}$ and location $z_{eq}$, where it
meets the other branch coming towards it.

$\displaystyle\Pi_{b}=\begin{cases}\Pi_{b_{L}}=\begin{cases}\mathcal{H}(\theta-
S_{B})\mathcal{H}(C-x)\mathcal{H}(x_{d}-x)+\mathcal{H}(\theta-
S_{D})\mathcal{H}(x-C)\mathcal{H}(x_{d}-x),\quad 0<t<t_{{ref}_{1}}\\\ \\\
\mathcal{H}(\theta-
S_{wi})\mathcal{H}(C-x)\mathcal{H}(x_{d}-x)+\mathcal{H}(\theta-S^{\star\star}_{L})\mathcal{H}(x-C)\mathcal{H}(x_{d}-x),\quad
t_{{ref}_{1}}<t<t_{j}\end{cases}\\\ \\\
\Pi_{b_{R}}=\begin{cases}\mathcal{H}(\theta-
S_{wi})\mathcal{H}(x-C)\mathcal{H}(x-x_{d})+\mathcal{H}(\theta-
S_{D})\mathcal{H}(C-x)\mathcal{H}(x-x_{d}),\quad 0<t<t_{{ref}_{2}}\\\ \\\
\mathcal{H}(\theta-
S_{B})\mathcal{H}(x-C)\mathcal{H}(x-x_{d})+\mathcal{H}(\theta-S^{\star\star}_{R})\mathcal{H}(C-x)\mathcal{H}(x-x_{d}),\quad
t_{{ref}_{2}}<t<t_{j}\end{cases}\end{cases}$ (33)

Where $\Pi_{b_{R}}$ and $\Pi_{b_{L}}$ manifest development of characteristics
in the right and left branches, respectively. Also, we use
$\Pi(t=0)=\mathcal{H}(\theta-S_{wi}),\Pi(x=0)=\mathcal{H}(\theta-
S_{B})=\mathcal{H}(\theta-1+S_{oi}),\Pi(x=x_{d})=\mathcal{H}(\theta-S_{D})$.
We define
$C(\theta,x,t)=\bigintss_{0}^{t}v(\theta,x(t^{\prime}),t^{\prime})dt^{\prime}=\dfrac{1}{\phi(x)}\dfrac{\partial
f(\theta)}{\partial\theta}\bigintss_{0}^{t}q(t^{\prime})dt^{\prime}=\dfrac{\partial
f(\theta)}{\partial\theta}\bigintss_{0}^{t}v(x,t^{\prime})dt^{\prime}$.

We will present the general formulation of CDF of water saturation, with
uncertainty in both porosity field and injection flux encapsulated in
$U(x,t)=\dfrac{q(t)}{\phi(x)}$. For the cases with one shock, i.e. a
horizontal reservoir as well as an inclined reservoir with updip flooding, the
CDF of saturation is then formulated as below ([37]),

$\displaystyle
F_{s}(\theta;x,t)=\begin{cases}\bigintss_{0}^{\infty}\Pi_{a}\mathcal{H}(x-x_{f})P_{U}(u)du,\hfill
S_{wi}\leq\theta<S^{\star}\\\
\bigintss_{0}^{\infty}\Pi_{a}\mathcal{H}(x-x_{f})P_{U}(u)du+\Pi_{b}\mathcal{H}(x_{f}-x)P_{U}(u)du,\quad\hfill
S^{\star}\leq\theta\leq 1-S_{oi}\end{cases}$ (34)

We introduce $F_{U(x,t)}$, CDF of the input random velocity field $U(x,t)$,
and employ the notation $a\vee b$ and $a\wedge b$ for the $max(a,b)$ and
$min(a,b)$ respectively. We expand this relation using a similar approach in
[14] as follows,

$\displaystyle F_{s_{a}}(\theta;x,t)$
$\displaystyle=\bigintsss_{0}^{\infty}\Pi_{a}\mathcal{H}(x-x_{f})P_{U}(u)du=\mathcal{H}(\theta-
S_{wi})\bigintsss_{0}^{\infty}\mathcal{H}\bigg{(}x-f^{\prime}(S^{\star})U(x,t)\bigg{)}P_{U}(u)du$
$\displaystyle=\mathcal{H}(\theta-
S_{wi})\bigintsss_{0}^{\infty}\mathcal{H}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star})}-u\bigg{)}P_{U}(u)du=\mathcal{H}(\theta-
S_{wi})F_{U(t)}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star})}\bigg{)}$ (35)

$\displaystyle
F_{s_{b}}(\theta;x,t)=F_{s_{a}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Pi_{b}\mathcal{H}(x_{f}-x)P_{U}(u)du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
S_{wi})\mathcal{H}(x-C)\mathcal{H}(x_{f}-x)+\mathcal{H}(\theta-
S_{B})\mathcal{H}(C-x)\mathcal{H}(x_{f}-x)\Bigg{]}P_{U}(u)du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
s_{wi})\mathcal{H}\bigg{(}x-f^{\prime}(\theta)U(x,t)\bigg{)}\mathcal{H}\bigg{(}f^{\prime}(S^{\star})U(x,t)-x\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\mathcal{H}(\theta-
S_{B})\mathcal{H}\bigg{(}f^{\prime}(\theta)U(x,t)-x\bigg{)}\mathcal{H}\bigg{(}f^{\prime}(S^{\star})U(x,t)-x\bigg{)}\Bigg{]}P_{U}du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
S_{wi})\mathcal{H}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}-u\bigg{)}\mathcal{H}\bigg{(}u-\dfrac{x}{f^{\prime}(S^{\star})}\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\mathcal{H}(\theta-
S_{B})\mathcal{H}\bigg{(}u-\dfrac{x}{f^{\prime}(\theta)}\bigg{)}\mathcal{H}\bigg{(}u-\dfrac{x}{f^{\prime}(S^{\star})}\bigg{)}\bigg{]}P_{U}du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\mathcal{H}(\theta-
S_{wi})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star})}\bigg{)}\Bigg{)}+\mathcal{H}(\theta-
S_{B})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S^{\star})}\bigg{)}\Bigg{)}$
(36)

For the downdip case as well as the pure gravity segregation case before
reflection of the waves from boundaries, the CDF equation reads,

$\displaystyle
F_{s}(\theta;x,t)=\begin{cases}\bigintss_{0}^{\infty}\Pi_{a}\mathcal{H}(x-x_{f_{R}})P_{U}(u)du,\hfill
S_{wi}\leq\theta<S_{R}^{\star}\\\
\bigintss_{0}^{\infty}\big{[}\Pi_{a}\mathcal{H}(x-x_{f_{R}})+\Pi_{b_{R}}\mathcal{H}(x_{f_{R}}-x)\big{]}P_{U}(u)du,\hfill
S_{R}^{\star}\leq\theta\leq S_{D}\\\
\bigintss_{0}^{\infty}\big{[}\Pi_{b_{L}}\mathcal{H}(x-x_{f_{L}})\big{]}P_{U}(u)du,\hfill
S_{D}\leq\theta\leq S_{L}^{\star}\\\
\bigintss_{0}^{\infty}\big{[}\Pi_{b_{L}}\mathcal{H}(x-x_{f_{L}})+\Pi_{c}\mathcal{H}(x_{f_{L}}-x)\big{]}P_{U}(u)du,\quad\hfill
S_{L}^{\star}<\theta\leq 1-S_{oi}\end{cases}$ (37)

Leveraging $F_{U(x,t)}$, this relation expands as follows. For $x>x_{D}$,
$S_{wi}<\theta<S^{\star}_{R}$ we have,

$\displaystyle F_{s_{a}}(\theta;x,t)$
$\displaystyle=\bigintsss_{0}^{\infty}\Pi_{a}\mathcal{H}(x-x_{f_{R}})P_{U}(u)du=\mathcal{H}(\theta-
S_{wi})\bigintsss_{0}^{\infty}\mathcal{H}\bigg{(}x-f^{\prime}(S^{\star}_{R})U(x,t)\bigg{)}P_{U}(u)du$
$\displaystyle=\mathcal{H}(\theta-
s_{wi})\bigintsss_{0}^{\infty}\mathcal{H}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star}_{R})}-u\bigg{)}P_{U}(u)du=\mathcal{H}(\theta-
S_{wi})F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star}_{R})}\bigg{)}$ (38)

For $x>x_{D}$, $S^{\star}_{R}<\theta<S_{D}$ we have,

$\displaystyle
F_{s_{b_{R}}}(\theta;x,t)=F_{s_{a}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Pi_{b}\mathcal{H}(x_{f_{R}}-x)P_{U}(u)du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
S_{wi})\mathcal{H}(x-C)\mathcal{H}(x_{f_{R}}-x)+\mathcal{H}(\theta-
S_{D})\mathcal{H}(C-x)\mathcal{H}(x_{f_{R}}-x)\Bigg{]}P_{U}(u)du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\mathcal{H}(\theta-
S_{wi})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star}_{R})}\bigg{)}\Bigg{)}+\mathcal{H}(\theta-
S_{D})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S^{\star}_{R})}\bigg{)}\Bigg{)}$
(39)

And for $x<x_{D}$, $S_{D}<\theta<S^{\star}_{L}$ we use,

$\displaystyle F_{s_{b_{L}}}(\theta;x,t)$
$\displaystyle=\bigintsss_{0}^{\infty}\Pi_{b}\mathcal{H}(x-x_{f_{L}})P_{U}(u)du$
$\displaystyle=\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
S_{B})\mathcal{H}(C-x)\mathcal{H}(x-x_{f_{L}})+\mathcal{H}(\theta-
S_{D})\mathcal{H}(x-C)\mathcal{H}(x-x_{f_{L}})\Bigg{]}P_{U}(u)du$
$\displaystyle=\mathcal{H}(\theta-
S_{B})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star}_{L})}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}\Bigg{)}+\mathcal{H}(\theta-
S_{D})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S^{\star}_{L})}\bigg{)}\Bigg{)}$
(40)

for $x<x_{D}$, $S^{\star}_{L}<\theta<S_{B}$ we leverage,

$\displaystyle F_{s_{c}}(\theta;x,t)$
$\displaystyle=F_{s_{b_{L}}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Pi_{c}\mathcal{H}(x_{f_{L}}-x)P_{U}(u)du$
$\displaystyle=F_{s_{b_{L}}}(\theta;x,t)+\mathcal{H}(S_{B}-\theta)\bigintsss_{0}^{\infty}\mathcal{H}\bigg{(}f^{\prime}(S^{\star}_{L})U(x,t)-x\bigg{)}P_{U}(u)du$
$\displaystyle=F_{s_{b_{L}}}(\theta;x,t)+\mathcal{H}(S_{B}-\theta)\bigintsss_{0}^{\infty}\mathcal{H}\bigg{(}u-\dfrac{x}{f^{\prime}(S^{\star}_{L})}\bigg{)}P_{U}(u)du$
$\displaystyle=F_{s_{b_{L}}}(\theta;x,t)+\mathcal{H}(S_{B}-\theta)\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S^{\star}_{L})}\bigg{)}\Bigg{)}$
(41)

Whereas the gravity segregation scenario after the reflection of waves follows
the CDF relation below,

$\displaystyle
F_{s}(\theta;x,t)=\begin{cases}\bigintss_{0}^{\infty}\Pi_{b_{L}}\mathcal{H}(x_{f_{L}}-x)P_{U}(u)du,\quad
S_{wi}\leq\theta<S_{L}^{\star\star}\\\
\bigintss_{0}^{\infty}\big{[}\Pi_{b_{L}}\mathcal{H}(x_{f_{L}}-x)+\Pi_{b}\mathcal{H}(x_{f_{R}}-x)\mathcal{H}(x-x_{f_{L}})\big{]}P_{U}(u)du,\quad
S_{R}^{\star}\leq\theta\leq S_{L}^{\star}\\\
\bigintss_{0}^{\infty}\big{[}\Pi_{b_{L}}\mathcal{H}(x_{f_{L}}-x)+\Pi_{b}\mathcal{H}(x_{f_{R}}-x)\mathcal{H}(x-x_{f_{L}})+\Pi_{b_{R}}\mathcal{H}(x-x_{f_{R}})\big{]}P_{U}(u)du,\\\
\hfill S_{R}^{\star\star}<\theta\leq 1-S_{oi}\end{cases}$ (42)

This relation is consequently expanded as below. For
$x<x_{S_{L}^{\star\star}}$, $S_{wi}<\theta<S^{\star\star}_{L}$ we have,

$\displaystyle
F_{s_{a}}(\theta;x,t)=\bigintsss_{0}^{\infty}\Pi_{b_{L}}\mathcal{H}(x_{f_{L}}-x)P_{U}(u)du$
$\displaystyle=\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
S_{L}^{\star\star})\mathcal{H}(x-C)\mathcal{H}(x_{f_{L}}-x)+\mathcal{H}(\theta-
S_{wi})\mathcal{H}(C-x)\mathcal{H}(x_{f_{L}}-x)\Bigg{]}P_{U}(u)du$
$\displaystyle=\mathcal{H}(\theta-
S_{L}^{\star\star})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S_{L}^{\star\star})}\bigg{)}\Bigg{)}+\mathcal{H}(\theta-
S_{wi})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S_{L}^{\star\star})}\bigg{)}\Bigg{)}$
(43)

Analogously to the previous derivations, for
$x_{S_{L}^{\star\star}}<x<x_{S_{R}^{\star\star}}$,
$S^{\star}_{R}<\theta<S^{\star}_{L}$ we have,

$\displaystyle F_{s_{b_{R}}}(\theta;x,t)$
$\displaystyle=F_{s_{a}}(\theta;x,t)+\mathcal{H}(\theta-
S_{wi})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S_{3})}\bigg{)}\Bigg{)}$
$\displaystyle+\mathcal{H}(\theta-
S_{D})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S_{3})}\bigg{)}\Bigg{)}$
(44) $\displaystyle F_{s_{b_{L}}}(\theta;x,t)$
$\displaystyle=F_{s_{a}}+F_{s_{b_{R}}}+\mathcal{H}(\theta-
S_{B})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S_{4})}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}\Bigg{)}$
$\displaystyle+\mathcal{H}(\theta-
S_{D})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S_{4})}\bigg{)}\Bigg{)}$
(45)

For $x>x_{S_{R}^{\star\star}}$, which corresponds to
$S^{\star\star}_{R}<\theta<1-S_{oi}$ we leverage,

$\displaystyle
F_{s_{c}}(\theta;x,t)=F_{s_{a}}(\theta;x,t)+F_{s_{b}}(\theta;x,t)+\bigintsss_{0}^{\infty}\Pi_{b_{R}}\mathcal{H}(x-x_{f_{R}})P_{U}(u)du=F_{s_{a}}(\theta;x,t)+F_{s_{b}}(\theta;x,t)$
$\displaystyle+\bigintsss_{0}^{\infty}\Bigg{[}\mathcal{H}(\theta-
S_{B})\mathcal{H}(x-C)\mathcal{H}(x-x_{f_{R}})+\mathcal{H}(\theta-
S_{R}^{\star\star})\mathcal{H}(C-x)\mathcal{H}(x-x_{f_{R}})\Bigg{]}P_{U}(u)du$
$\displaystyle=F_{s_{a}}(\theta;x,t)+F_{s_{b_{R}}}(\theta;x,t)+F_{s_{b_{L}}}(\theta;x,t)$
$\displaystyle+\mathcal{H}(\theta-
S_{B})\Bigg{(}1-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\vee\dfrac{x}{f^{\prime}(S_{R}^{\star\star})}\bigg{)}\Bigg{)}+\mathcal{H}(\theta-
S_{R}^{\star\star})\Bigg{(}F_{U}\bigg{(}\dfrac{x}{f^{\prime}(S_{R}^{\star\star})}\bigg{)}-F_{U}\bigg{(}\dfrac{x}{f^{\prime}(\theta)}\bigg{)}\Bigg{)}$
(46)

Where $F_{s_{b}}(\theta;x,t)$ should be found from the solutions of the before
reflection time in Eq. 3 and Eq.3. However, those solutions were developed for
$(S_{R}^{\star},S_{L}^{\star})$ rarefaction zone, whereas, the middle
rarefaction region here has a different saturation range $(S_{3},S_{4})$
(Fig.4), which keeps changing and becoming smaller over time, until it totally
fades away at $t_{j}$ where one shock forms in $z_{eq}$. That is, after
reflection from boundaries, the solution is no longer self-similar, and the
formula developed in Eq.3 and Eq.3 using the method of distributions, even
though following similar logic to the previous equations, fail to accurately
capture the $F_{s}(s;x,t)$ predicted from Monte Carlo scheme. Fig. 9
represents a reasonable prediction for CDF of saturation at $t_{1},t_{2}$
(before reflection), while at $t_{3},t_{4}$ (after reflection) the CDF of
saturation obtained from Eq.3 and Eq.3, don’t capture the solution of Monte
Carlo scheme.

Figure 4: Saturation profile for the inverted gravity column, before and after
being reflected from the boundaries, solid lines are the Monte Carlo results,
while dashed lines represent the corresponding deterministic solutions.

Note that CDF is a cumulative quantity, and hence we need to add the CDF of
all the previous saturation intervals, starting from the lowest saturation
range onward. Note that $F_{U(x,t)}(x,t)$ is only a function of time as
$F_{q(t)}(t)$ when the injection flux $q(t)$ is the sole source of
uncertainty, and is only a function of space $F_{\phi(x)}(x)$ when porosity
field $\phi(x)$ is the sole source of uncertainty. It is a space-time
dependent function for cases when both fields are random. We use Kernel
Density Estimation (KDE) to find $F_{U(x,t)}$.

We will validate and verify these assumptions by comparing the
$F_{s}(\theta;x,t)$ obtained from Eqs. (34)-(42) with their corresponding
$F_{s}(S,x,t)$ obtained from kernel density estimation applied on the Monte
Carlo simulation of Eq. (2) and Eq. (2). It’s worth mentioning that in general
$U(x,t)=\dfrac{q(t)}{\phi(x)}$, and hence when the injection flux is random,
we generate $F_{U}(t)(t)$ by applying KDE on 10000 samples of $U(t)=q(t)$
while keeping the porosity constant. On the other hand, when porosity is the
only random variable, we generate $F_{U}(x)(x)$ by applying KDE on
$U(x)=\dfrac{1}{\phi(x)}$ by generating 10000 samples of $\phi(x)$ while
keeping the injection flux constant. Similarly, when both fields are random,
10000 samples of each random variable is generated and then KDE is applied on
$U(x,t)=\dfrac{q(t)}{\phi(x)}$. In all three cases the bandwidth is selected
to be $n=10^{10}$ in order to guarantee the accuracy of matching with the
corresponding Monte Carlo solution of the finite volume Godunov scheme.

Figure 5: Spatio-temporal evolution of the mean saturation field represents
development of the shock(s), (top) horizontal flooding with one shock,
(middle) downdip flooding with two shocks, (bottom) gravity column with
reflection of two shocks from boundaries. Initial discontinuity for the
downdip and gravity column cases is located at $Z_{d}=0.3$ and all the plots
are scaled to the same spatial range for better comparison.

##  4 Numerical Experiments

###  4.1 Setting

For the Monte Carlo simulation of the nonlinear hyperbolic equation for
saturation Eq. (2) with complicated initial conditions, one must generally
resort to numerical schemes like finite element, difference, or volume
methods. Accordingly, we have employed finite volume discretization of Godunov
scheme to solve the dimensionless form of the nonlinear hyperbolic
conservation Eq. (2) in one-dimension. While previous studies by [13] , [14],
[15], have employed a Lagrangian streamline simulation methodology for
numerical treatment of the horizontal displacement, we have applied a finite
volume strategy to enable numerical manipulation of the multi-shock problems
like the downdip displacement in our study.

The CDF and PDF of water saturation $F_{s}(S_{w})$, and $f_{s}(S_{w})$ are
generated using kernel density estimation [2] with 10,000 samples of each
random variable, while the bandwidth parameter is taken to be $n=2^{10}$.
Intuitively, we would like to choose the smallest bandwidth parameter that
data will allow. After experimenting with the influence of this parameter on
the resulting CDF and PDF plots, we realized that values smaller than
$n=2^{7}-2^{12}$ result in an under-smoothed curve with too much spurious data
artifacts, while values greater than n yield an over-smoothed curve in which
much of the underlying structure is obscured. CDF plots usually need
$n=2^{10}$ to capture the transitions from shock to rarefaction zones
accurately, while such a large bandwidth leads to overfitting in the PDF
plots. Hence, we use $n=2^{8}$ to generate the PDF plots.

We will present the influence of parametric uncertainty on the saturation
uncertainty by inspecting the mean, variance, temporal evolution of CDF and
PDF of water saturation at a specific location, spatial evolution of CDF and
PDF of water saturation at a specific point in time and the point-wise average
of the CDFs generated at all locations along the domain. The aformentioned
profiles are all computed using both Monte Carlo sampling of the Buckley
Leverett equation, and the method of distributions. We will study distinct
cases for only $\phi(x)$ random, only $q(t)$ random, and when both are
stochastic fields. Random porosity field has been restricted to the range
$[0.3,0.6]$, and stochastic injection flux has established to vary between
$[0.3,0.8]$ as depicted in Fig.6. All cases have employed the Brooks and Corey
expressions for the relative permeabilities as a function of saturation. We
will demonstrate all the experiments for $m=0.5$, which was formerly studied
in test cases of [13], [14], [7]. The irreducible wetting and non-wetting
saturations are taken to be $S_{wi}=S_{oi}=0.1$. The injection, if present,
always happens at the inlet of domain $x=a$. Final simulation time is $T=1$
for all cases except in the gravity column simulation with $T=10$.

A convergence study at the end of results section will demonstrate the
sensitivity of solution to the spatial and temporal grid sizes, as well as
number of Monte Carlo trials used to find the results. Also, we will quantify
the error between the CDF obtained from Monte Carlo sampling of Godunov
equations and the method of distribution. Given the CDF of water saturation,
the first two statistical moments of saturation read,

$\displaystyle\langle
S_{w}\rangle=1-\bigintssss_{0}^{1}F_{S_{w}}(\zeta)d\zeta$ (47)
$\displaystyle\sigma^{2}_{S_{w}}=1-\bigintssss_{0}^{1}F_{S_{w}}(\sqrt{\zeta})d\zeta-\langle
S_{w}\rangle^{2}$ (48)

Fig. (5) represents spatio-temporal evolution of the average saturation fields
and development of the shocks for each case study using the Monte Carlo scheme
(in all cases the spatial domain is scaled to (0,1) range for easier
comparison).

Figure 6: Multiple realizations of the random porosity field (left) and random
injection flux (right) with the parameters defined in section 4.1.3

####  4.1.1 Numerical setting for stochastic porosity field $\phi(x)$

In the current scenario, we consider the injection flux $q_{T}(t)$ at the
boundary to be a deterministic constant as $q_{inlet}=0.3$. Hence, the
uncertainty in the porosity field is the only source of randomness. The
porosity field $\phi(x)$ is characterized by its know mean, variance,
covariance structure, and dimensionless correlation length.

Similarly to [13], in order to enforce the physical constraint $0<\phi(x)<1$,
we first define a random Gaussian field
$\gamma(x)\sim\mathcal{N}(\mu_{\gamma},C_{\gamma}(x))$, with mean $\mu_{x}$
and exponential covariance structure
$C_{\gamma}(x)=\sigma_{\gamma}^{2}exp(-\dfrac{x}{\lambda_{\gamma}})$, where
$\sigma_{\gamma}^{2}$ is the variance of $\gamma$, and $\lambda_{\gamma}$
represents the correlation length. Subsequently, we define porosity field as,

$\displaystyle\phi(x)=0.01\dfrac{2}{\pi}\lvert\text{arctan}(\gamma(x))|+0.5\bigg{(}1-\dfrac{2}{\pi}|\text{arctan}(\gamma(x))|\bigg{)}$
(49)

In consequence, we can choose $\mu_{\gamma}$ so that
$\langle\phi(x)\rangle=\mu_{\phi}$.

For a space discretization $\Delta x$ along the domain $[a,b]$, while defining
the spatial distance $x=(n-1)\Delta x$, we define the Gaussian vector
$(\gamma)^{n}=(\gamma(j\Delta x))_{0\leq j\leq n-1}$ by
$(\gamma)^{n}\sim\mathcal{N}(\mu_{\gamma}\mathds{1}_{n},(C_{\gamma})^{n})$,
where $(C_{\gamma})^{n}=(C_{\gamma}(|i-j|\Delta x))_{1\leq i,j\leq n}$, and
$\mathds{1}_{n}$ is a unit vector of size n.

Throughout our experiments with a stochastic porosity field, unless otherwise
stated, we use $\mu_{\phi}=0.3$, $\sigma_{\gamma}^{2}=0.5$, and
$\lambda_{\gamma}=0.5L$ where $L=1$ and $L=1$ are the dimensionless length of
the domain for horizontal and vertical cases, respectively. Nevertheless, in
the lower order approximation (SME) section, we keep modifying these
parameters and monitoring how the resulting first and second moments of water
saturation from Monte Carlo and CDF methods match with those resulting from
solving statistical moment equations for the different $\sigma_{\phi}^{2}$ and
$\lambda_{\phi}$ variables. Figures 7, 8, 9 show the results for random
porosity field as the only source of uncertainty. The mean saturation profiles
are plotted along with the deterministic solutions, and the smoother solutions
of MC and MD methods compared to the sharp solutions of the deterministic case
represent the results of averaged solution of a non-zero variance for the
underlying stochastic field. Note that the vertical shock regions in the
saturation profiles correspond to the horizontal regions in the CDF plots. As
these plots represent, our proposed method of distributions is able to predict
the behavior of Monte Carlo for mean, standard deviation, and CDF of
saturation. The only case for which our method shows limitation is in the
gravity column case and only after reflection of waves from boundaries (Fig.
9, third row). Unlike to the other cases, this mismatch occurs due to the non-
self-similarity of solution in the output state (saturation) domain. That is,
as we showed in Eq. 3 and Eq. 3, we have two non-constant specific saturation
values $S_{3}$ and $S_{4}$ (Fig. 4) needed to compute $F_{s_{b_{R}}}$ and
$F_{s_{b_{L}}}$ respectively, and in general it’s hard to keep track of these
constantly-changing saturation values. Consequently, solution of the CDF
equations is over-complicated and in general does not follow Eq. 3 and Eq. 3.
Therefore, even though the proposed method perfectly captures the two-wave
behavior in all the previous examples, it fails only for this specific case.

Figure 7: (First row) The first two moments of saturation, (second row)
spatial and temporal evolution of the CDF of saturation (dashed lines show the
results of method of distributions, while solid lines correspond to Monte
Carlo), (third row) point-wise average of the CDF and PDF at all points along
the domain, obtained from Monte Carlo and CDF methods, (fourth row) the first
two moments obtained using Low-order SME method, vs those of MC and MD. All
plots correspond to the horizontal flooding with $\phi(x)$ random,
$\mu_{\phi}=0.3$, $\sigma^{2}_{\phi}=0.5$, $\lambda_{\phi}=0.5L$, $q=0.3$, at
two dimensionless times $t_{1}=0.2$, $t_{2}=0.3$. CDF and PDF from Monte Carlo
are plotted using a bandwidth of $n=2^{12}$ and $n=2^{7}$ respectively for the
KDE post-processing.

Figure 8: (First row) The first two moments of saturation, (second row)
spatial and temporal evolution of the CDF of saturation (dashed lines show the
results of method of distributions, while solid lines correspond to Monte
Carlo), (third row) point-wise average of the CDF and PDF at all points along
the domain, obtained from Monte Carlo and CDF methods. All plots correspond to
the downdip flooding with $\phi(x)$ random, $\mu_{\phi}=0.3$,
$\sigma^{2}_{\phi}=0.5$, $\lambda_{\phi}=0.5L$, $q=0.3$, at two dimensionless
times $t_{1}=0.2$, $t_{2}=0.3$. CDF and PDF from Monte Carlo are plotted using
a bandwidth of $n=2^{12}$ and $n=2^{7}$ respectively for the KDE post-
processing.

Figure 9: (First row) The first two moments of saturation (solid lines are the
result of MC simulation, while markers show the results of method of
distributions), (second/third row) spatial evolution of the CDF of saturation
at four fixed dimensionless times $t_{1}=0.6,t_{2}=1.2,t_{3}=4.5,t_{4}=7.5$
and four dimensionless spatial positions along the domain, temporal evolution
of the CDF of saturation at a fixed dimensionless spatial positions along the
domain (dashed lines show the results of method of distributions, while solid
lines correspond to Monte Carlo), (fourth row) temporal evolution of the CDF
of saturation at a fixed dimensionless position $x=0.7$ and two dimensionless
times, as well as point-wise spatial average of CDF of saturation along the
domain at four dimensionless times. All simulations represent an inverted
gravity column with $\phi(x)$ random, $\mu_{\phi}=0.3$,
$\sigma^{2}_{\phi}=0.5$, $\lambda_{\phi}=0.6L$, $L=3$, $q=0.3$. Except in the
saturation profile, solid lines show the results of MC, while dashed lines are
resulted from CDF method.

####  4.1.2 Numerical setting for stochastic injection flux $q(t)$

Enforcing the general continuity constraint for incompressible flow,
$\nabla\cdot q_{T}(x,t)=0$ makes the total flux $q_{T}$ to be constant in x,
and is therefore equal to the injection flow rate at the left boundary of
domain $q(x,t)=q(t)=q_{0}(t)$. We are interested in finding the solutions to
Eq. (2) while treating $q_{0}(t)$ as an stochastic field in time with a
prescribed CDF. Henceforward, we choose to keep porosity as a deterministic
constant $\phi=0.3$ for this scenario, thereby only experimenting with the
effect of randomness in the injection flux. In order to find the distribution
for the random injection flux at the inlet, similarly to the previous case, we
assume the mean $\mu_{q}$, variance $\sigma^{2}_{q}$, covariance structure
$C_{q}(t)$, and correlation time $\tau_{q}$ are all known parameters. To this
end, similarly to [13], we first define a Gaussian random field
$\rho_{q}(t)\sim\mathcal{N}(\mu_{\rho},C_{\rho}(t))$, with an exponential
covariance structure $C_{\rho}(t)=\sigma^{2}_{\rho}exp(-\dfrac{t}{\tau_{q}})$,
where $\tau_{q}$ is the correlation parameter for events that are closer in
time. Eventually, we define the stochastic injection flux field as,

$\displaystyle q_{T}(t)=0.05+0.1(\rho_{q}(t))^{2}$ (50)

Such a definition allows us to first make sure that $q_{T}(t)>0$ and second by
controlling $\mu_{\rho}$ we can satisfy the criteria $\langle
q(t)\rangle=\mu_{q}$. Any other arbitrary non-negative distribution could be
used alternatively.

Having a time discretization of size $\Delta t$, for any point $t$ in the time
interval $[0,T]$, where $t=(n-1)\Delta t$, we define the Gaussian vector
$(\rho_{q})^{n}=(\rho_{q}(j\Delta t))_{0\leq j\leq n-1}$ by
$(\rho_{q})^{n}\sim\mathcal{N}(\mu_{\rho}\mathds{1}_{n},(C_{\rho})^{n})$,
where $C_{\rho}^{n}=(C_{q}(|i-j|\Delta t))_{1\leq i,j\leq n}$ and
$\mathds{1}_{n}$ is the vector of all ones which is of size n. Unless
otherwise stated, whenever experimenting with uncertainty in $q(t)$, we will
use $\mu_{q}=0.3$, $\sigma^{2}_{q}=0.5$, $\tau_{q}=0.5T$, $T=1$, $m=0.5$ .
Figures 10 and 11 shows the results for the injection flux as the sole source
of randomness.

Figure 10: (First row) The first two moments of saturation, (second row)
spatial and temporal evolution of the CDF of saturation (dashed lines show the
results of method of distributions, while solid lines correspond to Monte
Carlo), (third row) point-wise spatial average of the CDF and PDF at all
points along the domain, obtained from Monte Carlo and CDF methods, (fourth
row) the first two moments obtained using Low-order SME method, vs. those of
MC and MD. All simulations represent horizontal flooding, with $q(t)$ random,
$\mu_{q}=0.3$, $\sigma^{2}_{q}=0.25$, $\tau_{q}=0.1T$, $T=1$, $\phi=0.3$, at
two dimensionless times $t_{1}=0.2$, $t_{2}=0.3$. CDF and PDF from Monte Carlo
are plotted using a bandwidth of $n=2^{12}$ and $n=2^{7}$ respectively for the
KDE post-processing.

Figure 11: (First row) The first two moments of saturation, (second row)
spatial and temporal evolution of the CDF of saturation (dashed lines show the
results of method of distributions, while solid lines correspond to Monte
Carlo), (third row) point-wise spatial average of the CDF and PDF at all
points along the domain, obtained from Monte Carlo and CDF methods. All plots
correspond to the downdip flooding with $q(t)$ random, $\mu_{q}=0.3$,
$\sigma_{q}=0.5$, $\tau_{q}=0.5T$, at two dimensionless times $t_{1}=0.2$,
$t_{2}=0.3$.

####  4.1.3 Numerical setting for uncertainty in both $\phi(x)$ and $q(t)$

In this section, we will illustrate the joint effects of uncertainty in the
porosity field as well as the injection flux. We will resort to the same setup
as the two previous sections for the construction of $\phi(x)$ and $q_{T}(t)$
fields. As for the numerical parameters, we will employ $m=0.5$,
$\mu_{q}=0.3$, $\sigma^{2}_{q}=0.5$, $\tau_{q}=0.5T$, $T=1$, $\mu_{\phi}=0.3$,
$\sigma^{2}_{\phi}=0.5$, $\lambda_{\phi}=0.5L$. Fig. 6 illustrates four random
sampling of porosity field as well as injection flux field plotted using the
aforementioned values for mean, variance, and correlation length/time for each
random variable. In this section, we will illustrate the joint effects of
uncertainty in the porosity field as well as the injection flux. We will
encapsulate both uncertainties in the structure of the. stochastic velocity
field by representing it as $\textbf{v}_{T}=\textbf{v}_{T}(x,t)$. Fig. 12
shows the results of having both fields as the source of randomness

Figure 12: (First row) The first two moments of saturation, (second row)
spatial and temporal evolution of the CDF of saturation (dashed lines show the
results of method of distributions, while solid lines correspond to Monte
Carlo), (third row) point-wise average of the CDF and PDF at all points along
the domain, obtained from Monte Carlo and CDF methods, (fourth row) the first
two moments obtained using Low-order SME method, vs. those of MC and MD. All
simulations represent horizontal case with $\phi(x)$ random, $\mu_{\phi}=0.3$,
$\sigma^{2}_{\phi}=0.5$, $\lambda_{\phi}=0.5L$, $L=1$, and $q(t)$ random,
$\mu_{q}=0.3$, $\sigma^{2}_{q}=0.5$, $\tau_{q}=0.5T$, $T=1$. Solid lines are
the result of MC simulation, while dashed lines represent the results from CDF
method.

##  5 Low-Order Approximation (LOA)

So far we have demonstrated that estimations of the mean and standard
deviation, as well as distribution functions of saturation field from CDF
method exhibit a promising match with those obtained from Monte Carlo averaged
of a large number of realizations. That is, given an exhaustive number of
trials for the random porosity or injection flux followed by solving an
upwinding-type problem leads to the solution for saturation field which
matches well with the solution from method of distributions.

The first two moments of water saturation from Monte Carlo and method of
distributions were demonstrated to match well. More specifically, in Monte
Carlo and MD methods, the saturation moments are obtained through a post-
processing step, i.e. after finding the saturation field itself.
Alternatively, statistical moment equations (SME) directly solve equations for
the first two statistical moments of water saturation, i.e. mean and standard
deviation. In this section, we investigate the impact of variance and
correlation length of the uncertain input parameters on the first two
statistical moments of saturation field. This methodology also referred to as
low-order approximation (LOA) involves mathematical manipulations started with
the perturbation expansion of the stochastic parameters [7] and [13].

###  5.1 Formulation

We will present the effect of random input on the first two moments of water
saturation for different levels of standard deviation and correlation length
of the input stochastic field. First we assume that injection flux $q_{T}$ is
a random constant, while porosity $\phi$ is a deterministic constant. In the
second scenario, we will assume the opposite setting holds, i.e. we use the
equations resulting from perturbation expansion of the porosity field
$\phi(x)$, while $q_{T}$ is a deterministic constant.

In the first scenario, we assume a deterministic constant value of
$\phi=\phi_{0}=0.3$, and a lognormal distribution for $q$, with $\mu_{q}=0.3$,
while we will experiment with multiple values of variance $\sigma_{q}^{2}$ and
track the resulting first two moments of the state field $S_{w}$,

$\displaystyle q\sim
log\mathcal{N}\bigg{(}\langle\text{ln}(q)\rangle,\sigma_{\text{ln}(q)}\bigg{)}\qquad\langle
ln(q)\rangle=ln\bigg{(}\dfrac{\mu_{q}^{2}}{\sqrt{\sigma^{2}_{q}+\mu_{q}^{2}}}\bigg{)},\qquad\sigma_{ln(q)}=\sqrt{ln\bigg{(}\dfrac{\sigma_{q}^{2}}{\mu_{q}^{2}}+1\bigg{)}}$
(51)

Similarly to [7], we consider the one-dimensional version of Eq. (2), and
using the idea of displacement along a streamline ([7]), we define
$\tau(x;x_{0})$, the travel time (time of flight) of a particle from $x_{0}=0$
to $x$ in the total velocity field $\textbf{v}_{T}$ as,

$\displaystyle\dfrac{d\tau}{dx}=\dfrac{1}{\textbf{v}_{T}}=\dfrac{\phi_{0}}{q}$
(52)

Which has a solution in the form of
$\tau(x;x_{0})=(x-x_{0})\dfrac{\phi_{0}}{q}$. Since $\textbf{v}_{T}$ is
random, so is $\tau$, and hence, its statistical moments depend on those of
$\textbf{v}_{T}$. By applying a Reynolds decomposition on $q$, we can estimate
the random variable using its ensemble mean (expected value) plus the zero-
mean random perturbations (around the mean). Thereupon, Eq. (52) will be
recast in the following form,

$\displaystyle\tau(x;x_{0})=(x-x_{0})\dfrac{\phi_{0}}{(\mu_{q}+q^{\prime})}=(x-x_{0})\dfrac{\phi_{0}}{\mu_{q}}\bigg{[}1-\dfrac{q^{\prime}}{\mu_{q}}+\bigg{(}\dfrac{q^{\prime}}{\mu_{q}}\bigg{)}^{2}+\cdots\bigg{]}$
(53)

Correspondingly, LOA results in the first-order approximation of the first two
moments of time of flight as a function of $\mu_{q}$ and $\sigma^{2}_{q}$,

$\displaystyle\langle\tau(x;x_{0})\rangle\simeq(x-x_{0})\dfrac{\phi_{0}}{\mu_{q}},\qquad\qquad\qquad\sigma^{2}_{\tau}\simeq(x-x_{0})^{2}\phi_{0}^{2}\dfrac{\sigma_{q}^{2}}{\mu_{q}^{4}}$
(54)

Furthermore, by prescribing a probability distribution to $\tau$ for the given
mean and standard deviation, e.g. lognormal distribution $\tau\sim
log\mathcal{N}\bigg{(}\langle\text{ln}(\tau)\rangle,\sigma_{\text{ln}(\tau)}\bigg{)}$,
where the corresponding moments are defined as ([7]),

$\displaystyle\big{\langle}\text{ln}(\tau)\big{\rangle}=2\text{ln}(\langle\tau\rangle)-\dfrac{1}{2}\text{ln}\bigg{[}\langle\tau\rangle^{2}+\sigma^{2}_{\tau}\bigg{]},\qquad\qquad\sigma^{2}_{\text{ln}(\tau)}=\text{ln}\bigg{[}\langle\tau\rangle^{2}+\sigma^{2}_{\tau}\bigg{]}-2\text{ln}(\langle\tau\rangle)$
(55)

Therefore, distribution of $\tau$ is fully delineated as a lognormal
probability density function, along with the mean $\mu_{\tau}$ and standard
deviation $\sigma_{\tau}^{2}$, which are directly computable using $\mu_{q}$
and $\sigma_{q}^{2}$ in Eq. (54). Therefore, by controlling $\sigma_{q}^{2}$
parameter, we can find a low-order estimate of the ensemble mean (expected
value) and standard deviation of water saturation using the expressions below,

$\displaystyle\langle
S_{w}(x,t)\rangle=\bigintsss_{0}^{+\infty}\tilde{S}(\tau,t)p_{\tau}(\tau;x,x_{0})d\tau,$
$\displaystyle\sigma^{2}_{S}(x,t)=\bigintsss_{0}^{+\infty}\tilde{S}^{2}(\tau,t)p_{\tau}(\tau;x,x_{0})d\tau-\langle
S(x,t)\rangle^{2}$ (56)

Where where $p_{\tau}(\tau;x,x_{0})$ is the probability density function of
the (one-particle) travel time $\tau$, and $\tilde{S}(\tau,t)$ can be found by
combining Eq. (2) and Eq. (52) resulted in a transformed version of the
transport equation for water saturation as,

$\displaystyle\dfrac{\partial\tilde{S}}{\partial
t}+f_{w}^{\prime}(S)\dfrac{\partial\tilde{S}(\tau,t)}{d\tau}=0$ (57)

as described in [7]. It should be emphasized that even though [7] has studies
low-order approximation for a horizontal reservoir, we are indeed taking the
gravitational effects into account in the flux function while solving Eq.
(57).

In the second scenario, we assume a constant deterministic injection flux
$q=q_{0}=0.3$ and a random stationary Gaussian porosity field $\phi(x)$
characterized by its mean value of $\mu_{\phi}=0.3$, and covariance structure
$C_{\phi}(x)=\sigma^{2}_{\phi}exp(-\dfrac{|x|}{\lambda_{\phi}})$. Given such a
definition for the random porosity field, we can leverage the first-order
approximation of the first two moments of $\tau$ derived in [7] as following,

$\displaystyle\big{\langle}\tau(x;x_{0})\big{\rangle}\simeq(x-x_{0})\dfrac{\mu_{\phi}}{q_{0}},$
$\displaystyle\sigma^{2}_{\tau}(x;x_{0})\simeq\dfrac{2\sigma^{2}_{\phi}}{q_{0}^{2}}\bigg{[}\lambda_{\phi}(x-x_{0})-\lambda_{\phi}^{2}\bigg{(}1-exp\bigg{(}-\dfrac{|x-x_{0}|}{\lambda_{\phi}}\bigg{)}\bigg{)}\bigg{]}$
(58)

Furthermore, $\tau$ is assumed log-normally distributed and the same
definitions as Eq. (55) hold for $\langle\text{ln}(\tau)\rangle$ and
$\sigma^{2}_{\text{ln}(\tau)}$. To this end, by controlling $\mu_{\phi}$ and
$\lambda_{\phi}$, we will present comparisons of the low-order approximation
for the first two moments of water saturation using Eq. (5.1) against those
obtained from Monte Carlo and CDF methods.

###  5.2 Illustrative Examples and Numerical Results

While the the results from CDF method stay in agreement with those of MC up to
large variance and correlation lengths, the accuracy of LOA immediately starts
deteriorating even for small variances. That is because the statistical
moments of travel time and velocity are derived as first-order approximations.
Therefore, the first-order approximations require the variance of log
permeability to be (much) smaller than 1. As depicted in the plots showing
sensitivity to the increasing variance, by increasing the variance of input
parameters, the saturation variance is increasing, too. Consequently, the
saturation profile becomes smoother and deviates from a sharp deterministic
profile in the shock region. As the results of sensitivity to the correlation
length represents in Fig. 16, initially the LOA results are far from MC and
CDF results. As we increase the correlation length while keeping $\Delta x$
constant, the approximations for both mean and standard deviation become
closer to those of MC. This is because more grid blocks fall within one
correlation length, and therefore the underlying domain becomes more
homogeneous to resolve. Also, as Fig. 14 represents, for the cases where there
are two shocks, at early times LOA can not discern the two shocks and naively
predicts one combined shock, except for very small variances. Whereas, the CDF
method accurately captures the solution of MC at both early and later times.

Figure 13: Approximations of the first two moments from LOA, MC and CDF method
for the horizontal flooding with $q(t)$ random,
$\sigma_{q}=0.01,0.05,0.1,0.25$, for all cases $\mu_{q}=0.3$ and
$\tau_{q}=0.5T$, $T=1$. As variance increases, LOA approximations of the first
two moments deteriorate, whereas CDF results stay in agreement with the
corresponding Monte Carlo solutions.

Figure 14: Approximations of the first two moments from LOA, MC and CDF method
for the downdip flooding with $q(t)$ random,
$\sigma_{q}=0.005,0.01,0.02,0.05,0.1$

Figure 15: Approximations of the first two moments from LOA, MC and CDF method
for the horizontal flooding with $\phi(x)$ random at dimensionless times
$t_{1}=0.15,t_{2}=0.25$. For all cases, correlation length is
$\lambda_{\phi}=0.5L,L=1$ and variance is increasing as
$\sigma_{\phi}^{2}=10^{-5},10^{-4},10^{-3},10^{-2}$, respectively. By
increasing the variance, LOA approximations deviate from those of Monte Carlo,
whereas approximations of the CDF method stay in agreement with those of Monte
Carlo.

Figure 16: Approximations of the first two moments from LOA, MC and CDF
methods for the horizontal flooding with $\phi(x)$ random at dimensionless
times $t_{1}=0.15,t_{2}=0.25$. For all cases, $\sigma_{\phi}^{2}=10^{-3}$,
$\Delta x=0.005,\Delta t=0.01$, while the correlation length is increasing as
$\lambda_{\phi}=0.01,0.05,0.5,0.7$, respectively.

Figure 17: Low-order approximation of downdip flooding, $\phi(x)$ random, (1)
$\sigma^{2}_{\phi}=0.0001$, $\lambda_{\phi}=0.01L$, (2)
$\sigma^{2}_{\phi}=0.001$, $\lambda_{\phi}=0.01L$,(3)
$\sigma^{2}_{\phi}=0.01$, $\lambda_{\phi}=0.01L$, (4)
$\sigma^{2}_{\phi}=0.01$,$\lambda_{\phi}=0.1L$, (5) $\sigma^{2}_{\phi}=0.01$,
$\lambda_{\phi}=0.5L$, $t_{1}=0.02,t_{2}=0.1$

##  6 Convergence

For the MC simulations, $N_{MC}=3000-5000$ were required for the mean and
standard deviation profiles to converge, with a grid size of $\Delta
x=0.0005,\Delta t=0.01$. Also, a large number of realizations were essential
for the CDF plots to be smooth, and to be a perfectly horizontal line in the
shock region. Consequently, both shock and rarefaction regions would be
matching with those predicted by the method of characteristics. Regarding the
KDE estimations of the CDF of random fields in the MD method, $N=1000-3000$
realization are critical for the resulting $F_{s}(s;x,t)$ to be matching with
the corresponding of the MC benchmark. As Fig. 18 represents, by increasing
the number of MC samplings, the error will keep declining up to $N_{M}C=5500$,
after which increasing the number of trails was not observed to decrease the
error. For the sensitivity of solutions to the spatial and temporal grid
sizes, usually for $\Delta x<0.001$ and $\Delta t<0.001$ the lowest error was
achieved for most cases. The number of realizations required for MC to
converge increases with order of the moment. We experimented with skewness and
kurtosis, and even though the results are not included here, they usually need
up to $N_{MC}=10^{5}$ realizations in order to converge.

Figure 18: Grid convergence study for temporal grid size at constant $\Delta
x=0.001$ (top left), spatial grid size at constant $\Delta t=0.01$ (top
right), and number of MC realizations required for the moments to converge,
and for the spatial and temporal evolution of the CDF plots to be smooth.

##  7 Accuracy and Efficiency of the Method of Distributions

Figure 19: $L_{2}$ norm error of $F_{s}(S_{w})$ computed at two specific
dimensionless times, using Monte Carlo and CDF method for $q(t)$ as the sole
source of uncertainty. Sensitivity to the correlation length and variance of
the input random variable is represented.

Figure 20: $L_{2}$ norm error of $F_{s}(S_{w})$ computed at a specific time
using Monte Carlo and CDF method for the test cases with both $\phi(x)$ and
$q(t)$ random. Sensitivity to the correlation length/time and variances of
both input random variables quantifies accuracy of our CDF method.

The $L_{2}$ norm of error measures the discrepancy between average
$F_{s}(S,x,t)$ along the domain, obtained from Monte Carlo and the
corresponding results from the method of distributions. As shown in Fig. 19
and Fig. 20 . In Fig. 20, while experimenting with the sensitivity of the
$F_{s}(s)$ to the correlation length/time of the input stochastic fields, for
the $\lambda_{\phi(x)}>0.5L$ less error is observed. Also, initial values of
correlation time, i.e. $\tau_{q}(t)<0.5T$ incur a high error compared to
setting the correlation time higher than this threshold. Overall, the errors
are usually bounded by $0.1$ for the grid sizes $\Delta x=0.001,\Delta t=0.01$
that we have used in the accuracy experiments. Less error could be obtained if
the grid sizes, specially $\Delta x$ was selected to be smaller, however, the
MC simulations would take a long time as well. Also, while experimenting with
the values of variances of both input random fields, as Fig. 20 represents,
the highest error is observed when both variances are high.

In terms of the efficiency, our method is usually two orders of magnitude
faster than the corresponding MC experiment. The only time consuming part in
our CDF method, is generating the CDF of input random field, i.e.
$F_{q(t)}(t)$, $F_{\phi(x)}(x)$, $F_{q(t)(t),\phi(x)}(x,t)$ using KDE.
However, this process is usually two orders of magnitude faster than the Monte
Carlo process of generating 5000 samples, solving the nonlinear saturation
equation using Godunov method, and then KDE post-processing of the resulting
saturation field to find $F_{s}(s)$. In our CDF method, once the CDF of
underlying random field is generated, then a negligible computational cost is
associated with the ”semi-analytical” computation of the $F_{s}(s)$ using the
formula we derived before. To this end, the CDF method is always considerably
faster than corresponding Monte Carlo scheme.

##  8 Conclusion

We developed a semi-analytical formula for vertical heterogeneous reservoirs.
The test cases included horizontal flooding, updip flooding, downdip flooding,
and pure gravity segregation in an inverted gravity column, with sealed top
and bottom boundaries. For the updip case, the results were very similar to
those of the horizontal flooding, except that the rarefaction region is
shorter (Fig. (1), and hence the results are not included here for the sake of
brevity). We developed our CDF methodology in one dimension by extending a
previous study ([37]) for only horizontal reservoirs with a random injection
flux as the sole source of uncertainty, where the analysis was done in the
limit of $\sigma_{q}=0$, i.e. the deterministic case. Our method takes the
gravitational forces into account, and hence deals with the more complicated
spatio-temporal characteristics when two shocks are present in the problem.
Also, we covered all cases of random correlated injection flux and/or random
correlated porosity field, for a wide range of correlation lengths/times and
variances, with a reasonable accuracy for the CDF $F_{s}(s)$. The above-
mentioned uncertain fields lead to the uncertainty in the state variable
$s(x,t)$ described by the nonlinear Buckley Leverett problem, with a random
velocity field. The CDF method propagates uncertainty from input to the output
state variable by converting the original stochastic nonlinear problem to a
deterministic problem for the CDF of saturation, where the uncertainty is
encapsulated in the coefficients and initial/boundary conditions in a linear
fashion, making the problem more tractable. Having found the CDF of
saturation, the full probabilistic description of saturation is available and
hence we can easily find the higher order moments, i.e. skewness and kurtosis
which provide critical information about tails of the distribution. To the
best of our knowledge, this is the first work that converts the nonlinear
stochastic multi-shock hyperbolic conservation laws to a linear deterministic
equation for the CDF of state variable, and hence the methodology is able to
quantify the uncertainty for the more complicated physical scenarios than
single-shock physics in horizontal reservoirs. Additionally, we experimented
with both uniform and non-uniform initial conditions in our test cases, and
the CDF method was shown to match with the Monte Carlo counterparts for the
downdip case and gravity column with non-uniform boundary conditions. Also, no
closure approximations were required as the CDF method has its boundary
conditions straightforwardly defined. Moreover, we compared our results
against the so-called low-order approximation (LOA) of SME method, and the CDF
method was presented to match the results of the benchmark Monte Carlo,
whereas the SME results were valid only for very small random parameters, i.e.
variance and correlation length/time. Our LOA study extends the previous
studies done for the horizontal cases to the more general scenarios of
gravitational forces included. We presented the results for the imbibition
process, however, we can easily extend this methodology to the drainage
process by leveraging the descriptions in the flux function and Fig. 1. In
that case, the updip flooding needs treatment of two shocks, while downdip
flooding will have only one shock and one rarefaction zone, contrary to the
current study. Extension of this work to higher dimensions can be done by
manipulating the interaction of waves in multiple dimensions using a kinetic
defect term, which satisfies the entropy condition in all dimensions. A
comprehensive error analysis confirmed the robustness of the developed method
in this work.

## Acknowledgments

Financial support from Reservoir Simulation Industrial Affiliates Consortium
at Stanford University (SUPRI-B) is gratefully acknowledged.

## References

  * Aziz, [1979] Aziz, K. (1979).  Petroleum reservoir simulation.  Applied Science Publishers, 476.
  * Botev et al., [2010] Botev, Z. I., Grotowski, J. F., Kroese, D. P., et al. (2010).  Kernel density estimation via diffusion.  The annals of Statistics, 38(5):2916–2957.
  * Brenier and Jaffré, [1991] Brenier, Y. and Jaffré, J. (1991).  Upstream differencing for multiphase flow in reservoir simulation.  SIAM journal on numerical analysis, 28(3):685–696.
  * Brooks and Corey, [1964] Brooks, R. H. and Corey, A. T. (1964).  Hydraulic properties of porous media.  Hydrology papers (Colorado State University); no. 3.
  * Buckley et al., [1942] Buckley, S. E., Leverett, M., et al. (1942).  Mechanism of fluid displacement in sands.  Transactions of the AIME, 146(01):107–116.
  * Cheng et al., [2019] Cheng, M., Narayan, A., Qin, Y., Wang, P., Zhong, X., and Zhu, X. (2019).  An efficient solver for cumulative density function-based solutions of uncertain kinematic wave models.  Journal of Computational Physics, 382:138–151.
  * [7] Dongxiao, Z., Tchelepi, H., et al. (1999a).  Stochastic analysis of immiscible two-phase flow in heterogeneous media.  SPE journal, 4(04):380–388.
  * [8] Dongxiao, Z., Winter, C. L., et al. (1999b).  Moment-equation approach to single phase fluid flow in heterogeneous reservoirs.  SPE Journal, 4(02):118–127.
  * Fan et al., [2012] Fan, Y., Durlofsky, L. J., and Tchelepi, H. A. (2012).  A fully-coupled flow-reactive-transport formulation based on element conservation, with application to co2 storage simulations.  Advances in Water Resources, 42:47–61.
  * Fuks et al., [2019] Fuks, O., Ibrahima, F., Tomin, P., and Tchelepi, H. A. (2019).  Analysis of travel time distributions for uncertainty propagation in channelized porous systems.  Transport in Porous Media, 126(1):115–137.
  * Gelhar and Gelhar, [1993] Gelhar, L. W. and Gelhar, L. (1993).  Stochastic subsurface hydrology, volume 390.  Prentice-Hall Englewood Cliffs, NJ.
  * Graham and McLaughlin, [1989] Graham, W. and McLaughlin, D. (1989).  Stochastic analysis of nonstationary subsurface solute transport: 2. conditional moments.  Water Resources Research, 25(11):2331–2355.
  * Ibrahima et al., [2015] Ibrahima, F., Meyer, D. W., and Tchelepi, H. A. (2015).  Distribution functions of saturation for stochastic nonlinear two-phase flow.  Transport in porous media, 109(1):81–107.
  * Ibrahima and Tchelepi, [2017] Ibrahima, F. and Tchelepi, H. A. (2017).  Multipoint distribution of saturation for stochastic nonlinear two-phase transport.  SIAM/ASA Journal on Uncertainty Quantification, 5(1):353–377.
  * Ibrahima et al., [2018] Ibrahima, F., Tchelepi, H. A., and Meyer, D. W. (2018).  An efficient distribution method for nonlinear two-phase flow in highly heterogeneous multidimensional stochastic porous media.  Computational Geosciences, 22(1):389–412.
  * Jarman and Russell, [2003] Jarman, K. D. and Russell, T. F. (2003).  Eulerian moment equations for 2-d stochastic immiscible flow.  Multiscale modeling & simulation, 1(4):598–608.
  * Kitanidis, [1988] Kitanidis, P. K. (1988).  Prediction by the method of moments of transport in a heterogeneous formation.  Journal of Hydrology, 102(1-4):453–473.
  * Kwok and Tchelepi, [2008] Kwok, F. and Tchelepi, H. A. (2008).  Convergence of implicit monotone schemes with applications in multiphase flow in porous media.  SIAM Journal on Numerical Analysis, 46(5):2662–2687.
  * LeVeque et al., [2002] LeVeque, R. J. et al. (2002).  Finite volume methods for hyperbolic problems, volume 31.  Cambridge university press.
  * Li and Tchelepi, [2014] Li, B. and Tchelepi, H. A. (2014).  Unconditionally convergent nonlinear solver for multiphase flow in porous media under viscous force, buoyancy, and capillarity.  Energy Procedia, 59:404–411.
  * Li and Tchelepi, [2015] Li, B. and Tchelepi, H. A. (2015).  Nonlinear analysis of multiphase transport in porous media in the presence of viscous, buoyancy, and capillary forces.  Journal of Computational Physics, 297:104–131.
  * Li et al., [2013] Li, B., Tchelepi, H. A., and Benson, S. M. (2013).  Influence of capillary-pressure models on co2 solubility trapping.  Advances in water resources, 62:488–498.
  * Li et al., [2009] Li, H., Zhang, D., et al. (2009).  Efficient and accurate quantification of uncertainty for multiphase flow with the probabilistic collocation method.  Spe Journal, 14(04):665–679.
  * Li et al., [2005] Li, L., Tchelepi, H. A., et al. (2005).  Conditional statistical moment equations for dynamic data integration in heterogeneous reservoirs.  In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
  * Lie, [2019] Lie, K.-A. (2019).  An introduction to reservoir simulation using MATLAB/GNU Octave.  Cambridge University Press.
  * Likanapaisal, [2010] Likanapaisal, P. (2010).  Statistical moment equations for forward and inverse modeling of multiphase flow in porous media.  PhD thesis, Stanford University.
  * Meyer et al., [2010] Meyer, D. W., Jenny, P., and Tchelepi, H. A. (2010).  A joint velocity-concentration pdf method for tracer flow in heterogeneous porous media.  Water Resources Research, 46(12).
  * Meyer and Tchelepi, [2010] Meyer, D. W. and Tchelepi, H. A. (2010).  Particle-based transport model with markovian velocity processes for tracer dispersion in highly heterogeneous porous media.  Water Resources Research, 46(11).
  * Meyer et al., [2013] Meyer, D. W., Tchelepi, H. A., and Jenny, P. (2013).  A fast simulation method for uncertainty quantification of subsurface flow and transport.  Water Resources Research, 49(5):2359–2379.
  * Müller et al., [2013] Müller, F., Jenny, P., and Meyer, D. W. (2013).  Multilevel monte carlo for two phase flow and buckley–leverett transport in random heterogeneous porous media.  Journal of Computational Physics, 250:685–702.
  * Orr et al., [2007] Orr, F. M. et al. (2007).  Theory of gas injection processes, volume 5.  Tie-Line Publications Copenhagen.
  * Pettersson and Tchelepi, [2014] Pettersson, P. and Tchelepi, H. (2014).  Stochastic galerkin method for the buckley-leverett problem in heterogeneous formations.  In ECMOR XIV-14th European Conference on the Mathematics of Oil Recovery.
  * Pope, [1985] Pope, S. B. (1985).  Pdf methods for turbulent reactive flows.  Progress in energy and combustion science, 11(2):119–192.
  * Rubin, [2003] Rubin, Y. (2003).  Applied stochastic hydrogeology.  Oxford University Press.
  * Tchelepi, [1995] Tchelepi, H. A. (1995).  Viscous fingering, gravity segregation and permeability heterogeneity in two-dimensional and three-dimensional flows.
  * Van Genuchten, [1980] Van Genuchten, M. T. (1980).  A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.  Soil science society of America journal, 44(5):892–898.
  * Wang et al., [2013] Wang, P., Tartakovsky, D. M., Jarman Jr, K., and Tartakovsky, A. M. (2013).  Cdf solutions of buckley–leverett equation with uncertain parameters.  Multiscale Modeling & Simulation, 11(1):118–133.
  * Wang and Tchelepi, [2013] Wang, X. and Tchelepi, H. A. (2013).  Trust-region based solver for nonlinear transport in heterogeneous porous media.  Journal of Computational Physics, 253:114–137.
  * Winter et al., [2003] Winter, C. L., Tartakovsky, D., and Guadagnini, A. (2003).  Moment differential equations for flow in highly heterogeneous porous media.  Surveys in Geophysics, 24(1):81–106.
  * Yang et al., [2019] Yang, H. J., Boso, F., Tchelepi, H. A., and Tartakovsky, D. M. (2019).  Probabilistic forecast of single-phase flow in porous media with uncertain properties.  Water Resources Research.
  * Yousefzadeh, [2020] Yousefzadeh, M. (2020).  Numerical Simulation of Fluid-Mineral Interaction and Reactive Transport in Porous and Fractured Media.  PhD thesis, Stanford University.
  * Yousefzadeh and Battiato, [2017] Yousefzadeh, M. and Battiato, I. (2017).  Physics-based hybrid method for multiscale transport in porous media.  Journal of Computational Physics, 344:320–338.
  * Zaleski and Panfilov, [2017] Zaleski, S. and Panfilov, M. (2017).  Model of kinematic waves for gas–liquid segregation with phase transition in porous media.  Journal of Fluid Mechanics, 829:659–680.
  * Zhang et al., [2000] Zhang, D., Li, L., Tchelepi, H., et al. (2000).  Stochastic formulation for uncertainty analysis of two-phase flow in heterogeneous reservoirs.  Spe Journal, 5(01):60–70.