# Conditional Stochastic Simulations of Flow and  Transport with Karhunen-Loève Expansions, Stochastic  Collocation, and Sequential Gaussian Simulation

Authors: Mina E. Ossiander; Malgorzata Peszynska; Veronika S. Vasylkivska
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We derive a new method of conditional Karhunen-Loève (KL) expansions for stochastic coefficients in models of flow and transport in the subsurface, and in particular for the heterogeneous random permeability field. Exact values of this field are never known, and thus one must evaluate uncertainty of solutions to the flow and transport models. This is typically done by constructing independent realizations of the permeability field followed by numerical simulations of flow and transport for each realization and assembling statistical estimates of moments of desired quantities of interest. We follow the well-known framework of KL expansions and derive a new method that incorporates known values of the permeability at given locations so that the realizations of the permeability field honor this data exactly. Our method relies on projections to an appropriate subspace of random weights applied to the eigenfunctions of the covariance operator. We use the permeability realizations constructed with our stochastic simulation method in simulations of flow and transport and compare the results to those obtained when realizations are constructed with sequential Gaussian simulation (SGS). We also compare efficiency and stochastic convergence with that of stochastic collocation.


---

## Body

## 1. Introduction

Computational modeling of flow and transport in the subsurface requires detailed knowledge of coefficients of partial differential equations (PDEs), in particular of permeabilitiesK and porosities Φ. These are heterogeneous; that is, K=K(x) where x∈D and D⊂Rd is the domain of flow. However, the actual values of K(x) are never known exactly in a real subsurface reservoir; thus the values of K(x) can only be inferred from partial measurements at well locations, say x*=(x1*,…,xNm*)∈D, or from some other information such as seismic and other geological analyses. Therefore the PDEs are stochastic and a typical set of simulations of flow and transport involves geostatistical realizations of K(x) to account for the associated uncertainty of this data, the solutions to the PDE, and any associated quantities of interest Q. The randomness (uncertainty) of K(x) is denoted, as usual, by K(x,ω) where ω is a random element of a probability space.The use of Karhunen-Loève (KL) expansions, a close relative of Principal Component Analysis (PCA), for parametrization of a random field, is well known in computational mathematics and engineering community [1] and has been also applied in mathematical geosciences [2]. Assume that K(x)=K(x,ω) is a Gaussian log-normal field with known covariance function C(x,y). By design, KL expansions account for spatial variability of K(x) through a sequence of eigenfunctions of C smooth in x. Since these correspond to rapidly decreasing eigenvalues, one typically truncates the KL expansion to N terms in KN(x,ω) which capture the desired proportion of the variance of the field. Finally, KL expansions are independent of the resolution intended for x in the sense that a realization of KN(x,ω) can be constructed with arbitrary resolution.The framework of KL expansions discussed in [1, 2] does not however take into account existing point measurement data at locations x*. Though the variability, or “noise,” around a given mean K(x)¯ expected of a reservoir [2] is indicated, observed data are not honored.In this paper we construct (truncated) conditional Karhunen-Loève (KL) expansionsKN(x)~ for the permeability field K(x) with a given analytical spatial covariance model C(x,y) and observations K(x*). Our method uses the point data K(x*) in the construction of random weights applied to eigenfunctions. This gives a mean background field K(x)¯¯ which is essentially a kriged field based on K(x*). Also, the expansions KN(x)~-K(x)¯¯ provide random fluctuations which are zero at the observed locations. Consequently, all realizations KN(x)~ agree exactly with the point data: KN(x*)~=K(x*).Before KL expansions became popular, other geostatistics methods, for example, sequential Gaussian simulation (SGS), were used to construct conditional or unconditional realizations ofK(x). Thus we compare our expansions to those obtained with SGS implemented in a well known geostatistics package [3]. SGS uses the same covariance model and conditioning points as our method and outputs a set of independent realizations of K(x,ω) of desired size and at a prescribed resolution. We use next the permeability realizations obtained with our method and with SGS as data for a deterministic saturated Darcy flow and single component advection-diffusion transport solver. While the flow results have been published by many authors for unconditional stochastic simulations, the transport results coupled to the flow are less frequently considered and bring interesting insights.Our comparison of KL and SGS is of theoretical and practical value and we discuss various aspects of these two entirely different family of methods. Finally, we consider stochastic collocation, a highly accurate method of computing statistical moments of quantities of interest, and discuss its use in the context of conditional representations.Our approach gives a new method of using KL expansions to simulate conditioned random fields. An alternative approach to conditioning ofK(x) was considered in [4]. In order to honor the observed values, they compute the eigenfunction and eigenvalue pairs of the conditioned process. This results in an ensemble of eigenpairs that depend on the particular locations. In this paper we project the random coefficients of the KL expansion onto the appropriate subspace, so that the resulting realizations depend on the particular locations. This allows us to work with the eigenfunctions and eigenvalues derived for the unconditioned process. It also gives a simple decomposition of the simulated process in terms of the conditional mean (kriged values) and the fluctuations around that mean. Other stochastic numerical methods for flow and transport equations include the contributions in [5–9]; these consider very different avenues from direct stochastic parametrizations and will not be reviewed here.The analytical covariance models used in this paper are the ubiquitous exponential model and the less commonly considered Gaussian model. These stationary models depend only on two-point correlations. We do not address nonstationary field, experimental covariance, multipoint geostatistics or physical models beyond linear flow and transport. These important extensions will be discussed elsewhere, while we refer to [10] for some work in these directions.The plan of the paper is as follows. In Section2 we recall the flow and transport models which use field realizations of the random field K(x) as their input. In Section 3 we show how realizations of the field Y(x)=ln⁡K(x) are generated with a known covariance function CY(x,y) together with point measurement data given at x*. Section 3.1 develops the main theoretical result, the conditional KL expansions for a generic random field Y(x). We also recall in Section 3.2 how SGS creates realizations based on the same information. In Section 4 we present simulation results of the numerical model using realizations of K(x) obtained with KL and with SGS, paying particular attention to the use of stochastic collocation method with data based on prior measurements.


## 2. Stochastic Flow and Transport Models

In this section we make precise the flow and transport models and their numerical discretizations. These are well known [11, 12] but are provided for completeness together with appropriate boundary and initial conditions. The corresponding numerical discretization of these is also standard [12, 13]; details for the formally defined weak formulations of flow and transport are provided in Appendix D. We follow the structure of stochastic models established, for example, in [1, 2, 14].LetD∈Rd be the open bounded domain in Rd where the flow and transport take place. In principle, one can consider d=1,2,3 but our examples involve d=2. The time variable is denoted by t. The boundary ∂D of D is assumed to be smooth so that, in particular, the outer unit normal n(x) to ∂D is well defined. The boundary ∂D is split twofold into disjoint subboundaries depending on the prescribed boundary conditions for the flow ∂D=∂DD∪∂DN and transport ∂D=∂DI∪∂DO. These are made precise below.The regionD has associated permeability field K and porosity field ϕ. In this paper we assume that K=K(x,ω) is a random heterogeneous field and ϕ is constant. For simplicity we also assume that K(x) is isotropic. Furthermore, the fluid and the medium have an associated diffusion-dispersion coefficient D, but we ignore the dispersion by setting D equal to molecular diffusivity; that is, D=const=dmI. Since dm correlates with porosity, it makes sense to assume D is constant. The case of random ϕ, D will not be considered in this paper. In comparisons, we also consider no diffusion; that is, D=0.In the flow problem we seek fluid velocityu and its pressure p. Since K is random, so are p and u, as given by Doob-Dynkin Lemma [15]. The stochastic model of flow combines momentum (Darcy’s law) with mass conservation as follows:
(1a)u(x,ω)=-K(x,ω)∇p(x,ω),x∈D,(1b)∇·u(x,ω)=f(x),x∈D.
This flow model is complemented by boundary conditions of Dirichlet and Neumann type, respectively, prescribed on ∂DD and ∂DN as follows:
(1c)p(x,ω)=gD(x),x∈∂DD,(1d)u(x,ω)·n(x)=0,x∈∂DN.In (1a) above f(x) are the fluid sources and sinks due to, for example, the presence of wells, henceforth assumed absent.Once the flow model is solved, the pressuresp and fluxes u are known. Now the solute is subject to advection with velocity u and to diffusion. The transport model is solved for the solute concentration c and is a mass conservation combined with Fick’s law, with q denoting the solute source/sink term. Consider
(2a)ϕ∂c(x,t,ω)∂t+∇·[c(x,t,ω)u(x,ω)]-∇·[D∇c(x,t,ω)]=q(x),x∈D,t>0,
subject to the following initial conditions
(2b)c(x,0,ω)=c0(x),x∈D.
The boundary conditions are imposed on the inflow and outflow boundaries with ∂DI denoting the inflow boundary where u·n<0 and ∂DO the outflow boundary where u·n≥0. These are prescribed as
(2c)c(x,t,ω)=cI(x,t),(x,t)∈∂DI×(0,T],(2d)D(x)∇c(x,t,ω)·n(x)=0,x∈∂DO,t>0.An important quantity of interest is the average breakthrough curve
(3)BTC(t,ω)=∫∂DO‍c(x,t,ω)[u(x,ω)·n(x)]dx,
which represents the total amount of the substance leaving the region D. It is best plotted against another time dependent quantity, the pore volume injected which in the case of fully saturated flow grows simply linearly with time
(4)PVI(t,ω)=t∫∂DI‍[u(x,ω)·n(x)]dx.
Finally, the breakthrough times BTT are
(5)BTT(ω;b)≔sup⁡{t:BTC(t,ω)=b},
where b is some desired value of the breakthrough curve.


### 2.1. Numerical Solution for Flow and Transport

GivenK(x,ω) and boundary and initial data we can solve the flow ((1a), (1b), (1c), (1d)) and the transport models ((2a), (2b), (2c), (2d)) numerically. Let Th be a partition of the domain D into nonoverlapping rectangular elements Ei, where h denotes the largest diameter of elements Ei from Th. Also let 0=t0<t1<⋯<tn=T be a given partition of the time interval [0,T] and Δtn=tn-tn-1. (We also use the subscript n for random sequences in Section 3, but as they appear in separate contexts there is no risk of confusion.) We seek numerically the approximations ph(x,ω),uh(x,ω) to p(x,ω),u(x,ω), and Chn(x,ω) to c(x,tn,ω).To define these approximations, we use standard numerical methods for which natural stochastic extensions have either been already established [1, 2], or are straightforward. For the flow problem we employ cell-centered finite differences (CCFD). These are equivalent to mixed finite elements RT[0] of the lowest Raviart-Thomas order on rectangles [16]. For the transport we employ CCFD and upwinding, as well as implicit Backward Euler time discretization. These numerical methods are well known and are first order accurate in time and space. They are also superconvergent in some norms, and can be modified to obtain higher accuracy, but we do not discuss numerical error. See [2] for a recent formal extension of CCFD to stochastic methods and error estimates for flow. See also Appendix D for the formal weak setting of these methods for flow and transport, extending [13] to the stochastic case.See Figure4 for typical results of flow and transport simulations.


## 2.1. Numerical Solution for Flow and Transport

GivenK(x,ω) and boundary and initial data we can solve the flow ((1a), (1b), (1c), (1d)) and the transport models ((2a), (2b), (2c), (2d)) numerically. Let Th be a partition of the domain D into nonoverlapping rectangular elements Ei, where h denotes the largest diameter of elements Ei from Th. Also let 0=t0<t1<⋯<tn=T be a given partition of the time interval [0,T] and Δtn=tn-tn-1. (We also use the subscript n for random sequences in Section 3, but as they appear in separate contexts there is no risk of confusion.) We seek numerically the approximations ph(x,ω),uh(x,ω) to p(x,ω),u(x,ω), and Chn(x,ω) to c(x,tn,ω).To define these approximations, we use standard numerical methods for which natural stochastic extensions have either been already established [1, 2], or are straightforward. For the flow problem we employ cell-centered finite differences (CCFD). These are equivalent to mixed finite elements RT[0] of the lowest Raviart-Thomas order on rectangles [16]. For the transport we employ CCFD and upwinding, as well as implicit Backward Euler time discretization. These numerical methods are well known and are first order accurate in time and space. They are also superconvergent in some norms, and can be modified to obtain higher accuracy, but we do not discuss numerical error. See [2] for a recent formal extension of CCFD to stochastic methods and error estimates for flow. See also Appendix D for the formal weak setting of these methods for flow and transport, extending [13] to the stochastic case.See Figure4 for typical results of flow and transport simulations.


