DTIC ADA445711: Eulerian Moment Equations for 2-D Stochastic Immiscible Flow

19 Pages·0.31 MB·English

by  Defense Technical Information Center

EULERIAN MOMENT EQUATIONS FOR 2-D STOCHASTIC
IMMISCIBLE FLOW
KENNETH D. JARMAN∗ AND THOMAS F. RUSSELL†
Abstract. We solve statistical moment differential equations (MDEs) for immiscible flow in
porous media in the limit of zero capillary pressure, with application to secondary oil recovery.
ClosureisachievedbyTaylorexpansionofthefractionalflowfunctionandaperturbationargument.
Previousresultsin1-Dareextendedto2-D,inwhichabimodalprofileislessevident.
Mean and variance of (water) saturation exhibit a bimodal character; two shocks replace the
singleshockfrontevidentintheclassicalBuckley–Leverettsaturationprofile. ComparisontoMonte
Carlo simulations (MCS) shows that the MDE approach gives a good approximation to total oil
production. For such integrated or averaged quantities, or where a rough approximation of the
locationandmagnitudeofuncertaintyissufficient, MDEsmaybesubstantiallymoreefficientthan
MCS.
Key words. porousmedia,stochastic,randomfields,momentequations
AMS subject classifications. 76T99,76S05,35L60,35L65
1. Introduction. Stochastic representations of subsurface geologic properties
have become commonplace due to the difficulty in complete and certain characteriza-
tionoftheseproperties. Thisleadstouncertaintyinflowprofilesinsuchporousmedia,
so that statistical description of outcomes is appropriate. Additionally, from the per-
spective of practical macroscopic field-scale models, microscopic heterogeneities and
flows can be viewed as random processes to be upscaled. In two-phase flow on a field
scale, we are primarily interested in mean behavior and a measure of uncertainty in
thesaturation ofeach phase, which canbea measureof “macrodispersion” incertain
upscaling contexts.
A “zeroth-order” model of mean flow with averages of geologic properties ignores
correlations between flow variables. Monte Carlo simulations (MCS) of many real-
izations of geologic properties to estimate moments require much computation time
and careful sampling techniques [10,11,29]. Macrodispersion theory in contaminant
transportcapturesanapproximateeffectoffluctuationsbymodelingcovariancefunc-
tionsinanequationforthemeanconcentration[5,12]. Thistheoryhasalonghistory
in subsurface contaminant transport and is closely related to eddy diffusion models
of turbulence [15, p. 358 ff].
Rather than approximate covariance functions as macrodispersive terms, which
requires neglecting fluctuations in second moments, we solve PDEs for the covari-
ance functions and the mean of saturation, extending previous results to 2-D. These
moment differential equations (MDEs) are derived using a modified perturbation ex-
pansionasdescribedin[17]. TheMDEsallowdirectapproximationofthelocalmean
and covariance functions for general boundary conditions and general, nonstationary
stochastic geology [33].
∗Computational Sciences and Mathematics Division, Fundamental Science Directorate, Pa-
cific Northwest National Laboratory, P.O. Box 999 / MS K1-85, Richland, Washington 99352
([email protected]). Supported by NSF Graduate Traineeship in Applied Mathematics (DMS-9208685)
andDOELaboratoryDirectedResearchandDevelopment.
†DepartmentofMathematics,UniversityofColoradoatDenver,P.O.Box173364,CampusBox
170,Denver,Colorado80217-3364([email protected]). SupportedinpartbyNSFGrant
Nos.DMS-9706866,DMS-0084438andDMS-0222300,andAROGrantNo.DAAG55-97-1-0171.
1
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
2006 2. REPORT TYPE 00-00-2006 to 00-00-2006
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Eulerian Moment Equations for 2-D Stochastic Immiscible Flow
5b. GRANT NUMBER
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
Pacifica Northwest National Laboratory,Computational Sciences and REPORT NUMBER
Mathematics Division,PO Box 999,Richland ,WA,99352
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
see report
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 18
unclassified unclassified unclassified
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
2 K.D.JARMANANDT.F.RUSSELL
1.1. Applications. A statistical description of subsurface flow is of particular
interestforsecondaryoilrecovery. Intheadvectionequation,theprincipaldifficultyis
a non-convex nonlinear flux function that leads to discontinuous solutions. Standard
pressure-saturationequationsfor2-Dhorizontalflowoftwoincompressibleimmiscible
fluids in porous media in the limit of vanishing capillary pressure (pc =0) are
φv(x)=−K(x)∇h(x), ∇·v(x)=0, (1.1)
∂ts(x,t)+∇·[f(s(x,t))v(x)]=0. (1.2)
These are considered valid from laboratory (centimeters) to field scales of reservoir
depth (10–100 meters) and length (100–10,000 meters). Hydraulic conductivity K
may be an anisotropic tensor; here, for simplicity, it will be an isotropic scalar K.
Assume also that K depends weakly on (water) saturation s [1,3]. Apply (1.1)–(1.2)
to the flow of oil and water, for arbitrary fluid mobilities. Denote the total velocity,
a scaled total volumetric flux of both fluids, by v, hydraulic water head by h, and
porosity (assumed constant) by φ. The fractional flow function f(s), for pc = 0,
represents the fraction of v due to water. It is typically S-shaped as shown in Fig.
1.1; in §4.1 we use a form arising from quadratic relative permeabilities (see [1]).
However, our method does not depend on any such specific choice of f(s).
1
s)
f(
0
0 1
Saturation (s)
Fig. 1.1. The nonconvex fractional flow function.
Capillary pressure regularizes sharp fronts caused by the nonlinear advection
term. To obtain a linear approximation to this effect, add (cid:3)D∇2s(x,t) with (cid:3)D > 0
totherightsideof (1.2). Letting(cid:3)D →0definesthevanishing-viscositysolution[27],
which is the one we seek.
As is standard in subsurface applications, let Y = lnK be a random field with
prescribedmeanandcovariancefunctions;e.g.,itisoftenclaimedthatY ismultivari-
ateGaussian, basedonempiricalobservations[8](againourmethoddoesnotdepend
on such a limiting assumption). Through (1.1)–(1.2), v and s are thus random fields.
No other underlying sources of uncertainty are considered in this study.
Under the assumptions stated above, with steady boundary conditions, and ne-
glecting the weak dependence of K on s, a steady v can be determined from (1.1).
WeevolvesfromthestochasticPDE(1.2),assumingknownstatisticsofv. Moments
of v and h can be estimated from established theory ([31–33,36] use MDEs). We
MOMENTEQUATIONSFOR2-DIMMISCIBLEFLOW 3
combineanalyticalandnumericaltechniquestomodelthepropagationofuncertainty
from an underlying random field Y(x), through v(x) to the solution s(x,t).
1.2. Previouswork. ExistingworkonMDEsfocusesmostlyonadvectionequa-
tions with linear flux functions [13] and some nonlinear subsurface flow equations of
a form different from (1.2) [2,28,31,32,36].
LangloandEspedal[19,20]presentamacrodispersionapproachforthestochastic
versionof (1.2). ThefluxfunctionisexpandedinaTaylorseries,andhigh-orderterms
are neglected; then standard techniques represent macrodispersivity as a function of
flowvelocitycovariance. Zhang,Tchelepi,andLitakeadvantageofthesteadyvelocity
field, and transform 2-D flow to 1-D Lagrangian flow along streamlines [34,35]. They
formulateintegralequationsformomentsfromensembleaveragesoverthestreamlines
rather than a system of MDEs. Variations on the streamline approach are currently
popular in stochastic subsurface hydrology [4,6,7,25,26].
An Eulerian MDE approach has been successful for single- and multiphase pres-
sure and velocity equations [31], and a natural next step is to extend the theory from
flux equations to transport equations (as in [13]) and (1.2). The Eulerian framework
differs from streamlines not only in formulation but also in that the MDEs need no
velocity-distributionassumptionandextensiontotransientvelocityfieldsisrelatively
straightforward. The appeal of MDEs relative to macrodispersion approaches is that
covariances are computed directly, so that it is not necessary to approximate them
by solving an approximate fluctuation equation that neglects second-moment fluctu-
ations. In [17] and [18], we derived and solved second-order MDEs for (1.2) in 1-D,
using three different perturbation approaches. One of these approaches, a modified
expansion described in [17], is applied here in 2-D.
Equations are presented in §2. In §3, we discuss classification of MDEs and the
fundamental difference between 2-D and 1-D MDEs. Numerical solutions are given
in §4 and are compared to Monte Carlo simulations.
2. Moment Equations. We use a modified perturbation expansion, described
in detail in [17], to derive statistical MDEs. Closure implicitly depends upon the
assumption that random fluctuations in Y =lnK are small (σY (cid:3)1). Moments of h
and v are assumed known from (1.1) and moments of Y.
The saturation equation (1.2), with initial data s(x,0)=g(x), is
∂ts+∂xi(f(s)vi)=0 (2.1)
(Einsteinsummationconventionassumed). Weassumethatgisknownwithcertainty.
Solutions are defined in terms of vanishing viscosity as in §1.1; henceforth, this is
tacitly understood.
For examples of commonly used methods for deriving moment equations with
applications to subsurface flow and transport, see [5,9,12,13,23,31]. These meth-
ods originated in the correlation equations of turbulence models. Here we apply a
modification to a standard asymptotic expansion.
Let (cid:4)·(cid:5) denote the expectation operator, defined by
(cid:1)
(cid:4)ψ(cid:5)≡ ψ(ω)dP(ω) (2.2)
Ω
foranyintegrablefunctionψ :Ω→RonthesamplespaceΩwithprobabilitymeasure
P. We omit reference to ω in what follows.
4 K.D.JARMANANDT.F.RUSSELL
The random field Y is decomposed into deterministic mean plus random fluctu-
ation: Y = (cid:4)Y(cid:5)+δY. Each field dependent on Y is represented by a formal power
expansion in an as-yet unknown parameter (cid:3):
(cid:2)∞ (cid:2)∞ (cid:2)∞
n n n
h= (cid:3) hn(x), v= (cid:3) vn(x), s= (cid:3) sn(x,t). (2.3)
n=0 n=0 n=0
The expansion parameter (cid:3) = σY is shown to be appropriate within the context of
the velocity and head equations (1.1) [5, pp. 184–190], [33]. For example, for single-
phase,stationaryuniformmeanflowonaninfinite(cid:3)do(cid:4)mainin1-D,v0isadeterministic
scalar and σ2 is shown to be approximated by (cid:3)2 v2 =v2σ2. In 2-D, we denote the
v 1 0 Y
components of vn by (vi)n for i=1,2.
In the following, the term “order” applies to the power of (cid:3) rather than the order
of the statistical moment. In special cases, we may use the term “order” in the
latter sense and clearly differentiate from order in (cid:3). Thus, a “second moment” is a
covariance, but“second-ordermoment”isanymomentthatisapproximatedtoorder
(cid:3)2. Recall that we only need the decompositions of v and s here. The fractional flow
functionisexpandedinaTaylorseriesarounds˜=s0+(cid:3)(cid:4)s1(cid:5)+(cid:3)2(cid:4)s2(cid:5)(asecond-order
approximation to the mean saturation; we see below that (cid:4)s0(cid:5)=s0):
1
f(s)=f(s˜)+f(cid:2)(s˜)(s−s˜)+ f(cid:2)(cid:2)(s˜)(s−s˜)2+··· , (2.4)
2
where
(cid:2)∞
s−s˜=δs0+(cid:3)δs1+(cid:3)2δs2+ (cid:3)nsn(x,t).
n=3
The expansion center s˜is our approximation to the mean saturation, so we retain s˜
as the argument of f in the following modification to the standard equations of order
0,1, and 2. The modified equations for s0 and s1, obtained by substituting (2.3) and
(2.4) into (2.1) and collecting the powers (cid:3)0 and (cid:3)1, are given by
∂ts0+∂xi[f(s˜)(vi)0]=0, (2.5)
and
(cid:5) (cid:6)
(cid:2)
∂ts1+∂xi f(s˜)(vi)1+f (s˜)δs1(vi)0 =0. (2.6)
Note that (cid:4)s0(cid:5)=s0, because all other terms in the equation for s0 are deterministic.
We show that (cid:4)s1(cid:5) = 0. By definition, (cid:4)δs1(cid:5) = 0. Also, (cid:4)v1(cid:5) = 0 by a standard
result from the theory for moment equations of multidimensional single-phase flow
[36]. Applying the expectation operator (cid:4)·(cid:5) to the first-order equation (2.6) gives
(cid:5) (cid:6)
0=∂t(cid:4)s1(cid:5)+∂xi f(s˜)(cid:4)(vi)1(cid:5)+f(cid:2)(s˜)(cid:4)δs1(cid:5)(vi)0 =∂t(cid:4)s1(cid:5). (2.7)
The initial data are deterministic; therefore, (cid:4)s1(cid:5) is initially zero, and (2.7) implies
that (cid:4)s1(cid:5) remains zero for all time. The mean approximation is then s˜=s0+(cid:3)2(cid:4)s2(cid:5).
The correction term s2 satisfies
(cid:5) (cid:6)
1
(cid:2) (cid:2) (cid:2)(cid:2) 2
∂ts2+∂xi f(s˜)(vi)2+f (s˜)δs1(vi)1+f (s˜)δs2(vi)0+ 2f (s˜)δs1 (vi)0 =0.
MOMENTEQUATIONSFOR2-DIMMISCIBLEFLOW 5
(cid:3) (cid:4) (cid:3) (cid:4)
Apply the expectation operator and note that δs12 = (cid:4)s1s1(cid:5)−(cid:4)s1(cid:5)(cid:4)s1(cid:5) = s12 ,
(cid:4)δs2(cid:5) = 0, and v1 = δv1 implies (cid:4)δs1(vi)1(cid:5) = (cid:4)δs1(δvi)1(cid:5) = (cid:4)s1(vi)1(cid:5)−(cid:4)s1(cid:5)(cid:4)(vi)1(cid:5) =
(cid:4)s1(vi)1(cid:5), to obtain the following correction to (2.5):
(cid:5) (cid:3) (cid:4) (cid:6)
1
∂t(cid:4)s2(cid:5)+∂xi f(s˜)(cid:4)(vi)2(cid:5)+f(cid:2)(s˜)(cid:4)s1(vi)1(cid:5)+ 2f(cid:2)(cid:2)(s˜) s12 (vi)0 =0. (2.8)
Define v˜ = v0+(cid:3)2(cid:4)v2(cid:5), which is consistent with the second-order theory for single-
phase flow. Adding the zeroth-order equation (2.5) to (cid:3)2×(2.8) gives an equation for
the mean saturation:
(cid:5) (cid:3)2 (cid:3) (cid:4) (cid:6)
∂ts˜+∂xi f(s˜)v˜i+(cid:3)2f(cid:2)(s˜)(cid:4)s1(vi)1(cid:5)+ 2f(cid:2)(cid:2)(s˜) s12 (vi)0 =0. (2.9)
(cid:3) (cid:4)
We derive equations for the unknown covariance functions (cid:4)s1(vi)1(cid:5) and s21
usingthefollowingadditionalnotation. Theindependentvariablesarexandtexcept
where noted, and ·|y denotes the replacement of x by some y different from x. It is
convenient and useful to derive equations for the more general two-point covariances
(cid:4)s1(vi)1|y(cid:5) and (cid:4)s1s1|y(cid:5) rather than for one-point covariances.
To obtain an equation for (cid:4)s1(vi)1|y(cid:5), multiply (2.6) by (vi)1(y), and apply (cid:4)·(cid:5).
This results in
(cid:5) (cid:6)
∂t(cid:4)s1(vi)1|y(cid:5)+∂xi f(s˜)(cid:4)(vi)1(vi)1|y(cid:5)+f(cid:2)(s˜)(vi)0(cid:4)s1(vi)1|y(cid:5) =0. (2.10)
Similarly, multiply (2.6) by s1(y,t) and use the identity1 s1|y∂ts1 =∂t(s1s1|y)−
s1∂ts1|y, then use (2.6) with y in place of x, yielding this equation for the two-point
saturation covariance:
(cid:5) (cid:6)
∂t(cid:4)s1s1|y(cid:5)+∂xi f(s˜)(cid:4)s1|y(vi)1(cid:5)+f(cid:2)(s˜)(vi)0(cid:4)s1s1|y(cid:5)
(cid:5) (cid:6)
+∂yi f(s˜|y)(cid:4)s1(vi)1|y(cid:5)+f(cid:2)(s˜|y)(vi)0|y(cid:4)s1|ys1(cid:5) =0. (2.11)
Define csvi(x,y,t) = (cid:3)2(cid:4)s1(vi)1|y(cid:5), cs = (cid:3)2(cid:4)s1s1|y(cid:5), cvivj = (cid:3)2(cid:4)(vi)1(vj)1|y(cid:5),
(cid:4)s(cid:5)|y = (cid:4)s(cid:5)(y,t), σsvi(x,t) = csvi(x,x,t), and σs2(x,t) = cs(x,x,t). Let a caret over
(cid:7)
avariabledenotethemappingψ(x,y,t)=ψ(y,x,t)foranyfunctionψ. Theresulting
system is
(cid:5) (cid:6)
1
(cid:2) (cid:2)(cid:2) 2
∂ts˜+∂xi f(s˜)v˜(cid:5)i+f (s˜)σsvi + 2f (s˜)(vi)0(cid:6)σs =0, (2.12a)
(cid:2)
∂tcsvj +∂xi f(s˜)cvivj +f (s˜)(vi)0csvj =0, (2.12b)
(cid:5) (cid:6) (cid:5) (cid:8) (cid:9) (cid:8) (cid:9) (cid:6)
∂tcs+∂xi f(s˜)(cid:7)csvi +f(cid:2)(s˜)(vi)0cs +∂yi f (cid:7)s˜ csvi +f(cid:2) (cid:7)s˜ (vi)0|ycs =0; (2.12c)
i,j =1,2.
Recall that velocity moments (vi)0, v˜i, and cvivj are assumed known. The flux in
(2.12a) consists of a nonlinear advective mean flux term and two covariance terms.
In turbulence applications, terms such as csvj often are referred to as transport by
fluctuations [22]. Both (2.12b) and (2.12c) have advective flux terms, are coupled to
1Theidentityisnotvalidinastrongsensefordiscontinuoussolutions. Recall,however,thatwe
definesolutionsintermsofthe(smooth)viscoussolution,inthelimit(cid:1)D →0.
6 K.D.JARMANANDT.F.RUSSELL
the mean equation (2.12a), and are first-order in (cid:3)2 =σ2. This is consistent with the
Y
approximations˜,whichissecond-orderinσY. However,duetothechoiceofs˜,higher-
order effects are included in (2.12). We demonstrate this fact for the case f(s) = s,
where the mean advective flux term in (2.12a) is ∂xi[s˜v˜i], or ∂xi[(s0+(cid:3)2(cid:4)s2(cid:5))((vi)0+
(cid:3)2(cid:4)(vi)2(cid:5))]. The term (cid:3)4∂xi[(cid:4)s2(cid:5)(cid:4)(vi)2(cid:5)] in this expression is a partial fourth-order
correction.
We refer to this problem as 2-D although the covariance functions involve four
spatial variables and the saturation covariance has fluxes in all directions. The four
variables represent two different points in the same 2-D domain.
To formally define the vanishing-viscosity solution, first add the deterministic
dteirffmusi(cid:3)oDn∂xt2eisrmtos (cid:3)(2D.∂1x)2ias˜,nd(cid:3)Dc∂axr2riycsvoju,tatnhde (cid:3)eDxp∂ax2nicssiontoatnhdedreigrihvta-thiaonndabsiodvees.oTfh(i2s.1a2dad)s,
(2.12b),and(2.12c),respectively. Thenfindsolutionsinthelimit(cid:3)D →0. Inpractice
we solve the system (2.12) directly.
3. Classification. Classification is of interest for several reasons. Macrodisper-
siontheoryproducesaparabolicequationforthemean,withamacrodispersivitythat
depends on cvivj and (cid:4)s(cid:5). In contrast, we have shown that 1-D MDEs are hyperbolic
[18]. To achieve this result, we exploited special structure to reduce the equations.
Namely,1-Dvelocityhasaninfinitecorrelationlengthbecauseitisaspatiallyconstant
random variable. This led to the covariance relationship cs cv = csv (cid:7)csv. Velocity
correlationshaveafinitelengthscalein2-D.Wearguethatthisfactprohibitsclassifi-
cationandsuggestsatransitiontoparabolicbehaviorasdescribedbymacrodispersion
theory.
We expect the equations to be nearly hyperbolic since the original deterministic
PDE (1.2) is hyperbolic. However, σsvi and σs2, rather than csvi and cs, appear
in the mean equation. This mix of one-point and two-point covariance functions
prevents us from immediately treating the MDEs as classically hyperbolic. Also,
independent variables x and y are permuted in s˜ and csvj in (2.12c). In general,
s˜|y = s˜(y,t) (cid:7)= s˜(x,t), and (cid:7)csvj (cid:7)= csvj. The system (2.12) cannot even be written
in conservation-law form. The inconsistency in independent variables is dealt with
directly by writing a redundant set of equations for the functions(cid:7)s˜and (cid:7)csvj.
Transposing x and y in (2.12a) and (2.12b) results in the expanded system
(cid:5) (cid:6)
1
(cid:2) (cid:2)(cid:2) 2
∂ts˜+∂xi f(s˜)v˜(cid:5)i+f (s˜)σsvi + 2f (s˜)(vi)(cid:6)0σs =0,
(cid:2)
∂tcsvj +∂xi f(s˜)cvivj +f (s˜)(vi)0csvj =0,
(cid:5) (cid:6) (cid:5) (cid:8) (cid:9) (cid:8) (cid:9) (cid:6)
∂tcs+∂xi f(s˜)(cid:7)csvi +f(cid:2)(s˜)(vi)0cs +∂yi f (cid:7)s˜ csvi +f(cid:2) (cid:7)s˜ (vi)0|ycs =0, (3.1)
(cid:5) (cid:8) (cid:9) (cid:8) (cid:9) (cid:6)
∂t(cid:7)csvj +∂yi f (cid:7)s˜ cvjvi +f(cid:2) (cid:7)s˜ (vi)0|y(cid:7)csvj =0,
(cid:5) (cid:8) (cid:9) (cid:8) (cid:9) (cid:8) (cid:9) (cid:6)
∂t(cid:7)s˜+∂yi f (cid:7)s˜ v˜i|y+f(cid:2) (cid:7)s˜ σ(cid:7)svi + 21f(cid:2)(cid:2) (cid:7)s˜ (vi)0|yσ(cid:7)s2 =0;
with i,j = 1,2, where σ(cid:7)svi(y,t) = csvi(y,y,t), σ(cid:7)s2(y,t) = cs(y,y,t). The system is
(cid:7)
symmetric under the permutation ψ(x,y,t) = ψ(y,x,t). To be consistent with the
originalsystem(2.12a)–(2.12c),theinitialconditionsandthenumericalmethodmust
obey this symmetry.
4. Solution. Wenumericallysolve(3.1)withafirst-orderupwindschemewithin
the framework of CLAWPACK [21]. However, we do not take full advantage of the
MOMENTEQUATIONSFOR2-DIMMISCIBLEFLOW 7
functionality of CLAWPACK due to the difficulty in defining a Riemann solver for
the MDEs. The solutions are shown to be bimodal as in 1-D (see Fig. 4.1, and see
[18] for 1-D details). Solutions are also compared to MCS moment estimates.
1
n
o 0.8
ati
ur 0.6
at
n s 0.4
a
e 0.2
M
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x
0.4
v.
e 0.3
d
d.
st 0.2
n
atio 0.1
ur
at 0
S
−0.1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
x
Fig.4.1. Meanandstandarddeviationofsaturation(1-D).Thesolidcurveisasemi-analytical
solution described in [18]; dashed curve is obtained from fine-grid numerical PDE scheme.
4.1. MDEs. Our example problem is on a rectangular domain R = [0,L1]×
[0,L2] = [0,2]×[0,1], x = (x1,x2) ∈ R, with ∂x2h(x) = 0 on x2 = 0 and x2 = 1
and constant h(x) at x1 = 0 and x1 = 2. The entire domain R is initially oil-
saturated (s = 0), with water (s = 1) pumped in along x1 = 0. We choose a form
of the fractional flow function f(s) that arises from quadratic relative permeabilities:
f(s) = s2/(s2 +m(1−s)2). We set the viscosity ratio m to 0.5 and porosity φ to
0.2. Statistical parameters for MDEs are (cid:4)Y(cid:5) = 0 and σ2 =0.25, and we use the
Y
exponential covariance function given by
(cid:10) (cid:11)
C(r)=σ2 exp −|r1| − |r2| , (4.1)
Y
λ1 λ2
where r=(r1, r2). Weset correlation lengths λ1 =λ2 =0.2. This gives 5 correlation
lengths in the x2 direction and 10 correlation lengths in the x1 direction. Hydraulic
headandvelocitymomentswerecomputedusingthesemi-analyticalflowapproachof
ZhangandWinter[30,36]. TheauthorsderivedMDEsusingperturbationexpansions
in a manner similar to ours. Due to the mathematical complexity caused by the
presenceofboundariesandnonstationarity,theseMDEsaresolvednumerically. Finite
velocity correlation length is evident in examples of velocity covariance shown in Fig.
4.2.
We use a very coarse grid (40 by 20) due to the fact that two-point covariance
functions make the 2-D MDEs essentially 4-D. The 2-D solutions for the mean s˜and
variance σ2 of saturation are shown at an early and late time in Fig. 4.3. These
s
indicate a pair of saturation fronts moving at distinct speeds. Saturation-velocity
8 K.D.JARMANANDT.F.RUSSELL
0.8
0.6
0.4
1
v
1
v
c
0.2
0
−0.2
0 0.2 0.4 0.6 0.8 1
x/L
i i
Fig. 4.2. Velocity covariance cv1v1, showing finite correlation length. The solid line shows
covariance between v1 at the center of the domain and at all points on x2=L2/2. The dotted line
shows covariance between v1 at the center and at all points along x1=L1/2.
one-point covariances σsvi, i=1,2 are shown in Fig. 4.4. The σsv1 surface is similar
to that of σ2, and is very close to the 1-D profile.
s
The results also are shown for s˜, σsv1, and σs along a slice at x2 = L2/4 on the
longitudinal axis in Fig. 4.5. In each case, a bimodal solution is evident. On the left,
a rarefaction zone trails a slow front. Physically this represents a smoothly varying
mixture of two phases. The uncertainty as measured by σ2 is relatively small in this
s
region. Beyondtheslowfront,intheintervalofuncertainty,thereisanearlyconstant
ratio of phases. Finally, at the fast front, mean saturation and covariances drop to
zero; beyond the front, only one phase (oil) is present.
Acrudetestnowcanbeperformedtoseewhetheracovariancerelationshipholds
similar to that in 1-D. The flow is nearly 1-D; transverse velocities are much smaller
than longitudinal velocities, and σsv2 is much smaller than σsv1. So it is reasonable
to expect that the products σs2σv1v1 and σs2v1 are similar, because they are identical
in1-D.However, observethedipinσsv1 betweentwopeaks(Fig.4.5), ontheinterval
ofuncertainty. Thesaturationvariancedoesnothavesuchaprofile,andneitherdoes
σv1v1 (seeAppendixA). Thusthereisadifferenceinshapebetweenthetwoproducts
σs2σv1v1 and σs2v1, so that the 1-D identity σs2σv1v1 =σs2v1 cannot hold.
The 2-D solutions are comparable to solutions in the limit of 1-D flow. This
limiting solution is obtained by letting the longitudinal velocity componentsv˜1, (v1)0
be constant and setting transverse velocities to zero. The limiting solution matches
a 1-D solution computed on the coarse grid. The result along a longitudinal slice
is shown in Fig. 4.6. (The 1-D result is diffused compared to Fig. 4.1 because of
the coarse grid.) The nearly uniform 2-D flow scenario chosen gives results nearly
identical to thosein1-Dinspiteofthelocalized velocity correlation structurein2-D.
Theadditionofsmalllineardiffusionterms,whichareemployedinthenumerical
schemetostabilizethesolution,doesnoteliminatebimodality. Thefactthatwesolve
on a coarse grid, introducing large numerical diffusion, makes it even more striking
that bimodal solutions are clearly evident in the solution. Grid refinement and more
accurate numerical methods can be expected to emphasize the wave separation even
more.
MOMENTEQUATIONSFOR2-DIMMISCIBLEFLOW 9
(a) Mean saturation (b) Variance
0.07
0.8
0.06
0.8 0.8
0.05
0.6
0.6 0.6
0.04
0.4 0.03
0.4 0.4
0.02
0.2
0.2 0.2
0.01
0 0
0.5 1 1.5 0.5 1 1.5
(c) Mean saturation (d) Variance
0.07
0.8
0.06
0.8 0.8
0.7
0.05
0.6
0.6 0.6
0.5 0.04
0.4 0.03
0.4 0.4
0.3
0.02
0.2
0.2 0.2
0.01
0.1
0
0.5 1 1.5 0.5 1 1.5
Fig. 4.3. Saturation mean and variance contours at (a,b) early and (c,d) late times. Lighter
shades indicate a higher mean water saturation or higher variance.
The separation between wave fronts increases with time. In many specific appli-
cationstherelevanttimescalemaybesuchthataseparationisnotevident; however,
even when diffusion dominates initially, bimodality will appear over sufficiently long
times. This property of the solution must be accounted for in any use of MDEs for
advection and diffusion, for both single-phase and multiphase applications.
4.2. Comparison to MCS. MonteCarlosimulationsarefrequentlyusedinap-
plications to subsurface flow and transport and for reservoir modeling. The primary
drawback is that MCS is computationally costly. Separate and very important issues
include the problem of fair sampling of the true sample space and the standard prob-
lemofautocorrelationinrandomnumbergenerators. However,itisawell-established
method and serves as a standard for comparison.
For a bound on the error |s¯(x,t)−(cid:4)s(cid:5)(x,t)| of just 10−2, we estimate that about
4000simulationsareneeded. Thisisasignificantburdenontimeandspaceresources,
so we present results of only 350 simulations. This gives an estimated absolute error
bound of about 0.032, which is sufficient for a rough comparison (recall that mean
saturation values are restricted to [0,1]).
Local correlation structure is again evident in MCS velocity moments in Fig. 4.7
(compare to Fig. 4.2). Sample mean andan envelope ofonestandard deviation along
the midline x2 =L2/2 are shown in Fig. 4.8. Contours of sample moments are shown
in Fig. 4.9. The roughness of the contours highlights the fact that a greater number
of simulations is required to obtain high accuracy. The mean front is near x1 = 1.3,