Nested sparse grid collocation method with delay and transformation for subsurface flow and transport problems

Accepted Manuscript

Nested sparse grid collocation method with delay and transformation for subsurface flow and transport problems Qinzhuo Liao , Dongxiao Zhang , Hamdi Tchelepi PII: DOI: Reference:

S0309-1708(16)30183-X 10.1016/j.advwatres.2017.03.020 ADWR 2811

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

24 June 2016 22 March 2017 23 March 2017

Please cite this article as: Qinzhuo Liao , Dongxiao Zhang , Hamdi Tchelepi , Nested sparse grid collocation method with delay and transformation for subsurface flow and transport problems, Advances in Water Resources (2017), doi: 10.1016/j.advwatres.2017.03.020

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights 

A delayed Kronrod-Patterson-Hermite nested quadrature rule is proposed, which has the fastest convergence rate compared to other quadrature rules.



A transformed approach is introduced to improve accuracy in case of low regularity or unsmoothness.



The sparse grid collocation method combining the delayed and transformed

CR IP T

processes is tested in various cases including three-dimensional three-phase

AC

CE

PT

ED

M

AN US

flow problems, and provides accurate results in an efficient manner.

1

ACCEPTED MANUSCRIPT

Nested sparse grid collocation method with delay and transformation for subsurface flow and transport problems Qinzhuo Liaoa, Dongxiao Zhanga and Hamdi Tchelepib a

BIC-ESAT, ERE, and SKLTCS, College of Engineering, Peking University, Beijing, China

b

Department of Energy Resources Engineering, Stanford University, Stanford, California, USA

CR IP T

Abstract

In numerical modeling of subsurface flow and transport problems, formation properties may not be deterministically characterized, which leads to uncertainty in simulation results. In this study, we propose a sparse grid collocation method, which

AN US

adopts nested quadrature rules with delay and transformation to quantify the uncertainty of model solutions. We show that the nested Kronrod-Patterson-Hermite quadrature is more efficient than the unnested Gauss-Hermite quadrature. We compare the convergence rates of various quadrature rules including the domain

M

truncation and domain mapping approaches. To further improve accuracy and efficiency, we present a delayed process in selecting quadrature nodes and a

ED

transformed process for approximating unsmooth or discontinuous solutions. The proposed method is tested by an analytical function and in one-dimensional

PT

single-phase and two-phase flow problems with different spatial variances and correlation lengths. An additional example is given to demonstrate its applicability to

CE

three-dimensional black-oil models. It is found from these examples that the proposed method provides a promising approach for obtaining satisfactory estimation of the

AC

solution statistics and is much more efficient than the Monte-Carlo simulations.

1. Introduction

Uncertainty quantification of subsurface flow and transport models has become a subject of wide interest in the past few decades [13, 31, 44, 50, 61]. In general, this approach recognizes that the fluid flow properties such as hydraulic conductivity or permeability are unknown and may be regarded as Gaussian random fields with 2

ACCEPTED MANUSCRIPT

spatial variability. This randomness leads to uncertain model solutions such as hydraulic head and phase saturation. The uncertainty quantification that aims to estimate the statistical moments and probability density functions (PDFs) of the solution usually involves multidimensional integrals and interpolations. The commonly used Monte Carlo (MC) method is straightforward and robust, but could be computationally expensive [3]. As an efficient alternative, the collocation

CR IP T

methods have drawn a great deal of interests [34, 58]. Such approaches either construct Lagrange interpolating polynomials [2, 42], or use generalized polynomial chaos while determining the coefficients by spectral projection [34, 45, 57] or matrix inversion [36, 53]. In either approach, one needs to select a proper set of collocation

AN US

points, the choice of which is crucial for accuracy and efficiency. For univariate problems, a convenient way is to use appropriate quadrature nodes according to the distribution of the random parameter. For multivariate problems, one may use tensor products [23] or Smalyak sparse grids [57] based on univariate formulas.

M

In practice, one usually first evaluates the integral or constructs the interpolant with a low level quadrature formula, and then successively increases the level to refine the

ED

approximation. Ideally, the lower-level quadrature nodes could be re-used in the higher-level scheme, to minimize the integrand evaluations. This is called a nested

PT

quadrature rule. For uniform

distribution, popular

nested

quadrature rules

include the

CE

Clenshaw–Curtis rule [12] and Fejer’s rules [20], both of which are based on Chebyshev polynomials. Another choice is the Kronrod and Patterson extension of

AC

Gauss-Legendre rule [32, 48], which is referred to as Kronrod-Patterson-Legendre (KPL) rule in this study. To further enhance the efficiency of the KPL rule, Petras [49] presented a delayed process in which the number of nodes grows slowly in the sequence of formulas. It has been proved that the delayed process requires the asymptotically minimal number of nodes for a given degree of polynomial exactness. As for normal distribution, one may extend Kronrod and Patterson’s idea to the unnested

Gauss-Hermite

(GH)

rule

[17,

22],

which

is

referred

to

as

Kronrod-Patterson-Hermite (KPH) rule in this study. Although it has been applied to 3

ACCEPTED MANUSCRIPT

biology [27], finance [28, 29], electrical engineering [30] and solid mechanics [60], it has not been tested in subsurface flow and transport problems, and its performance has not been carefully investigated compared to other quadrature rules. In this work, we aim for the optimality of Smolyak sparse grids and propose a delayed KPH rule, which is more efficient than the standard KPH rule. It results from a repetition of quadrature formulas in the sequence, similar to the delayed process for

CR IP T

KPL in [49]. We investigate degree of exactness and stability factor for a series of GH-type unnested formulas and the KPH-type nested formulas. We also compare various quadrature rules including domain truncation and domain mapping approaches. To improve the performance of the collocation method when facing

AN US

discontinuous model responses, we introduce an additional transformed process [37], in which we approximate the arrival time that has higher regularity instead of the model response. We analyze the performance of these approaches in single-phase and multi-phase flow problems from one-dimensional to three-dimensional physical

M

spaces.

The paper is organized as follows. Section 2 introduces governing equations of

ED

single-phase and multi-phase flow problems. Section 3 describes the general framework of stochastic numerical methods. Section 4 presents the sparse grid

PT

collocation method using the KPH quadrature rule with delay and transformation. Section 5 shows numerical examples. Sections 6 and 7 provide discussions and

CE

conclusions.

AC

2. Governing equations

2.1. Single-phase flow

The steady-state, single-phase fluid flow satisfies the following equation [5]:

   K (x)h(x)  q(x)

(1)

where x is location, K (x) is the hydraulic conductivity, h(x) is the hydraulic head, 4

ACCEPTED MANUSCRIPT

and q(x) is the source/sink term.

2.2. Multiphase flow and transport model

The three-phase black-oil model is expressed by the following equations as [1]:

CR IP T

 k (x)kro     ( x) S o   (po  o gz )      qo  o Bo  t  Bo   k (x)krw     ( x) S w   (pw   w gz )      qw   w Bw  t  Bw 

 k (x)kro  k (x)krg    Rs (po  o gz )  (pg   g gz )  o Bo  g Bg    S S g     (x)  Rs o      Rs qo  qg  t  B B o g   

AN US



(2)

where o, w and g denotes oil, water and gas, respectively; i , pi , i , Si , Bi is the

