Probabilistic forecast of multiphase transport under viscous and buoyancy forces in heterogeneous porous media

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.


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

Problem formulation
We consider the nonlinear, incompressible, two-phase (Darcy) flow in a heterogeneous reservoir with intrinsic heterogeneous permeability K and spatially-varying porosity φ(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 nonwetting 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 (µ w and µ o ), parameters describing rock properties, such as porosity and permeability (φ(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 φ(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 ϕ (φ) and p q (Q) to them, fully characterized by their known mean (µ), standard deviation (σ 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 α = 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, Where capillary pressure relating the pressures of two phases is P c (S w ) = P o − P w . In the above-mentioned equations, u α (α = {w, o}) are the phase Darcy velocities of the wetting (water) and nonwetting (oil) phases. k rα , µ α and ρ α (α = {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), 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 ∇ · u T = 0 for total Darcy flux. Subsequently, we obtain an equation in terms of pressures, Where the relative mobilities are defined as λ α = kk rα µ α (α = {w, o}) for each phase and the total mobility reads λ T = λ o + λ 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, 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 θ, where 0 ≤ θ ≤ 90, ranging from 0 for the horizontal reservoirs to 1 for the vertical case studies, resulting in ∇h = sinθ. The dimensionless gravity number N g is defined as the ratio of buoyancy to viscous forces, Subsequently, defining the viscosity ratio m as m = µ w µ o , while employing the definitions for mobilities, the fractional flow of water can be recast in the following form, 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, Where v T (x, t) is defined as the interstitial velocity (total seepage ve- . 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, In this study, we assume the total velocity field 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, ∇f w (S w ) = f w (S w )∇S w . Thereupon, supplemented with the following initial/boundary conditions, Buckley Leverett equation defines an IBVP for water saturation which reads, Where S wi is the initial water saturation in the heterogeneous reservoir domain Ω. 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 Γ i , i.e. x = 0, where the wetting phase is being injected. S oi is the irreducible oil saturation. Γ 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, 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 θ, • N g sinθ = 0 turns Eq. (7) to that of a horizontal reservoir, meaning gravity plays no role and viscous forces are the only driving mechanism.
• N g sinθ > 0 represents an updip flooding, in which gravity retards the water flow, thereby reducing velocity of the front movement. It is noteworthy that sinθ > 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.
• N g sinθ < 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θ < 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.
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 = 1, and the corresponding flux function is asymmetric, and hence, the two waves move asynchronously, i.e. one may hit its boundary much earlier  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.
than the other, depending on m. Once reflected from both boundaries, the waves follow a very small rarefaction region. Afterwards, at S R and S 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 u T = 0, we have u o = 0, u w = 0, where the continuity condition results in u o + u w = constant. Assuming boundary conditions u o = u w = 0 at the top and bottom of the domain, we conclude the constant is zero, and hence u o = −u w . Consequently, 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 u o = −u w in the absence of capillarity, we end up with, Where the first fraction is assumed to be a characteristic velocity resulted from gravity 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 ) = u w u c , and is equal to, Therefore, the transport equation is recast as, 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, A closer inspection of Eq. (9) immediately reveals the fact that hyperbolic Buckley Leverett equation is indeed in the form of an advection equa- . 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), Where the second term drops for the horizontal cases with a gravity number of (N g = 0). Similarly, for the pure segregation scenario, A hyperbolic conservation law with piece-wise constant initial condition 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. (9) 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 v T (x, t), can be expresses as following ( [1], [31]), 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, , 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), Where the continuity condition imposed on the speed of the shocks deter- For the special case of horizontal reservoir domain (with fractional flow supplemented by Brooks and Corey relative permeabilities), S has an analytical expression (Eq. 4.7 in [14]). Nevertheless, we resort to a general numerical root finding scheme which handles N g = 0 cases as well.

Nondimensionalization
We will investigate the impact of uncertainty in the inputs on the output saturation field, first by assuming that φ(x) is random and q is a deterministic constant, second by assuming that q T (t) is random and φ 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 µ q and correlation length τ q , we can make the transport Eq. (9) dimensionless by defining, This rescaling will recast Eq. (9) as the dimensionless form in the onedimensional domain Ω, For the cases in which q is a deterministic constant, we can simply use the dimensionless quantities below to make Eq. (9) dimensionless, This definitions lead to the dimensionless form of Eq. (9), In order to make transport Eq. (13) 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, Where H is the total length of the vertical domain and τ = µ w H g(ρ o − ρ w )k is the characteristic time scale for gravity segregation. That is, if τ = 1, the particle has enough time to travel along the whole H. Such definitions allow us to work with the dimensionless form of Eq. (13) as,

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 ; x, t) for more general cases of gravitational forces taken into consideration. This approach overcomes the limitations of solving the nonlinear hyperbolic Eq. (9), by introducing a linear hyperbolic equation for the so-called random "raw" CDF function of water saturation, where H(.) is the Heaviside step function and θ 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. (9), the linear stochastic hyperbolic equation for the raw CDF reads [37], 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. (15) for the first three cases, and it follows Eq. (16) for pure gravity segregation.
As mentioned before, Eq.  27), into two sub-regions similarly to [37], 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 , It should be noted that Π a and Π 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 π 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 , The ensemble average of Π 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, Therefore, by solving a linear equation for π, followed by finding only the first ensemble moment of π, 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. (9) 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 (θ = S wi ; x, t) = 0 and F s (θ = 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. (27), 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], In order to define a family of characteristics x = x(t; x 0 ) along which the original hyperbolic PDE becomes an ODE dπ dt = 0, we start by using x(t = 0) = x 0 , where x 0 specifies where the characteristic line has originated from. Thus x = t 0 v(θ, x(t ), t )dt + 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 nonsynchronism 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 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 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.
Where Π b R and Π b L manifest development of characteristics in the right and left branches, respectively. Also, we use Π(t = 0) = H(θ − S wi ), Π(x = We will present the general formulation of CDF of water saturation, with uncertainty in both porosity field and injection flux encapsulated in . 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]), We introduce F U (x,t) , CDF of the input random velocity field U (x, t), and employ the notation a ∨ b and a ∧ b for the max(a, b) and min(a, b) respectively. We expand this relation using a similar approach in [14] as follows, For the downdip case as well as the pure gravity segregation case before reflection of the waves from boundaries, the CDF equation reads, Leveraging F U (x,t) , this relation expands as follows. For x > x D , S wi < θ < S R we have, And for x < x D , S D < θ < S L we use, Whereas the gravity segregation scenario after the reflection of waves follows the CDF relation below, This relation is consequently expanded as below. For x < x S L , S wi < θ < S L we have, Analogously to the previous derivations, for x S L < x < x S R , S R < θ < S L we have, For x > x S R , which corresponds to S R < θ < 1 − S oi we leverage, Where F s b (θ; x, t) should be found from the solutions of the before reflection time in Eq. 39 and Eq.40. However, those solutions were developed for (S R , S L ) 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.43 and Eq.46 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 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 φ(x) (x) when porosity field φ(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 (θ; x, t) obtained from Eqs.
. 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.

Setting
For the Monte Carlo simulation of the nonlinear hyperbolic equation for saturation Eq. (9) 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. (9) 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 pointwise 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 φ(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 nonwetting 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, 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).

Numerical setting for stochastic porosity field φ(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 φ(x) is characterized by its know mean, variance, covariance structure, and dimensionless correlation length.
Throughout our experiments with a stochastic porosity field, unless otherwise stated, we use µ φ = 0.3, σ 2 γ = 0.5, and λ γ = 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 σ 2 φ and λ φ 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. 45 and Eq. 44, 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. 45 and Eq. 44. Therefore, even though the proposed method perfectly captures the two-wave behavior in all the previous examples, it fails only for this specific case.

Numerical setting for stochastic injection flux q(t)
Enforcing the general continuity constraint for incompressible flow, ∇ · 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. (9) 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 φ = 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 µ q , variance σ 2 q , covariance structure C q (t), and correlation time τ q are all known parameters. To this end, similarly to [13], we first define a Gaussian random field ρ q (t) ∼ N (µ ρ , C ρ (t)), with an exponential covariance structure C ρ (t) = where τ q is the correlation parameter for events that are closer in time. Eventually, we define the stochastic injection flux field as, q T (t) = 0.05 + 0.1(ρ q (t)) 2 Such a definition allows us to first make sure that q T (t) > 0 and second by controlling µ ρ we can satisfy the criteria q(t) = µ q . Any other arbitrary non-negative distribution could be used alternatively.
Having a time discretization of size ∆t, for any point t in the time interval [0, T ], where t = (n − 1)∆t, we define the Gaussian vector (ρ q ) n = (ρ q (j∆t)) 0≤j≤n−1 by (ρ q ) n ∼ N (µ ρ 1 n , (C ρ ) n ), where C n ρ = (C q (|i−j|∆t)) 1≤i,j≤n and 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 µ q = 0.3, σ 2 q = 0.5, τ 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.

Numerical setting for uncertainty in both φ(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 φ(x) and q T (t) fields. As for the numerical parameters, we will employ m = 0.5, µ q = 0.3, σ 2 q = 0.5, τ q = 0.5T , T = 1, µ φ = 0.3, σ 2 φ = 0.5, λ φ = 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 v T = v T (x, t). Fig. 12 shows the results of having both fields as the source of randomness

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

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 φ 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 φ(x), while q T is a deterministic constant.
In the first scenario, we assume a deterministic constant value of φ = φ 0 = 0.3, and a lognormal distribution for q, with µ q = 0.3, while we will experiment with multiple values of variance σ 2 q and track the resulting first two moments of the state field S w , Similarly to [7], we consider the one-dimensional version of Eq. (9), and using the idea of displacement along a streamline ( [7]), we define τ (x; x 0 ), the travel time (time of flight) of a particle from x 0 = 0 to x in the total velocity field v T as, Which has a solution in the form of τ (x; Since v T is random, so is τ , and hence, its statistical moments depend on those of 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, Correspondingly, LOA results in the first-order approximation of the first two moments of time of flight as a function of µ q and σ 2 q , Furthermore, by prescribing a probability distribution to τ for the given mean and standard deviation, e.g. lognormal distribution τ ∼ logN ln(τ ) , σ ln(τ ) , where the corresponding moments are defined as ( [7]), Therefore, distribution of τ is fully delineated as a lognormal probability density function, along with the mean µ τ and standard deviation σ 2 τ , which are directly computable using µ q and σ 2 q in Eq. (54). Therefore, by controlling σ 2 q parameter, we can find a low-order estimate of the ensemble mean (expected value) and standard deviation of water saturation using the expressions below, Where where p τ (τ ; x, x 0 ) is the probability density function of the (oneparticle) travel time τ , andS(τ, t) can be found by combining Eq. (9) and Eq. (52) resulted in a transformed version of the transport equation for water saturation as, 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 φ(x) characterized by its mean value of µ φ = 0.3, and covariance structure Given such a definition for the random porosity field, we can leverage the first-order approximation of the first two moments of τ derived in [7] as following, Furthermore, τ is assumed log-normally distributed and the same definitions as Eq. (55) hold for ln(τ ) and σ 2 ln(τ ) . To this end, by controlling µ φ and λ φ , we will present comparisons of the low-order approximation for the first two moments of water saturation using Eq. (56) against those obtained from Monte Carlo and CDF methods.

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

Convergence
For the MC simulations, N M C = 3000 − 5000 were required for the mean and standard deviation profiles to converge, with a grid size of ∆x = 0.0005, ∆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 ∆x < 0.001 and ∆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 M C = 10 5 realizations in order to converge.     Figure 18: Grid convergence study for temporal grid size at constant ∆x = 0.001 (top left), spatial grid size at constant ∆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. Figure 19: L2 norm error of Fs(Sw) 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: L2 norm error of Fs(Sw) computed at a specific time using Monte Carlo and CDF method for the test cases with both φ(x) and q(t) random. Sensitivity to the correlation length/time and variances of both input random variables quantifies accuracy of our CDF method.

Accuracy and Efficiency of the Method of Distributions
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 λ φ(x) > 0.5L less error is observed. Also, initial values of correlation time, i.e. τ 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 ∆x = 0.001, ∆t = 0.01 that we have used in the accuracy experiments. Less error could be obtained if the grid sizes, specially ∆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 φ(x) (x), F q(t)(t),φ(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.

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