A stochastic dimension reduction multiscale finite element method for groundwater flow problems in heterogeneous random porous media

Journal of Hydrology 478 (2013) 77–88

Contents lists available at SciVerse ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

A stochastic dimension reduction multiscale finite element method for groundwater flow problems in heterogeneous random porous media Xinguang He a, Lijian Jiang b,⇑, J. David Moulton b a b

College of Resource and Environment Science, Hunan Normal University, Changsha 410081, China Applied Mathematics and Plasma Physics, Los Alamos National Laboratory, NM 87545, USA

a r t i c l e

i n f o

Article history: Received 13 July 2012 Received in revised form 13 November 2012 Accepted 21 November 2012 Available online 29 November 2012 This manuscript was handled by Peter K. Kitanidis, Editor-in-Chief, with the assistance of Jian Luo, Associate Editor Keywords: High-dimensional model representation Sparse grid stochastic collocation method Multiscale finite element method Groundwater flow Heterogeneous random porous media

s u m m a r y In this paper we present a stochastic dimension reduction multiscale finite element method for solving groundwater flow problems in heterogeneous random porous media. The stochastic conductivity is parameterized in a stochastic space using a truncated Karhunen–Loève expansion with finite random variables. This often results in high stochastic dimensionality due to small correlation length and large variance of the inherent covariance function. In order to treat the high-dimensional stochastic problems efficiently, a truncated high-dimensional model representation (HDMR) technique is proposed that decomposes the high-dimensional stochastic problem into a moderate-dimensional stochastic problem and a few one-dimensional stochastic problems. The derived low-dimensional stochastic problems are solved by incorporating a sparse grid stochastic collocation method (SGSCM) with the proposed HDMR. This leads to a set of uncoupled deterministic problems at the sparse grid collocation points. To capture the small-scale heterogeneities in physical space, a multiscale finite element method (MsFEM) is utilized to solve these deterministic problems. Since the SGSCM and MsFEM are implemented in the different spaces, the integrated approach inherits the advantages from both of the SGSCM and the MsFEM, in which each low-dimensional stochastic problem is decoupled each other at small number of stochastic collocation points and is solved on a coarse spatial mesh. To demonstrate the efficiency and accuracy of the proposed stochastic MsFEM with the HDMR technique, a few numerical experiments are carried out for the transient water flow problems in heterogeneous random media with different correlation lengths and various degrees of spatial variabilities. Published by Elsevier B.V.

1. Introduction Medium properties for groundwater models, such as hydraulic conductivity, are often highly heterogeneous and may vary at all length scales, from pore scale to field scale. Computer simulations of such models that use a very fine grid to capture the heterogeneity are computationally very expensive and possibly infeasible. However, disregarding the heterogeneity can lead to large errors. In order to achieve reasonable computational performance, the simulation of groundwater flows is often performed on a coarse grid using multiscale methods. Because of measurement errors and the incomplete knowledge of medium properties, modeling of groundwater flow often involves large uncertainties in the model’s inputs (e.g., hydraulic conductivity). The heterogeneities and uncertainties can be integrated in a stochastic multiscale model. The model’s output can be predicted accurately by simulating the ⇑ Corresponding author. Address: Applied Mathematics and Plasma Physics Theoretical Division, Los Alamos National Laboratory Los Alamos, NM 87545, USA. Tel.: +1 505 665 4497; fax: +1 505 665 5757. E-mail address: [email protected] (L. Jiang). 0022-1694/$ - see front matter Published by Elsevier B.V. http://dx.doi.org/10.1016/j.jhydrol.2012.11.052

stochastic multiscale model. The existence of heterogeneity at multiple scales combined with its uncertainty brings significant challenges to the development of efficient algorithms for the simulation of the dynamic processes in random heterogeneous porous media. Therefore, the interest in developing stochastic multiscale methods for stochastic multiscale models has steadily grown in recent years (Dostert et al., 2008; Ganapathysubramanian and Zabaras, 2009; Ganis et al., 2011; Nguyen, 2008; Shi et al., 2010). One of the challenges in the simulation of transport through heterogeneous random porous media is the multiscale nature of the media properties. The direct full-scale spatial and temporal resolution simulation of such problems is difficult and impractical due to the cost required for resolving the smallest scales. Multiscale methods therefore have been used, including the multiscale finite element method (MsFEM) (Hou and Wu, 1997), variational multiscale method (VMM) (Hughes et al., 1998) and heterogeneous multiscale method (HMM) (Engquist, 2003), to upscale complex fine-scale equations onto coarse-scale equations that capture the fine-scale affect on coarse-scale solutions. Among these multiscale methods, MsFEM has become a class of important simulation methods for multiscale systems (Efendiev and Hou, 2009). The

78

X. He et al. / Journal of Hydrology 478 (2013) 77–88

main idea of MsFEM is to incorporate the small-scale information into the finite element basis functions and capture their effect on the large scale via finite element computations. In many cases, the multiscale basis functions can be pre-computed and used repeatedly in subsequent computations with different source terms and boundary conditions. This leads to a large computational saving in upscaling saturated flows where the pressure head needs to be solved many times dynamically. To efficiently treat spatial heterogeneities, we use the MsFEM for spatial discretization in this paper. In addition to multiscale features, the uncertainties associated with porous media also pose significant challenges for simulations. One way to deal with this difficulty is to view the medium property variations as a random field that satisfies certain statistical correlations. This naturally results in describing the flow problem using stochastic partial differential equations (SPDEs). Many numerical methods have been developed for solving SPDEs (e.g., Babus˘ka et al., 2007; Matthies and Keese, 2005; Xiu and Hesthaven, 2005; Zhang et al., 2000). Among the existing stochastic numerical methods, the stochastic collocation method has been paid much attention in the research community because it shares the fast convergence property of the stochastic finite element method while having the decoupled nature of Monte Carlo sampling (Ma and Zabaras, 2011). The stochastic collocation method represents the solution of stochastic PDEs as a polynomial approximation in stochastic space. The interpolation operator is constructed via independent functions, which call for deterministic problem solvers at different collocation points which are selected based on special rules. The different choice of collocation points leads to the different collocation methods including the full-tensor product method (Babus˘ka et al., 2007) and Smolyak sparse grid method (Ganapathysubramanian and Zabaras, 2007; Nobile et al., 2008). Sparse grid stochastic collocation is known to have the same asymptotic accuracy as full-tensor product collocation but uses significantly fewer collocation points. However, the number of collocation points required in the sparse grid method increases drastically for high-dimensional stochastic problems. Thus, the method is still limited to the problems with a moderate number of random variables. However, the uncertainties of porous media are usually parameterized by a large number of random variables in the transport through the random media due to the small correlation length and large variability of the covariance structure. Consequently, the model’s input is defined in a high-dimensional random parameter space. Sampling in high-dimensional random space is very difficult. If the sampling is conducted in full random space, then the number of samples increases exponentially with respect to the dimension of the random space. This is the notorious curse of dimensionality, and is critical challenge for stochastic approximation in the high-dimensional stochastic space. Dimension reduction techniques such as high-dimensional model representation (HDMR) (Rabitz et al., 1999), proper orthogonal decomposition and principal component analysis, can partially alleviate this problem. Among these methods, HDMR is one of the most efficient stochastic dimension reduction techniques and has drawn the attention of researchers in recent years (Li et al., 2001; Ma and Zabaras, 2010, 2011; Nguyen, 2008; Yang et al., 2011). The HDMR approach was originally developed in the application of chemical models by Rabitz et al. (1999) and has become a general set of quantitative assessment and analysis tools for capturing the high-dimensional relationship between model inputs and model outputs. HDMR has been used for improving the efficiency of deducing high-dimensional input–output system behavior. Instead of dealing with the full high-dimensional random space, the HDMR technique can decompose a high-dimensional model into a set of low-dimensional models, each of which requires less computa-