## 3. Methods for Generating Realizations of K(x,ω) and Computing Moments of Quantities of Interest Q(·,ω)

In this section we describe how to generate independent realizations of the random heterogeneous permeability fieldK(x,ω),x∈D. Here ω is a random element in an underlying probability space. Most of the crucial information about K comes from measurements and its spatial structure with ω indexing the uncertainty of the information. In physical simulations of flow and transport we consider M independent realizations K(x,ω1),…,K(x,ωM) with which we assess statistical moments of the simulation results; that is, some quantities of interest Q which depend on the (ωj)’s and possibly on x,t. For example, Q could be pointwise values of pressures or some averages of fluxes or concentrations. In what follows we will denote this dependence by Q(·,ωj) where the “·” accounts for any relevant nonstochastic variables or parameters.We assume the following probabilistic information aboutK(x,ω): it is a stationary field, its logarithm
(6)Y(x,ω)=ln⁡K(x,ω),
has a normal distribution, with known bounded (analytical) covariance function
(7)CY(x,y)=Cov[Y(x),Y(y)].We also assume thatNm measurements Y1,Y2,…,YNm of Y are given at Nm points x*≔(x1*,x2*,…,xNm*) so that
(8)Y(xi*)=Yi,i=1,…,Nm,
is known. If no measurements are available, this is denoted by Nm=0.Both KL expansions and SGS require a model for (7) and can work with either Nm=0 or Nm>0. In the examples reported in this paper we assume that the covariance is stationary, so CY(x,y)=CY(∥x-y∥). We also assume CY is a quickly decaying function, exponential or Gaussian, given, respectively, as
(9)CY(x,y)=σ2exp⁡(-∥x-y∥η),(10)CY(x,y)=σ2exp⁡(-∥x-y∥2η2),
where the fixed parameters σ2 and η, representing respectively the stationary variance (or sill) and the correlation length, give the assumed second-order spatial properties of the field.KL expansions additionally require that the mean field,(11)Y(x)¯=E[Y](x),
is known. In practice, it appears that (11) provides a highly variable background and a realistic looking field to which the spatially distributed “noise” constructed by unconditional KL expansions is added. However, this does not allow the realizations of the field K(x) to satisfy the prescribed given values exactly. Thus the unconditional expansions investigated in [1, 2] could not be compared to conditional Gaussian fields.Below we describe how to getK(x,ω1),…,K(x,ωM) using KL and SGS. Next we show how the moments of Q(·,ω) are calculated.


### 3.1. Karhunen-Loève (KL) Representation and Realizations ofK(x,ω)

First we consider the well known case whenNm=0, next we develop Nm>0. By Mercer’s theorem [17] the spectral decomposition of the positive definite covariance matrix is
(12)CY(x,y)=∑n=1∞‍λnϕn(x)ϕn(y),
where λn are positive eigenvalues and ϕn(x) are mutually orthogonal eigenfunctions, that is, solutions to the Fredholm integral equation of the second kind [14]
(13)λϕ(x)=∫D‍CY(x,y)ϕ(y)dy,
which satisfy ∫D‍ϕn(x)ϕm(x)dx=δn,m with δn,m being the Kronecker-delta.In what follows we denote byΛ the infinitely dimensional diagonal matrix of eigenvalues λn, n≥0. By convention, the eigenvalues are arranged in Λ from the largest to smallest. Also, we denote by Φ(x) the vector of eigenfunctions evaluated at x. One can thus rewrite (12) as CY(x,y)=ΦT(y)ΛΦ(x). We denote the upper left diagonal block of Λ of size N×N by ΛN and also define ΦN(x) as the finite length vector corresponding to the first N eigenfunctions.


#### 3.1.1. Unconditional Representation WhenNm=0

The KL expansion [18] of Y is well known
(14)Y(x,ω)=Y(x)¯+∑n=1∞‍λnϕn(x)ξn(ω),
where {ξn(ω)}n=1∞ are independent standard Gaussian (N(0,1)) random variables. In other words, (14) represents the fluctuation
(15)Y(x,ω)-Y(x)¯=ΦT(x)ΛΞ(ω),
where Ξ(ω)={ξn(ω)}n=1∞.In practice (14) is truncated to the first N terms
(16)Y(x,ω)≈YN(x,ω)=Y(x)¯+∑n=1N‍λnϕn(x)ξn(ω),
that is, YN(x,ω)-Y(x)¯=ΦNT(x)ΛNΞN. N is chosen so that most of the “energy” or “volatility” represented by the sum of eigenvalues in ΛN is captured to the desired accuracy; see Section 3.1.3.


#### 3.1.2. Conditional Representation WhenNm>0

In order forY(x) to have a given measurement value Y(xi*) at every point xi* in x*, we consider the associated conditional process Y(x)~. Our main idea is that in the representation of Y(x)~ we use the original eigenfunctions listed in Φ(x) but we project the random coefficients in Ξ onto the appropriate subspace. In other words, we replace the original ξn by ξ~n in (14). The variable {ξ~n} has the same distribution as the original random sequence {ξn} after being conditioned on X*≔σ({Y(x*)}), the σ-algebra generated by the Nm measurements Y(x*).DefineΣ∈RNm×Nm to be the covariance matrix corresponding to the observed locations
(17)Σi,j≔CY(xi*,xj*)=∑n=1∞‍λnϕn(xi*)ϕn(xj*),i,j=1,…,Nm,
and R to be the matrix with nth row given by the values of the nth eigenfunction at the observed locations x*(18)ϕn(x*)=(ϕn(x1*),ϕn(x2*),…,ϕn(xNm*)).
This gives Σ=(Σij)=RTΛR. We assume that Σ is full rank. Notice that for every n and xi*,
(19)Cov[ξn,Y(xi*)]=λnϕn(xi*).
Using the usual formula for conditional means and covariances of Gaussian random variables [19], we compute the mean μ~ and covariance M=(mn,k) of the sequence of random variables {ξn}n=1∞ conditioned on X*(20)μ~n=E[ξn∣X*]=λnϕn(x*)Σ-1(Y(x*)-Y(x*)¯),(21)mn,k=Cov[ξn,ξk∣X*]=δn,k-λkϕk(x*)Σ-1ϕn(x*)λn.
Thus M=I-Λ1/2RΣ-1RTΛ1/2, and it is easy to show M=M2; that is, M is a projection matrix. (See Appendix A for details.) Finally, the sequence {ξn}n=1∞ conditioned on X* has the same distribution as
(22)ξ~=μ~+Mη,
where {ηk}k=1∞ is a sequence of i.i.d. N(0,1) random variables. In particular M projects η onto the subspace that gives the Y~ process conditional variance 0 at the locations x1*,…,xNm*. We can thus state the conditional representation Y(x)~ of Y(x)(23)Y(x,ω)~=Y(x)¯+∑n=1∞‍λnϕn(x)ξ~n(ω),x∈D,
which modifies (14).The truncation atN finite of the expansion (23) written componentwise is
(24)YN(x,ω)~=Y(x)¯+∑n=1N‍λnϕn(x)ξ~N,n(ω),x∈D,
where the {ξ~N,n}n=1N are the random coefficients for the truncated expansion constructed from the analogues μN~ and MN of μ~ and M. In particular Σ is replaced by ΣN, the covariance matrix (also assumed full rank) of YN(x*). (For more details, see Appendix A.)It should be noted thatMN projects the N dimensional vector ΞN onto the N-Nm dimensional subspace corresponding to YN(x*)~=Y(x*). In order to have a sufficient number of degrees of freedom for the projection, we must then have
(25)N>Nm.Obviously, the elements of the sequence{ξ~n}n=1N are not independent. Hence their joint probability density function, which is needed later in moment calculations, is not a product. While this can be handled as in [1], it is much more convenient to work directly with the ηj’s which are independent and consequently have a product density. Combining (24) with the truncated version of (22) we get
(26)YN(x)~=Y(x)¯+∑n=1N‍λnϕn(x)(μ~N,n+(MNη)n)=(Y(x)¯+∑n=1N‍λnϕn(x)μ~N,n)+∑n=1N‍λnϕn(x)(∑j=1N‍(MN)n,jηj).
The practical expression (26) is used in moment calculations using stochastic collocation discussed in the sequel. We also see that the first term in (26)
(27)YN(x)¯¯≔Y(x)¯+∑n=1N‍λnϕn(x)μ~N,n=E[YN(x)∣X*],
is essentially deterministic, depending on the Yi but not on the ηj’s. It can be readily interpreted as the kriged mean of Y(x) which honors the known data. The second term in (26) is associated with MNη and forces the fluctuations of YN to be zero at the x* locations.


#### 3.1.3. Truncation Error and Conditioning

The truncationYN of Y reduces total variability as follows:
(28)∫D‍Var⁡[YN(x)]dx=∫D‍∑n≤N‍λnϕn2(x)dx=∑n≤N‍λn≤∑n≥1‍λn=∫D‍Var⁡[Y(x)]dx.
The conditioned process has, conditionally, reduced total variability in the following sense:
(29)∫D‍Var⁡[YN(x)~∣X*]dx=∑n≤N‍λnm¯n,
where m¯n=(MN)n,n∈[0,1]. (See Appendix A for computations.) It should be stressed however that the totalunconditional variability of the truncated conditioned process matches that of the truncated YN process; in particular for any N≥1,
(30)∫D‍Var⁡[YN(x)~]dx=∑n≤N‍λn.
As a consequence, the conditioned process replicates the irregularity and other distributional properties of the unconditioned process; compare [17]. Again, see Appendix A for details.


#### 3.1.4. Example of Conditional Representation with KL (24)

A given covariance functionCY(x,y), in principle, uniquely determines the spectra (λn,ϕn)n=1∞ of CY(x,y) via (13). The analytical solutions of (13) are known only for a few selected stationary separable models of CY(x,y) [20, 21], and even these require some numerical calculations. In general, a numerical solution of (13) must be obtained via, for example, collocation or Galerkin methods [22]. We use the collocation method to obtain the spectra for covariances corresponding to a few variograms common in geostatistics [23, Chapter 4] and in particular, for the exponential, Gaussian, spherical, and hole variograms common in geostatistics [3]. Examples of the first two are given in Figure 1 along with the plot of their eigenvalues.(a) plot of exponential (9) and Gaussian (10) covariance models as function of the lag ∥x-y∥, with η=0.3 and σ=1. (b) first 15 eigenvalues of their KL representations.
(a)
(b)We provide now an example illustrating conditional realizations using KL ind=1. Let D=(0,1) and let the Gaussian covariance CY(x,y) be given by (10) with η=0.3, σ=1. We consider the truncated conditional series expansion (24) with N=100 terms and Nm=3 data points x*=(0.2,0.6,0.8) and Y(x*)=(1,0,-1). We consider the background mean Y(x)¯=0. In Figure 2 we plot the realizations and illustrate how they agree with Y(x*), and how they are decomposed into the conditioned mean Y(x)¯¯ and the fluctuations Y(x)~-Y(x)¯¯.Figure 2
Illustration of a conditional representationY(x)~ of a random one-dimensional field on [0,1] with Gaussian covariance function (10) with η=0.3 and σ=1. We use N=100 terms and Nm=3 data points x*=(0.2,0.6,0.8) with Y(x*)=(1,0,-1) Left: first 10 realizations of Y(x,ωk)~, k=1,…,10. Middle: the conditional mean Y(x)¯¯ as in (27). Right: the fluctuations of individual realizations which are added to Y(x)¯¯ to produce Y(x)~.


### 3.2. Sequential Gaussian Simulation SGS

Geostatistical Software Library (GSLIB) is an open source library of routines which allows simulation of realizations of porous media using different variograms (or covariance functions) as a proxy for assumed physical properties [3]. GSLIB is well known in the geostatistical community and SGS is a relatively simple but quite robust simulation method. The realizations of Y(x)~ honor input data at their locations, and the global histogram and variogram are reproduced within ergodic fluctuations. SGS follows a sequence of well defined steps including transformations to and from normal scores and establishing a random path through the locations where the new data is needed.The major differences betweenY simulated with SGS and with KL are that (i) KL realizations are grid independent while realizations with SGS are grid dependent, and (ii) the KL fields are, by construction, a linear combination of basis functions; thus they are smooth if the basis functions are smooth. On the other hand, SGS realizations are not smooth. Smoothness is important in history matching and reservoir optimization [10].


### 3.3. Examples of ConditionalK(x) Constructed with SGS and KL Expansions