M

viscosity, pressure, density, saturation and formation volume factor of phase i, respectively; k (x) is the absolute permeability; kri is the relative permeability of

ED

phase i, which is a function of Si;  (x) is porosity; Rs is the solution gas-oil ratio; qi

PT

is the source or sink term; g is the gravitational acceleration; and z is depth. Eq. (2) is

CE

coupled with:

So  S w  S g  1 pcow ( S w , S g )  po  pw

(3)

pcog ( S w , S g )  pg  po

AC

where pcow and pcog are the capillary pressures, which are functions of Sw and Sg. In this study, the hydraulic conductivity and permeability are treated as random

processes (space functions), thus the governing equations become stochastic differential equations, whose solutions are described by statistical moments and distributions.

3. Stochastic numerical methods 5

ACCEPTED MANUSCRIPT

3.1. Representation of input random field

To solve stochastic equations, one needs to represent the input random field properly. Let the input m be a random function of space coordinate x, time coordinate

CR IP T

t and random event θ, and denoted by m(x, t; )  m(x, t )  m(x, t; ) , where m(x, t ) is mean and m(x, t; ) is fluctuation. Assume the covariance is described by

Cm (x1 , x2 , t )  m(x1 , t; )m(x2 , t; ) , where



denotes expectation. It can be

decomposed as [24]: 

AN US

Cm (x1 , x 2 , t )   n f n (x1 , t ) f n (x 2 , t )

(4)

n 1

where n and f n (x, t ) are eigenvalues and eigenfunctions, respectively. Thus, the

M

random field m(x, t; ) can be expressed by the Karhunen-Loeve expansion (KLE): 

m(x, t; )  m(x, t )   n f n (x, t ) n ( )

(5)

ED

n 1

where  n are the independent random variables. We sort the eigenvalues n in the

N



n

 0.9 n , where N is the number of random variables retained [38]. n 1

CE

n 1



PT

non-increasing order and truncate the KLE according to a certain criterion, e.g.,

AC

3.2. Estimation of solution statistics

Consider a stochastic partial differential equation: s(x, t; ξ)  f (x, t; ξ)  0

where

(6)

is a differential operator, f is a known function, s is the solution,

ξ  (1 , 2 ,...,  N )T is the random parameter. For input ξ and output s(ξ ) , we study quadrature formulas: 6

ACCEPTED MANUSCRIPT

P

Q( s)   wk s(ξ k ), wk  , ξ k 

(7)

k 1

for N-dimensional integrals:

I ( s)   s(ξ)  (ξ)dξ

(8)



where ξ k are nodes/abscissas in sample space   1 

 N with intervals

CR IP T

 j , 1  j  N , wk are associated weights, P is the number of nodes, and

 (ξ)  1 (1 )    N ( N ) is the joint PDF. In this study, we focus on normal distribution, as the log-conductivity or log-permeability is usually assumed to be a

AN US

Gaussian random field.

3.2.1. Monte Carlo method

The Monte Carlo (MC) method is a simple ―brute-force‖ yet robust simulation technique. The typical procedure for the MC method is: first randomly generating P

M

realizations of random field, then solving the deterministic problems for each realization and finally estimating the solution statistics as (e.g., mean and variance): 1 P 1 P 2 2 s ( ξ ),    s(ξ k )  s    k s P k 1 P  1 k 1

ED s 

(9)

PT

In the traditional MC simulations, the mean converges asymptotically as P−0.5, which

CE

is independent of the number of random variables.

3.2.2. Collocation method

AC

The general procedure in the collocation method is similar to that of MC, but with

different nodes and weights: selecting P collocation points based on a certain quadrature rule, solving the deterministic problems on each collocation point and finally estimating the mean and variance as: P

P

k 1

k 1

s   s(ξ k )wk ,  s2   s 2 (ξ k ) wk  s 2

(10)

Skewness, kurtosis and higher moments can also be calculated. In this study, we 7

ACCEPTED MANUSCRIPT

compute the mean and variance to compare the convergence rates. To obtain the PDFs of the solution, we construct polynomial interpolations in the random space. This is usually done by a Lagrange interpolation with the collocation points [57]: P

sˆ(ξ)   s(ξ k ) Lk (ξ)

(11)

k 1

CR IP T

where L j (ξi )   ij ,1  i, j  Nc are Lagrange polynomials. Readers may refer to [2] for more aspects including convergence analysis and anisotropy approximation.

4. Sparse grid collocation with delay and transformation

AN US

4.1. Multivariate quadrature rules

The selection of collocation points is an important component of the collocation method, especially for multivariate problems. In this study, we focus on the problems

, define a sequence of univariate quadrature rules with level i:

ED

s : (, ) N 

M

where the input random fields follow normal distributions. Consider functions

ni

U i ( s)   s( X ki )  aki , i 

(12)

k 1

, X ni i } and weights

, ani i } . We denote mi as the degree of exactness, i.e., Eq. (12) is exact for all

CE