tional effort. The curse of dimensionality can be considerably suppressed using HDMR. However, there exist some drawbacks in traditional HDMR (Li et al., 2001; Rabitz et al., 1999) and adaptive HDMR (Ma and Zabaras, 2010; Yang et al., 2011). For example, a large number of low-dimensional models still need to be computed, and this is a disadvantage when combining HDMR with MsFEM because a large number of multiscale basis functions need to be computed. To overcome the above mentioned drawbacks, we propose a new truncated HDMR for stochastic dimension reduction in this paper. Then, we incorporate the proposed HDMR with a sparse grid stochastic collocation method (SGSCM) and a MsFEM to develop a stochastic dimension reduction multiscale finite element method. We investigate the applications of the proposed method to the groundwater flow in heterogeneous random porous media. The goal of the proposed HDMR is to alleviate the curse of dimensionality in high-dimensional stochastic spaces whereas the MsFEM is utilized to coarsen the flow equations in the spatial domain. The proposed HDMR decomposes a high-dimensional model into a moderate-dimensional model and a few one-dimensional models. Then each low-dimensional random model derived by the proposed HDMR is solved utilizing SGSCM and MsFEM. Hence the high-dimensional stochastic and fine-scale spatial model is significantly reduced and the model’s output can be efficiently predicted. The mean and variance of the output can be computed efficiently. We note that the proposed HDMR provides a good trade-off between computational complexity and dimension reduction error. We will discuss the performance of the stochastic dimension reduction MsFEM for saturated transient flow models with different correlation lengths and various degrees of spatial variabilities. Extensive numerical studies are made to confirm the accuracy, efficiency and computation complexity for the presented method. The rest of the paper is organized as follows. In Section 2, we briefly introduce the groundwater flow problem. Section 3 is devoted to presenting stochastic variational formulation and a sparse grid collocation method. The proposed truncated HDMR is introduced in Section 4. In Section 5, the multiscale finite element formulation is presented for the flow equations at the resulting collocation points. In Section 6, a few numerical examples are presented to demonstrate the performance of the stochastic dimension reduction MsFEM. Some conclusions and closing remarks are made finally. 2. Flow problem Let D denote a two-dimensional bounded physical domain, and ðX; F ; PÞ denote a complete probability space with x 2 X denoting a sampling point, where F is the complete r-algebra over the subsets of X and P is the probability measure. The saturated flow in heterogeneous random porous media with a spatially varying conductivity can be described by the following stochastic parabolic partial differential equation:

@wðx; t; xÞ  r  ½Kðx; xÞrwðx; t; xÞ @t ¼ Rðx; tÞ in D  ð0; T  X

Ss ðxÞ

ð1Þ

subject to the following initial and boundary conditions:

wðx; t; xÞjt¼0 ¼ w0 ðxÞ x 2 D;

ð2Þ

wðx; t; xÞ ¼ wC ðx; tÞ x 2 CD ;

ð3Þ

 Kðx; xÞrwðx; t; xÞ  n ¼ qðx; tÞ x 2 CN ;

ð4Þ

where x 2 D is the spatial coordinate, K is the hydraulic conductivity, which is often highly variable in the physical space D, Ss is the specific storage coefficient, w is the hydraulic head, R is the

79

X. He et al. / Journal of Hydrology 478 (2013) 77–88

source/sink term, T is the study time, D  ð0; T  X is the usual Cartesian product, w0 is the initial head in the domain D, wC is the fixed head on the Dirichlet boundary CD , q is the flux across the Neumann boundary CN, n is the unit outer normal vector of T CN, and CD CN ¼ £. In the paper, we consider Y(x, x) = lnK(x, x) and is a Gaussian field with second-order moment. Without loss of generality, we may assume that Ss is a deterministic function in the domain D. We note that in any case, Ss varies less than lnK (Dagan, 1982). 3. Stochastic variational formulation and sparse grid collocation method In this section, we present a stochastic variational formulation for flow problem (1) and describe the sparse grid stochastic collocation method (SGSCM) for solving the stochastic flow problems in random heterogeneous porous media by following the recent works developed by Babus˘ka et al. (2007), Nobile and Tempone (2009) and Nobile et al. (2008). We first introduce the Karhunen– Loève expansion (KLE) of random fields to parameterize stochastic conductivity fields. 3.1. Karhunen–Loève expansion The Karhunen–Loève expansion of the random field Y(x, x) is based on the spectral expansion of its covariance function cov[Y](x1, x2), where x1 and x2 denote two different spatial points. The covariance function is real, symmetric, and positive-definite. All its eigenfunctions are mutually orthogonal and form a complete set spanning the function space to which Y(x, x) belongs (Loève, 1977). Then the random field Y(x, x) can be approximately represented by the following truncated KLE:

Yðx; xÞ  E½Yðx; xÞ þ

N pffiffiffiffi X ki bi ðxÞhi ðxÞ;

ð5Þ

i¼1

where fhi ðxÞgNi¼1 are uncorrelated random variables. If Y(x, x) is a Gaussian field, then the random variables in (5) are standard identically independent Gaussian random variables. Also, fki gNi¼1 and fbi ðxÞgNi¼1 are the eigenvalues and eigenfunctions of the correlation function, respectively. They are the solutions of the following integral equation:

Z

cov ½Yðx1 ; x2 Þbðx2 Þdx2 ¼ kbðx1 Þ:

ð6Þ

Nobile and Tempone, 2009): find stochastic function wðx; t; xÞ : D  ð0; T  X ! R and almost everywhere (a.e.) in (0, T] such that

d dt ¼

Z Z

Z

PðdxÞðSs w; v Þ þ

X

PðdxÞðK rw; rv Þ þ

Z

X

PðdxÞðR; v Þ;

X

PðdxÞ

Z

qv ds

CN

X

8v 2 H10 ðDÞ  L2P ðXÞ;

ð7Þ

where the notation (u, v) stands for the L2 inner product of functions R u, v over the domain D, that is ðu; v Þ ¼ D uv dx. According to Eq. (5), K(x, x) can be represented by KN(x, H(x)) = KN(x, h1(x), . . . , hN(x)). We denote with ci  hi(X), the image of hi ; !N ¼ Ni¼1 ci , and assume that the random variables H have a joint probability density function q : !N ! Rþ , with q 2 L1(!N). Then the solution w(x, t, x) of the stochastic problem (7) can be approximated by w(x, t, H) = w(x, t, h1(x), . . . , hN(x)). Thus, our goal is to approximate the function w = w(x, t, H), where x 2 D and H 2 !N. Accordingly, the finite dimensional version of the formulation (7) reads as following: find wðx; t; HÞ : D  ð0; T !N ! R and a.e. in (0, T] such that

d dt

Z !N

¼

ðSs w; v ÞqdH þ

Z

!N

ðR; v ÞqdH;

Z !N

ðK N rw; rv ÞqdH þ

Z !N

Z

qqv dsdH

CN

8v 2 H10 ðDÞ  L2q ð!N Þ:

ð8Þ

The stochastic problem (7) now becomes a deterministic problem with an N-dimensional parameter. Then, we can utilize any stochastic discretization method in the random space and the resulting equations become a set of deterministic equations. We use the notation w(H):¼w(x, t, H) when we emphasize the dependence of the parameter H. We use a similar notation for the random input coefficient KN. Then it can be shown that the problem (7) is equivalent to (cf., Nobile et al., 2008)

d ðSs wðHÞ; /Þ þ ðK N ðHÞrwðHÞ; r/Þ þ dt

8/ 2 H10 ðDÞ;

for a:s: H 2 !N :

Z

q/ds ¼ ðR; /Þ;

CN

ð9Þ

We can use the stochastic collocation method (Babus˘ka et al., 2007) to discretize the stochastic space !N. The basic idea is to approximate the multi-dimensional stochastic space !N using interpolating Np functions associated to a set of collocation points fHk gk¼1 2 !N , and use a finite element approximation for the spatial domain D.

D