Now we show examples of realizations inD=[0,1]2 constructed with SGS and with KL expansions on a 50×50 grid. For both we use the Gaussian covariance model (10) with η=0.2, σ=1.We select as the “true” data a sample field of the same spatial structure generated without conditioning using thesgsim routine in GSLIB 2.0 [3]. See Figure 3 and Appendix B for the GSLIB parameter file. Next we select Nm=30 data locations x* to be used for conditioning, as shown in Figure 3.Figure 3
Far left: Permeability field used as the “true field” (and a source of data) for conditional simulations and the location ofNm=30 points x* (left). Right and far right: realizations of a conditional permeability field based on Nm=30 data points: GSLIB (right) and KL expansion (far right).Figure 4
Typical results of flow and transport simulations for one realizationω1. From left to right: permeability field K(x,ω1), pressure solution p(x,ω1), and concentrations c(x,t,ω1) at different time steps.Simulations with SGS. We generate M=6400 realizations of Y(x)~ conditioned on the “true" data at the selected locations using SGS and follow up with (6) to obtain the M realizations K(·,ω1),…,K(·,ωM) of the permeability field.Simulations with KL. To get KL expansions, we calculate the eigenfunctions and eigenvalues solving (13). The first N=250 eigenvalues contribute 97% of the variance of the logarithm of the permeability field. We further generate M=6400 conditional realizations using the KL expansion (24) with N=250 terms and follow with (6) to obtain the M realizations of the permeability field KN(·,ω1),…,KN(·,ωM).Comparison of SGS and KL. Figure 3 shows x* and two typical conditional realizations produced, respectively, with SGS and KL expansion. As expected, the latter looks substantially smoother than that obtained with SGS.


### 3.4. Computing Moments of Quantities of Interest

The solutions to physical models of flow and transport (described below) use the realizationsK(x,ω) as their data. These physical models compute some quantities of interest Q(·,ω), which may also depend on x, t, denoted by the dot (·). For brevity we skip (·) in Q(·) when no confusion results.The mean and variance ofQ(·) can be computed using the probability density ρ(y) of Q(·,y) via
(31)E[Q(·)]=∫Γ‍Q(·,y)ρ(y)dy,(32)Var⁡[Q(·)]=∫Γ‍(Q(·,y)-E[Q(·)])2ρ(y)dy,
if ρ(y) is known. Alternatively, we can approximate (31) and (32) with the Monte Carlo method. This is especially helpful when the vector y∈Γ is multidimensional. The approximation process is the same for both simulation methods given in Section 3. However, the choice of the method affects the accuracy of the approximation. We provide details below, focusing on the expectation (31).


#### 3.4.1. Monte Carlo with SGS (SGS-MC)

If SGS is used, the realizations,(33)K(x,ω1),…,K(x,ωM),
are independent realizations of the random field K(x), as in the Monte Carlo method. Thus the moments E[Q], Var⁡[Q] are estimated by the sample mean and variance
(34)E[Q(·)]≈EM[Q(·)]=1M∑k=1M‍Q(·,ωk),(35)Var⁡[Q(·)]≈VarM[Q(·)]=1M-1∑k=1M‍(Q(·,ωk)-EM[Q(·)])2.


#### 3.4.2. KL Expansions with Monte Carlo (KL-MC)

For KL expansions, the randomness inK comes from the vector Ξ in (24). When we draw M independent realizations of Ξ (that is of Y) we treat this process as a Monte Carlo simulation. Since the representation of K(x,ω) is truncated to KN(x,ω), we generate
(36)KN(x,ω1),…,KN(x,ωM).
Thus it is appropriate to denote the corresponding quantity of interest as QN. We estimate (31) by the sample mean
(37)E[Q(·)]≈EN,M[Q(·)]≔EM[QN(·)],
where the latter is defined in (34). Var⁡[Q(·)] is analogously estimated by VarN,M[Q(·)].The following point is important. LetΓn=ξn(Ω) denote the image of ξn and Γ=∏n=1N‍Γn⊂RN. When Nm=0, the joint probability density function ρ=∏n=1N‍ρn. In this paper we assume that each ξn is N(0,1); thus ρ is a multivariate normal density. However if Nm>0 the ξ~n’s are not independent and their joint density is not a product of the individual densities. Via the shift expressed in (26) though, we can work with the product density of the independent ηn’s.


#### 3.4.3. KL with Stochastic Collocation (KL-SC)

An approach yielding more accuracy than that of Monte Carlo is the 2-step stochastic collocation (SC) method as given in [1, 14]. First Q(y) is approximated by a (multivariate) tensor product polynomial Qm of degree m≔(m1,…,mN) where mn is the degree of the component polynomial in yn, the stochastic dimension n. This is followed by a highly accurate numerical integration method which is optimal for the polynomial approximation. The polynomial degrees in m and the total number of degrees of freedom Mm≔∏n=1N‍(mn+1) of Qm(y) are chosen to yield the desired degree of accuracy. In practice we choose mk=m uniform in k so that
(38)Mm=(m+1)N.
The details are somewhat involved but standard [1, 2, 14]. We provide some here for the sake of direct comparison with (34); the remaining details are given in Appendix C.The polynomialQm(y) is constructed in a special way intended to make (31) easy to evaluate, while maintaining accuracy. As in [14], Qm is an interpolating polynomial of Q constructed with collocation points {yk}k=1Mm⊂Γ as follows. Let {Lkm(y)}k=1Mm be the Lagrange basis with Lkm(yj)=δk,j, j, k=1,…,Mm. Then
(39)Qm(y)=∑k=1Mm‍Q(yk)Lkm(y),
and we approximate the integral
(40)E[Q]≈E[Qm]=∫Γ‍∑k=1Mm‍Q(yk)Lkm(y)ρ(y)dy=∑k=1Mm‍wkQ(yk),
where wk=∫Γ‍Lkm(y)ρ(y)dy. Details on the polynomials Lk, weights wk, and points yk are in Appendix C. In practice, only the weights wk and points yk are needed for (40), since the approximation (39) is never formed explicitly.


### 3.5. Summary: Stochastic Simulations of the Model Problems

Consider the formulas (34), (37), and (40). We see that the three corresponding methods, respectively, SGS-MC, KL-MC, and KL-SC, require the knowledge of Q(·) evaluated for each realization ωk or at each collocation point yk. While the implementation of all three methods is very similar, their accuracy differs. For the estimation of moments in (31) or (32), SGS-MC and KL-MC use similar uniform weights 1/M≈1/(M-1) since all realizations are equally probable, while KL-SC uses weights wk optimal for accuracy and optimal choice of realizations.We summarize the simulations steps for each of the three methods.


#### 3.5.1. SGS-MC Outline

(1)
Fixh and the grid Th.
(2)
Define a desired numberM∈N+ of realizations of the permeability field.
(3)
GenerateM independent realizations K(·,ωk), k=1,…,M with SGS; these automatically honor the data K(x*).
(4)
For each realization ofK(·,ωk), k=1,…,M, find the approximate solutions to the flow equations ((1a), (1b), (1c), (1d)), and at each time step find the approximate solutions to the transport model ((2a), (2b), (2c), (2d)). Calculate Q(·,ωk) for each k=1,…M.
(5)
ApproximateE[Q] by the sample average (34) and Var⁡[Q] by the sample variance (35).All the realizations generated in Step(4) of SGS-MC are equally probable and oscillate around the mean.


#### 3.5.2. KL-MC Outline

(1)
Decide on the finite numberN of terms in KL expansion (24).
(2)
ObtainΦ,Λ or ΦN,ΛN analytically or numerically. If finding these numerically, choose a spatial grid on D at least as fine as Th chosen below. If found analytically, there is no grid dependence.
(3)
ComputeΣN,RN for the given x* and Φ.
(4)
Fixh and the grid Th.
(5)
Fix a desired numberM∈N+.
(6)
GenerateM realizations of Ξ~(ωk). For each vector Ξ~(ωk) build YN(x,ωk)~ with (24) and calculate KN(x,ωk) with (6).
(7)
(Same as Step(4) from SGS-MC): Calculate QN(·,ωk) for each k=1,…,M.
(8)
ApproximateE[QN] by the sample average (37) and Var⁡[QN] analogously.In contrast to SGS-MC, the realizations of KL-MC in Step(6), while equally probable, can capture the desired stochastic accuracy by varying N.


#### 3.5.3. KL-SC Outline

(1)
Follow steps  (1)–(4) from KL-MC.
(5)
Define a desired order of polynomials with multi-indexm as in Section 3.4.3. (For simplicity, consider (38) with m a small integer.) Set up the collocation points yk, k=1,…,Mm, corresponding to m.
(6)
For each collocation pointyk, build Y(x,yk)~ with (24), and calculate K(x,yk) with (6).
(7)
For each realization ofK(·,yk), k=1,…,Mm, find the solutions to the flow equations ((1a), (1b), (1c), (1d)), and for each tn find the solutions to the transport model ((2a), (2b), (2c), (2d)). Calculate Q(·,yk) for each k=1,…,Mm.
(8)
ApproximateE[Q] by (40) and variance analogously.The realizations of KL-SC in Step(6) are actually not random, but the moments calculated with KL-MC capture a desired stochastic accuracy by varying N and m.


## 3.1. Karhunen-Loève (KL) Representation and Realizations ofK(x,ω)

First we consider the well known case whenNm=0, next we develop Nm>0. By Mercer’s theorem [17] the spectral decomposition of the positive definite covariance matrix is
(12)CY(x,y)=∑n=1∞‍λnϕn(x)ϕn(y),
where λn are positive eigenvalues and ϕn(x) are mutually orthogonal eigenfunctions, that is, solutions to the Fredholm integral equation of the second kind [14]
(13)λϕ(x)=∫D‍CY(x,y)ϕ(y)dy,
which satisfy ∫D‍ϕn(x)ϕm(x)dx=δn,m with δn,m being the Kronecker-delta.In what follows we denote byΛ the infinitely dimensional diagonal matrix of eigenvalues λn, n≥0. By convention, the eigenvalues are arranged in Λ from the largest to smallest. Also, we denote by Φ(x) the vector of eigenfunctions evaluated at x. One can thus rewrite (12) as CY(x,y)=ΦT(y)ΛΦ(x). We denote the upper left diagonal block of Λ of size N×N by ΛN and also define ΦN(x) as the finite length vector corresponding to the first N eigenfunctions.


### 3.1.1. Unconditional Representation WhenNm=0

The KL expansion [18] of Y is well known
(14)Y(x,ω)=Y(x)¯+∑n=1∞‍λnϕn(x)ξn(ω),
where {ξn(ω)}n=1∞ are independent standard Gaussian (N(0,1)) random variables. In other words, (14) represents the fluctuation
(15)Y(x,ω)-Y(x)¯=ΦT(x)ΛΞ(ω),
where Ξ(ω)={ξn(ω)}n=1∞.In practice (14) is truncated to the first N terms
(16)Y(x,ω)≈YN(x,ω)=Y(x)¯+∑n=1N‍λnϕn(x)ξn(ω),
that is, YN(x,ω)-Y(x)¯=ΦNT(x)ΛNΞN. N is chosen so that most of the “energy” or “volatility” represented by the sum of eigenvalues in ΛN is captured to the desired accuracy; see Section 3.1.3.


### 3.1.2. Conditional Representation WhenNm>0