{a1i ,

PT

where ni is the number of nodes with nodal sets i  { X1i ,

polynomials of degree at most mi.

AC

A natural choice for an N-dimensional integral is the tensor product of univariate

integrals:



Q( s )  U 

which needs P  ni1

i1

U

iN

ni1

niN

 (s)    s( X k1 1

k N 1

i1 k1

,

, X kiNN )  aki11

akiNN

(13)

niN nodes and grows exponentially as N increases (i.e., the

―curse of dimensionality‖). To reduce the computational burden, the sparse grid strategy based on the Smolyak algorithm can be employed by a linear combination of 8

ACCEPTED MANUSCRIPT

product formulas as [23]: Q( s)  A(q, N ) 

  1

q i

q  N 1 i  q

 N  1 i   U 1  q  i 



 U iN



(14)

where the sparseness parameter q determines the order of the formula, and

i  i1  i2 

 iN . Set k  q  N as the ―level‖ of the Smolyak construction. To

(q, N ) 

q  N 1 i  q





 iN



AN US

4.2. Univariate quadrature rules

i1

CR IP T

compute A(q, N ) , one only needs to evaluate the function on the sparse grids: (15)

In this section, we discuss a variety of approaches for dealing with quadrature in the infinite domain (, ) as: 

I ( s)   s( )  ( )d 

(16)

M

These approaches include [10]: (1) domain truncation (approximating (, ) by

ED

[ L, L] with L  1); (2) domain mapping (changing coordinate from (, ) to

PT

[1, 1] ); (3) using basic functions intrinsic to an infinite interval (e.g., Hermite

function). Note that the domain truncation and domain mapping approaches are based

CE

on quadrature rules in a finite domain, we first introduce some well-known quadrature

AC

rules in the finite domain [1, 1] .

4.2.1. Quadrature rules in finite domain Gauss quadrature may be one of the most widely used quadrature rules. The idea is

to choose the weighting coefficients and the abscissas at which the function is to be evaluated in order to achieve a maximal degree of exactness. For example, the Gauss–Legendre (GL) rule yields a degree of exactness of mi  2ni  1 , for a quadrature in [1, 1] with uniform weight distribution [14]. If endpoints are used as 9

ACCEPTED MANUSCRIPT

abscissas, the Gauss-Lobatto-Legendre (GLL) rule can be considered [14]. However, Gauss quadrature suffers from the drawback of unnestedness, i.e., i  i1 . To address this issue in uniform distribution, Kronrod [32] added n  1 points to an n-point Gauss-Legendre quadrature formula such that the resulting rule has a degree of exactness of 3n  1 . Patterson [48] later extended Kronrod’s scheme, which is

CR IP T

referred to as Kronrod-Patterson-Legendre (KPL) in this work, through augmenting an n-point quadrature formula by p-points and reached an optimal degree of exactness of n  2 p  1 . Considering that the number of nodes in Patterson’s scheme grows too fast as the level increases, Petras [49] suggested a delayed sequence, in which the

AN US

quadrature formulas are repeated for some levels. He showed that the number of nodes for Smolyak algorithm is expected to be regularly asymptotically minimal, if the degree of exactness in a certain quadrature rule grows as mi  2i  1 . Here, regularly means the degree of exactness for the sequence of quadrature formula starts

M

from zero.

Another series of formulas are based on Chebyshev polynomials [54]. Fejer [20]

ED

suggested using the zeros of Chebyshev polynomial of first or second kind as interpolation points for quadrature rules. Fejer’s first rule (Fejer1) uses the nodes:

PT

xk  cos  k , k 

(2k  1) , k  1: n 2n

(17)

CE

We set the nodal sets as ni  3i 1 so that the sets are nested. Fejer’s second rule

AC

(Fejer2) uses the nodes:

xk  cos  k ,  k 

k , k  1: n  1 n

(18)

We set the nodal sets as ni  2i  1 . Both Fejer’s rules are open quadrature rules, i.e. they do not use the end points of the interval [1, 1] . Clenshaw-Curtis (CC) rule [12] is a closed version of Fejer’s second rule using the nodes:

xk  cos  k , k 

k , k  0:n n

(19)

We set the nodal sets as n1  1 and ni  2i 1  1, i  1 . The explicit formulas for the 10

ACCEPTED MANUSCRIPT

associated weights can be found, e.g., in [14].

4.2.2. Domain truncation Domain truncation is a simple strategy for solving problems in an infinite domain because it requires no modifications of strategies in a finite domain as: 

L



L

I (s)   s( )  ( )d   s( )  ( )d

CR IP T

(20)

Then we can apply a quadrature rule in a finite domain to approximate the integral in the infinite domain. If the function s( )  ( ) decays exponentially as    , the error in approximating the infinite domain (, ) by the finite domain [ L, L]

AN US

will decrease exponentially with the domain size L [10]. However, the error could be sensitive to the choice of the domain size, which will be discussed in Section 5.1.

4.2.3. Domain mapping

M

By using a transformation that maps an infinite interval into a finite domain, it is possible to calculate the integral as: 

1

I ( s)   s( )  ( )d   s( g (t ))  ( g (t )) g '(t )dt

ED



1

(21)

  g (t )   ,   , t  g ( )  [1,1] 1

PT

We remark that the accuracy depends on the regularity of s( g (t ))  ( g (t )) g '(t ) as a function of t . An infinite variety of maps is possible, but we can introduce some

CE

commonly used choices by considering four particular types:

AC

(1) Type 1 mapping, which is essentially the normal score transform [34, 42]:

 t 1      L   , t  2    1,  2  L 1

d L         dt 2   L  

1

(22)

where  is the standard normal probability distribution function,  is the standard normal cumulative distribution function. (2) Type 2 mapping, which is referred to as logarithmic mapping [10]: L    d   L arc tanh  t  , t  tanh   ,  dt 1  t 2 L

(23) 11

ACCEPTED MANUSCRIPT

(3) Type 3 mapping, which is referred to as algebraic mapping [10]:



Lt 1 t

2

, t

  L 2

2

,

d L  dt (1  t 2 )3 2

(24)

(4) Type 4 mapping, which is used in the Chebfun software [15]:

Lt 2  , t , 2 1 t L  L2  4 2

d L(1  t 2 )  dt (1  t 2 )2

(25)

CR IP T

In these types of mapping, L is a user-specified constant, whose optimal choice is problem-dependent. The above mapping methods can be combined with any quadrature rule in a finite domain. We will compare the performance of these

AN US

mapping methods in section 5.1.

4.2.4. Hermite-type quadrature rules

For normal distribution, the commonly used univariate quadrature rule is i Gauss-Hermite (GH) quadrature U GH with ni  i and mi  2i  1 [33]. Recently,

M

people considered some modifications on the standard GH rule by changing the quadrature sequences. Eldred [16] presented a linearly growing rule (GHlinear)

ED

i 2i 1 U GH  U GH with ni  {1,3,5,7,...} , so that the nodes are weakly nested (i.e., linear 2 i / 2 1

with

PT

i  U GH repeated at the origin). A similar slowly growing rule (GHslow) U GH slow

ni  {1,3,3,5,5,...} can also be considered, where  denotes the floor function. Lin

CE

ni i and Tartakovsky [38] used an exponentially growing rule (GHexp+1) U GH  U GH exp+1

AC

with n1  1 and ni  2i 1  1, i  1 . Eldred [16] presented another exponentially i 2 1  U GH growing rule (GHexp−1) U GH with ni  2i  1 . These rules can be referred exp1 i

to as the GH-type rules. A generalization of Patterson’s idea to Gauss-Hermite quadrature for normal distribution, which is called Kronrod-Patterson-Hermite (KPH) quadrature in this i work, are presented in [17] and [22]. They reported a sub-optimal choice U KPH

where 12

ACCEPTED MANUSCRIPT

ni  {1,3,9,19,35}, mi  {1,5,15, 29,51}, i  {1, 2,3, 4,5}

(26)

The nodes and weights can be found, e.g., at http://www.sparse-grids.de/. Recently, Bourquin [6] extended this sequence after an exhaustive search. Following the idea in Petras [49], we suggest a delayed Kronrod-Patterson-Hermite (KPHdelay) quadrature rule as: 9 15 4 16  26 5 U KPH  U KPH , U KPH  U KPH , delay delay

(27)

CR IP T

1 1 2 3 2 4 8 3 U KPH  U KPH , U KPH  U KPH , U KPH  U KPH , delay delay delay

with sequences ni  {1,3,3,9,9,9,9,9,19,19,...}, i  {1, 2,...} . For example, for level i  4 , the optimal degree of exactness is mi  2i  1  7 , and the smallest mi that

AN US

satisfy mi  7 is mi  15 , then the number of nodes is ni  9 .

Actually, Genz and Keister [22] presented a similar delayed sequence, where

ni  {1,3,3,7,9,9,9,9,17,19,19,...}, i  {1, 2,...} . For level i  4 , they used a number

M

of nodes ni  7 instead of ni  9 , based on the fact that ni  7 is enough for the degree of exactness mi  7 . Such an adjustment was also made for ni  17 when

ED

level i  9 . In some cases, its performance is close to that of KPHdelay. However, in other cases, it could be unstable and worse than the KPHdelay, as will be illustrated in

rules.

PT

Section 5.2. We refer to the KPH, KPHdelay and Genz-Keister rules as the KPH-type

CE

We point out that for GH-type rules, when the quadrature is computed twice using different number of points to check for accuracy, the model runs for the first

AC

computation are almost useless for the second computation. The reason is that for GH-type rules, the abscissas for a given number of points never equal the abscissas of a different number (except for the origin). This is unfortunate when the model run is computationally expensive. In contrast, for KPH-type rules, previous model runs can be reused when the level of sparse grids increases. Fig. 1 shows the degrees of exactness and corresponding nodes for different levels of univariate quadrature rules. The GH rule has a relatively slow (which is good) 13

ACCEPTED MANUSCRIPT

increasing speed in exactness, but it is unnested. The KPH rule is nested, but the exactness increases too fast. The KPHdelay rule combines the advantages of nestedness and slow increase in the exactness. Fig. 2 shows the degrees of exactness and corresponding nodes for multivariate (here is two) quadrature rules. To reach a total degree of exactness of 11, the GH, KPH and KPHdelay rules need 89, 65 and 45 nodes,

PT

ED

M

AN US

CR IP T

respectively.

CE

Fig. 1. Quadrature sequences (first row) and corresponding nodes (second row) for different

quadrature

rules:

(a)

and

(d)

Gauss-Hermite,

(b)

and

(e)

AC

Kronrod-Patterson-Hermite, and (c) and (f) delayed Kronrod-Patterson-Hermite.

14

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 2. Exactnesses (first row) and corresponding sparse grids (second row) for different

quadrature

rules:

(a)

and

(d)

Gauss-Hermite,

(b)

and

(e)

M

Kronrod-Patterson-Hermite, and (c) and (f) delayed Kronrod-Patterson-Hermite. To

ED

reach a total degree of exactness of 11, the KPHdelay requires the least number of

PT

nodes.

CE

4.3. Degree of exactness and stability factor

For different quadrature rules, we analyze the degree of exactness and the stability

AC

factor, whereas some other issues (e.g., error bound estimation and adaptive construction) require further investigation. Novak and Ritter [46] proved that if a univariate quadrature rule satisfies m1  1 ,

m2  3 , mi 1  mi  mi  mi 1 , defining m0  1 and assuming    q / N  and

  q mod N , then Smolyak algorithm A(q, N ) has degree of exactness: l (q, N )  (m 1  1)  N  (  1)  (m  1)(  1)  1

(28)

From Eq. (28), we get the following results for different choices of univariate 15

ACCEPTED MANUSCRIPT

quadrature rules:

l (q, N )GH  2(q  N )  1. l (q, N )GHslow  2(q  N )  1.

l (q, N )Genz-Keister  2(q  N )  1.

(29)

CR IP T

q  2N , 2(q  N )  1, l (q, N )GHlinear   4(q  N )  2 N  3, otherwise. 2(q  N )  1, q  2N ,  l (q, N )GHexp+1  4(q  N )  2 N  3, 2 N  q  3N , 2 1 ( N    1)  2 N  1, otherwise.  q  2N , 2(q  N )  1, l (q, N )GHexp1    2 ( N    1)  2 N  1, otherwise. l (q, N ) KPHdelay  2(q  N )  1.

AN US

We remark that in reality, l (q, N )GHslow , l (q, N )KPHdelay and l (q, N )Genz-Keister are expected to be better than 2(q  N )  1, since the degree of exactness mi  2i  1 for some level i. There is no explicit formula for l (q, N )KPH , since the number of nodes

ED

directly using Eq. (28).

M

increases irregularly in the sequence; nevertheless, we can calculate l (q, N )KPH

To compare the efficiency of different quadrature rules, we observe the number of

(q, N ) required for a given degree of exactness l (q, N ) with N = 2, as

PT

points

shown in Fig. 3a. We can see that when l (q, N ) is small, there is no significant

CE

difference between these rules. As l (q, N ) increases, GH needs the largest number

AC

of points; GHslow, GHlinear, GHexp+1 and GHexp−1 need moderate numbers of points; and KPH, KPHdelay and Genz-Keister need the least number of points. The result of GHlinear is essentially a smooth form of that of GHslow, which is not a surprise.

16

CR IP T

ACCEPTED MANUSCRIPT

Fig. 3. (a) Degree of exactness l (q, N ) vs. number of points

(q, N ) , and (b)

AN US

degree of exactness l (q, N ) vs. stability factor C (q, N ) , with N = 2.

Another issue that arises when using an integration rule is stability [22]. It is usually measured by the sum of the absolute values of weights normalized by the sum

M

of the weights, which is round-off error magnification factor in the worst case. For the quadrature rule in Eq. (10), we denote the stability factor as: P

ED

w k 1 P

k

w

PT

C

k 1

1

(30)

k

CE

The larger the stability factor, the more unstable the rule is expected to be. In Fig. 3b, we show the stability factors C (q, N ) for different rules to reach a certain degree of

AC

exactness l (q, N ) with N  2 . The KPH-type rules clearly outperform other rules with much smaller stability factor values, because the univariate quadrature nodes are nested and the multivariate weights are partially canceled out. We remark that the stability may also be affected by numerical error or discretization error, e.g., in some cases the stability issue of sparse grid method is more pronounced when the governing equations are solved with lower accuracy [35].

17

ACCEPTED MANUSCRIPT

4.4. Collocation method with transformation

The collocation method relies on the regularity of the model responses with regard to the parameters; thus, its accuracy deteriorates quickly as the regularity decreases, and its convergence is not guaranteed in systems with unsmooth dependence on input parameters. To address this issue, Liao and Zhang [37] proposed a time-based

CR IP T

transformed collocation method, by introducing a transformed process between the response and the arrival time. Specifically, for any fixed location x, instead of directly approximating the response s at time t, they first determine the arrival time when the response takes a possible value (forward transform from the response to the arrival

AN US

time): t (x, s; ξ k )  t s(x, t; ξ k )  s

(31)

then approximate the arrival time by the Lagrange interpolation polynomials: Nc

tˆ(x, s; ξ)   t (x, s; ξ k ) Lk (ξ)

(32)

M

k 1

and finally determine the response when the arrival time takes a possible value

ED

(backward transform from the arrival time to the response):





(33)

PT

s(x, t; ξ)  s tˆ(x, s; ξ)  t

The forward and backward transforms are essentially one-dimensional interpolation,

CE

which can be done by lookup tables. This approach significantly increases accuracy since the arrival times are much easier to approximate by polynomials than are the

AC

responses, especially for advection-dominated problems. We remark that the transformation is helpful not only for estimating saturation in multiphase flow problems as shown in this study, but also for estimating concentration [37], density [43] and other variables.

5. Numerical examples 5.1. Analytical function

18

ACCEPTED MANUSCRIPT

Since the log-conductivity or log-permeability is usually assumed to be a Gaussian random field, we consider an analytical function with normally distributed parameters: 



N

s(ξ)  exp   aii  , i ~ N (0,1)  i 1

(34)





s  







N



exp   aii   (1 )  i 1



 ( N )d1

CR IP T

and estimate its mean: d N ,  (i ) 

1 i2 /2 e 2

(35)

We set the coefficients to be N  3 , and ai  1, i  1,2,3 , and calculate Eq. (35) using Smolyak sparse grid algorithm with different univariate quadrature rules.

AN US

In the domain truncation approach, we select a proper domain size L by analyzing its effect on accuracy. Fig. 4a shows the estimation error by the domain truncation approach using the KPL rule for different values of L. When L  3 , the error ceases to decrease after about 1,000 points, indicating L is not large enough to cover

M

(, ) and the total error mainly comes from truncation. When L  5 , the error

ED

only starts to decrease after about 1,000 points, indicating L is too large so that the points are dispersed and the total error mainly comes from quadrature in [ L, L] . If

PT

we aim for an accurate result using O(103) points, L  4 is a better option than L  3 or L  5 . To find the optimal L for a fixed number of points, we plot error vs.

CE

L as shown in Fig. 4b. If we use 6,079 points (corresponding to the level 8 sparse grids using the KPL rule with three random variables), the optimal L is around 4.5.

AC

We can see that the optimal L is related to both the quadrature rule and the number of points. Furthermore, it is problem-dependent and usually not known beforehand. In this test, we choose a near-optimal L for each quadrature rule with O(103) points. The domain mapping approach with the user-defined constant L also has this issue, which is treated similarly.

19

CR IP T

ACCEPTED MANUSCRIPT

Fig. 4. (a) The convergence rate of the domain truncation approach is affected by the

AN US

domain size L; (b) for a fixed number of points, there exists an optimal L.

Fig. 5 shows the estimation errors using different approaches. Note that Fig. 5f uses a different scale in y-axis from those in Figs. 5a to 5e, we can see that the Hermite-type approach (Fig. 5f) converges faster than the domain truncation approach

M

(Fig. 5a) or the domain mapping approach (Figs. 5b to 5e).

ED

In the domain truncation and mapping approaches, the GL rule converges the slowest of all since it is unnested. The KPL, GLL, Fejer1 and Fejer2 rules have

PT

similar performance and are slightly better than the CC rule. These findings are consistent with [23] who reported that the KPL rule performs best in comparison to

CE

Gauss-Legendre, Clenshaw-Curtis or trapezoidal formulas. We can see that the type 2 mapping is more attractive than other types of mapping. Among all possible

AC

combinations, the type 2 mapping using the KPL rule and the Fejer1 rule seems to be the most efficient, thus we will perform more tests on these two rules in the following examples. In the Hermite-type approach, the KPH-type rules generally converge faster than the GH-type rules since the former ones are nested whereas the latter ones are not. In the KPH-type rules, the KPHdelay rule is the most efficient one followed by the Genz-Keister rule. Both are more efficient than the standard KPH rule owing to the delayed process. In the GH-type rules, the standard GH rule is the least efficient one 20

ACCEPTED MANUSCRIPT

as it is completely unnested, whereas other modified rules are more efficient thanks to weak nestedness (repeated at the origin). The GHlinear and GHslow rules are slightly more efficient than the GHexp+1 and GHexp-1 rules because the former ones are kinds of

AC

CE

PT

ED

M

AN US

CR IP T

delayed sequences (see Eq. (29)) as suggested by Petras [49].

21

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 5. Estimation error using different quadrature rules: (a) domain truncation; (b) type 1 domain mapping; (c) type 2 domain mapping; (d) type 3 domain mapping; (e) type 4 domain mapping; (f) Hermite-type approach. The domain truncation and domain mapping approaches are combined with various quadrature rules on a finite 22

ACCEPTED MANUSCRIPT

interval (e.g., GL, KPL, etc.).

5.2. One-dimensional single-phase flow example

In this section, we consider a single-phase flow example. We aim for the convergence rates of different quadrature rules, and thus choose a relatively small

CR IP T

number of random variables to make a high level of sparse grids (with a large number of model runs) affordable. The log-conductivity ln K (in m/day) is treated as a second-order stationary Gaussian random field with a mean of 0, and a covariance of Cln K   ln2 K exp( x1  x2

2

 2 ) , where  ln2 K  1 and the correlation length is

AN US

  0.5L , L  50 m is the domain length. The first two random variables in the KLE are retained, due to the fast decay of eigenvalues. We assume following boundary conditions:

dh  q0 , dx x 0

h xL  0

(36)

M

K

ED

where q0  1 m/day .

Fig. 6 shows the solution profiles of mean and variance of the hydraulic head. The

PT

results by the MC method with 100 and 1,000 model runs are obviously different from those by the MC method with 10,000 or more model runs. The collocation

CE

method using the KPHdelay rule with 45 model runs is as accurate as the MC method

AC

with over 10,000 model runs.

23

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 6. (a) Mean head and (b) head variance, with   0.5L and  ln2 K  1 .

To further compare the accuracy of different quadrature rules, we calculate the moment errors as:

1/2

(37)

M

 L  y ( x)  y ( x) 2 dx      0 L   0 dx  

ED

where y is the approximated value, and y is the reference value, which is obtained from an extremely high level of sparse grids (level 25 sparse grids using GHlinear with

PT

18,497 points in this test). Fig. 7a shows the error of the mean head. We illustrate ten sets of MC runs, the convergence rates of which are approximately P−0.5 as expected.

CE

We can see that KPH-type rules converge faster than the GH-type rules, and the KPHdelay rule converges fastest of all. These findings are consistent with Fig. 3. The

AC

mapped KPL rule is slightly better than the mapped Fejer1 rule, but both are less accurate than the Hermite-type rules, indicating the domain mapping approach is not an appropriate choice. This can be partially explained by that a function with high regularity on the infinite interval may turn into a function with low regularity on the finite interval after mapping. Fig. 7b reveals similar phenomena for the error of the head variance. The variance is essentially more difficult to estimate than is the mean, as the variance requires a polynomial approximation with a higher degree of exactness. We can see that in this case, the advantage of KPH-type rules over other quadrature 24

ACCEPTED MANUSCRIPT

AN US

CR IP T

rules is more obvious.

Fig. 7. Error of (a) mean head and (b) head variance, with   0.5L and  ln2 K  1 .

M

We next test a second case with a larger input variability of  ln2 K  4 . Fig. 8a shows the error of the mean head. The errors of the KPH-type rules quickly reach the

ED

machine epsilon using about 300 model runs, whereas the GH-type rules need over 1,000 model runs. Fig. 8b shows the error of the head variance. To reach the machine

PT

epsilon, the KPH-type rules use about 1,000 model runs, whereas GH-type rules need over 3000 model runs. It is noted that the Genz-Keister error may increase as the

AC

CE

number of nodes grows, which is clearly unfavorable.

25

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 8. Error of (a) mean head and (b) head variance, with   0.5L and  ln2 K  4 .

We then test a third case with a smaller correlation length of   0.2L and retain the first five random variables in the KLE. Fig. 9 shows that the KPHdelay error

M

decreases not only fastest of all but also in a smooth way. It is two to six orders of

ED

magnitude more accurate than the GH-type rules for the mean estimation, and is two to four orders of magnitude more accurate than the GH-type rules for the variance

PT

estimation. The Genz-Keister error also decreases fast but is not as smooth as the

AC

CE

KPHdelay error, especially for the variance.

26

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 9. Error of (a) mean head and (b) head variance, with   0.2L and  ln2 K  4 .

Ideally, a proper univariate quadrature rule for multivariate integration using the Smolyak sparse grid algorithm should better have four characteristics: (1) the

M

estimation error decreases fast as the number of nodes increases, i.e., fast convergence rate; (2) the nodes are nested; (3) the number of nodes increases relatively slowly as

ED

the level increases; and (4) the error decreases monotonically. The first two characteristics are directly related to accuracy and efficiency. The third one helps to

PT

increase the number of model runs gradually, so as not to waste too many function evaluations on the unnecessary nodes for a given accuracy threshold. The last

CE

characteristic makes it reasonable to terminate the calculation if the estimations from two successive levels of sparse grid construction are close enough. From the above

AC

tests, we see that the KPHdelay rule is superior to other quadrature rules.

5.3. One-dimensional multiphase flow example

In this section, we illustrate the performance of the transformed process through a two-phase flow example. We consider a one-dimensional domain (length L = 500 m) initially saturated with oil, which is displaced by water. There is a water injector at the 27

ACCEPTED MANUSCRIPT

left end and a producer at the right end, subject to the following boundary condition (BC) and initial condition (IC): BC: p x 0  p1 , S w

x 0

 1; p x  L  p2 , S w

xL

0

IC: p t 0  p2 , S w t 0  0

(38)

where p1  180 bar and p2  100 bar. We assume the porosity is 0.2 throughout the

CR IP T

reservoir and the log-permeability (in mD) is a second-order stationary Gaussian random field with mean equal to 3.0, and covariance as Cln k   ln2 k exp( x1  x2  ) , where  ln k  0.5 and   200 m. The relative permeabilities follow the Corey-type model in quadratic forms as krw  Sw2 , kro  (1  Sw )2 . The capillary pressure is

AN US

neglected.

In the MC method, to obtain relatively stable values of the sample moments, we carried out 100,000 realizations by direct model runs as a reference. In the collocation method, we retain the first six random variables in the KLE and use 73 collocation

M

points based on the KPHdelay quadrature rule.

Fig. 10 shows the mean and variance of pressure. The results from the MC method

ED

show that 100 or 1,000 realizations are not enough, whereas 10,000 realizations are visually sufficient. The collocation method with only 73 realizations produces more

PT

accurate results than the MC method with 1,000 realizations. Since saturation is less smoother than pressure, the former is more difficult to estimate than the latter by the

CE

collocation method. Fig. 11 shows the mean and variance of water saturation. Although we used the efficient KPHdelay quadrature rule, the results are still deviated

AC

from the reference. This deviation is successfully corrected after introducing the transformed process without additional collocation points. Fig. 12 shows the PDFs of pressure and saturation at location x  300 m. The reference is generated by 10,000 samples from forward simulations. In the collocation method, the PDF is estimated by 10,000 samples using the surrogate polynomials. The reference pressure PDF is unimodal and is approximated well by the collocation method using the KPHdelay quadrature rule (Fig. 12a). The reference saturation PDF is bimodal and is 28

ACCEPTED MANUSCRIPT

approximated well by the collocation method using the KPHdelay quadrature rule with transformation (Fig. 12b). There is a single peak point at Sw  0 since the saturation

AN US

CR IP T

values are zero in a large number of realizations.

AC

CE

PT

ED

M

Fig. 10. (a) Mean pressure and (b) pressure variance.

Fig. 11. (a) Mean saturation and (b) saturation variance.

29

CR IP T

ACCEPTED MANUSCRIPT

Fig. 12. PDFs of (a) pressure and (b) saturation, at location x = 300 m.

AN US

5.4. Three-dimensional multiphase flow example

We proceed to demonstrate the capability of the proposed method to handle three-dimensional three-phase flow problems. We consider the PUNQ-S3 test case,

M

which is a black-oil model set up by the production forecasting with uncertainty quantification (PUNQ) project and widely used for performing benchmarks [26]. The

ED

model contains 19×28×5 gridblocks, 1761 of which are active. The field is bounded to the east and south by a fault, and links to the north and west to a fairly strong aquifer.

PT

It contains 6 production wells, and no injection wells. A full description of this case can

be

found

on

the

webpage:

CE

http://www.ese.imperial.ac.uk/earth-science/research/research-groups/perm/standardmodels/. In this study, we set bottom-hole pressures as 150 bars for all wells and

AC

observe the production data up to 10 years. As the true permeability and porosity are given, we perturb the permeability to analyze how it affects the production data. We assume the permeability for each layer is the true permeability multiplied by independent coefficients i , i  1,...,5 , where  i follows log-normal distribution as

ln i ~ N (0,1) . In the MC method, we carried out 100,000 realizations as a reference. Fig. 13 shows the mean and variance of field oil production total (FOPT) and field gas 30

ACCEPTED MANUSCRIPT

production total (FGPT). The MC method with 10,000 realizations starts to produce results close to the reference. The collocation method uses only 51 realizations and matches the reference well, thanks to the efficient KPHdelay quadrature rule. Fig. 14 shows the results at time t  10 years as the number of realizations increases. We can

PT

ED

M

AN US

CR IP T

see that the collocation method converges much faster than the MC method.

Fig. 13. Statistical moments of production data compared between different methods:

AC

CE

(a) mean FOPT, (b) FOPT variance, (c) mean FGPT, and (d) FGPT variance.

31

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 14. Convergence check for (a) mean FOPT, (b) FOPT variance, (c) mean FGPT,

M

and (d) FGPT variance, at time t  10 years.

ED

We also analyze the results for individual wells. Fig. 15 shows the production data, including well oil production rate (WOPR), well gas production rate (WGPR) and

PT

well water cut (WWCT), for well PRO-12. We can see that both the MC method and the collocation method provide accurate estimations for the WOPR and WGPR. As

CE

for the WWCT, the MC method with over 1,000 realizations is able to provide visually accurate estimations, whereas the collocation method without transformation

AC

deviates from the reference due to the discontinuity caused by water breakthrough. This issue is successfully addressed after considering the transformed process. Fig. 16 illustrates the PDFs of the well production data at time t  10 years. We can see that the proposed collocation method using the KPHdelay quadrature nodes (with transformation for advection-dominated variables) captures different types of PDFs quite well.

32

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 15. Statistical moments of production data compared between different methods

CE

for well PRO-12: (a) mean WOPR, (b) WOPR variance, (c) mean WGPR, (d) WGPR

AC

variance, (e) mean WWCT, and (f) WWCT variance.

33

CR IP T

ACCEPTED MANUSCRIPT

Fig. 16. PDFs of (a) WOPR, (b) WGPR and (c) WWCT, at time t = 10 years.

AN US

6. Discussions 6.1. Finite and infinite intervals

A large amount of work has been devoted to collocation methods in finite intervals

M

with uniform distribution over the past decades [e.g., 23, 39, 40, 55-59]. However, in geostatistics, the parameters (e.g., log-conductivity) are usually considered as

ED

Gaussian random fields that are related to normal distribution in infinite intervals. Unfortunately, the unboundedness of interval results in some penalties. Boyd [10]

PT

showed that the rate of convergence for quadrature in an infinite interval is usually sub-geometric rather than geometric as in a finite interval. Gonnet [25] pointed out the

CE

error estimation is more difficult to provide for an infinite interval than a finite one. Most previous studies used GH-type rules for normal distribution [e.g., 2, 9, 36, 38],

AC

whereas some others used domain mapping approach combined with quadrature rules for uniform distribution [e.g., 34, 42]. In this study, we compare a variety of quadrature approaches and find that these approaches are not as efficient as the KPH-type rules. Here, we focus mainly on the computational aspects and less the theoretical derivation.

6.2. Adaptive quadrature 34

ACCEPTED MANUSCRIPT

Adaptive quadrature refers to the process of calculating the integral of a given function by adaptively subdividing the integration interval into smaller subintervals. This subdivision process is carried out recursively until the required accuracy is achieved. Many general purpose integration programs are adaptive, since such a strategy can be successful over a very wide range of integrands. For example, for univariate quadrature, the function quad, quadl and quadgk in MATLAB are based on

CR IP T

adaptive recursive Simpson rule, Gauss-Lobatto rule and Gauss-Kronrod rule, respectively. Gander [21] proposed adaptsim and adaptlob with improved efficiency and stability using different termination criteria and artificial limitation for recursion levels. Espelid [18] developed two doubly adaptive codes coteda, and coteglob, which

AN US

turn out to be superior to the previous ones. Espelid [19] extended the code coteglob and developed several new globally adaptive codes, in which the overall preferred code is da2glob. As for multivariate quadrature, Bungartz and Dirnstorfer [11] studied the potential of adaptive sparse grids for multivariate numerical quadrature in the

M

moderate or high dimensional case. Wan and Karniadakis [55, 56] introduced the multi-element generalized polynomial chaos method to address discontinuities in the

ED

random space. Ma and Zabaras [39, 40] proposed an adaptive hierarchical sparse grid collocation and utilized a high-dimensional model representation.

PT

However, the high accuracy from adaptivity is at the cost of more function evaluations near singularity or discontinuity. Thus, the computational burden could be

CE

increased, especially in high-dimensional problems. This issue is still an open question for adaptive collocation methods. Here, we suggest using transformation by

AC

approximating the arrival time instead of the model response. By doing so, the discontinuity or low-regularity issue can be substantially mitigated without extra nodes.

6.3. Integration and interpolation

For polynomial interpolation, there is a well developed theory that quantifies the convergence by the Lebesgue constant, which should be minimized for good point 35

ACCEPTED MANUSCRIPT

sets [54]. The Fekete and Leja points, which maximize the determinant of the Vandermonde matrix, are near-optimal points for polynomial interpolation that are defined for a compact set [7-9, 51]. Bos et al. [9] reported that tensor-product Gauss-Lobatto points are Fekete points for the cube. Sommariva et al. [52] concluded that the nontensorial cubature formulas at the Padua points concerning Lissajous curves turn out to be competitive with tensor-product Clenshaw-Curtis and

CR IP T

Gauss-Legendre formulas.

However, the Fekete points or Padua points are, at least now, only for a compact set, e.g., a cube [1, 1]N . The extension of these rules to infinite intervals, e.g.,

AN US

(, ) N , is a possible improvement in future research. We here limit our attention to find the best quadrature nodes for integration, and show that these nodes can also be used to construct surrogate polynomials to get PDFs. Further investigation is required for the interplay between interpolation and integration, which is beyond the scope of

M

this work. Interested readers are referred to [54, 59] for more details.

ED

7. Conclusions

PT

In this study, a new sparse grid collocation method has been developed for uncertainty quantification of flow and transport problems in heterogeneous porous

CE

media:

1. The proposed method is based on the Kronrod-Patterson-Hermite nested

AC

quadrature rule with a delayed process (KPHdelay), which has high degree of exactness and good stability. It is compared to the original Kronrod-Patterson-Hermite rule, a series of Gauss-Hermite-type unnested quadrature rules, domain truncation and domain mapping approaches combined with various quadrature rules for a finite interval

(including

Gauss-Legendre,

Gauss-Lobatto-Legendre,

Kronrod-Patterson-Legendre, Clenshaw-Curtis, and Fejer’s first and second rules), and Monte Carlo (MC) simulations. We find that the KPHdelay converges faster than other quadrature rules, and can be more efficient than the MC method by a few orders 36

ACCEPTED MANUSCRIPT

of magnitude. The Genz-Keister rule has a convergence rate similar to the KPHdelay rule, but it is not as stable as the latter. 2. For unsmooth model responses that have low regularity, such as saturation or water cut, a transformed process is introduced to improve accuracy. Essentially, the arrival time of the model response is approximated by the interpolating polynomials. Since the regularity of the arrival time is usually higher than that of the model

CR IP T

response, especially in advection-dominated problems, the arrival time is easier to approximate than the model response.

3. Combining delay and transformation, the proposed method provides satisfactory estimations for statistical moments and probability distributions in various numerical

three-dimensional

physical

AN US

tests including single-phase and multiphase flow problems in one-dimensional to spaces.

Meanwhile,

the

obtained

polynomial

approximation may be combined with reduced-order modeling (e.g., [47]) or serve as an accurate proxy for data assimilation and Bayesian inference (e.g., [4, 41]), which

M

will be considered in future work.

ED

Acknowledgements

This work is partially funded by the National Natural Science Foundation of China

PT

(Grant No. U1663208 and 51520105005), the National Key Special Science and Technology Program of China (Grant No. 2016ZX05037-003 and 2017ZX05009005),

CE

the Platform Construction Project for Research on Water and Ecology in the Ordos Plateau (Grant No. 201311076) and the China Postdoctoral Science Foundation

AC

(Grant No. 2016M590021). We also thank the Reservoir Simulation Research Consortium at Stanford University (SUPRI-B) for supporting this work.

37

ACCEPTED MANUSCRIPT

References

AC

CE

PT

ED

M

AN US

CR IP T

[1] Aziz K, Settari A. Petroleum reservoir simulation. Chapman & Hall; 1979. [2] Babuška I, Nobile F, Tempone R. A stochastic collocation method for elliptic partial differential equations with random input data. Siam Review 2010;52:317–355. [3] Ballio F, Guadagnini A. Convergence assessment of numerical Monte Carlo simulations in groundwater hydrology. Water Resources Research 2004;40. [4] Bazargan H, Christie M, Elsheikh AH, Ahmadi M. Surrogate accelerated sampling of reservoir models with complex structures using sparse polynomial chaos expansion. Advances in Water Resources 2015;86, Part B:385–99. [5] Bear J. Dynamics of fluids in porous media. Eisevier, New York, 764p 1972. [6] Bourquin R. Exhaustive search for higher-order Kronrod-Patterson Extensions. Research Report; 2015. [7] Bos L, De Marchi S, Sommariva A, Vianello M. Computing multivariate Fekete and Leja points by numerical linear algebra. SIAM Journal on Numerical Analysis 2010;48:1984–1999. [8] Bos L, Levenberg N. On the Approximate Calculation of Fekete Points: the Univariate Case. Electronic Transactions on Numerical Analysis 2008;30:377–397. [9] Bos L, Taylor MA, Wingate BA. Tensor product Gauss-Lobatto points are Fekete points for the cube. Mathematics of Computation 2001;70:1543–1547. [10] Boyd JP. Chebyshev and Fourier Spectral Methods. Dover Publications; 2001. [11] Bungartz HJ, Dirnstorfer S. Multivariate quadrature on adaptive sparse grids. Computing 2003;71:89–114. [12] Clenshaw CW, Curtis AR. A method for numerical integration on an automatic computer. Numerische Mathematik 1960;2:197–205. [13] Dagan G. Flow and transport in porous formations. Springer; 1989. [14] Dahlquist G, Björck Å. Numerical methods in scientific computing, volume I. SIAM, Philadelphia; 2008. [15] Driscoll TA, Hale N, Trefethen LN. Chebfun Guide. Pafnuty Publications, Oxford; 2014. [16] Eldred MS. Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design. AIAA Paper 2009;2274:37. [17] Elhay S, Kautsky J. Generalized Kronrod Patterson type imbedded quadratures. Applications of Mathematics 1992;37:81–103. [18] Espelid TO. Doubly adaptive quadrature routines based on newton–cotes rules. BIT Numerical Mathematics 2003;43:319–337. [19] Espelid TO. Algorithm 868: globally doubly adaptive quadrature—reliable matlab codes. ACM Transactions on Mathematical Software 2007;33:1–21. [20] Fejér L. Mechanische quadraturen mit positiven cotesschen zahlen. Mathematische Zeitschrift 1933;37:287–309. [21] Gander W, Gautschi W. Adaptive quadrature—revisited. BIT Numerical Mathematics 2000;40:84–101. 38

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[22] Genz A, Keister BD. Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight. Journal of Computational and Applied Mathematics 1996;71:299–309. [23] Gerstner T, Griebel M. Numerical integration using sparse grids. Numerical Algorithms 1998;18: 209–232. [24] Ghanem RG, Spanos PD. Stochastic finite elements: a spectral approach. Courier Corporation; 2003. [25] Gonnet P. A review of error estimation in adaptive quadrature. ACM Computing Surveys 2010;44:1173–1184. [26] Gu Y, Oliver DS. History matching of the PUNQ-S3 reservoir model using the ensemble Kalman filter. SPE Journal 2005;10:217–24. [27] Guedj J, Thiébaut R, Commenges D. Maximum likelihood estimation in dynamical models of HIV. Biometrics 2007;63:1198–206. [28] Heiss F, Winschel V. Likelihood approximation by numerical integration on sparse grids. Journal of Econometrics 2008;144:62–80. [29] Holtz M. Sparse grid quadrature in high dimensions with applications in finance and insurance. vol. 77. Springer Science & Business Media; 2010. [30] Horwood JT, Poore AB. Adaptive Gaussian sum filters for space surveillance. Automatic Control, IEEE Transactions on 2011;56:1777–90. [31] Ibrahima F, Meyer DW, Tchelepi HA. Distribution Functions of Saturation for Stochastic Nonlinear Two-Phase Flow. Transport in Porous Media 2015;109:81–107. [32] Kronrod AS. Nodes and weights of quadrature formulas: sixteen-place tables. Not Avail; 1965. [33] Kythe PK, Schäferkotter MR. Handbook of computational methods for integration. CRC Press; 2004. [34] Le Matre OP, Knio OM. Spectral methods for uncertainty quantification: with applications to computational fluid dynamics. 2010. [35] Lei H, Yang X, Zheng B, Lin G, Baker NA. Constructing surrogate models of complex systems with enhanced sparsity: Quantifying the influence of conformational uncertainty in biomolecular solvation. SIAM Journal on Multiscale Modeling and Simulation 2015;13:1327–1353. [36] Li H, Zhang D. Probabilistic collocation method for flow in porous media: Comparisons with other stochastic methods. Water Resources Research 2007;43. [37] Liao Q, Zhang D. Probabilistic collocation method for strongly nonlinear problems: 3. Transform by time. Water Resources Research 2016;52:2366–75. [38] Lin G, Tartakovsky AM. An efficient, high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media. Advances in Water Resources 2009;32:712–722. [39] Ma X, Zabaras N. An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations. Journal of Computational Physics 2009;228:3084–3113. [40] Ma X, Zabaras N. An adaptive high-dimensional stochastic model representation 39

ACCEPTED MANUSCRIPT

[46] [47]

[48] [49] [50]

[51]

CE

[52]

CR IP T

[45]

AN US

[44]

M

[43]

ED

[42]

PT

[41]

technique for the solution of stochastic partial differential equations. Journal of Computational Physics 2010;229:3884–3915. Man J, Li W, Zeng L, Wu L. Data assimilation for unsaturated flow models with restart adaptive probabilistic collocation based Kalman filter. Advances in Water Resources 2016;92:258–270. Mathelin L, Hussaini MY. A stochastic collocation algorithm for uncertainty analysis 2003. Mazzia A, Putti M. Three-dimensional mixed finite element-finite volume approach for the solution of density-dependent flow in porous media. Journal of Computational and Applied Mathematics 2006;185:347–359. Meyer DW, Tchelepi HA, Jenny P. A fast simulation method for uncertainty quantification of subsurface flow and transport. Water Resources Research 2013;49:2359–79. Müller F, Jenny P, Meyer DW. Probabilistic collocation and Lagrangian sampling for advective tracer transport in randomly heterogeneous porous media. Advances in Water Resources 2011;34:1527–1538. Novak E, Ritter K. Simple Cubature Formulas with High Polynomial Exactness. Constructive Approximation 1999;15:499–522. Pasetto D, Guadagnini A, Putti M. A reduced-order model for monte carlo simulations of stochastic groundwater flow. Computational Geosciences 2014;18:157–169. Patterson T cL. The optimum addition of points to quadrature formulae. Mathematics of Computation 1968;22:847–856. Petras K. Smolyak cubature of given polynomial degree with few nodes for increasing dimension. Numerische Mathematik 2003;93:729–753. Riva M, Guadagnini A, Dell’Oca A. Probabilistic assessment of seawater intrusion under multiple sources of uncertainty. Advances in Water Resources 2015;75:93–104. Sloan IH, Womersley RS. Extremal systems of points and numerical integration on the sphere. Advances in Computational Mathematics 2004;21:107125. Sommariva A, Vianello M, Zanovello R. Nontensorial clenshaw–curtis cubature. Numerical Algorithms 2008;49:409–427. Sun AY, Zeidouni M, Nicot J-P, Lu Z, Zhang D. Assessing leakage detectability at geologic CO2 sequestration sites using the probabilistic collocation method. Advances in Water Resources 2013;56:49–60. doi:10.1016/j.advwatres.2012.11.017. Trefethen LN. Approximation Theory and Approximation Practice. Society for Industrial and Applied Mathematics; 2013. Wan X, Karniadakis GE, An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. Journal of Computational Physics 2005;209:617–642. Wan X, Karniadakis GE, Multi-element generalized polynomial chaos for arbitrary probability measures. SIAM Journal on Scientific Computing 2006;28:901–928.

AC

[53]

[54] [55]

[56]

40

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

[57] Xiu D, Hesthaven JS. High-order collocation methods for differential equations with random inputs. SIAM Journal on Scientific Computing 2005;27:1118–1139. [58] Xiu D. Efficient collocational approach for parametric uncertainty analysis. Communications in Computational Physics 2007;2:293–309. [59] Xiu D. Numerical methods for stochastic computations: a spectral method approach. Princeton University Press; 2010. [60] Xu H, Rahman S. A generalized dimension-reduction method for multidimensional integration in stochastic mechanics. International Journal for Numerical Methods in Engineering 2004;61:1992–2019. [61] Zhang D. Stochastic methods for flow in porous media: coping with uncertainties. Academic press; 2002.

41