The finite number N of terms retained depends on the spectrum enP ergy measured by Ni¼1 ki =ðr2 jDjÞ, where r2 is the variance of Y(x, x) and jDj is the area of the spatial domain D. Chang and Zhang (2009) provided some guidance on the choice of N under various correlation lengths. From Eq. (5), we represent Y(x, H(x))  YN(x, h1 (x), . . . , hN(x)) in the paper, where H = (h1, . . . , hN). As a result, K(x, x) is approximated by the finite dimensional parameter representation KN(x, H(x)) = exp (YN(x, H(x))). 3.2. Stochastic variational formulation In this subsection, we briefly present the stochastic variational formulation of the stochastic flow problem (1) using the finite dimensional representation of random input coefficient. For presentation, we introduce some basic notations. Denote the Hilbert R space L2P ðXÞ with the inner product ðw; v ÞX ¼ X PðdxÞwðxÞv ðxÞ, where, as before, P is the probability measure. Denote by H1 ðDÞ the usual Sobolev space. Let the tensor product space H10 ðDÞ  L2P ðXÞ ¼ fwjw 2 H1 ðDÞ  L2P ðXÞ; w ¼ 0 on CD g: Then, multiplying the flow Eq. (1) by appropriate test functions and using integration by parts leads to the following weak formulation (cf.,

3.3. Sparse grid stochastic collocation In this subsection, we describe sparse grid stochastic collocation techniques for Eq. (8). We seek an approximate solution to Eq. (8) in a suitable finite dimensional subspace associated to a set of collocation points {Hk}. We solve Eq. (9) at each collocation point Hk, and build the approximate solution w(x, t, H) of Eq. (8) by interpolating the collocated solutions. There are several options to generate the sets of collocation points {Hk}. In general, we use the roots of an orthogonal polynomial to find the collocation points. Regardless of the choice of collocation points, the corresponding interpolation can be constructed by using either full tensor product polynomials or sparse polynomials. The Smolyak sparse grid collocation method is known to have the same asymptotic accuracy as full-tensor product collocation method, while requiring many fewer collocation points as the parameter dimension increases. Below, we briefly introduce Smolyak sparse grid collocation method. Let Ii be a one-dimensional Lagrange interpolation operator, Ii wðhÞ ¼ Pmi ki ¼1 wðhi;ki Þ  li;ki ðhÞ, where li;ki ðhÞ is the interpolation polynomial of degree (mi  1) and fhi;ki g is a set of mi points for the onedimensional variable h. Let I Ni be an N-dimensional interpolation

80

X. He et al. / Journal of Hydrology 478 (2013) 77–88

operator. We denote by i the multi-index i ¼ ði1 ; i2 ; . . . ; iN Þ 2 NNþ . Then, for each w(H) 2 C(!N), the full tensor product interpolant of w can be written as

I Ni wðHÞ ¼ ðIi1      IiN ÞwðHÞ ¼

mi 1 X



k1 ¼1

mi N X

wðhi1 ;k1 ; . . . ; hiN ;kN Þ  ðli1 ;k1      liN ;kN Þ:

wðHÞ ¼ w0 þ

ð1Þqjij

qNþ16jij6q



N1 q  jij



 ðIi1      IiN Þ;

ð11Þ

where jij = i1 +    + iN for multi-index i 2 NNþ . Equivalently, formula (11) can be written as

X

Aðq; NÞðwÞ ¼

ð1Þqjij

qNþ16jij6q



N1 q  jij

 X  wðHk Þ  li;k ;

ð12Þ

k

where

Hk ¼ ðhi1 ;k1 ; . . . ; hiN ;kN Þ 2 !N ; k ¼ ðk1 ; k2 ; . . . ; kN Þ 2 NNþ ; kj ðj ¼ 1; . . . ; NÞ denotes the location of the given support node Hk in the jth dimension, and li;k ¼ li1 ;k1      liN ;kN denotes the N-dimensional Lagrange polynomial. From the above sparse grid interpolation formulas, we can see that to compute the approximate solution w ¼ Aðq; NÞðwÞ of the problem (8), we only need to know function values on the sparse grid

Hðq; NÞ ¼

[

ð#i1      #iN Þ;

N L X X wi ðhi Þ þ i¼1

Clearly, the above product requires function values j¼1 mij wðHk Þ ¼ wðhi1 ;k1 ; . . . ; hiN ;kN Þ sampled on random parameter space. Because each of the function values w(Hk) is obtained by solving Eq. (9), full-tensor product collocation is prohibitively expensive. Smolyak sparse grid collocation method can be used to alleviate this difficulty (see, Barthelmann et al., 2000; Dostert et al., 2008). The Smolyak interpolant Aðq; NÞ (q > N) is simply a linear combinations of tensor product interpolations with the following property: only products with a relatively small number of nodes are used, and the linear combination is chosen in such a way that an interpolation property for N = 1 is preserved for N > 1 (cf., Smolyak, 1963; Barthelmann et al., 2000). The Smolyak algorithm constructs the sparse interpolant Aðq; NÞ by using the formulas (10) as the building blocks. Aðq; NÞ is given by (Barthelmann et al., 2000)

X

Let w(H) = w(h1, h2, . . . , hN) be a real-valued function representing the map between the input variables h1, h2, . . ., hN and the output. For L N, the traditional truncated HDMR of w(H) is represented in the form

ð10Þ

kN ¼1

QN

Aðq; NÞ ¼

4.1. Traditional truncated HDMR technique

ð13Þ

qNþ16jij6q

where #i ¼ fhi;1 ; . . . ; hi;mi g denotes the set of points used in 1D interpolation operators Ii. According to suggestion by Barthelmann et al. (2000) and Nobile et al. (2008), in this work we use Clenshaw-Curtis abscissas, which are the extrema of Chebyshev polynomials, to conr struct the Smolyak interpolant. This leads to nðr þ N; NÞ  2r!  N r collocation points used by AðN þ r; NÞ (see, Barthelmann et al., 2000), where r usually refers to the interpolation level and defines the ‘‘size’’ of the sparse grid. We note that AðN þ r; NÞ is exact for all polynomials of degree r using the Smolyak interpolant. 4. High dimensional model representation technique If the computation of the stochastic problem (8) is directly conducted in full stochastic space !N based on the Smolyak formulas (12), then the number of collocation points for the stochastic simulation increases drastically with respect to the dimension of the stochastic space. This poses great difficulties for the stochastic simulation in a high-dimensional stochastic space. To enhance the efficiency of the simulation, we develop a new truncated HDMR for stochastic dimension reduction.

X

wi1 im ðhi1 ; . . . ; him Þ:

ð14Þ

m¼216i1 <
Here w0 is the zeroth-order component denoting the mean effect of w(H). The first-order component function wi(hi) represents the individual contribution of the input hi, and the second-order component function wi1 i2 ðhi1 ; hi2 Þ represents the cooperative effects of hi1 and hi2 and so on. It is argued in the work by Rabitz et al. (1999) that the low-order component functions (e.g., up to the third-order) usually have the main contribution on the output in the high-dimensional models for most realistic physical systems. HDMR is roughly classified into two categories: ANOVA-HDMR and Cut-HDMR (Li et al., 2001; Ma and Zabaras, 2010). Since the computation of the components in Cut-HDMR is cheap, we adopt Cut-HDMR in the paper. For Cut-HDMR, we define a cut-point H ¼ ð h1 ; . . . ;  hn Þ. In the paper, H is the mean vector of the random input H. Then the component functions of Cut-HDMR are explicitly given as follows (cf., Ma and Zabaras, 2010):

b iÞ  w ; wi ðhi Þ ¼ wð H 0 b i;j Þ  w  w  w ; . . . ; wij ðhi ; hj Þ ¼ wð H i j 0

w0 ¼ wðHÞ;

ð15Þ

b i ¼ ð where H h1 ; . . .  hi1 ; hi ;  hiþ1 ; . . . ;  hn Þ, and b      H i;j ¼ ðh1 ; . . . hi1 ; hi ; hiþ1 ; . . . ; hj1 ; hj ; hjþ1 ; . . . ; hn Þ: Although the traditional truncated HDMR in Eq. (14) can reduce dimensionality of stochastic space, they have at least two drawbacks: (a) There are many terms in the truncated HDMR expansion   P 50 (e.g., 1 þ 2i¼1 ¼ 1276 terms if N = 50 and truncating up to i 2nd order), the total computation for all the components is very expensive. (b) For the truncated Cut-HDMR, the total variance of w(H) may be not equal to the summation of variances of all terms in right-hand side of (14). This leads to difficulty in computing variance of w(H). To overcome these drawbacks, we propose a new truncated HDMR technique and combine it with MsFEM for model reduction. 4.2. A new truncated HDMR formulation In this paper, we focus on the case of random conductivity field K(x, H) = exp (Y(x, H)) and Y(x, H) is generated by truncated KLE (cf., Eq. (5)). From Eq. (5), we can know that the eigenvalue kj is nonincreasing as j increases. Moreover, kj characterizes the energy of input information hj in the jth term of the KLE. We know that the larger the magnitude of kj is, the stronger the impact of the corresponding hj on the output is. Therefore, we can set a threshold constant with 0 < f < 1, and find the smallest M such that M N X X kj = kj P f: j¼1

ð16Þ

j¼1

M Then we define fhj gM j¼1 corresponding to the eigenvalue set fkj gj¼1 to be the most M active dimensions. The criterion defined in Eq. (16) gives a straightforward approach to find the most active random parameters. We note that some different approaches were proposed in Ma and Zabaras (2010) and Yang et al. (2011) to find the most active dimensions for more general situation. Their basic idea is using Fourier amplitude sensitivity test (Cukier et al., 1975) and some extra computation efforts are required to pick up most active dimensions from the component of the high-dimensional random vector

81

X. He et al. / Journal of Hydrology 478 (2013) 77–88

H. According to practical situation, we can also adjust the value of f in Eq. (16) such that M is much less than N. We define an index set M ¼ f1; 2; . . . ; Mg. Then, we define a new truncated HDMR as following,

b MÞ þ w wðHÞ ¼ wð H Mþ1 ðhMþ1 Þ þ    þ wN ðhN Þ;

ð17Þ

b M ¼ ðh1 ; h2 ; . . . ; hM ;  where H hMþ1 ; . . . ;  hN Þ is an M-dimensional random vector, and wi(hi)(i = M + 1, . . . , N) are defined in the formulas b M Þ depends on the intermediate-dimensional (15). The term wð H b M , which contains all of the most active random random vector H b M Þ has main contribution to the output variables. Therefore, wð H

w(H) of the underlying stochastic model. Compared with the traditional HDMR truncated to 2nd order terms, the proposed HDMR defined in Eq. (17) has the following advantages: (a) Since the proposed HDMR inherently includes the all cooperative contribution from the most active dimensions, better approximation accuracy may be achieved in the proposed HDMR than the traditional truncated HDMR. (b) The proposed HDMR consists of many fewer terms than the traditional truncated HDMR, the total computational effort of the proposed HDMR may be much less than the traditional truncated HDMR. (c) Since all terms in the expansion of the proposed HDMR are uncorrelated, the total variance of the proposed HDMR can be represented as the summation of the variances of all terms in the right-hand side of (17). This is desirable for calculating variance. We note that the truncated HDMR and MsFEM are independent of each other. The truncated HDMR is to reduce the complexity from a highdimensional random parameter space and MsFEM is to coarsen the underlying fine-scale deterministic models at given realizations. If high level (e.g., level 3 and above) sparse grid collocation is reb MÞ quired to approximate the intermediate-dimensional term wð H in (17), then we can further use the traditional truncated HDMR to b M Þ to reduce computation complexity. If approximate the term wð H b M Þ is smooth with respect to H b M , then level 2 sparse grid colwð H location method can provide an accurate approximation and the proposed HDMR (17) can be directly utilized. Many problems fall into this case. If the number M of active dimensions is much larger than N/2, we may use adaptive HDMR (Ma and Zabaras, 2010; Yang et al., 2011) instead. However, we find that the number M of active dimensions is usually less than N/2 for most practical problems as N is large.

As stated in Section 3, the sparse grid stochastic collocation method reduces the stochastic parabolic problem (1) into a series of deterministic collocation equations (see, Eq. (12)). However, it suffers from curse of dimensionality with increasing the number of dimensions. Integrating HDMR and SGSCM is a practical way to overcome this difficulty. For the simplicity of the notation, the formula (12) is written as

X X s wðHk Þ  li;k ;

ð18Þ

jij6Nþr k

 N1  li;k , and r is the interpolation N þ r  jij level. Combining the proposed HDMR formula (17) and (18) we can define an interpolation formula for the approximate solution of Eq. (8) which is given by s

where li;k ¼ ð1ÞNþrjij



b M ÞÞ þ wðHÞ ¼ AðwðHÞÞ ¼ Aðwð H

N X

j¼Mþ1jij j61þr kj

jiM j6Mþr kM



s lij ;kj

 ðN  MÞw0 ;

ð20Þ

b M;k ¼ ðhi ;k ; hi ;k ; . . . ; hi ;k ; h Mþ1 ; . . . ; h N Þ denote the sparse where H M M M 1 1 2 2 grid collocation points in the M-dimensional stochastic space !M, b j;k ¼ ð and H h1 ; . . .  hj1 ; hi ;k ;  hjþ1 ; . . . ;  hn Þ denote the sparse grid colloj

j

j

cation points in the one-dimensional stochastic space cj. Thus, an Ndimensional sparse grid stochastic collocation problem is split into to a lower M-dimensional and a few one-dimensional problems using the proposed HDMR. We would like to note that the number of collocation points in the expansion (20) is defined as the sum of the number of collocation points for each sub-problem. By Eq. (20), to get solution w(H), we only need to compute b M;k Þ and wð H b j;k Þ at collocation points in the stochastic subwð H M

j

spaces !M and cj(j = M + 1, . . . , N), respectively. Let Hl denote a colb M;k 2 !M or H b j;k 2 cj , which is used in Eq. (20). location point H M

j

Let Ns be the total number of collocation points in !M and cj(j = M + 1, . . . , N). Then Eq. (20) is reduced to: for each hydraulic conductivity Kl: = K(x, Hl) at collocation point Hl, l = 1, 2, . . ., Ns, we solve the deterministic problem: find wl:¼w(x, t, Hl) such that for l = 1, 2, . . ., Ns,

d ðSs wl ; /Þ þ ðK l rwl ; r/Þ þ dt

Z

q/ds ¼ ðR; /Þ;

CN

8/ 2 H10 ðDÞ: ð21Þ

This set of deterministic equations can be solved by deterministic solvers. However, at each set of collocation points, the highresolution simulation is still needed to resolve the fine-scale heterogeneity. It would be computationally expensive if we directly carry out the numerical simulation on a fine mesh for the flow problems. To enhance efficiency, we use MsFEM to solve the deterministic equation (21) in the paper. Meanwhile, Eq. (20) also provides an efficient approach to compute mean and variance of w(x, t, H). The idea is to use the sparse grid based on Gaussian quadrature rule (Xiu and Hesthaven, 2005) to integrate the component functions of the expansion (20). The head mean is obtained directly through the weighted sum of the function values at the collocation points. That is, the mean of w(x, t, H) is given by

X X b M;k Þ wiM ;kM wðx; t; H M

 E½wðx; tÞ  w

jiM j6Mþr kM

4.3. Integration of HDMR and SGSCM

AðwÞ :¼ AðN þ r; NÞðwÞ ¼

N X X X X X b M;k Þ  ls b j;k Þ wð H wð H iM ;kM þ M j

wðHÞ ¼

Aðwj ðhj ÞÞ:

ð19Þ

j¼Mþ1

Substituting the first-order component functions in Eq. (15) into the first-order components on the right-hand side of Eq. (19), we have

N X X X b j;k Þ  ðN  MÞw ; ð22Þ wij ;kj wðx; t; H 0 j

þ

j¼Mþ1jij j61þr kj

R

s

where wiM ;kM ¼ !M liM ;kM qðHM ÞdHM and HM = (h1, . . . , hM), and R s wij ;kj ¼ cj lij ;kj qðhj Þdhj denote the weights, and q is the (joint) probability density functions. To compute the variance of the head, we can similarly construct an approximation for w2 and use the formula Var[w](x, t) = E[w2](x, t)  (E[w](x, t))2 for each of the terms in right-hand side of (17). The weights in the formula (22) are evaluated by applying Gaussian quadrature rule in the multiple dimensions. Clearly, this is a pre-processing procedure which can be done completely parallel. In fact, we had built up the database of the weights for the different dimensions in advance. These weights can be repeatedly used when we make the integration with respect to the random variables with the same distribution. 5. Multiscale finite element formulation The key ingredient of MsFEM is the construction of an appropriate multiscale finite dimensional space in which the solution is sought. In particular, the fine-scale heterogeneity in the hydraulic

82

X. He et al. / Journal of Hydrology 478 (2013) 77–88

conductivity realization Kl will be imbedded in this finite dimensional space. Then, this information is incorporated into the coarse-scale formulation through the coarse-scale stiffness matrix. 5.1. Basis functions in MsFEM For MsFEM, a coarse partition is first performed in the domain D. Let TH be a quasi-uniform partition of D and K be a representative coarse-scale element with diamðKÞ ¼ HK . Define H ¼ maxfHK ; K 2 TH g. Following the work by Hou and Wu (1997), we define multiscale basis functions by /K;i for nodes i = 1, . . ., d of coarse element K for the problem (21), which satisfy

 r  ðK l r/K;i Þ ¼ 0 in K;

ð23Þ

/K;i ¼ lK;i

ð24Þ

on @K;

wH ðx; tÞ ¼

Nc X wj ðtÞ/j ;

ð27Þ

j¼1

where Nc is the total number of nodes over the coarse-scale partition TH ; wj ðtÞ is the nodal point value of w(x, t) at node j, and /j 2 XH is the multiscale base function at node j. Substituting Eq. (27) into Eq. (26), we get Nc Nc X dwj ðtÞ X ð/i ; Ss /j Þ ðr/i ; K r/j Þwj ðtÞ þ dt j¼1 j¼1 Z ¼ ð/i ; RÞ  /i qds;

ð28Þ

CN

where i = 1, . . ., Nc. This gives a linear ordinary differential system about

c fwj ðtÞgNj¼1 .

Let

W ¼ ðw1 ; w2 ; . . . ; wNc ÞT ; cij ¼ ð/i ; Ss /j Þ; R

where lK;i is the boundary condition associated with node i along the boundary @K of the coarse element K, the subscript l denotes a realization of hydraulic conductivity. There exist some options for the boundary condition lK;i (see, He and Ren, 2006; Hou and Wu, 1997; Jiang et al., 2007). Eq. (23) defines basis functions for local MsFEM if lK;i using local information. Incorporating global information into lK;i leads to global MsFEM (Jiang et al., 2007). Although there are many different choices for the boundary condition lK;i , in the paper we only uses the oscillatory boundary condition proposed by Hou and Wu (1997). This is because the MsFEM basis functions with oscillatory boundary condition can accurately span the solutions of the flow models discussed in Section 6. For compressible problem, it may be required to dynamically update multiscale basis functions using an adaptive strategy (He and Ren, 2009). We now introduce the multiscale finite element space for the MsFEM as follows:

aij ¼ ðr/i ; K r/j Þ, and bi ¼ ð/i ; RÞ 

V H ¼ spanf/K;i : i ¼ 1; . . . ; d;

where Wm and bm denote the values of W and b at tm, respectively. At each time level m + 1, the solution Wm+1 of the system of algebraic Eq. (30) incorporating the boundary conditions is obtained using the value of Wm by the iterative process, where m = 0, 1, . . . It is sometimes of interest to recover the fine-scale solution in some regions. Here we explain how to reconstruct the fine-scale solution from MsFEM solution. The idea is to use the constructed multiscale basis functions to project the MsFEM solution onto fine-grid (cf., Eq. (27)). In our numerical experiments, we use the reconstructed MsFEM solution to measure the error of the stochastic multiscale models.

for all K 2 TH g:

ð25Þ

In general, the base functions /K;i have no analytical solutions. Their numerical solutions are obtained by solving the local elliptic problem (23) with specific boundary condition on a partition of the coarse element K with a mesh size resolving the small-scale details using appropriate numerical method such as the standard FEM. We pre-compute the multiscale basis functions for each realization of hydraulic conductivity. These MsFE basis functions are repeatedly used at all times for the deterministic parabolic problem (21) with a given hydraulic conductivity realization Kl. We note that the MsFE space VH depends on the collocation point Hl (l = 1, 2, . . . , Ns) for Eq. (21). Once we obtain the MsFE basis functions, we follow the standard finite element procedure to compute the MsFEM solution of Eq. (21) for each sample Hl. 5.2. Numerical formulation of MsFEM The formulation of MsFEM shares the same procedure as the standard finite element method except for the basis functions. Here, we present MsFE formulation for the deterministic Eq. (21). For simplicity of presentation, we omit the subscript l and assume the deterministic equation is satisfied for an arbitrary conductivity realization in the stochastic space. We first introduce the semidiscrete approximate solution wH(x, t) of Eq. (21), obtained by projecting (21) onto the subspace V H H1 ðDÞ, that is

d ðSs wH ; /i Þ þ ðK rwH ; r/i Þ þ dt

8/i 2 V H ;

Z

q/i ds ¼ ðR; /i Þ;

CN

ð26Þ

where /i is a multiscale basis function. The approximated solution wH(x, t) is in the form of

CN

/i qds, where ( )T refers to

transposition, and i = 1, . . ., Nc, j = 1, . . ., Nc. In matrix notation, Eq. (28) can be rewritten as

dW þ AW ¼ b; dt

C

0 < t 6 T;

ð29Þ

where C is the matrix with the (i, j)th element cij, A is the matrix with the (i, j)th element aij, and b is the vector with the ith element bi. The Euler method or the others can be used to discretize the semidiscrete Eq. (29). In this study, we utilize the Crank-Nicolson (CN) method to construct the fully discretized version of Eq. (29). Let Dt be time step and tm = mDt (m = 0, 1, . . .). Then the CN discretization of Eq. (29) is



   mþ1 m C A C A b þb Wmþ1 ¼ Wm þ þ  ; Dt 2 Dt 2 2

ð30Þ

6. Numerical experiments In this section, we consider the transient flow problem (1) in the absence of sources and sinks and present a few numerical experiments for random hydraulic conductivities with different correlation lengths l and variances. The numerical results will demonstrate the performance of the proposed stochastic dimension reduction MsFEM. We assume that the specific storage and the thickness of the aquifer are 8  104 m1 and 10 m, respectively. Let Y(x, H) be a stochastic field and the conductivity field K(x, H) = exp (Y(x, H)). For the numerical experiments, we use a two point exponential covariance function for the stochastic field Y,

cov ½Yðx1 ; y1 ; x2 ; y2 Þ ¼ r2 exp 

jx1  x2 j2 2

2lx



jy1  y2 j2 2

2ly

! ;

ð31Þ

where (x, y) is the spatial coordinate, r2 is the variance of stochastic field Y, and lx and ly denote the correlation length in the x-and y-direction, respectively. Given the covariance, Y is numerically generated by using truncated KLE. For all the examples, we consider the same spatial domain with different random conductivity fields.

X. He et al. / Journal of Hydrology 478 (2013) 77–88

The spatial domain D is a square with size 200 m  200 m. A uniform finite element mesh is constructed by dividing the domain D into Nx  Ny sub-rectangles, and is referred as a Nx  Ny mesh. MsFEM is performed on a 16  16 coarse mesh in the spatial domain unless otherwise stated, while the fine-scale solution is computed on a 128  128 fine mesh in all the examples. That is to say, each coarse-scale element contains 8  8 fine sub-rectangular elements. The stochastic MsFEM solution using the proposed HDMR technique is compared with the stochastic fine-scale FEM solution solved in a full stochastic space. We employ the sparse grid stochastic collocation method in the full stochastic space to solve the flow Eq. (1) on the fine mesh. This solution serves as the reference solution in the paper. To simplify the presentation, we use the following abbreviations: MsFEM-HDMR, SGSCM-F, FEM-HDMR and MsFEM-Full, where MsFEM-HDMR stands for the stochastic dimension reduction multiscale finite element method with the proposed HDMR technique, SGSCM-F for the sparse grid stochastic collocation method in the full stochastic space and FEM on fine mesh, FEM-HDMR for the stochastic FEM on fine mesh with the proposed HDMR and sparse grid stochastic collocation technique, and MsFEM-Full for the stochastic multiscale finite element method in the full stochastic space. The numerical simulations are implemented on a PC with Intel Core i7-2600 CPU. In all the examples below, we take the relative L2 norm to measure the error between MsFEM-HDMR and SGSCM-F. The relative L2 error of the mean or variance between the head by stochastic MsFEM-HDMR and the reference head by SGSCM-F is defined by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PNf  2  i¼1 ðwms ðxi Þ  wr ðxi ÞÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; L2 ¼ PNf  2 i¼1 wr ðxi Þ

ð32Þ

where Nf is the total number of fine-grid nodes in the spatial do ms and w  r denote the mean or variance of the head by stomain, w chastic MsFEM-HDMR and SGSCM-F, respectively. The error of FEM-HDMR and MsFEM-Full solutions is similarly measured by using (32).

time. Meanwhile, we also note that the results of FEM-HDMR are the most accurate. By Fig. 1, we find that the independent effect of the proposed truncated HDMR is much smaller than the independent effect of MsFEM. MsFEM provides the dominant contribution to the total error between MsFEM-HDMR and SGSCM-F. To assess the computational efficiency of different methods, we record the CPU time for the different methods on the same computer. The CPU time is 0.64, 2.52, 16.3 and 58.3 h, respectively for MsFEM-HDMR, MsFEM-Full, FEM-HDMR and SGSCM-F to achieve the mean and variance of the head at time t = 20 h (Dt = 5 min). This demonstrates that the proposed MsFEM-HDMR provides a great computational efficiency in this example. The CPU time for the stochastic coarse grid model (16  16) with the proposed HDMR technique is only about 1.1% of that for the fine grid model (128  128) using sparse grid collocation in full stochastic space. We note that the total number of the sparse grid collocation points is 5101 for the full stochastic space N = 50. However, when M = 25 is taken in the proposed HDMR, the total number of the collocation points is only 1426. This shows that the HDMR offers a considerable saving in computing time due to the decreasing of the collocation points in stochastic space. The coarsening of mesh in MsFEM also significantly enhances the computational efficiency. By comparison, it is apparent that the new truncated HDMR reduces the CPU time of the full method by a factor of 3.6. Similarly, a reduction in CPU time by a factor of 23.1 is achieved by employing MsFEM in the full dimensional space. Finally, combining these approaches reduces the CPU time by nearly two orders of magnitude, with the dominant contribution from the use of MsFEM. Using the MsFEM-HDMR approach will lead to even greater gains in efficiency for larger spatial scale and higher-dimensional stochastic problems. In addition, the distributions of the head mean and head variance at time t = 11 h on the spatial domain obtained from MsFEM-HDMR and SGSCM-F are depicted in Figs. 2 and 3, respectively. As we can see from these two figures, the head mean obtained by stochastic MsFEM-HDMR matches well the reference fine-scale mean, and the head variance by MsFEM-HDMR also provides good agreement with the reference head variance. 6.2. Effect of active dimensions

6.1. Responding to a sudden change In this subsection, we consider the aquifer’s response to a sudden change in reservoir level. Specifically, we simulate changes in head through time if, at t = 0, we suddenly drop the water level in the reservoir on @D from 15 m to 12 m. Then we impose the following initial condition and boundary condition for the transient flow problem:

wðx; 0Þ ¼ 15 m in D; and wðx; tÞ ¼ 12 m on ð0; T  @D:

83

ð33Þ

We fix a time step size Dt = 5 min for the temporal discretization, and interpolation level r = 2 for the Smolyak sparse grid collocation. We first present a set of simulation results for a stochastic conductivity field Y with the variance r2 = 1.5. Assume the mean of the stochastic field Y (or lnK) equals 6.215. The correlation structure of the conductivity field is anisotropic with lx = 25 m and ly = 30 m. We truncate N = 50 terms in KLE to get random field Y, which represents 99% of the total variance of the exponential covariance P function, i.e. the spectrum energy ratio Ni¼1 ki =ðr2 jDjÞ  0:99. To determine the active dimensions in the proposed HDMR, we take the threshold f = 0.9 in Eq. (16) and the number of active dimensions is M = 25. Fig. 1 shows the relative L2 error for the mean and variance of the hydraulic head from different methods at times t = 8, 11, 14, 17, and 20 h. From Fig. 1, we can observe that all methods give a reasonable accuracy for the mean and variance of the head compared with the reference SGSCM-F solution, and the relative error of the mean remains less than 0.3% through all the

In this subsection, we concentrate on the effect of the active dimensions in the proposed stochastic MsFEM-HDMR. For the simulation, we consider a stochastic conductivity field Y with variance r2 = 2.0. The mean of the stochastic field Y (or lnK) equals 6.215, and its correlation structure is anisotropic with lx = 15 m and ly = 20 m. To ensure that the truncated KLE can represent an approximate 98% of the total variance of the exponential covariance function, we truncate N = 100 terms in KLE to numerically get random field Y. Then, in order to investigate the effect of the active dimensions on the accuracy and efficiency of the stochastic MsFEM-HDMR, we choose a variety of threshold fs to obtain M in the new truncated HDMR expansion. In general, with increasing f, more active dimensions are used. In particular, all of the random dimensions are included when f = 1. By extensive numerical studies for f  0.89, we found that the relative L2 error of the head mean using the stochastic MsFEM-HDMR is usually less than 0.33%, and the relative L2 error of the variance less than 7.5%. Below, we present a set of computational results for different fs. We choose five different threshold values f = 0.73, 0.83, 0.89, 0.93 and 1.0. Consequently, the number of the active dimensions in the proposed HDMR is M = 30, 40, 50, 60 and 100, respectively. Fig. 4 depicts the relative L2 error for the mean and variance of the head at time t = 11 h by the stochastic MsFEM-HDMR with the five different numbers of active dimensions. From Fig. 4, we can see that the relative error of the mean and variance of the head first decreases monotonically and then plateaus as M increases further.

84

X. He et al. / Journal of Hydrology 478 (2013) 77–88 0.2 FEM−HDMR MsFEM−Full MsFEM−HDMR

0.009 0.008

FEM−HDMR MsFEM−Full MsFEM−HDMR

0.18 0.16

0.007

0.14

0.006

0.12

0.005

0.1

0.004 0.08 0.003 0.06

0.002

0.04

0.001 8

10

12

14

16

18

20

0.02

8

10

12

14

16

18

20

Fig. 1. The relative errors for the head mean (left) and variance (right) from different methods at different times.

12.1 12.4 12.6

200 180

13.2

13

13.1

140

13.2

12.4 12.6

180

12.8

160

200

12.8 13

13.1

160 140

12.8

12.8 120

120

100

12.6

80

12.4

13 12.7 12.5

40 20

12.2

13.2 12.4

60

13 12.2 12.7

40

12.2 50

100

150

200

12

0

12.2

12.5

20

0 0

12.6

80

13.2

60

100

0

50

100

150

200

12

Fig. 2. Comparison of the mean between the stochastic MsFEM-HDMR and reference solution at time t = 11 h. Case N = 50, M = 25.

200

0.241 0.402 0.4

0.482

180

0.563 0.724

160

200 0.8

0.595 0.68 0.51

180 0.7

160

0.804 0.724 0.563

140

0.6 0.5

0.804 0.402 100

50

0.

100

100

150

0.51

60 0.2

0.482

0

0.7 0.6

200

0.1 0

0.5

0.425 0.34

80

0.482 0.643 0.322 0.161

20

0.85

120

0.3

0.804

40

0

0.4

0.8040.482 0.643

60

0.8

0.595

120

80

140

0.255 0.425 0.5 0.595 0.6 0.51 0.765 0.85 0.765

0.51 0.68 0.85

0.4

0.85

0.3

0.68 0.17 0.34

0.51 0.595 0.68

40

0.51 0.59 0.6

20 0

0

50

100

150

200

0.2 0.1 0

Fig. 3. Comparison of the variance between the stochastic MsFEM-HDMR and reference solution at time t = 11 h. Case N = 50, M = 25.

85

X. He et al. / Journal of Hydrology 478 (2013) 77–88 0.2 Mean of head Variance of head

0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0

30

40

50

60

70

80

90

100

Fig. 4. Relative errors of the head mean and variance at time t = 11 h from stochastic MsFEM-HDMR with different active dimensions.

140

120

100

80

60

40

20

0

M=30

M=40

M=50

M=100 SGSCM−F

M=60

Fig. 5. Comparison of CPU time among the stochastic MsFEM-HDMR with different active dimensions M and reference fine-scale model.

200

6.3. Effect of correlation length In the subsection, we aim to study the impact of correlation length on the accuracy of the stochastic MsFEM-HDMR. We fix 200

12.4 12.7

180

Also, we observe that the relative error of the head mean remains very small through all the five M values of the active dimensions, and the relative error of the variance rapidly decrease for the first three M values of the active dimensions. This indicates that there exists an optimal dimension number M for the HDMR dimension reduction. To assess the effect of the active dimensions on the computational efficiency of the stochastic MsFEM-HDMR, we record the computation time for the different methods on the same computer. Fig. 5 shows the comparison of CPU time among the stochastic MsFEM-HDMR with different active dimensions M and reference fine-scale model. It is apparent that all the stochastic MsFEM-HDMR with different Ms save considerable CPU time compared with SGSCM-F, where the CPU time is 120.83 h to obtain the fine-scale reference solution at time t = 11 h (Dt = 5 min), while the CPU time is 0.83, 1.35, 2.01, 3.13 and 7.7 h, respectively for the stochastic MsFEM-HDMR solutions at time t = 11 h with the different active dimensions M = 30, 40, 50, 60 and 100. In addition, we note that the total number of Smolyak sparse grid collocation points with level r = 2 is 2211, 3581, 5351, 7521, and 20,201, respectively, for the stochastic MsFEM-HDMR with M = 30, 40, 50, 60 and 100. This affirms that the computational cost will significantly increase as M increases. Thus the advantage of using the new truncated HDMR with a moderate M is obvious. Based on our numerical experiments, the stochastic MsFEM-HDMR produces very good results with significantly less computational cost when the threshold f 2 (0.88, 0.90). Figs. 6 and 7 compare the head mean and head variance at time t = 11 h on the spatial domain obtained from the stochastic MsFEM-HDMR with f = 0.89 (i.e., M = 50) with the reference results by SGSCM-F. From these two figures, we observe that the head mean for stochastic MsFEM-HDMR is in excellent agreement with the reference fine-scale mean, and the MsFEMHDMR variance is a good approximation of the reference variance. Finally, we would like to emphasize that in this case (M = 50), the CPU time needed to achieve the stochastic MsFEM-HDMR results at time t = 11 h is only about 1.66% of that for the reference results by SGSCM-F.

13.2

12.4 12.7

180

160

13.2

12.9

12.9 160

13.2

13.2

13

13

140

140 12.8

120

12.8

120 100

100

13.3 12.6

12.6 80 60

12.4

13.1 12.8

40

0

13.1 12.8 12.5 12.3

60

12.4

40

12.5 12.3

20 0

80

13.3

12.2

12.2 20

100

200

12

0

0

100

200

12

Fig. 6. Comparison of the mean between the stochastic MsFEM-HDMR and reference solution at time t = 11 h. Case N = 100, M = 50.

86

X. He et al. / Journal of Hydrology 478 (2013) 77–88

200

0.553

180

0.553 160

0.8

0.237 0.395 0.47 0.

0.395 0.711 0.39 0.474

140

160

0.632 0.79 0.79

0.3

0.474 0.632 0.474 0.395 0.316 0.158 0.553

60 40

0

100

0.6

120

80

0.2

0.1

200

0.4

0.811 0.486 0.649 0.568 0.568 0.162 0.486 0.4

40

0.4

0

0.5

0.811 0.73 0.649

100

60

20

0.7

140

0.4

80

0.8

0.6

0.79 0.711 0.79

100

0

180 0.7

0.5

120

0. 0.243 0.405 0.568 0.649 0.6 0.486 0.486 0.73 0.568 0.5

200

0.324 0.649

0.3

0.2

0. 0.1

20 0

0

100

200

Fig. 7. Comparison of the variance between the stochastic MsFEM-HDMR and reference solution at time t = 11 h. Case N = 100, M = 50.

the variance of the conductivity field Y by r2 = 1.5, and consider three different correlation lengths: l = 15, 20, and 25 m. In the three cases, the correlation structure of the stochastic field is isotropic with lx = ly = l. In addition, all other related parameters with respect to the flow problem are the same as in Subsection 6.1. Since the variance of the covariance function is fixed, the change of correlation length alters the energy ratio of each random dimension in the truncated KLE (Huang et al., 2001). By using a fixed spectrum enPN 2 ergy ratio i¼1 ki =ðr jDjÞ  0:98, we truncate N = 120, 85 and 60 terms in the KLE, respectively for l = 15, 20, and 25 m. To investigate the impact of the correlation lengths on the performance of stochastic MsFEM-HDMR, we fix the threshold f  0.89 to obtain M in the proposed HDMR expansion for the three different correlation lengths. The number of active dimensions is M = 60, 42 and 30, respectively for l = 15, 20, and 25 m. Fig. 8 shows the relative L2 error for the mean and variance of the hydraulic head by the stochastic MsFEM-HDMR at time t = 11 h for the three different correlation lengths. From Fig. 8, we observe that the stochastic dimension reduction MsFEM gives very good accuracy for the mean and variance of the head in these cases. This figure also shows that although three correlation lengths result in different truncated stochastic dimensions from KLE and different active dimensions from the fixed threshold value, the stochastic MsFEM-HDMR gives nearly identical relative errors of the head mean and head variance for the three different correlation lengths. This demonstrates that the proposed stochastic dimension reduction MsFEM is efficient and robust. Fig. 9 depicts the comparison of the head mean and head variance obtained from the stochastic MsFEM-HDMR with those from reference fine-scale results at time t = 11 h in section y = 100 m for the three different correlation lengths. We see that the head means obtained by stochastic dimension reduction model on the coarse-scale mesh are able to match the reference fine-scale means in all the cases, and that the coarse-scale variances from stochastic MsFEM-HDMR can effectively capture the large-scale behavior of the fine-scale variances. 6.4. Effect of spatial variability In this subsection, we fix the correlation length of the conductivity field Y by lx = 25 m and ly = 30 m, and then investigate the

0.2 Mean of head Variance of head

0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 15

16

17

18

19

20

21

22

23

24

25

Fig. 8. Relative errors of the head mean and head variance from stochastic MsFEMHDMR at time t = 11 h for three different correlation lengths.

impact of the spatial variability r2 to the accuracy of the proposed stochastic MsFEM-HDMR. We consider three different variances r2 = 1.5, 2.0 and 2.5, and assume all other related parameters are the same as in the subsection 6.1 for flow problem. We note that the increase of r2 results in the larger variability in hydraulic head. The fixed correlation length implies that almost the same energy ratio is kept for each stochastic dimension in the truncated KLE for the three different variances (Huang et al., 2001). In all the three cases, the number of stochastic dimensions is N = 50 based P on a fixed spectrum energy ratio Ni¼1 ki =ðr2 jDjÞ  0:99 in the truncated KLE. The number of active dimensions is M = 25 in new truncated HDMR expansions for three different spatial variabilities by fixing the threshold f  0.90. In such a situation, it is expected that the accuracy of the stochastic MsFEM-HDMR depends on the ability of deterministic MsFEM capturing the spatial variability on the coarse mesh. Fig. 10 depicts the relative L2 error for the mean and variance of the head at time t = 11 h by stochastic MsFEM-HDMR on two different coarse meshes: 16  16 and 32  32. We observe that the stochastic MsFEM-HDMR gives a good approximation in the head mean and head variance. As expected, the errors in the head mean and head variance increase with increasing

87

X. He et al. / Journal of Hydrology 478 (2013) 77–88 14.5

2.5

14

2

13.5

1.5

13

1

12.5

0.5

12

0 0

50

100

150

200

0

50

100

150

200

Fig. 9. Comparison of the head mean and head variance obtained from the stochastic MsFEM-HDMR with those from reference fine-scale model at time t = 11 h in section y = 100 m for three different correlation lengths.

0.2 Mean of head, 16 × 16 mesh Variance of head, 16 × 16 mesh Mean of head, 32 × 32 mesh Variance of head, 32 × 32 mesh

0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 1.5

1.6

1.7

1.8

1.9

2

2.1

2.2

2.3

2.4

2.5

Fig. 10. Relative errors of the head mean and head variance from stochastic MsFEM-HDMR at time t = 11 h for three different spatial variabilities.

14.5

r2. Although the stochastic MsFEM-HDMR with 16  16 mesh maintains quite good accuracy for the head mean, the relative error for the head variance grows as r2 increases. However, its performance is still reasonable even with r2 = 2.5, which corresponds to very high variability in the physical space. Specifically, we observe that many conductivity realizations vary by seven orders of magnitude when r2 = 2.5. This poses a significant challenge for MsFEM to capture the influence of fine-scale heterogeneities on a coarse mesh. As illustrated in Fig. 10, the improvement of the results in the head variance is obvious as the coarse mesh is refined to 32  32 mesh. This again shows that the proposed HDMR is efficient and robust for stochastic dimension reduction. Here, we would like to note that although the coarse mesh is refined, the stochastic MsFEM-HDMR with 32  32 mesh can still save 97.5% CPU time compared with the reference SGSCM-F at the end of the simulation (t = 11 h). Fig. 11 depicts the comparison of the head mean and head variance by the stochastic MsFEM-HDMR with 16  16 mesh with those from reference fine-scale model at time 2.5

σ2= 1.5, fine−scale

σ2= 1.5, fine−scale

2

σ = 1.5, coarse−scale σ2= 2.0, fine−scale

14

2

σ = 2.0, coarse−scale

2

σ = 1.5, coarse−scale σ2= 2.0, fine−scale

2

σ2= 2.0, coarse−scale

σ2= 2.5, fine−scale

2

σ = 2.5, fine−scale

2

σ = 2.5, coarse−scale

σ2= 2.5, coarse−scale

13.5

1.5

13

1

12.5

0.5

12

0

0

50

100

150

200

0

50

100

150

200

Fig. 11. Comparison of the head mean and variance obtained from the stochastic MsFEM-HDMR with those from reference fine-scale model at time t = 11 h in section y = 100 m for three different spatial variabilities.

88

X. He et al. / Journal of Hydrology 478 (2013) 77–88

t = 11 h in section y = 100 m for the three different variances. From Fig. 11, we see that the head mean and head variance obtained with the stochastic dimension reduction model on the coarse mesh approximate the reference fine-scale results in the case r2 = 1.5 quite well. We also observe that the coarse-scale mean and variance deviate gradually from the fine-scale results as r2 increases from 2.0 to 2.5. However, the stochastic MsFEM-HDMR still captures the large-scale behavior of the fine-scale results even when r2 = 2.5. 7. Conclusions In this paper, a stochastic dimension reduction multiscale finite element method for the groundwater flow problems in random heterogeneous porous media has been developed, implemented and analyzed. To reduce high-dimensional stochastic complexity, a new truncated HDMR technique has been proposed and decomposes a high-dimensional stochastic problem into a moderatedimensional stochastic problem and a few one-dimensional stochastic problems. Each derived stochastic sub-problem is then reduced into a set of deterministic problems by applying SGSCM. These deterministic problems are uncoupled and can be solved by any appropriate numerical method. However, at each collocation point, a high-resolution simulation is still needed to resolve the fine-scale heterogeneity of porous media. It may be computationally expensive and even infeasible to capture the heterogeneity on a very fine grid in the large-scale spatial domain. Thus, we combined the SGSCM with MsFEM to solve each low-dimensional stochastic sub-problem from the proposed HDMR efficiently. Using MsFEM, the fine-scale heterogeneity of the water flow can be effectively captured in coarse-scale computations. The numerical examples with different stochastic conductivity fields have been studied to demonstrate the accuracy and efficiency of the proposed method. The numerical results have shown that using a proper threshold, the proposed stochastic dimension reduction multiscale method can provide a good approximation to the reference finescale results and leads to a reduction of nearly two orders of magnitude in CPU time, compared with SGSCM-F. In conclusion, the proposed method can efficiently capture the large-scale behavior of the solution on a coarse grid without resolving all the smallscale details in full stochastic space, and hence, results in a substantial reduction in computational cost. Acknowledgments We thank the anonymous reviewers for many constructive comments which improve the paper. X. He acknowledges that the work is funded by the National Natural Science Foundation of China Grant 41272271 and Hunan Provincial Natural Science Foundation of China Grant 11JJ3057. L. Jiang acknowledges the Chinese NSF 10901050 for the work. L. Jiang and J.D. Moulton acknowledge the support by the Department of Energy at Los Alamos National Laboratory under Contracts DE-AC52-06NA25396 and the DOE Office of Science Advanced Computing Research (ASCR) program in Applied Mathematical Sciences. References Babus˘ka, I., Nobile, F., Tempone, R., 2007. A stochastic collocation method for elliptic partial differential equations with random input data. SIAM J. Numer. Anal. 45, 1005–1034.

Barthelmann, V., Novak, E., Ritter, K., 2000. High dimensional polynomial interpolation on sparse grids. Adv. Compuat. Math. 12, 273–288. Chang, H., Zhang, D., 2009. A comparative study of stochastic collocation methods for flow in porous media. Commun. Comput. Phys. 6 (3), 509–535. Cukier, R.I., Schaibly, J.H., Shuler, K.E., 1975. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients: III. Analysis of the approximations. J. Chem. Phys. 63, 1140–1149. Dagan, G., 1982. Analysis of flow through heterogeneous random aquifers 2. Unsteady flow in confined formations. Water Resour. Res. 18 (5), 1571–1585. Dostert, P., Efendiev, Y., Hou, T.Y., 2008. Multiscale finite element methods for stochastic porous media flow equations and application to uncertainty quantification. Comput. Methods Appl. Mech. Eng. 197, 3445–3455. E, W., Engquist, B., 2003. The heterogeneous multi-scale methods. Commun. Math. Sci. 1, 87–133. Efendiev, Y., Hou, T.Y., 2009. Multiscale Finite Element Methods: Theory and Applications. Springer, New York. Ganapathysubramanian, B., Zabaras, N., 2007. Sparse grid collocation schemes for stochastic natural convection problems. J. Comput. Phys. 225, 652–685. Ganapathysubramanian, B., Zabaras, N., 2009. A stochastic multiscale framework for modeling flow through random heterogeneous porous media. J. Comput. Phys. 228, 591–618. Ganis, B., Yotov, I., Zhong, M., 2011. A stochastic mortar mixed finite element method for flow in porous media with multiple rock types. SIAM J. Sci. Comput. 33, 1439–1474. He, X.G., Ren, L., 2006. A modified multiscale finite element method for well-driven flow problems in heterogeneous porous media. J. Hydrol. 329, 674–684. He, X.G., Ren, L., 2009. An adaptive multiscale finite element method for unsaturated flow problems in heterogeneous porous media. J. Hydrol. 374, 56–70. Hou, T.Y., Wu, X.H., 1997. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 134, 169–189. Huang, S.P., Quek, S.T., Phoon, K.K., 2001. Convergence study of the truncated Karhunen–Loève expansion for simulation of stochastic processes. Int. J. Numer. Methods Eng. 52, 1029–1043. Hughes, T.J.R., Feijoo, G.R., Mazzei, L., Quincy, J.B., 1998. The variational multiscale method – a paradigm for computational mechanics. Comput. Methods Appl. Mech. Eng. 166, 3–24. Jiang, L., Efendiev, Y., Ginting, V., 2007. Multiscale methods for parabolic equations with continuum spatial scales. DCDS Ser. B 8, 833–859. Li, G.Y., Wang, S., Rosenthal, C., Rabitz, H., 2001. High dimensional model representations generated from low dimensional data samples: I. mp-CutHDMR. J. Math. Chem. 30, 1–30. Loève, M., 1977. Probability Theory, fourth ed. Springer Verlag, Berlin. Ma, X., Zabaras, N., 2010. An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations. J. Comput. Phys. 229, 3884–3915. Ma, X., Zabaras, N., 2011. A stochastic mixed finite element heterogeneous multiscale method for flow in porous media. J. Comput. Phys. 230, 4696–4722. Matthies, H.G., Keese, A., 2005. Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. Comput. Methods Appl. Mech. Eng. 194, 1295–1331. Nguyen, N.C., 2008. A multiscale reduced-basis method for parameterized elliptic partial differentia equations with multiple scales. J. Comput. Phys. 227, 9807– 9822. Nobile, F., Tempone, R., 2009. Analysis and implementation issues for the numerical approximation of parabolic equations with random coefficients. Int. J. Numer. Methods Eng 80, 979–1006. Nobile, F., Tempone, R., Webster, C.G., 2008. An anisotropic sparse grid stochastic collocation method for elliptic partial differential equations with random input data. SIAM J. Numer. Anal. 46, 2411–2442. Rabitz, H., Allis, ö.F., Shorter, J., Shim, K., 1999. Efficient input-output model representation. Comput. Phys. Commun. 117, 11–20. Shi, L.S., Zhang, D.X., Lin, L., Yang, J.Z., 2010. A multiscale probabilistic collocation method for subsurface flow in heterogeneous media. Water Resour. Res. 46, W11562. http://dx.doi.org/10.1029/2010WR009066. Smolyak, S., 1963. Quadrature and interpolation formulas for tensor products of certain classes of functions. Soviet Math. Dokl. 4, 240–243. Xiu, D., Hesthaven, J.S., 2005. High-order collocation methods for differential equations with random inputs. SIAM J. Sci. Comput. 27, 1118–1139. Yang, X., Choi, M., Lin, G., Karniadakis, G.E., 2011. Adaptive ANOVA decomposition of stochastic incompressible and compressible flows. J. Comput. Phys. http:// dx.doi.org/10.1016/j.jcp.2011.10.028. Zhang, D., Li, L., Tchelepi, H., 2000. Stochastic formulation for uncertainty analysis of two-phase flow in heterogeneous reservoirs. SPE J. 5, 60–70.