In order forY(x) to have a given measurement value Y(xi*) at every point xi* in x*, we consider the associated conditional process Y(x)~. Our main idea is that in the representation of Y(x)~ we use the original eigenfunctions listed in Φ(x) but we project the random coefficients in Ξ onto the appropriate subspace. In other words, we replace the original ξn by ξ~n in (14). The variable {ξ~n} has the same distribution as the original random sequence {ξn} after being conditioned on X*≔σ({Y(x*)}), the σ-algebra generated by the Nm measurements Y(x*).DefineΣ∈RNm×Nm to be the covariance matrix corresponding to the observed locations
(17)Σi,j≔CY(xi*,xj*)=∑n=1∞‍λnϕn(xi*)ϕn(xj*),i,j=1,…,Nm,
and R to be the matrix with nth row given by the values of the nth eigenfunction at the observed locations x*(18)ϕn(x*)=(ϕn(x1*),ϕn(x2*),…,ϕn(xNm*)).
This gives Σ=(Σij)=RTΛR. We assume that Σ is full rank. Notice that for every n and xi*,
(19)Cov[ξn,Y(xi*)]=λnϕn(xi*).
Using the usual formula for conditional means and covariances of Gaussian random variables [19], we compute the mean μ~ and covariance M=(mn,k) of the sequence of random variables {ξn}n=1∞ conditioned on X*(20)μ~n=E[ξn∣X*]=λnϕn(x*)Σ-1(Y(x*)-Y(x*)¯),(21)mn,k=Cov[ξn,ξk∣X*]=δn,k-λkϕk(x*)Σ-1ϕn(x*)λn.
Thus M=I-Λ1/2RΣ-1RTΛ1/2, and it is easy to show M=M2; that is, M is a projection matrix. (See Appendix A for details.) Finally, the sequence {ξn}n=1∞ conditioned on X* has the same distribution as
(22)ξ~=μ~+Mη,
where {ηk}k=1∞ is a sequence of i.i.d. N(0,1) random variables. In particular M projects η onto the subspace that gives the Y~ process conditional variance 0 at the locations x1*,…,xNm*. We can thus state the conditional representation Y(x)~ of Y(x)(23)Y(x,ω)~=Y(x)¯+∑n=1∞‍λnϕn(x)ξ~n(ω),x∈D,
which modifies (14).The truncation atN finite of the expansion (23) written componentwise is
(24)YN(x,ω)~=Y(x)¯+∑n=1N‍λnϕn(x)ξ~N,n(ω),x∈D,
where the {ξ~N,n}n=1N are the random coefficients for the truncated expansion constructed from the analogues μN~ and MN of μ~ and M. In particular Σ is replaced by ΣN, the covariance matrix (also assumed full rank) of YN(x*). (For more details, see Appendix A.)It should be noted thatMN projects the N dimensional vector ΞN onto the N-Nm dimensional subspace corresponding to YN(x*)~=Y(x*). In order to have a sufficient number of degrees of freedom for the projection, we must then have
(25)N>Nm.Obviously, the elements of the sequence{ξ~n}n=1N are not independent. Hence their joint probability density function, which is needed later in moment calculations, is not a product. While this can be handled as in [1], it is much more convenient to work directly with the ηj’s which are independent and consequently have a product density. Combining (24) with the truncated version of (22) we get
(26)YN(x)~=Y(x)¯+∑n=1N‍λnϕn(x)(μ~N,n+(MNη)n)=(Y(x)¯+∑n=1N‍λnϕn(x)μ~N,n)+∑n=1N‍λnϕn(x)(∑j=1N‍(MN)n,jηj).
The practical expression (26) is used in moment calculations using stochastic collocation discussed in the sequel. We also see that the first term in (26)
(27)YN(x)¯¯≔Y(x)¯+∑n=1N‍λnϕn(x)μ~N,n=E[YN(x)∣X*],
is essentially deterministic, depending on the Yi but not on the ηj’s. It can be readily interpreted as the kriged mean of Y(x) which honors the known data. The second term in (26) is associated with MNη and forces the fluctuations of YN to be zero at the x* locations.


### 3.1.3. Truncation Error and Conditioning

The truncationYN of Y reduces total variability as follows:
(28)∫D‍Var⁡[YN(x)]dx=∫D‍∑n≤N‍λnϕn2(x)dx=∑n≤N‍λn≤∑n≥1‍λn=∫D‍Var⁡[Y(x)]dx.
The conditioned process has, conditionally, reduced total variability in the following sense:
(29)∫D‍Var⁡[YN(x)~∣X*]dx=∑n≤N‍λnm¯n,
where m¯n=(MN)n,n∈[0,1]. (See Appendix A for computations.) It should be stressed however that the totalunconditional variability of the truncated conditioned process matches that of the truncated YN process; in particular for any N≥1,
(30)∫D‍Var⁡[YN(x)~]dx=∑n≤N‍λn.
As a consequence, the conditioned process replicates the irregularity and other distributional properties of the unconditioned process; compare [17]. Again, see Appendix A for details.


### 3.1.4. Example of Conditional Representation with KL (24)

A given covariance functionCY(x,y), in principle, uniquely determines the spectra (λn,ϕn)n=1∞ of CY(x,y) via (13). The analytical solutions of (13) are known only for a few selected stationary separable models of CY(x,y) [20, 21], and even these require some numerical calculations. In general, a numerical solution of (13) must be obtained via, for example, collocation or Galerkin methods [22]. We use the collocation method to obtain the spectra for covariances corresponding to a few variograms common in geostatistics [23, Chapter 4] and in particular, for the exponential, Gaussian, spherical, and hole variograms common in geostatistics [3]. Examples of the first two are given in Figure 1 along with the plot of their eigenvalues.(a) plot of exponential (9) and Gaussian (10) covariance models as function of the lag ∥x-y∥, with η=0.3 and σ=1. (b) first 15 eigenvalues of their KL representations.
(a)
(b)We provide now an example illustrating conditional realizations using KL ind=1. Let D=(0,1) and let the Gaussian covariance CY(x,y) be given by (10) with η=0.3, σ=1. We consider the truncated conditional series expansion (24) with N=100 terms and Nm=3 data points x*=(0.2,0.6,0.8) and Y(x*)=(1,0,-1). We consider the background mean Y(x)¯=0. In Figure 2 we plot the realizations and illustrate how they agree with Y(x*), and how they are decomposed into the conditioned mean Y(x)¯¯ and the fluctuations Y(x)~-Y(x)¯¯.Figure 2
Illustration of a conditional representationY(x)~ of a random one-dimensional field on [0,1] with Gaussian covariance function (10) with η=0.3 and σ=1. We use N=100 terms and Nm=3 data points x*=(0.2,0.6,0.8) with Y(x*)=(1,0,-1) Left: first 10 realizations of Y(x,ωk)~, k=1,…,10. Middle: the conditional mean Y(x)¯¯ as in (27). Right: the fluctuations of individual realizations which are added to Y(x)¯¯ to produce Y(x)~.


## 3.1.1. Unconditional Representation WhenNm=0

The KL expansion [18] of Y is well known
(14)Y(x,ω)=Y(x)¯+∑n=1∞‍λnϕn(x)ξn(ω),
where {ξn(ω)}n=1∞ are independent standard Gaussian (N(0,1)) random variables. In other words, (14) represents the fluctuation
(15)Y(x,ω)-Y(x)¯=ΦT(x)ΛΞ(ω),
where Ξ(ω)={ξn(ω)}n=1∞.In practice (14) is truncated to the first N terms
(16)Y(x,ω)≈YN(x,ω)=Y(x)¯+∑n=1N‍λnϕn(x)ξn(ω),
that is, YN(x,ω)-Y(x)¯=ΦNT(x)ΛNΞN. N is chosen so that most of the “energy” or “volatility” represented by the sum of eigenvalues in ΛN is captured to the desired accuracy; see Section 3.1.3.


## 3.1.2. Conditional Representation WhenNm>0

In order forY(x) to have a given measurement value Y(xi*) at every point xi* in x*, we consider the associated conditional process Y(x)~. Our main idea is that in the representation of Y(x)~ we use the original eigenfunctions listed in Φ(x) but we project the random coefficients in Ξ onto the appropriate subspace. In other words, we replace the original ξn by ξ~n in (14). The variable {ξ~n} has the same distribution as the original random sequence {ξn} after being conditioned on X*≔σ({Y(x*)}), the σ-algebra generated by the Nm measurements Y(x*).DefineΣ∈RNm×Nm to be the covariance matrix corresponding to the observed locations
(17)Σi,j≔CY(xi*,xj*)=∑n=1∞‍λnϕn(xi*)ϕn(xj*),i,j=1,…,Nm,
and R to be the matrix with nth row given by the values of the nth eigenfunction at the observed locations x*(18)ϕn(x*)=(ϕn(x1*),ϕn(x2*),…,ϕn(xNm*)).
This gives Σ=(Σij)=RTΛR. We assume that Σ is full rank. Notice that for every n and xi*,
(19)Cov[ξn,Y(xi*)]=λnϕn(xi*).
Using the usual formula for conditional means and covariances of Gaussian random variables [19], we compute the mean μ~ and covariance M=(mn,k) of the sequence of random variables {ξn}n=1∞ conditioned on X*(20)μ~n=E[ξn∣X*]=λnϕn(x*)Σ-1(Y(x*)-Y(x*)¯),(21)mn,k=Cov[ξn,ξk∣X*]=δn,k-λkϕk(x*)Σ-1ϕn(x*)λn.
Thus M=I-Λ1/2RΣ-1RTΛ1/2, and it is easy to show M=M2; that is, M is a projection matrix. (See Appendix A for details.) Finally, the sequence {ξn}n=1∞ conditioned on X* has the same distribution as
(22)ξ~=μ~+Mη,
where {ηk}k=1∞ is a sequence of i.i.d. N(0,1) random variables. In particular M projects η onto the subspace that gives the Y~ process conditional variance 0 at the locations x1*,…,xNm*. We can thus state the conditional representation Y(x)~ of Y(x)(23)Y(x,ω)~=Y(x)¯+∑n=1∞‍λnϕn(x)ξ~n(ω),x∈D,
which modifies (14).The truncation atN finite of the expansion (23) written componentwise is
(24)YN(x,ω)~=Y(x)¯+∑n=1N‍λnϕn(x)ξ~N,n(ω),x∈D,
where the {ξ~N,n}n=1N are the random coefficients for the truncated expansion constructed from the analogues μN~ and MN of μ~ and M. In particular Σ is replaced by ΣN, the covariance matrix (also assumed full rank) of YN(x*). (For more details, see Appendix A.)It should be noted thatMN projects the N dimensional vector ΞN onto the N-Nm dimensional subspace corresponding to YN(x*)~=Y(x*). In order to have a sufficient number of degrees of freedom for the projection, we must then have
(25)N>Nm.Obviously, the elements of the sequence{ξ~n}n=1N are not independent. Hence their joint probability density function, which is needed later in moment calculations, is not a product. While this can be handled as in [1], it is much more convenient to work directly with the ηj’s which are independent and consequently have a product density. Combining (24) with the truncated version of (22) we get
(26)YN(x)~=Y(x)¯+∑n=1N‍λnϕn(x)(μ~N,n+(MNη)n)=(Y(x)¯+∑n=1N‍λnϕn(x)μ~N,n)+∑n=1N‍λnϕn(x)(∑j=1N‍(MN)n,jηj).
The practical expression (26) is used in moment calculations using stochastic collocation discussed in the sequel. We also see that the first term in (26)
(27)YN(x)¯¯≔Y(x)¯+∑n=1N‍λnϕn(x)μ~N,n=E[YN(x)∣X*],
is essentially deterministic, depending on the Yi but not on the ηj’s. It can be readily interpreted as the kriged mean of Y(x) which honors the known data. The second term in (26) is associated with MNη and forces the fluctuations of YN to be zero at the x* locations.


## 3.1.3. Truncation Error and Conditioning

The truncationYN of Y reduces total variability as follows:
(28)∫D‍Var⁡[YN(x)]dx=∫D‍∑n≤N‍λnϕn2(x)dx=∑n≤N‍λn≤∑n≥1‍λn=∫D‍Var⁡[Y(x)]dx.
The conditioned process has, conditionally, reduced total variability in the following sense:
(29)∫D‍Var⁡[YN(x)~∣X*]dx=∑n≤N‍λnm¯n,
where m¯n=(MN)n,n∈[0,1]. (See Appendix A for computations.) It should be stressed however that the totalunconditional variability of the truncated conditioned process matches that of the truncated YN process; in particular for any N≥1,
(30)∫D‍Var⁡[YN(x)~]dx=∑n≤N‍λn.
As a consequence, the conditioned process replicates the irregularity and other distributional properties of the unconditioned process; compare [17]. Again, see Appendix A for details.


## 3.1.4. Example of Conditional Representation with KL (24)

A given covariance functionCY(x,y), in principle, uniquely determines the spectra (λn,ϕn)n=1∞ of CY(x,y) via (13). The analytical solutions of (13) are known only for a few selected stationary separable models of CY(x,y) [20, 21], and even these require some numerical calculations. In general, a numerical solution of (13) must be obtained via, for example, collocation or Galerkin methods [22]. We use the collocation method to obtain the spectra for covariances corresponding to a few variograms common in geostatistics [23, Chapter 4] and in particular, for the exponential, Gaussian, spherical, and hole variograms common in geostatistics [3]. Examples of the first two are given in Figure 1 along with the plot of their eigenvalues.(a) plot of exponential (9) and Gaussian (10) covariance models as function of the lag ∥x-y∥, with η=0.3 and σ=1. (b) first 15 eigenvalues of their KL representations.
(a)
(b)We provide now an example illustrating conditional realizations using KL ind=1. Let D=(0,1) and let the Gaussian covariance CY(x,y) be given by (10) with η=0.3, σ=1. We consider the truncated conditional series expansion (24) with N=100 terms and Nm=3 data points x*=(0.2,0.6,0.8) and Y(x*)=(1,0,-1). We consider the background mean Y(x)¯=0. In Figure 2 we plot the realizations and illustrate how they agree with Y(x*), and how they are decomposed into the conditioned mean Y(x)¯¯ and the fluctuations Y(x)~-Y(x)¯¯.Figure 2
Illustration of a conditional representationY(x)~ of a random one-dimensional field on [0,1] with Gaussian covariance function (10) with η=0.3 and σ=1. We use N=100 terms and Nm=3 data points x*=(0.2,0.6,0.8) with Y(x*)=(1,0,-1) Left: first 10 realizations of Y(x,ωk)~, k=1,…,10. Middle: the conditional mean Y(x)¯¯ as in (27). Right: the fluctuations of individual realizations which are added to Y(x)¯¯ to produce Y(x)~.


## 3.2. Sequential Gaussian Simulation SGS

Geostatistical Software Library (GSLIB) is an open source library of routines which allows simulation of realizations of porous media using different variograms (or covariance functions) as a proxy for assumed physical properties [3]. GSLIB is well known in the geostatistical community and SGS is a relatively simple but quite robust simulation method. The realizations of Y(x)~ honor input data at their locations, and the global histogram and variogram are reproduced within ergodic fluctuations. SGS follows a sequence of well defined steps including transformations to and from normal scores and establishing a random path through the locations where the new data is needed.The major differences betweenY simulated with SGS and with KL are that (i) KL realizations are grid independent while realizations with SGS are grid dependent, and (ii) the KL fields are, by construction, a linear combination of basis functions; thus they are smooth if the basis functions are smooth. On the other hand, SGS realizations are not smooth. Smoothness is important in history matching and reservoir optimization [10].


## 3.3. Examples of ConditionalK(x) Constructed with SGS and KL Expansions

Now we show examples of realizations inD=[0,1]2 constructed with SGS and with KL expansions on a 50×50 grid. For both we use the Gaussian covariance model (10) with η=0.2, σ=1.We select as the “true” data a sample field of the same spatial structure generated without conditioning using thesgsim routine in GSLIB 2.0 [3]. See Figure 3 and Appendix B for the GSLIB parameter file. Next we select Nm=30 data locations x* to be used for conditioning, as shown in Figure 3.Figure 3
Far left: Permeability field used as the “true field” (and a source of data) for conditional simulations and the location ofNm=30 points x* (left). Right and far right: realizations of a conditional permeability field based on Nm=30 data points: GSLIB (right) and KL expansion (far right).Figure 4
Typical results of flow and transport simulations for one realizationω1. From left to right: permeability field K(x,ω1), pressure solution p(x,ω1), and concentrations c(x,t,ω1) at different time steps.Simulations with SGS. We generate M=6400 realizations of Y(x)~ conditioned on the “true" data at the selected locations using SGS and follow up with (6) to obtain the M realizations K(·,ω1),…,K(·,ωM) of the permeability field.Simulations with KL. To get KL expansions, we calculate the eigenfunctions and eigenvalues solving (13). The first N=250 eigenvalues contribute 97% of the variance of the logarithm of the permeability field. We further generate M=6400 conditional realizations using the KL expansion (24) with N=250 terms and follow with (6) to obtain the M realizations of the permeability field KN(·,ω1),…,KN(·,ωM).Comparison of SGS and KL. Figure 3 shows x* and two typical conditional realizations produced, respectively, with SGS and KL expansion. As expected, the latter looks substantially smoother than that obtained with SGS.


## 3.4. Computing Moments of Quantities of Interest

The solutions to physical models of flow and transport (described below) use the realizationsK(x,ω) as their data. These physical models compute some quantities of interest Q(·,ω), which may also depend on x, t, denoted by the dot (·). For brevity we skip (·) in Q(·) when no confusion results.The mean and variance ofQ(·) can be computed using the probability density ρ(y) of Q(·,y) via
(31)E[Q(·)]=∫Γ‍Q(·,y)ρ(y)dy,(32)Var⁡[Q(·)]=∫Γ‍(Q(·,y)-E[Q(·)])2ρ(y)dy,
if ρ(y) is known. Alternatively, we can approximate (31) and (32) with the Monte Carlo method. This is especially helpful when the vector y∈Γ is multidimensional. The approximation process is the same for both simulation methods given in Section 3. However, the choice of the method affects the accuracy of the approximation. We provide details below, focusing on the expectation (31).


### 3.4.1. Monte Carlo with SGS (SGS-MC)

If SGS is used, the realizations,(33)K(x,ω1),…,K(x,ωM),
are independent realizations of the random field K(x), as in the Monte Carlo method. Thus the moments E[Q], Var⁡[Q] are estimated by the sample mean and variance
(34)E[Q(·)]≈EM[Q(·)]=1M∑k=1M‍Q(·,ωk),(35)Var⁡[Q(·)]≈VarM[Q(·)]=1M-1∑k=1M‍(Q(·,ωk)-EM[Q(·)])2.


### 3.4.2. KL Expansions with Monte Carlo (KL-MC)

For KL expansions, the randomness inK comes from the vector Ξ in (24). When we draw M independent realizations of Ξ (that is of Y) we treat this process as a Monte Carlo simulation. Since the representation of K(x,ω) is truncated to KN(x,ω), we generate
(36)KN(x,ω1),…,KN(x,ωM).
Thus it is appropriate to denote the corresponding quantity of interest as QN. We estimate (31) by the sample mean
(37)E[Q(·)]≈EN,M[Q(·)]≔EM[QN(·)],
where the latter is defined in (34). Var⁡[Q(·)] is analogously estimated by VarN,M[Q(·)].The following point is important. LetΓn=ξn(Ω) denote the image of ξn and Γ=∏n=1N‍Γn⊂RN. When Nm=0, the joint probability density function ρ=∏n=1N‍ρn. In this paper we assume that each ξn is N(0,1); thus ρ is a multivariate normal density. However if Nm>0 the ξ~n’s are not independent and their joint density is not a product of the individual densities. Via the shift expressed in (26) though, we can work with the product density of the independent ηn’s.


### 3.4.3. KL with Stochastic Collocation (KL-SC)

An approach yielding more accuracy than that of Monte Carlo is the 2-step stochastic collocation (SC) method as given in [1, 14]. First Q(y) is approximated by a (multivariate) tensor product polynomial Qm of degree m≔(m1,…,mN) where mn is the degree of the component polynomial in yn, the stochastic dimension n. This is followed by a highly accurate numerical integration method which is optimal for the polynomial approximation. The polynomial degrees in m and the total number of degrees of freedom Mm≔∏n=1N‍(mn+1) of Qm(y) are chosen to yield the desired degree of accuracy. In practice we choose mk=m uniform in k so that
(38)Mm=(m+1)N.
The details are somewhat involved but standard [1, 2, 14]. We provide some here for the sake of direct comparison with (34); the remaining details are given in Appendix C.The polynomialQm(y) is constructed in a special way intended to make (31) easy to evaluate, while maintaining accuracy. As in [14], Qm is an interpolating polynomial of Q constructed with collocation points {yk}k=1Mm⊂Γ as follows. Let {Lkm(y)}k=1Mm be the Lagrange basis with Lkm(yj)=δk,j, j, k=1,…,Mm. Then
(39)Qm(y)=∑k=1Mm‍Q(yk)Lkm(y),
and we approximate the integral
(40)E[Q]≈E[Qm]=∫Γ‍∑k=1Mm‍Q(yk)Lkm(y)ρ(y)dy=∑k=1Mm‍wkQ(yk),
where wk=∫Γ‍Lkm(y)ρ(y)dy. Details on the polynomials Lk, weights wk, and points yk are in Appendix C. In practice, only the weights wk and points yk are needed for (40), since the approximation (39) is never formed explicitly.


## 3.4.1. Monte Carlo with SGS (SGS-MC)

If SGS is used, the realizations,(33)K(x,ω1),…,K(x,ωM),
are independent realizations of the random field K(x), as in the Monte Carlo method. Thus the moments E[Q], Var⁡[Q] are estimated by the sample mean and variance
(34)E[Q(·)]≈EM[Q(·)]=1M∑k=1M‍Q(·,ωk),(35)Var⁡[Q(·)]≈VarM[Q(·)]=1M-1∑k=1M‍(Q(·,ωk)-EM[Q(·)])2.


## 3.4.2. KL Expansions with Monte Carlo (KL-MC)

For KL expansions, the randomness inK comes from the vector Ξ in (24). When we draw M independent realizations of Ξ (that is of Y) we treat this process as a Monte Carlo simulation. Since the representation of K(x,ω) is truncated to KN(x,ω), we generate
(36)KN(x,ω1),…,KN(x,ωM).
Thus it is appropriate to denote the corresponding quantity of interest as QN. We estimate (31) by the sample mean
(37)E[Q(·)]≈EN,M[Q(·)]≔EM[QN(·)],
where the latter is defined in (34). Var⁡[Q(·)] is analogously estimated by VarN,M[Q(·)].The following point is important. LetΓn=ξn(Ω) denote the image of ξn and Γ=∏n=1N‍Γn⊂RN. When Nm=0, the joint probability density function ρ=∏n=1N‍ρn. In this paper we assume that each ξn is N(0,1); thus ρ is a multivariate normal density. However if Nm>0 the ξ~n’s are not independent and their joint density is not a product of the individual densities. Via the shift expressed in (26) though, we can work with the product density of the independent ηn’s.


## 3.4.3. KL with Stochastic Collocation (KL-SC)

An approach yielding more accuracy than that of Monte Carlo is the 2-step stochastic collocation (SC) method as given in [1, 14]. First Q(y) is approximated by a (multivariate) tensor product polynomial Qm of degree m≔(m1,…,mN) where mn is the degree of the component polynomial in yn, the stochastic dimension n. This is followed by a highly accurate numerical integration method which is optimal for the polynomial approximation. The polynomial degrees in m and the total number of degrees of freedom Mm≔∏n=1N‍(mn+1) of Qm(y) are chosen to yield the desired degree of accuracy. In practice we choose mk=m uniform in k so that
(38)Mm=(m+1)N.
The details are somewhat involved but standard [1, 2, 14]. We provide some here for the sake of direct comparison with (34); the remaining details are given in Appendix C.The polynomialQm(y) is constructed in a special way intended to make (31) easy to evaluate, while maintaining accuracy. As in [14], Qm is an interpolating polynomial of Q constructed with collocation points {yk}k=1Mm⊂Γ as follows. Let {Lkm(y)}k=1Mm be the Lagrange basis with Lkm(yj)=δk,j, j, k=1,…,Mm. Then
(39)Qm(y)=∑k=1Mm‍Q(yk)Lkm(y),
and we approximate the integral
(40)E[Q]≈E[Qm]=∫Γ‍∑k=1Mm‍Q(yk)Lkm(y)ρ(y)dy=∑k=1Mm‍wkQ(yk),
where wk=∫Γ‍Lkm(y)ρ(y)dy. Details on the polynomials Lk, weights wk, and points yk are in Appendix C. In practice, only the weights wk and points yk are needed for (40), since the approximation (39) is never formed explicitly.


## 3.5. Summary: Stochastic Simulations of the Model Problems

Consider the formulas (34), (37), and (40). We see that the three corresponding methods, respectively, SGS-MC, KL-MC, and KL-SC, require the knowledge of Q(·) evaluated for each realization ωk or at each collocation point yk. While the implementation of all three methods is very similar, their accuracy differs. For the estimation of moments in (31) or (32), SGS-MC and KL-MC use similar uniform weights 1/M≈1/(M-1) since all realizations are equally probable, while KL-SC uses weights wk optimal for accuracy and optimal choice of realizations.We summarize the simulations steps for each of the three methods.


### 3.5.1. SGS-MC Outline

(1)
Fixh and the grid Th.
(2)
Define a desired numberM∈N+ of realizations of the permeability field.
(3)
GenerateM independent realizations K(·,ωk), k=1,…,M with SGS; these automatically honor the data K(x*).
(4)
For each realization ofK(·,ωk), k=1,…,M, find the approximate solutions to the flow equations ((1a), (1b), (1c), (1d)), and at each time step find the approximate solutions to the transport model ((2a), (2b), (2c), (2d)). Calculate Q(·,ωk) for each k=1,…M.
(5)
ApproximateE[Q] by the sample average (34) and Var⁡[Q] by the sample variance (35).All the realizations generated in Step(4) of SGS-MC are equally probable and oscillate around the mean.


### 3.5.2. KL-MC Outline

(1)
Decide on the finite numberN of terms in KL expansion (24).
(2)
ObtainΦ,Λ or ΦN,ΛN analytically or numerically. If finding these numerically, choose a spatial grid on D at least as fine as Th chosen below. If found analytically, there is no grid dependence.
(3)
ComputeΣN,RN for the given x* and Φ.
(4)
Fixh and the grid Th.
(5)
Fix a desired numberM∈N+.
(6)
GenerateM realizations of Ξ~(ωk). For each vector Ξ~(ωk) build YN(x,ωk)~ with (24) and calculate KN(x,ωk) with (6).
(7)
(Same as Step(4) from SGS-MC): Calculate QN(·,ωk) for each k=1,…,M.
(8)
ApproximateE[QN] by the sample average (37) and Var⁡[QN] analogously.In contrast to SGS-MC, the realizations of KL-MC in Step(6), while equally probable, can capture the desired stochastic accuracy by varying N.


### 3.5.3. KL-SC Outline

(1)
Follow steps  (1)–(4) from KL-MC.
(5)
Define a desired order of polynomials with multi-indexm as in Section 3.4.3. (For simplicity, consider (38) with m a small integer.) Set up the collocation points yk, k=1,…,Mm, corresponding to m.
(6)
For each collocation pointyk, build Y(x,yk)~ with (24), and calculate K(x,yk) with (6).
(7)
For each realization ofK(·,yk), k=1,…,Mm, find the solutions to the flow equations ((1a), (1b), (1c), (1d)), and for each tn find the solutions to the transport model ((2a), (2b), (2c), (2d)). Calculate Q(·,yk) for each k=1,…,Mm.
(8)
ApproximateE[Q] by (40) and variance analogously.The realizations of KL-SC in Step(6) are actually not random, but the moments calculated with KL-MC capture a desired stochastic accuracy by varying N and m.


## 3.5.1. SGS-MC Outline

(1)
Fixh and the grid Th.
(2)
Define a desired numberM∈N+ of realizations of the permeability field.
(3)
GenerateM independent realizations K(·,ωk), k=1,…,M with SGS; these automatically honor the data K(x*).
(4)
For each realization ofK(·,ωk), k=1,…,M, find the approximate solutions to the flow equations ((1a), (1b), (1c), (1d)), and at each time step find the approximate solutions to the transport model ((2a), (2b), (2c), (2d)). Calculate Q(·,ωk) for each k=1,…M.
(5)
ApproximateE[Q] by the sample average (34) and Var⁡[Q] by the sample variance (35).All the realizations generated in Step(4) of SGS-MC are equally probable and oscillate around the mean.


## 3.5.2. KL-MC Outline

(1)
Decide on the finite numberN of terms in KL expansion (24).
(2)
ObtainΦ,Λ or ΦN,ΛN analytically or numerically. If finding these numerically, choose a spatial grid on D at least as fine as Th chosen below. If found analytically, there is no grid dependence.
(3)
ComputeΣN,RN for the given x* and Φ.
(4)
Fixh and the grid Th.
(5)
Fix a desired numberM∈N+.
(6)
GenerateM realizations of Ξ~(ωk). For each vector Ξ~(ωk) build YN(x,ωk)~ with (24) and calculate KN(x,ωk) with (6).
(7)
(Same as Step(4) from SGS-MC): Calculate QN(·,ωk) for each k=1,…,M.
(8)
ApproximateE[QN] by the sample average (37) and Var⁡[QN] analogously.In contrast to SGS-MC, the realizations of KL-MC in Step(6), while equally probable, can capture the desired stochastic accuracy by varying N.


## 3.5.3. KL-SC Outline

(1)
Follow steps  (1)–(4) from KL-MC.
(5)
Define a desired order of polynomials with multi-indexm as in Section 3.4.3. (For simplicity, consider (38) with m a small integer.) Set up the collocation points yk, k=1,…,Mm, corresponding to m.
(6)
For each collocation pointyk, build Y(x,yk)~ with (24), and calculate K(x,yk) with (6).
(7)
For each realization ofK(·,yk), k=1,…,Mm, find the solutions to the flow equations ((1a), (1b), (1c), (1d)), and for each tn find the solutions to the transport model ((2a), (2b), (2c), (2d)). Calculate Q(·,yk) for each k=1,…,Mm.
(8)
ApproximateE[Q] by (40) and variance analogously.The realizations of KL-SC in Step(6) are actually not random, but the moments calculated with KL-MC capture a desired stochastic accuracy by varying N and m.


## 4. Results of Stochastic Simulations of Flow and Transport

Now we present numerical experiments and discuss the properties of associated moments of the quantities of interestQ(·,ω).TakeD=[0,1]2 covered with a uniform 50×50 spatial grid Th using the Gaussian (10) and exponential (9) covariance models with K(x,ω) computed as discussed in Section 3. Different choices of the number Nm of conditioning points are used; we also vary the number of terms N in the KL expansions (24) and the number of realizations M, as well as the polynomial order m in KL-SC. In all simulations we assume for simplicity that the background Y(x)¯≡0. This is of course easy to modify for practical simulations.Boundary conditions for the flow and transport models are the same for all examples, and so are the homogeneous source/sink terms(41)f(x)=0,q(x)=0,x∈D,
and thus the flow and transport are driven only by boundary data.For the flow model (1a)–(1c) we use the conditions
(42)p(x)=1,x∈{(x1,x2)∈∂D:x1=0},(43)p(x)=0,x∈{(x1,x2)∈∂D:x1=1},(44)u(x)·n(x)=0,x∈{(x1,x2)∈∂D:x2=0orx2=1}.
These boundary conditions mean that the prevailing direction of the flow is from left to right. Thus the left boundary is also the inflow, and the right boundary is the outflow boundary, and the bottom and top boundaries are the no-flow boundaries.For the transport model (2a), (2b)–(2d) we use porosity ϕ(x)=1 and either diffusion D=0 or D=dm=0.0252I. We also define the following initial and boundary conditions
(45)c(x,0)=1,x∈D,(46)c(x,t)=0,(x,t)∈{(x1,x2)∈∂D:x1=0}×(0,T],(47)D(x)∇c(x,t)·n(x)=0,(x,t)∈{(x1,x2)∈∂D:x1≠0}×(0,T].
The condition (47) imposes the no-flow condition on top and bottom boundaries, as well as the usual numericaloutflow boundary for D≠0 on the right (outflow) boundary of D. If D=0, no boundary condition is imposed on the outflow boundary. Simulations are run for T sufficiently long so that BTC(T)≈0.An illustration of the typical behavior of the flow and transport solutions is given in Figure4. As we show below, some of the cumulative quantities of interest associated with the transport solutions have a somewhat different behavior than the pointwise quantities for the flow.Our emphasis in the examples below is threefold. We compare SGS-MC and KL-MC methods in Section4.1. In Section 4.2 we consider the impact of conditioning on the KL-MC method on transport results. In Section 4.3 we evaluate stochastic convergence for KL-MC versus that for KL-SC.In what follows we droph when referring to numerical results.


### 4.1. Comparison of SGS-MC and KL-MC Methods

This section gives comparisons of the flow and transport results usingM=6400 realizations of K(x,ω1),…,K(x,ωM), that is, Monte Carlo simulations corresponding to the example in Section 3.3. Our focus is on the quality of the approximations to the moments of the quantities of interest E[Q(·)], Var⁡[Q(·)]. For conditioning we use Nm=30 points; see examples in Figure 3. For Q(·) we choose either pointwise values of pressure p(x) or fluxes u(x) or breakthrough curves BTC(t). The approximations to their moments are EM[Q(·)] and EN,M[Q(·)] as defined in (34) and (37), respectively. See results in Figure 5.Conditional simulations with SGS-MC (a) and KL-MC (b) withM=6400, Nm=30.
(a)
(b)The first observation is that, as expected, the structure of results corresponding to SGS-MC and KL-MC is very similar. This is also true about unconditional realizations but we do not show these here (see [23]).Next the flow results are examined with a bit more detail; in particular the dependence of the quality of the approximations to moments asM increases is shown. Figures 6 and 7 show the statistical moments of solution of the flow problem along the profile x2=0.5 approximated by unconditional and conditional simulations with SGS-MC and KL-MC, respectively. We see that, as expected, conditional moments of the solutions converge faster than unconditional ones to the asymptotic values for large M. More significantly, the fluxes appear to settle earlier (around M=800) for KL-MC than for SGS-MC (around M=1200).Profiles obtained with SGS-MC of moments of pressures and fluxes for unconditional fieldNm=0 (a) and conditional field with Nm=30 (b).
(a)
(b)Profiles obtained with KL-MC of moments of pressures and fluxes for unconditional fieldNm=0 (a) and conditional field with Nm=30 (b).
(a)
(b)Next we discuss transport results shown in Figures8 and 9. The breakthrough time distributions are qualitatively similar for both GSLIB and KL based experiments. Unconditional breakthrough times have a more skewed distribution and a wider spread in comparison with conditional breakthrough times. As expected, the variances converge faster when Nm>0 than when Nm=0. However, in the averaged quantities of interest shown in Figure 8, we do not see a dramatic difference between SGS-MC and KL-MC, but the variance of BTT appears smaller for KL-MC than for SGS-MC. On the other hand, Figure 9 seems to suggest that both variance and expectation estimates settle down more quickly for KL-MC than for SGS-MC.Transport results with SGS-MC (a) and KL-MC (b) withM=1600. From left to right: moments of breakthrough curves BTC  (PVI), and histograms of breakthrough times BTT.
(a)
(b)Figure 9
Example from Section4.1. Convergence of expectations and variances of breakthrough times BTT depending on the number of realizations M obtained with SGS-MC and KL-MC. Each figure shows results for Nm=0 and Nm=30.In summary, KL-MC seems to require half as many realizations than SGS-MC for pointwise quantities of interest. This is not as strongly exhibited for those transient quantities of interest that are averaged over space such as breakthrough curves. More simulations are needed over various correlation lengths and other parameters to determine the relative benefits of KL-MC versus SGS-MC.


### 4.2. Stochastic Convergence of Moments Depending on the NumberNm of Conditioning Points

In this example we consider KL-MC only and assess the dependence onNm of the convergence of moments of various Q(·) associated with the flow and transport. This example is motivated by [24] where experiments of this type were performed with a nonlinear flow and transport model and SGS.We use a nonseparable exponential covariance model (9) with η=0.3, σ=1, and N=200 and generate M=Mmax⁡ realizations of K(x,ωj); j=1,…,M with KL-MC as described in Section 3.1. First we consider the unconditional field with Nm=0. For conditional simulations we use one arbitrarily chosen unconditional realization as the “true data” providing a source of measurements. We first select a set x* of size Nm=12 for conditioning locations and then a superset x* of size Nm=25 containing the smaller set. We generate Mmax⁡ conditional realizations for all choices of Nm and investigate the effect on the moments of Q(·).Table1 presents convergence of the moments of the flow solutions to ((1a), (1b), (1c), (1d)) as M increases. The convergence is assessed experimentally by comparison of the quantity of interest QN,M with that for QN,Mmax⁡; here we use Mmax⁡=25600. We compute the norm ∥QN,M(·)-QN,Mmax⁡(·)∥ appropriate for the quantity of interest. For pressures it is the discrete l2 norm at the cell centers; for fluxes it is the discrete l2 norm of the values at the midpoints of the edges. The results exhibit the usual decrease in the stochastic error which essentially follows the usual O(1/M) ratio. It is important to notice that the magnitude of the stochastic error seems to weakly decrease with the number Nm of conditioning points.Table 1
Stochastic convergence results for KL-MC simulations withMmax⁡=25600 in the pressures and fluxes. Convergence ratios are given in parentheses.
M
Pressure error
Fluxu1 error
Fluxu2 error
Results forNm=0
100
0.871395160535
0.065990195314
0.028573904583
400
0.364829041492 (2.39)
0.069916447041 (0.94)
0.018699389288 (1.52)
1600
0.055858569054 (6.53)
0.021215209498 (3.29)
0.006586753979 (2.83)
6400
0.054169889340 (1.03)
0.012411413073 (1.71)
0.003615613202 (1.82)
Results forNm=12
100
0.254380891374
0.049476177143
0.026290307316
400
0.112240817113 (2.27)
0.024602883474 (2.01)
0.011912256200 (2.21)
1600
0.064059879872 (1.75)
0.011673912035 (2.11)
0.004237191244 (2.81)
6400
0.050133900679 (1.28)
0.008381613088 (1.39)
0.003806688304 (1.11)
Results forNm=25
100
0.441427542393
0.047070576489
0.022927539745
400
0.204776531658 (2.16)
0.025904577051 (1.82)
0.014762469928 (1.55)
1600
0.093540699576 (2.19)
0.015060205886 (1.72)
0.006659195033 (2.22)
6400
0.035279177226 (2.65)
0.005514284413 (2.73)
0.002559067034 (2.60)The transport model ((2a), (2b), (2c), (2d)) has a considerably larger complexity than the flow model ((1a), (1b), (1c), (1d)) since the transient solutions for all time steps of ((2a), (2b), (2c), (2d)) must be obtained with the size of the time step chosen to satisfy CFL condition. In addition, any comparisons at fixed time steps may be difficult because time stepping varies from realization to realization due to a large variability of the fluxes. Therefore the use of large M is not feasible.Instead, we useM=1600 and present a comparison of cumulative rather than pointwise quantities of interest for Nm=0 and Nm=25. Figure 10 shows the breakthrough curves and breakthrough times for conditional and unconditional simulations. We see that the spread of the conditional breakthrough curves is considerably smaller than that for Nm=0. We also observe a significant reduction in the variances of the breakthrough curves after conditioning. As expected, we see that the presence of diffusion D≠0 when compared to the case D=0 changes the shape of the breakthrough curves and affects the distribution of breakthrough times.Example from Section4.2. In each row we show expectations (far left) and variances (left) approximated with M=1600 unconditional and conditional realizations and histograms of unconditional (right) and conditional (far right) breakthrough times. The rows correspond to the cases with D=0 (a) and D≠0 (b).
(a)
(b)


### 4.3. Comparison of KL-MC and KL-SC

It is well known that KL-SC has a faster convergence rate compared to KL-MC as shown in for example [1, 2] for unconditional simulations of stationary fields. We provide results showing the same behavior for conditional simulations of coupled flow and transport.Here the Gaussian covariance model (10) with parameters η=0.45, σ=1, and N=6 terms is used for KN(x). (In this case, the first N=6 eigenvalues of the covariance matrix CY contribute around 70% of the variance of Y). For conditional simulations we use one of the unconditional realizations as the “true data” and a source of measurements (see Figure 11). Then Nm=3 conditioning locations x* are selected, and realizations of the permeability field are generated based on this data at the collocation points yk, k=1,…,Mm. Note that for KL-SC we need to use a relatively small N and thus a small number of conditioning points Nm. Figure 11 shows the construction. Due to the small value of N it is only possible to simulate large scale variability. This is in contrast to the scale of variability illustrated in Figure 3.Figure 11
Example from Section4.3. Far left: permeability field used as a source of data for conditional simulations with KL-SC. Left: location of Nm=3 data points x*. Right and far right: two realizations of conditional permeability field for KL-SC.Next we perform flow and transport simulations and compute moments of quantities of interest. We provide quantitative analysis but omit plots of flow results due to the similarity with those given in Sections4.1 and 4.2. The results in Tables 2 and 3 demonstrate stochastic convergence for KL-MC and KL-SC methods, respectively. In parentheses we put the ratios between the errors on successive levels of refinement. These are to be understood as follows. For KL-MC we use as a “true solution” the mean calculated over Mmax⁡ realizations, and we calculate the error. Thus for KL-MC we see a decrease of the error as M increases, as to be expected, but the decrease is rather slow. On the other hand, for KL-SC, we choose as the “true solution” (or “reference solution”) that calculated with Mmax⁡ corresponding to 66=46656. Recall that here we are increasing the order of the polynomial used for approximations. We observe very fast decrease of the error with the polynomial order m and M≡mN.Table 2
Stochastic convergence for KL-MC method,Mmax⁡= 102,400. Convergence ratios are given in parentheses.
M
Pressure error
Fluxu1 error
Fluxu2 error
Results forNm=0
100
4.48455775
E
-
1
1.48350640
E
-
1
1.17914638
E
-
2
400
2.66856506
E
-
1 (1.68)
9.76989490
E
-
2 (1.52)
1.33512891
E
-
2 (0.88)
1600
6.38750772
E
-
2 (4.18)
2.85173883
E
-
2 (3.43)
3.17813228
E
-
3 (4.20)
6400
4.62308181
E
-
2 (1.38)
6.09263880
E
-
3 (4.68)
1.03454966
E
-
3 (3.07)
25600
1.58598120
E
-
2 (2.91)
3.18129278
E
-
3 (1.92)
7.24256686
E
-
4 (1.43)
Results forNm=3
100
3.17329908
E
-
1
9.52829319
E
-
2
8.51644272
E
-
3
400
2.01750784
E
-
1 (1.57)
6.54737794
E
-
2 (1.46)
1.70397643
E
-
2 (0.49)
1600
1.30910195
E
-
1 (1.54)
2.52870559
E
-
2 (2.59)
3.99852809
E
-
3 (4.26)
6400
1.85236630
E
-
2 (7.07)
6.47435184
E
-
3 (3.91)
1.84532439
E
-
3 (2.17)
25600
8.95884920
E
-
3 (2.07)
1.06851043
E
-
3 (6.06)
3.53995511
E
-
4 (5.21)Table 3
Stochastic convergence for KL-SC method,Mmax⁡=66= 46,656. Convergence ratios are given in parentheses.
M
m
Pressure error
Fluxu1 error
Fluxu2 error
Results forNm=0
2
6
=
64
1.73171578
E
-
2
1.53781262
E
-
2
1.19104228
E
-
3
3
6
=
729
1.35136013
E
-
3 (12.81)
2.64723907
E
-
4 (58.09)
2.63058974
E
-
5 (45.28)
4
6
=
4096
1.35036246
E
-
4 (10.01)
4.24765346
E
-
5 (6.23)
6.34172897
E
-
6 (4.15)
5
6
=
15625
2.38436555
E
-
5 (5.66)
7.24388553
E
-
6 (5.86)
1.10585528
E
-
6 (5.73)
Results forNm=3
2
6
=
64
1.14545082
E
-
2
2.01435814
E
-
3
2.58675228
E
-
4
3
6
=
729
6.84148724
E
-
4 (16.74)
9.72102804
E
-
5 (20.72)
1.53448331
E
-
5 (16.86)
4
6
=
4096
4.82400762
E
-
5 (14.18)
7.09604144
E
-
6 (13.69)
1.46771202
E
-
6 (10.45)
5
6
=
15625
4.40509288
E
-
6 (10.95)
7.67343049
E
-
7 (9.25)
2.27047888
E
-
7 (6.46)It is then important to compare KL-MC and KL-SC to evaluate the accuracy of the moments of quantities of interest corresponding to the same order of computational effort. For example, compare the errors reported forM=6400 and that for Mm=46, that is, fourth and third rows in Tables 2 and 3 for KL-MC and KL-SC, respectively. These correspond roughly to the same order of computational effort, but the error for KL-SC is about two orders of magnitude smaller.Next we discuss transport results illustrated in Figure12. In each plot we compare the moments for various selected values of M against those for Mmax⁡=6400. First, we see that the statistical moments approximated by KL-MC and KL-SC methods seem to agree visually with each other and appear close to each other for large M. Since the qualitative nature of expectations of breakthrough curves does not differ significantly from those reported in Section 4.2, we skip detailed discussion. We only mention that the reduction in variance is about fivefold between Nm=0 and Nm=3, for both KL-MC and KL-SC. In addition, KL-MC seems to overpredict the variance in BTC while KL-SC underpredicts it, at least for Nm=0.Moments of the breakthrough curves approximated by KL-MC method (a) and KL-SC (b) withN=6 and varying M and m, respectively.
(a)
(b)These qualitative observations are consistent with the stochastic error estimates shown in Table4. The superiority of KL-SC over KL-MC is confirmed, however it is not as dramatic as for pointwise quantities of interest for the stationary (flow) problems, as we seem to be gaining approximately only one order of magnitude accuracy between KL-MC and KL-SC.Table 4
Stochastic convergence for transport results unconditional KL-MC method. Errors are reported against the solution obtained withM=6400 MC realizations.
Results for KL-MC withMmax⁡=6400
M
BTC error,Nm=0
BTC error,Nm=3
100
3.22185276
E
-
5
2.34897907
E
-
5
400
8.45299874
E
-
5
4.45654307
E
-
5
1600
2.37960024
E
-
5
2.48353941
E
-
5
Results for KL-SC withMmax⁡=56=15625
M
m
BTC errorr,Nm=0
BTC error,Nm=3
2
6
=
64
2.80396825
E
-
5
3.89977320
E
-
6
3
6
=
729
2.07419462
E
-
6
7.03030755
E
-
7
4
6
=
4096
1.04005211
E
-
6
6.82138851
E
-
7


## 4.1. Comparison of SGS-MC and KL-MC Methods

This section gives comparisons of the flow and transport results usingM=6400 realizations of K(x,ω1),…,K(x,ωM), that is, Monte Carlo simulations corresponding to the example in Section 3.3. Our focus is on the quality of the approximations to the moments of the quantities of interest E[Q(·)], Var⁡[Q(·)]. For conditioning we use Nm=30 points; see examples in Figure 3. For Q(·) we choose either pointwise values of pressure p(x) or fluxes u(x) or breakthrough curves BTC(t). The approximations to their moments are EM[Q(·)] and EN,M[Q(·)] as defined in (34) and (37), respectively. See results in Figure 5.Conditional simulations with SGS-MC (a) and KL-MC (b) withM=6400, Nm=30.
(a)
(b)The first observation is that, as expected, the structure of results corresponding to SGS-MC and KL-MC is very similar. This is also true about unconditional realizations but we do not show these here (see [23]).Next the flow results are examined with a bit more detail; in particular the dependence of the quality of the approximations to moments asM increases is shown. Figures 6 and 7 show the statistical moments of solution of the flow problem along the profile x2=0.5 approximated by unconditional and conditional simulations with SGS-MC and KL-MC, respectively. We see that, as expected, conditional moments of the solutions converge faster than unconditional ones to the asymptotic values for large M. More significantly, the fluxes appear to settle earlier (around M=800) for KL-MC than for SGS-MC (around M=1200).Profiles obtained with SGS-MC of moments of pressures and fluxes for unconditional fieldNm=0 (a) and conditional field with Nm=30 (b).
(a)
(b)Profiles obtained with KL-MC of moments of pressures and fluxes for unconditional fieldNm=0 (a) and conditional field with Nm=30 (b).
(a)
(b)Next we discuss transport results shown in Figures8 and 9. The breakthrough time distributions are qualitatively similar for both GSLIB and KL based experiments. Unconditional breakthrough times have a more skewed distribution and a wider spread in comparison with conditional breakthrough times. As expected, the variances converge faster when Nm>0 than when Nm=0. However, in the averaged quantities of interest shown in Figure 8, we do not see a dramatic difference between SGS-MC and KL-MC, but the variance of BTT appears smaller for KL-MC than for SGS-MC. On the other hand, Figure 9 seems to suggest that both variance and expectation estimates settle down more quickly for KL-MC than for SGS-MC.Transport results with SGS-MC (a) and KL-MC (b) withM=1600. From left to right: moments of breakthrough curves BTC  (PVI), and histograms of breakthrough times BTT.
(a)
(b)Figure 9
Example from Section4.1. Convergence of expectations and variances of breakthrough times BTT depending on the number of realizations M obtained with SGS-MC and KL-MC. Each figure shows results for Nm=0 and Nm=30.In summary, KL-MC seems to require half as many realizations than SGS-MC for pointwise quantities of interest. This is not as strongly exhibited for those transient quantities of interest that are averaged over space such as breakthrough curves. More simulations are needed over various correlation lengths and other parameters to determine the relative benefits of KL-MC versus SGS-MC.


## 4.2. Stochastic Convergence of Moments Depending on the NumberNm of Conditioning Points

In this example we consider KL-MC only and assess the dependence onNm of the convergence of moments of various Q(·) associated with the flow and transport. This example is motivated by [24] where experiments of this type were performed with a nonlinear flow and transport model and SGS.We use a nonseparable exponential covariance model (9) with η=0.3, σ=1, and N=200 and generate M=Mmax⁡ realizations of K(x,ωj); j=1,…,M with KL-MC as described in Section 3.1. First we consider the unconditional field with Nm=0. For conditional simulations we use one arbitrarily chosen unconditional realization as the “true data” providing a source of measurements. We first select a set x* of size Nm=12 for conditioning locations and then a superset x* of size Nm=25 containing the smaller set. We generate Mmax⁡ conditional realizations for all choices of Nm and investigate the effect on the moments of Q(·).Table1 presents convergence of the moments of the flow solutions to ((1a), (1b), (1c), (1d)) as M increases. The convergence is assessed experimentally by comparison of the quantity of interest QN,M with that for QN,Mmax⁡; here we use Mmax⁡=25600. We compute the norm ∥QN,M(·)-QN,Mmax⁡(·)∥ appropriate for the quantity of interest. For pressures it is the discrete l2 norm at the cell centers; for fluxes it is the discrete l2 norm of the values at the midpoints of the edges. The results exhibit the usual decrease in the stochastic error which essentially follows the usual O(1/M) ratio. It is important to notice that the magnitude of the stochastic error seems to weakly decrease with the number Nm of conditioning points.Table 1
Stochastic convergence results for KL-MC simulations withMmax⁡=25600 in the pressures and fluxes. Convergence ratios are given in parentheses.
M
Pressure error
Fluxu1 error
Fluxu2 error
Results forNm=0
100
0.871395160535
0.065990195314
0.028573904583
400
0.364829041492 (2.39)
0.069916447041 (0.94)
0.018699389288 (1.52)
1600
0.055858569054 (6.53)
0.021215209498 (3.29)
0.006586753979 (2.83)
6400
0.054169889340 (1.03)
0.012411413073 (1.71)
0.003615613202 (1.82)
Results forNm=12
100
0.254380891374
0.049476177143
0.026290307316
400
0.112240817113 (2.27)
0.024602883474 (2.01)
0.011912256200 (2.21)
1600
0.064059879872 (1.75)
0.011673912035 (2.11)
0.004237191244 (2.81)
6400
0.050133900679 (1.28)
0.008381613088 (1.39)
0.003806688304 (1.11)
Results forNm=25
100
0.441427542393
0.047070576489
0.022927539745
400
0.204776531658 (2.16)
0.025904577051 (1.82)
0.014762469928 (1.55)
1600
0.093540699576 (2.19)
0.015060205886 (1.72)
0.006659195033 (2.22)
6400
0.035279177226 (2.65)
0.005514284413 (2.73)
0.002559067034 (2.60)The transport model ((2a), (2b), (2c), (2d)) has a considerably larger complexity than the flow model ((1a), (1b), (1c), (1d)) since the transient solutions for all time steps of ((2a), (2b), (2c), (2d)) must be obtained with the size of the time step chosen to satisfy CFL condition. In addition, any comparisons at fixed time steps may be difficult because time stepping varies from realization to realization due to a large variability of the fluxes. Therefore the use of large M is not feasible.Instead, we useM=1600 and present a comparison of cumulative rather than pointwise quantities of interest for Nm=0 and Nm=25. Figure 10 shows the breakthrough curves and breakthrough times for conditional and unconditional simulations. We see that the spread of the conditional breakthrough curves is considerably smaller than that for Nm=0. We also observe a significant reduction in the variances of the breakthrough curves after conditioning. As expected, we see that the presence of diffusion D≠0 when compared to the case D=0 changes the shape of the breakthrough curves and affects the distribution of breakthrough times.Example from Section4.2. In each row we show expectations (far left) and variances (left) approximated with M=1600 unconditional and conditional realizations and histograms of unconditional (right) and conditional (far right) breakthrough times. The rows correspond to the cases with D=0 (a) and D≠0 (b).
(a)
(b)


## 4.3. Comparison of KL-MC and KL-SC

It is well known that KL-SC has a faster convergence rate compared to KL-MC as shown in for example [1, 2] for unconditional simulations of stationary fields. We provide results showing the same behavior for conditional simulations of coupled flow and transport.Here the Gaussian covariance model (10) with parameters η=0.45, σ=1, and N=6 terms is used for KN(x). (In this case, the first N=6 eigenvalues of the covariance matrix CY contribute around 70% of the variance of Y). For conditional simulations we use one of the unconditional realizations as the “true data” and a source of measurements (see Figure 11). Then Nm=3 conditioning locations x* are selected, and realizations of the permeability field are generated based on this data at the collocation points yk, k=1,…,Mm. Note that for KL-SC we need to use a relatively small N and thus a small number of conditioning points Nm. Figure 11 shows the construction. Due to the small value of N it is only possible to simulate large scale variability. This is in contrast to the scale of variability illustrated in Figure 3.Figure 11
Example from Section4.3. Far left: permeability field used as a source of data for conditional simulations with KL-SC. Left: location of Nm=3 data points x*. Right and far right: two realizations of conditional permeability field for KL-SC.Next we perform flow and transport simulations and compute moments of quantities of interest. We provide quantitative analysis but omit plots of flow results due to the similarity with those given in Sections4.1 and 4.2. The results in Tables 2 and 3 demonstrate stochastic convergence for KL-MC and KL-SC methods, respectively. In parentheses we put the ratios between the errors on successive levels of refinement. These are to be understood as follows. For KL-MC we use as a “true solution” the mean calculated over Mmax⁡ realizations, and we calculate the error. Thus for KL-MC we see a decrease of the error as M increases, as to be expected, but the decrease is rather slow. On the other hand, for KL-SC, we choose as the “true solution” (or “reference solution”) that calculated with Mmax⁡ corresponding to 66=46656. Recall that here we are increasing the order of the polynomial used for approximations. We observe very fast decrease of the error with the polynomial order m and M≡mN.Table 2
Stochastic convergence for KL-MC method,Mmax⁡= 102,400. Convergence ratios are given in parentheses.
M
Pressure error
Fluxu1 error
Fluxu2 error
Results forNm=0
100
4.48455775
E
-
1
1.48350640
E
-
1
1.17914638
E
-
2
400
2.66856506
E
-
1 (1.68)
9.76989490
E
-
2 (1.52)
1.33512891
E
-
2 (0.88)
1600
6.38750772
E
-
2 (4.18)
2.85173883
E
-
2 (3.43)
3.17813228
E
-
3 (4.20)
6400
4.62308181
E
-
2 (1.38)
6.09263880
E
-
3 (4.68)
1.03454966
E
-
3 (3.07)
25600
1.58598120
E
-
2 (2.91)
3.18129278
E
-
3 (1.92)
7.24256686
E
-
4 (1.43)
Results forNm=3
100
3.17329908
E
-
1
9.52829319
E
-
2
8.51644272
E
-
3
400
2.01750784
E
-
1 (1.57)
6.54737794
E
-
2 (1.46)
1.70397643
E
-
2 (0.49)
1600
1.30910195
E
-
1 (1.54)
2.52870559
E
-
2 (2.59)
3.99852809
E
-
3 (4.26)
6400
1.85236630
E
-
2 (7.07)
6.47435184
E
-
3 (3.91)
1.84532439
E
-
3 (2.17)
25600
8.95884920
E
-
3 (2.07)
1.06851043
E
-
3 (6.06)
3.53995511
E
-
4 (5.21)Table 3
Stochastic convergence for KL-SC method,Mmax⁡=66= 46,656. Convergence ratios are given in parentheses.
M
m
Pressure error
Fluxu1 error
Fluxu2 error
Results forNm=0
2
6
=
64
1.73171578
E
-
2
1.53781262
E
-
2
1.19104228
E
-
3
3
6
=
729
1.35136013
E
-
3 (12.81)
2.64723907
E
-
4 (58.09)
2.63058974
E
-
5 (45.28)
4
6
=
4096
1.35036246
E
-
4 (10.01)
4.24765346
E
-
5 (6.23)
6.34172897
E
-
6 (4.15)
5
6
=
15625
2.38436555
E
-
5 (5.66)
7.24388553
E
-
6 (5.86)
1.10585528
E
-
6 (5.73)
Results forNm=3
2
6
=
64
1.14545082
E
-
2
2.01435814
E
-
3
2.58675228
E
-
4
3
6
=
729
6.84148724
E
-
4 (16.74)
9.72102804
E
-
5 (20.72)
1.53448331
E
-
5 (16.86)
4
6
=
4096
4.82400762
E
-
5 (14.18)
7.09604144
E
-
6 (13.69)
1.46771202
E
-
6 (10.45)
5
6
=
15625
4.40509288
E
-
6 (10.95)
7.67343049
E
-
7 (9.25)
2.27047888
E
-
7 (6.46)It is then important to compare KL-MC and KL-SC to evaluate the accuracy of the moments of quantities of interest corresponding to the same order of computational effort. For example, compare the errors reported forM=6400 and that for Mm=46, that is, fourth and third rows in Tables 2 and 3 for KL-MC and KL-SC, respectively. These correspond roughly to the same order of computational effort, but the error for KL-SC is about two orders of magnitude smaller.Next we discuss transport results illustrated in Figure12. In each plot we compare the moments for various selected values of M against those for Mmax⁡=6400. First, we see that the statistical moments approximated by KL-MC and KL-SC methods seem to agree visually with each other and appear close to each other for large M. Since the qualitative nature of expectations of breakthrough curves does not differ significantly from those reported in Section 4.2, we skip detailed discussion. We only mention that the reduction in variance is about fivefold between Nm=0 and Nm=3, for both KL-MC and KL-SC. In addition, KL-MC seems to overpredict the variance in BTC while KL-SC underpredicts it, at least for Nm=0.Moments of the breakthrough curves approximated by KL-MC method (a) and KL-SC (b) withN=6 and varying M and m, respectively.
(a)
(b)These qualitative observations are consistent with the stochastic error estimates shown in Table4. The superiority of KL-SC over KL-MC is confirmed, however it is not as dramatic as for pointwise quantities of interest for the stationary (flow) problems, as we seem to be gaining approximately only one order of magnitude accuracy between KL-MC and KL-SC.Table 4
Stochastic convergence for transport results unconditional KL-MC method. Errors are reported against the solution obtained withM=6400 MC realizations.
Results for KL-MC withMmax⁡=6400
M
BTC error,Nm=0
BTC error,Nm=3
100
3.22185276
E
-
5
2.34897907
E
-
5
400
8.45299874
E
-
5
4.45654307
E
-
5
1600
2.37960024
E
-
5
2.48353941
E
-
5
Results for KL-SC withMmax⁡=56=15625
M
m
BTC errorr,Nm=0
BTC error,Nm=3
2
6
=
64
2.80396825
E
-
5
3.89977320
E
-
6
3
6
=
729
2.07419462
E
-
6
7.03030755
E
-
7
4
6
=
4096
1.04005211
E
-
6
6.82138851
E
-
7


## 5. Summary

In this paper we presented a new method of constructing conditional Karhunen-Loève (KL) expansions for the permeability fieldK(x,ω) used as data in simulations of flow and transport to compute certain pointwise and cumulative quantities of interest Q(·,ω). Our new method is important due to its ability to honor given data. It has a natural decomposition giving (i) a kriged mean field agreeing with observed data and (ii) random fluctuations around the kriged mean which are zero at the observed locations; thus all realizations honor the point data.We give details of Monte Carlo implementation KL-MC of our method as well as of its use in stochastic collocation KL-SC. We compare the results of simulations with KL-MC and KL-SC to those of SGS-MC obtained with sequential Gaussian simulation (SGS) as implemented in GSLIB [3].All three methods rely on the same initial information: a covariance modelCY and the data K(x*) at Nm conditioning locations x*. All three can generate any desired number of realizations M; those in KL-SC have a particular structure.The comparison between SGS-MC, KL-MC, and KL-SC can be summarized as follows. When used for calculation of the mean and variance of various quantities of interest, KL-MC seems to converge about twice as fast compared to SGS-MC. In turn, KL-SC converges orders of magnitude faster than KL-MC. These observations depend on the quantities of interest, with the pointwise values seeming to be more sensitive to the method chosen than the cumulative ones. However, sinceN>Nm as in (25), in practice we can use tensor grid KL-SC implementation only with a few conditioning points due to the O(mN) complexity. At this time we thus see KL-MC as the most versatile and accurate technique.The conditional expansions proposed in our method come from simple algebraic modifications of the unconditional expansions. The additional cost associated with the conditioning has complexity equal to that of finding an inverse of a matrix of sizeNm×Nm. In contrast, the expansions in [4] require determining the spectra for each set of conditioning points. As shown in Section 3.5.2, our method can reuse unconditional expansions and vary the number and position x* of conditioning points. This is in constrast to SGS which works with a fixed set x* and a fixed spatial grid or resolution.In order to fully understand the relative merits of the three methods, more testing is needed with different covariance models and parameters and with a variety of conditioning data setsx*. Moreover, different physical models and quantities of interest need to be evaluated. To further improve on KL-MC as presented and the use of the more efficient KL-SC variant, future plans include exploration of nontensor grids in stochastic collocation. Other current work includes numerical error analysis and the use of multipoint statistics for nonstationary covariance models.



---

*Source: 652594-2014-10-15.xml*