{'fraction_non_alphanumeric': 0.042630230572160546, 'fraction_numerical': 0.011733561058923996, 'mean_word_length': 4.506536255055018, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A hybrid HDMR for mixed multiscale finite element method with application for flows in random porous media

Stochastic modeling has become a popular approach to quantify uncertainty in flows through heterogeneous porous media. The uncertainty in heterogeneous structure properties is often parameterized by a high-dimensional random variable. This leads to a deterministic problem in a high-dimensional parameter space and the numerical computation becomes very challengeable as the dimension of the parameter space increases. To efficiently tackle the high-dimensionality, we propose a hybrid high dimensional model representation (HDMR) technique, through which the high-dimensional stochastic model is decomposed into a moderate-dimensional stochastic model in a most active random space and a few one-dimensional stochastic models. The derived low-dimensional stochastic models are solved by incorporating sparse grid stochastic collocation method into the proposed hybrid HDMR. The porous media properties such as permeability are often heterogeneous. To treat the heterogeneity, we use a mixed multiscale finite element method (MMsFEM) to simulate each of derived stochastic models. To capture the non-local spatial features of the porous media and the important effects of random variables, we can hierarchically incorporate the global information individually from each of random parameters. This significantly enhances the accuracy of the multiscale simulation. The synergy of the hybrid HDMR and the MMsFEM reduces the stochastic model of flows in both stochastic space and physical space, and significantly decreases the computation complexity. We carefully analyze the proposed HDMR technique and the derived stochastic MMsFEM. A few numerical experiments are carried out for two-phase flows in random porous media and support the efficiency and accuracy of the MMsFEM based on the hybrid HDMR.

1. Introduction. The modeling of dynamic flow and transport processes in geologic porous media plays a significant role in the management of natural resources, such as oil reservoirs and water aquifers. These porous media are often created by complex geological processes and may contain materials with widely varying abilities to transmit fluids. Thus, multiscale phenomena are inherent in these applications and must be captured accurately in the model. In addition, due to measurement errors and limited knowledge of the material properties and external forcing, modeling of subsurface flow and transport often uses random fields to represent model inputs (e.g., permeability). Then the system's behavior can be accurately predicted by efficiently simulating the stochastic multiscale model. The existence of heterogeneity at multiple scales combined with this uncertainty brings significant challenges to the development of efficient algorithms for the simulation of the dynamic processes in random porous media. As a result, the interest in developing stochastic multiscale methods for stochastic subsurface models has steadily grown in recent years.
The existence of uncertainty in random porous media is an important challenge for simulations. One way to describe the uncertainty is to model the random porous media as a random field which satisfies certain statistical correlations. This naturally results in describing the flow and transport problem using stochastic partial differential equations (SPDE). Over the last few decades several numerical methods have been developed for solving SPDEs. The existing stochastic numerical approaches roughly fall into two classes: (1) non-intrusive schemes and (2) intrusive schemes. In non-intrusive schemes, the existing deterministic solvers are used without any modification to solve a (large) set of deterministic problems, which correspond to a set of samples from the random space. This leads to a set of outputs, which are used to recover the desired statistical quantities. Monte Carlo methods [14] and stochastic collocation methods fall into this category. In contrast, intrusive schemes treat the approximation of the probabilistic dependence directly in the SPDE. Thus, discretization leads to a coupled set of equations that are fundamentally different from discretizations of the original deterministic model, and hence, raises new challenges for nonlinear solvers. Typical examples from this class are stochastic Galerkin [19,13,34] and perturbation methods [48]. A broad survey of these methods can be found in [29]. Among these methods, stochastic collocation methods have gained significant attention in the research community because they share the fast convergence property of the stochastic finite element methods while having the decoupled nature of Monte Carlo methods. A stochastic collocation method consists of two components: a set of collocation points and an interpolation operator. The collocation points are selected based on specific objectives, and choosing different objectives or constraints leads to different methods. Two popular collocation methods are the full-tensor product method [5] and the Smolyak sparse-grid method [6,36,40]. Sparse-grid stochastic collocation is known to have the same asymptotic accuracy as full-tensor product collocation, but uses significantly fewer collocation points. Unfortunately, the number of collocation points required in the sparse-grid method still increases dramatically for high-dimensional stochastic problems. Hence, application of this method is still limited to problems in a moderate-dimensional stochastic space.
Due to small correlation lengths in the covariance structure, the uncertainty in characterizations of porous media are often parametrized by a large number of random variables, e.g., using a truncated Fourier type expansion for random fields. Consequently, the model's input is defined in a high-dimensional random parameter space. Sampling in high-dimensional random spaces is computationally demanding and very time-consuming. If the sampling of random space is conducted in the full random space through the stochastic collocation method, then the number of samples increases sharply with respect to the dimension of the random space. This is the notorious curse of dimensionality, and is a fundamental problem facing stochastic approximations in high-dimensional stochastic spaces. Dimension reduction techniques such as proper orthogonal decomposition, principal component analysis, reduced basis methods [39], novel random field expansion techniques [27,45,46] and high-dimensional model representations (HDMRs) [38], strive to address this problem. Among these methods, the HDMR is one of the most promising methods for efficient stochastic dimension reduction in subsurface applications, and has received considerable attention in recent years [32,33,35,47]. The HDMR was originally developed in the application of chemical models by [38], but was later cast as 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, and can be employed to reduce the computational effort. In practice, to avoid dealing with the full high-dimensional random space, a truncated HDMR technique is used to decompose a high-dimensional model into a set of low-dimensional models, each of which needs much less computational effort. The curse of dimensionality can be suppressed significantly by using this approach to HDMR. However, some drawbacks exist in the traditional truncated HDMR [9,38] and related methods such as the adaptive HDMR [32,47]. For example, often too many low-dimensional models need to be computed, and high-order cooperative effects from important random variables may be neglected. To overcome these drawbacks, we propose a hybrid HDMR, which implicitly uses the complete HDMR in a moderate-dimensional space and explicitly uses the first-order truncated HDMR for the remaining dimensions. Thus the hybrid HDMR decomposes a highdimensional model into a moderate-dimensional model and a few one-dimensional models. The moderate-dimensional space is spanned by the most active random dimensions. We use sensitivity analysis to identify the most active random dimensions. Then each low-dimensional stochastic model in the hybrid HDMR is solved by using sparse-grid stochastic collocation method. The proposed hybrid HDMR renders a good approximation in stochastic space and substantially improves the efficiency for high-dimensional stochastic models. Therefore, a good trade-off between computational complexity and dimension reduction error can be achieved with the hybrid HDMR technique.
Porous media often exhibit complex heterogeneous structures that are inherently hierarchical and multiscale in nature, and hence, pose a significant challenge for developing accurate and efficient numerical methods. Simulating flow in porous media using a very fine grid to resolve this heterogeneous structure is computationally very expensive and possibly infeasible. However, disregarding the heterogeneities can lead to large errors. Thus, many multiscale methods including the MMsFEM [7], variational multiscale method [23], two-scale conservative subgrid approaches [3] and heterogeneous multiscale method [10], multiscale finite volume method [24], spectral MsFEM [11] (see [12] for more complete references), have been developed over the last few decades to capture the influence of fine-and multi-scale heterogeneities in underresolved simulations. We note that these multiscale methods share some similarities [12]. In this paper, we focus on the MMsFEM. The main idea behind MMsFEM is to incorporate fine-scale information into the finite element velocity basis functions, and hence, capture the influence of these fine scales in a large-scale mixed finite element formulation. The MMsFEM retains local conservation of velocity flux and has been shown to be effective for solving flow and transport equations in heterogeneous porous media. In many cases, the overhead of computing multiscale basis functions can be amortized as they can be computed once for a particular medium, and then reused in subsequent computations with different source terms and boundary conditions. This leads to a large computational savings in simulating the flow and transport process where the flow equation needs to be solved many times dynamically. When porous media exhibit non-separable scales, some global information is needed to represent non-local effects (e.g., channel, fracture and shale barriers) and is used to construct the multiscale basis functions. Using global information can significantly improve the simulation accuracy in these important geological settings [1,2,25,26]. If the global information is incorporated into MMsFEM, we refer to the MMsFEM as global MMsFEM (G-MMsFEM). In the paper, we extend the central concept of the global information to MMsFEM based on the hybrid HDMR.
Combining different multiscale methods and stochastic numerical approaches yields various stochastic multiscale methods [15,16,20,26,27,33,44]. These methods established a general approach to solve problems in multiscale stochastic media. Specifically, dimension reduction techniques are applied, particularly focusing on the multiscale spatial dimension of these stochastic models, to reduce the overall computational cost. For example, within the stochastic Galerkin framework several earlier works (e.g., [4,35,39]) proposed solving an optimization problem to select a set of optimal basis functions to define a reduced model. Although, these approaches reduce the resolution of the spatial dimension, the stochastic dimension is not reduced. Therefore, the computational cost of these approaches may still be too expensive in high-dimensional stochastic spaces.
In this paper, we propose a hybrid HDMR framework, and use the sparse-grid stochastic collocation method with MMsFEM to develop a reduced multiscale stochastic model. Sensitivity analysis is used to control the reduction of the stochastic dimension. For the MMsFEM, we hierarchically utilize approximate global information that is aligned with the terms of the hybrid HDMR. Specifically, only the required basis functions are computed and to amortize their computational cost they are held constant in time. Thus, the proposed multiscale stochastic model reduction approach is able to reduce a high-dimensional stochastic multiscale model in both stochastic space, and the resolution of physical space. Compared with traditional truncated HDMR techniques, much better efficiency and very good accuracy are achieved with this hybrid HDMR. We analyze the proposed multiscale stochastic model reduction approach and investigate its application to flows in random and heterogeneous porous media. Important statistical properties (e.g., mean and variance) of the outputs of the stochastic flow models are computed and discussed.
The rest of the paper is organized as follows. In Section 2, we briefly introduce the background on the flow and transport models in random porous media. In Section 3, we present a general framework for HDMR and propose the hybrid HDMR technique. Some theoretical results and computational complexity are also addressed in this section. Section 4 is devoted to presenting MMsFEM integration with the hybrid HDMR. In Section 5, numerical examples using two-phase flow are presented to demonstrate the performance of the proposed the hybrid HDMR with MMsFEM. Finally, some conclusions and closing remarks are made.

Background and notations.
2.1. Two-phase flow system and its stochastic parametrization. Let D be a convex bounded domain in R d (d = 2, 3) and (Ω, F, P ) be a probability space, where Ω is the set of outcomes, F is the σ-algebra generated by Ω, and P is a probability measure.
We consider two-phase flow and transport in a random permeability field, k(x, ω). Here the two phases are referred to as water and oil, and designated by subscripts w and o, respectively. The equations of two-phase flow and transport (in the absence of gravity and capillary effects) can be written: where the total mobility λ(S) is given by λ(S) = λ w (S)+λ o (S) and q is a source term.
Here λ w (S) = k rw (S)/µ w and λ o (S) = k ro (S)/µ o , where µ o and µ w are viscosities of oil and water phases, respectively, and k rw (S) and k ro (S) are relative permeabilities of oil and water phases, respectively. Here f w (S) is the fractional flow of water and given by f w = λ w /(λ w + λ o ), Equation (2.1) is the flow equation governing the water pressure, and (2.2) is the transport (or saturation) equation. According to Darcy's law, the total velocity u in (2.2) is given by A random field may be parameterized by a Fourier type expansion, such as the expansion of Karhunen-Loève (ref. [31]), polynomial chaos or wavelet [29]. This often gives rise to an infinite-dimensional random space. For computation, we truncate such an expansion to approximate the random field. Then k(x, ω) can be formally described by For example, if a random field is characterized by a covariance structure, then the random field can be approximated by a finite sum of uncorrelated random variables through a truncated Karhunen-Loève expansion (KLE). To obtain an accurate approximation, a large number N of random parameters is required in (2.4). This leads to a family of deterministic models in a high-dimensional random parameter space.
To simplify the presentation, we make the following assumption, Let I N ⊂ R N be the image of Θ, i.e., I N = Θ(Ω), and ρ(Θ) = Π N i=1 ρ i (θ i ) be the joint probability function of (θ 1 , · · · , θ N ). By equations (2.1), (2.2), (2.3) and parametrization of permeability field, we formulate the following stochastic two-phase flow system: Find random fields p(x, Θ, t), v(x, Θ, t), S(x, Θ, t) : D × I N × (0, T ] −→ R such that they almost surely (a.s) satisfy the following equations subject to initial and boundary conditions, (2.5) To further simply the notation, we will suppress the spatial variable x and temporal variable t in the rest of paper when no ambiguity occurs.

Sparse grid stochastic collocation method.
For stochastic two-phase flow systems (2.5), the statistical properties (e.g., mean and variance) of solutions are of interest. These properties may be obtained by first sampling the parameter random space using, for example, a Monte Carlo method or sparse grid collocation method, then solving the deterministic problems for the samples and analyzing the corresponding results to obtain the desired statistical quantities. The convergence of Monte Carlo methods is slow and a very large number of samples may be required, which leads to high computational cost. Instead, we use the Smolyak sparse-grid collocation method [40], where the Smolyak interpolant A(N + ℓ, N ) (ℓ ≥ 1) is a linear combination of tensor product interpolants with the 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 [6].
In the notation A(N + ℓ, N ), the ℓ represents the interpolation level. The sparse-grid collocation method is known to have the same asymptotic accuracy as tensor product collocation method, while requiring many fewer collocation points as the parameter dimension increases [36].
Sparse grids have been successfully applied to stochastic collocation in many recent works (e.g., [32,33,26]). Based on Smolyak formula (ref. [6]), a set of collocation points {Θ (j) } Nc j=1 in I N are specially chosen, where N c is the number of collocation points. With these chosen collocation points and corresponding weights {w (j) } Nc j=1 , the statistical properties of the solutions can be obtained. At each of the collocation points, the deterministic system (2.5) is solved and the output, for example, S(x, Θ (j) , t) is obtained. Then the mean of S(x, Θ, t) can be estimated by Here the weights {w (j) } Nc j=1 are determined by the basis functions of A(N + ℓ, N ) and joint probability function ρ(Θ) (ref. [18]). Similarly, the variance of S(x, Θ (j) , t) can be obtained by (2.6) Let H(N + ℓ, N ) denote the number of collocation points for Smolyak sparse grid interpolation A(N + ℓ, N ). Then it follows that (see [6]) This implies that the number of collocation points algebraically increases with respect to the dimension N . We utilize Smolyak sparse grid collocation for the numerical computation. The stochastic approximation of the Smolyak sparse grid collocation method depends on the number of collocation points and the dimension N of the random parameter space. The convergence analysis in [36] implies that the convergence of Smolyak sparse grid collocation is exponential with respect to the number of Smolyak points, but depends on the parameter dimension N . This exponential convergence rate behaves algebraically for N ≫ 1.
3. High dimensional model representation. The truncated KLE leads to the family of two-phase flow deterministic models (2.5) with a high-dimensional parameter Θ. The most challenging part of solving such a high-dimensional stochastic system is to discretize the high-dimensional stochastic space. There exist a few methods for the discretization of the random space [4,6,5,13,14,16,19,29,33,35,36,42,48] (more references can be found therein). Among these methods, the sparsegrid collocation method has been widely used and generates completely decoupled systems, each of which is the same size as the deterministic system. This method is usually very efficient in moderate-dimensional spaces. However, when the dimension of the random parameter is large, a large number of collocation points are required (2.7), and the deterministic model (2.5) must be solved at each of these collocation points. Consequently, the efficiency of this collocation method will deteriorate in a high-dimensional space. To overcome this difficulty, we use a high-dimensional model representation (HDMR) to reduce the stochastic dimension and enhance the efficiency of the simulation. By truncating HDMR, the high-dimensional model can be decomposed into a set of low-dimensional models. Hence the computational effort can be significantly reduced [32,47].
In the following section, we present a general HDMR framework and propose a hybrid HDMR method.
3.1. A general HDMR framework. In this section, we adopt the decomposition of multivariate functions described in [28] to present a general HDMR framework in terms of operator theory. Let F be a linear space of real-valued functions f (Θ) defined on a cube I N and Θ := (θ 1 , θ 2 , · · · , θ N ) ∈ I N . In this discussion f (Θ) represents a relationship between the random vector model input, Θ, and a model output (e.g., saturation, water-cut).
To introduce HDMR approach, we define a set of projection operators as follows.
be a set of commuting projection operators on F satisfying the following property: Let u ⊆ {1, · · · , N } be a subset and I the identity operator, we define P u := Π j∈u P j , and P ∅ := I.
Associated with projection P u , we define P u := P {1,··· ,N }\u . Then P u f only depends on the variables with indices in u.
Let µ be a measure on Borel subsets of I N and a product measure with unit mass, i.e., The inner product ·, · on F induced by the measure µ is defined as follows: The norm · F on F is defined by f F := f, f 1/2 . Given the measure µ, we can specify the projection operator in (3.1). We assume that the functions in F are integrable with respect to µ. For any f ∈ F and j ∈ {1, · · · , N }, we define Consequently, P u f (Θ) = I N −|u| f (Θ)Π j / ∈u dµ j (θ j ) for any u ⊆ {1, · · · , N }. By the definitions of these projection operator, we can obtain a decomposition of F as following (see [37]), For any u ⊆ {1, · · · , N }, we define Then we can show the operators M u are commutative projection operators and mutually orthogonal, i.e., The equation (3.4) gives an abstract HDMR expansion of f by The equation (3.4) also implies that Any set of commutative projectors {M v } generate a distributive lattice whose elements are obtained by all possible combinations (addition and multiplication) of the projectors in the set. The lattice has a unique maximal projection operator M, which gives the algebraically best approximation to the functions in F [21]. The range of the maximal projection operator M for the lattice is the union of the ranges of {M v }.
Because the commutative projectors {M v } are mutually orthogonal here, the maximal projection operator M and the range F M of M have explicit expressions as follows: As more orthogonal projectors are retained in the set, the resulting approximation by the maximal projection operator M will become better. If the measure µ in (3.2) is taken to be the probability measure then the resulting HDMR defined in (3.5) is called ANOVA-HDMR [38] . Let us fix a point Θ := (θ 1 ,θ 2 , · · · ,θ N ). If the measure µ in (3.2) is taken as the Dirac measure located at the point Θ, i.e., then the resulting HDMR is Cut-HDMR. The point Θ is so called cut-point or anchor point. In ANOVA-HDMR, the P v f involves N − |v| dimensional integration and the computation of the components of f is expensive. In Cut-HDMR, the P v f only involves the evaluation in a N − |v| dimensional space and the computation is cheap and straightforward. Because of this reason, we use Cut-HDMR for the analysis and computation in the paper.

A hybrid HDMR.
In this subsection, we use a less abstract representation of the HDMR to motivate approximations suitable for practical computation. Specifically, we review the traditional approach to truncated HDMR, and then propose and analyze an alternative hybrid HDMR.
Expanding elements of (3.5), the HDMR of f (Θ) can be written in the form, Here f 0 is the zeroth-order component denoting the mean effect of f (Θ). The firstorder component f i (θ i ) represents the individual contribution of the input θ i and the second-order component f ij (θ i , θ j ) represents the cooperative effects of θ i and θ j and so on. We define q-th order truncated HDMR by For most realistic physical systems, low-order HDMR M q f (e.g., q ≤ 3) may give a good approximation [38]. The M q f can approximate f (Θ) through truncating the HDMR. The M q f consists of a large number of component terms for high-dimensional models, the computation of all the components may be costly. For many cases, some high-order cooperative effects cannot be neglected in the model's output, especially when a model strongly relies on a few dependent variables. Consequently, the traditional truncated HDMR in (3.7) may not yield a very good approximation. To eliminate or alleviate these drawbacks of traditional truncated HDMR, we propose a new truncated HDMR, which we refer to as hybrid HDMR. We use the Fourier amplitude sensitivity test [9] to select the most active dimensions from the components of the high-dimensional random vector Θ. Let {f j } N j=1 be the first-order components defined in Eq.(3.6) and let σ 2 (f j ) denote the variance of f j . We may assume that σ 2 (f 1 ) ≥ σ 2 (f 2 ) ≥ · · · ≥ σ 2 (f N ) as we can re-order the index set {j} N j=1 to satisfy this requirement. Alternatively, we can order the index set to produce monotonically decreasing expectations, and use this ordering to select the most active dimensions. In this paper, we only consider the variance sensitivity, and we calculate σ 2 (f i ) (i = 1, · · · N ) using the method described in (2.6). Because the first-order components {f i } N i=1 are defined in one-dimensional random parameter spaces, the computational cost for will be reused when we calculate the variance of outputs represented by hybrid HDMR. We set a threshold constant with 0 < ζ < 1 and find an optimal J such that Then we define {θ j } J j=1 to be the most J active dimensions. We note that σ 2 (f j ) provides information on the impact of θ j when it is acting alone on the output. It is clear that if a change of the random input θ j within its range leads to a significant change in the output, then σ 2 (f j ) is large. Therefore, (3.8) provides a reasonable criteria for identifying the most active dimensions, and for a given stochastic model, we can adjust the value of ζ in Eq. (3.8) such that J is much less than N . Defining the index set J = {1, 2, · · · , J}, the proposed hybrid HDMR may be written, Here the set of projectors lattice and its maximal projection operator is M J . By (3.9), the hybrid HDMR M J f consists of two parts: the first part v⊆J M v f is the complete HDMR on the most active dimensions indexed in J , the second part i∈{1,··· ,N }\J M i f is the firstorder truncated HDMR on the remaining dimensions {1, · · · , N } \ J . We note that the component f u of f defined in (3.5) can be recursively computed and explicitly computed, respectively by (ref. [28]) By (3.10), we can rewrite equation (3.9) by We note that the operator P J projects the N -variable function f to a function defined on the J the most active dimensions. The term P J f gives the dominant contribution to M J f . Since any first-order components are usually important, these components are retained in M J f . In Cut-HDMR, there is no error for the approximation M J f of f (Θ) whenever the point Θ is located the J-th dimensional subvolume {Θ J ,θ J+1 , · · · ,θ N } across the cut-point Θ.
In this paper, identification of the J most important dimensions is based on a global sensitivity analysis on the univariate terms of HDMR. This criteria has been shown to be reasonable for many cases and is often used [9,32,47,22]. However, it is important to note that this criteria may not effectively identify the most important dimensions in some cases. For example, in some situations the individual contribution of an input parameter may not be significant, even though its cooperative influence with other inputs is significant. Similarly, the identification of the most active dimensions for highly variable solutions may require different criteria. Recent work [41] presents a preliminary discussion on the choice of this critera. Since, the new hybrid HDMR naturally includes the cooperative effects of higher-order terms within the J most important dimensions, this selection criteria is key to addressing these challenging cases and is an important topic for furture work.
If we use traditional truncated HDMR to approximate the term P J f in (3.11), then we can obtain the adaptive HDMR developed in [32,47], To simplify the notation, we will suppress q in M ad J ,q in the paper when the truncation order q is not emphasized. The following proposition gives the relationship among f , (3.11) and (3.12), respectively. Then Proof. We note the equality Then the equation (3.13) follows immediately by combining equation (3.14) and (3.15).
The equation (3.13) implies that hybrid HDMR has better approximation properties than adaptive HDMR.
The mean of f can be computed by summing the mean of all HDMR components of f for both ANOVA-HDMR and Cut-HDMR. The variance of f is the sum of the variance of HDMR components of f for ANOVA-HDMR. However, the direct summation of variance of components of f may not equal to variance of f for Cut-HDMR. Consequently, we may want to have a truncated ANOVA-HDMR to approximate the variance of f for Cut-HDMR. If we directly derive a truncated ANOVA-HDMR, the computation is very costly and more expensive than the direct computation of the variance of f itself. This is because a truncated ANOVA-HDMR requires computing many high-dimensional integrations. To overcome the difficulty, we can use an efficient two-step approach to derive a hybrid ANOVA-HDMR through the hybrid Cut-HDMR.
Let Θ J be the J-dim variable with indices in J . Using the hybrid HDMR formulation (3.11), the hybrid Cut-HDMR f cut (θ) of f has the following form We use M J to act on f cut (Θ) to get a hybrid ANOVA-HDMR f anova (θ) of f , which has the following form Then we have the following theorem. Theorem 3.3. Let f cut (Θ) and f anova (Θ) be defined in (3.16) and (3.17), respectively. Then The proof of Theorem 3.3 is presented in Appendix A. Theorem 3.3 implies that we can compute the mean and variance of f cut by directly summing the means and variances of the components of f cut . This observation can significantly reduce the complexity of this computation. In addition, the derived hybrid ANOVA-HDMR f anova is easily obtained by using the two-step approach through the hybrid Cut-HDMR f cut . Our further calculation shows that the traditional truncated Cut-HDMR (3.7) and adaptive Cut-HDMR (3.12) do not have these properties. We use sparse-grid quadrature [18] to compute the mean and variance in the paper.
Compared with the traditional truncated HDMR and adaptive HDMR, the hybrid HDMR defined in Eq.(3.11) has the following advantages: (a) The hybrid HDMR consists of significantly fewer terms than the traditional truncated HDMR and adaptive HDMR, the total computational effort of the hybrid HDMR can be substantially reduced. (b) Since the hybrid HDMR inherently includes all the cooperative contributions from the most active dimensions, approximation accuracy is not worse (maybe better) in hybrid HDMR than the traditional truncated HDMR and adaptive HDMR; (c) According to Theorem 3.3, the computation of variance of hybrid HDMR is much more efficient.

Analysis of computational complexity.
In the previous subsection, we have addressed the accuracy of the hybrid HDMR and developed some comparisons with traditional truncated HDMR and adaptive HDMR. In this subsection, we discuss the computational efficiency for the various HDMR techniques when Smolyak sparsegrid collocation is used. We remind the reader that the HDMR refers to Cut-HDMR.
Let C(f, ℓ) be the number of sparse grid collocation points with level ℓ in full random parameter dimension space and C(M q f, ℓ) the total number of sparse grid collocation points with level ℓ in the traditional truncated HDMR M q f . We define C(M J f, ℓ) and C(M ad J ,q f, ℓ) in a similar way. We first consider the case ℓ = 2 and q = 2 and calculate the total number of sparse grid collocation points for the different approaches. By using (2.7), we have for N ≫ 1 This relationship tells us that adaptive HDMR M ad J ,q f may be the most efficient if higher-level collocation is required and J is large as well.
For a high-dimensional stochastic model, if the number J of the most active dimensions is large and high-level (e.g., level 3 and above) sparse-grid collocation is required to approximate the term P J f in (3.11), using the adaptive HDMR (3.12) may improve efficiency according to the bounds given in (3.20). Adaptive HDMR has been extensively considered in many recent papers [22,32,33,47]. However, in our experience, the number of the most active dimensions J is often less than N 2 as N is large. If P J f is smooth with respect to Θ J , low-level sparse grid collocation (e.g., level 2) will generally provide an accurate approximation. Many problems fall into this class, and we focus on them in this paper.

3.4.
Integrating HDMR and sparse-grid collocation. As stated in Subsection 2.2, the sparse grid stochastic collocation method can reduce to the stochastic two-phase flow system (2.5) into a set of deterministic two-phase flow systems. However, it suffers from curse of dimensionality with increasing stochastic dimension. Integrating HDMR and a sparse-grid collocation methods provides an approach to overcome this difficulty and may significantly enhance the efficiency.
Without loss of generality, we consider the saturation solution S(Θ) as an example to present the technique for the hybrid Cut-HDMR. Using (3.11), (3.10) and Smolyak sparse-grid interpolation, we have whereŜ(Θ J ) = P J S(Θ),Ŝ(θ i ) = P i S(Θ) and S 0 = P ∅ S(Θ) = S(Θ). Here Θ is the cut-point for a Cut-HDMR. The choice of the cut-point may affect the accuracy of the truncated Cut-HDMR approximation. The study in [17] argued that an optimal choice of the cut-point is the center point of a sparse-grid quadrature. In this paper, we will use such a cut-point for computation. Due to Theorem 3.3, the mean and variance of S can be approximated by  Smolyak sparse-grid interpolations. The computation of these interpolations is usually very expensive when N and J are large, and hence, computing the variance using adaptive HDMR is usually much more expensive than using hybrid HDMR. The numerical experiments in Section 5 demonstrate this performance advantage for hybrid HDMR.
4. Mixed multiscale finite element method. In Section 3.4, we have discussed integrating hybrid HDMR and sparse grid collocation method to reduce the computation complexity from high-dimensional stochastic spaces. The permeability field is often heterogeneous in porous media. It is necessary to use a numerical method to capture the heterogeneity. MMsFEM is one of such numerical methods and has been widely used in simulating flows in heterogeneous porous media [2,26]. To simulate the two-phase flow system (2.5), it is necessary to retain local conservation for velocity (or flux). To this end, we use MMsFEM to solve the flow equation and obtain locally conservative velocity. Using MMsFEM coarsens the multiscale model in spatial space and can significantly enhance the computation efficiency.
Let V h ⊂ H 0 (div, D) and Q h ⊂ L 2 (D)/R be the finite element spaces for velocity and pressure, respectively. Then the numerical mixed formulation of (4.2) is to find We use MMsFEM for (4.3). It means that mixed finite element approximation is performed on coarse grid, where the multiscale basis functions are defined. In MMsFEM, we use piecewise constant basis functions on coarse grid for pressure. For the velocity, we define multiscale velocity basis functions. The degree of freedom of multiscale velocity basis function is defined on interface of coarse grid. Let e K i be a generic edge or face of the coarse block K. The velocity multiscale basis equation associated with e K i is defined by For local mixed MsFEM [7], b K i = n, the normal vector. If the media demonstrate strong non-local features including channels, fracture and shale barriers, some limited global information is needed to define the boundary condition b K i to improve accuracy of approximation [2,25]. We will specify the boundary condition b K i for different parts in (4.1). Then ψ K i = −k∇w K i defines multiscale velocity basis function associated to e K i , and the multiscale finite dimensional space for velocity is defined by It is well-known that using approximate single-phase global velocity information can considerably improve accuracy for multiscale simulation of two-phase flows [1,2,25,26]. For the two-phase flow system (2.5), the single-phase global velocity u sg (Θ) solves the following equation  By using hybrid HDMR and sparse grid interpolation, the u sg (Θ) admits the following approximation whereû sg (θ j ) = P j u sg (Θ) andû sg (θ i ) = P i u sg (Θ). Here the interpolation levels ℓ and ℓ ′ can be different. Because {θ j } j∈J are the most active dimensions, it is usually desirable that ℓ ≥ ℓ ′ . Now we are ready to specifically describe the multiscale finite element space forû(Θ J ),û(θ i ) (i ∈ {1, . . . , N } \ J ) and u 0 .
• To construct the multiscale basis functions forû(Θ J ), we takek =k(Θ J ) in (4.4) and The multiscale finite element space forû(Θ J ) is defined by • To construct the multiscale basis functions forû(θ i ) (i ∈ {1, . . . , N } \ J ), we takek =k(θ i ) in (4.4) and The multiscale finite element space forû(θ i ) is defined by • To construct the multiscale basis functions for u 0 , we replacek by k(Θ) in (4.4) and b K i = u sg,0 | e K i . The multiscale finite element space for u 0 is defined by For different parts of the hybrid HDMR expansion M J u(Θ) defined in (4.1), we use different boundary conditions for multiscale basis equations. This increases the hierarchies of multiscale basis functions, which are in tune with the sensitivity of random parameter dimensions.
By (3.23), to compute the moments of outputs we only need to construct the multiscale finite element space V h (Θ J ) at the set of collocation point {Θ The convergence analysis for the mixed MsFEM using approximate global information can be found in [25,26]. From the above discussion, here the approximate global information is completely determined by stochastic single-phase velocity contributed from individual random dimensions θ i (i = 1, . . . , N ). According to the analysis in [25], the resonance error between the coarse mesh size h and the effect of individual contributions of θ i (i = 1, . . . , N ) is negligible. This error is a major source using local MMsFEM [26]. However, the resonance error between the coarse mesh size h and the effect of cooperative contributions of Θ is retained, but is usually small. This is our motivation to use limited global information determined in 1-dimensional random space. If the random permeability field has low variance and does not exhibit global features, we can use local mixed MsFEM for the stochastic simulation. Define For velocity, we actually compute u h (x, Θ), namely, the numerical solutionsû h (x, Θ J ) andû h (x, θ i ) are evaluated at their collocation points. We define u(x, Θ) in the same way. Now we analyze the total error between u(x, Θ) and u h (x, Θ). By triangle inequality, where E hdmr represents the error from hybrid HDMR, E col the error by sparse grid collocation and E ms the error of MMsFEM. The error E hdmr depends on the choice of the most active dimensions J , the error E col depends on the number of collocation points and random dimensions, and the error E ms relies on the coarse mesh size h and the approximation of the boundary conditions (4.6) and (4.8) toû sg (Θ J ) and u sg (θ i ), respectively. According to the works in [2,25], if the boundary conditions in (4.6) and (4.8) are replaced byû sg (Θ J ) andû sg (θ i ), respectively, then the error E ms only depends on coarse mesh size h. The error E col has been extensively discussed in many recent literatures, e.g., [5,6,42], and we focus on the errors E hdmr and E ms in the paper.

Numerical results.
In this section, we present numerical results for twophase flow in random porous media. The hybrid HDMR and MMsFEM are incorporated to enhance simulation efficiency. For random porous media with non-local spatial features and high variance, we use limited global information to capture the multiscale velocity information and achieve better accuracy. We combine MMsFEM with the two different truncated HDMR techniques techniques (hybrid HDMR M J and adaptive HDMR M ad J ) for the flow simulations. The accuracy and efficiency of the different methods will be compared.
For numerical simulation, we assume that k(x, Θ) is a logarithmic random field, i.e., k(x, Θ) := exp (a(x, Θ)). This assumption assures the positivity of k(x, Θ) and the well-posedness of the flow equation in (2.5). Here a(x, ω) is a stochastic field and can be characterized by its covariance function cov[a]: For the numerical experiments, we use a two point exponential covariance function for the stochastic field a, where (x i , y i ) (i = 1, 2) is the spatial coordinate in 2D, σ 2 is the variance of stochastic field a, and l x and l y denote the correlation length in the x− and y−direction, respectively. We will specify these parameters to define the covariance function in numerical examples. Using Karhunen-Loève expansion (KLE) the random field a(x, Θ) admits the following decomposition where the random vector Θ := (θ 1 , θ 2 , · · · , θ N ) ∈ R N and the random variables {θ i (ω)} N i=1 are mutually orthogonal and have zero mean and unit variance. The effects of permeability field a(x, Θ) with uniform, beta and Gaussian distributions on the mean and standard deviation of output were discussed in [30], where the numerical results showed similar peak values of standard deviation for the three different distributions. Therefore, in this section we assume for numerical convenience that each θ i is i.i.d. and uniform on [−1, 1]. Although, this assumption leads to potentially less physical fields, it generates variability suitable for evaluating the performance of the proposed methods.
When MMsFEM is used, the fine grid is coarsened to form a uniform coarse grid. We solve the pressure equation on the coarse grid using MMsFEM and then reconstruct the fine-scale velocity field as a superposition of the multiscale basis functions. The reconstructed velocity field is used to solve the saturation equation with a finite volume method on the fine grid. We solve the two-phase flow system (2.5) using the classical IMPES (Implicit Pressure Explicit Saturation). The temporal discretization is an implicit scheme, which is unconditionally stable but leads to a nonlinear system (Newton-Raphson iteration solves the nonlinear system). In all numerical simulations, mixed multiscale basis functions are constructed once at the beginning of computation. In the discussion, we refer to the grid where multiscale basis functions are constructed as the coarse grid. Limited global information is computed on the fine grid.
We compare the saturation fields and water-cut data as a function of pore volume injected (PVI). The water-cut is defined as the fraction of water in the produced fluid and is given by q w /q t , where q t = q o + q w , with q o and q w being the flow rates of oil and water at the production edge of the model. In particular, q w = ∂D out f (S)u ·nds, q t = ∂Ω out u · nds, where ∂D out is the out-flow boundary. Pore volume injected is defined as PVI = 1 Vp t 0 q t (τ )dτ , with V p being the total pore volume of the system, and provides the dimensionless time for the displacement. We consider a traditional quarter five-spot problem (e.g., [1]) on a square domain D, where the water is injected at left-top corner and oil is produced at the right-lower corner of the rectangular domain (in the 2D examples).

Linear transport.
In this subsection, we consider the case when the transport (saturation) equation is linear in the flow system (2.5). For this purpose, we consider the mobility of water and mobility of oil defined as following, Consequently, the total mobility λ(S) = 1 and the fractional flow of water f w (S) = S. Then the two-phase flow system (2.5) reduces to a linear flow model. Since λ(S) = 1, the velocity is not updated to compute saturation. We use the example to investigate the performance of the hybrid HDMR technique and MMsFEM for linear models.
The stochastic permeability field k(x, Θ) is given by k(x, Θ) = exp a(x, Θ) , where a(x, Θ) is characterized by the covariance function cov[a] defined in (5.1), whose parameters are defined as follows: σ 2 = 1, l x = l y = 0.1. We note that the deceay rate of the eigenvalues in the KLE (5.2) expansion depends inversely on the correlation length and is often used as a guideline to determine the number of terms, N [43]. Here the decay rate is relatively slow and we use 80 terms to represent the random field a(x, Θ) with a relative error of less than 5%. This implies that k(x, Θ) is defined in an 80-dimensional random parameter space, i.e., N = 80. Fig. 5.1 depicts a realization of the random field a(x, Θ). The stochastic field k(x, Θ) is defined in a 60 × 60 fine grid. We choose 6 × 6 coarse grid and apply MMsFEM to compute velocity. The time step is taken to be 0.02 PVI for discretizing temporal variable.
To demonstrate the efficacy of the proposed method to simulate the stochastic linear two-phase flow system, we compare the performance of several methods: finescale mixed finite element method (MFEM) with full random dimensions, MMsFEM on full random dimensions, MMsFEM based on adaptive HDMR, and MMsFEM based on hybrid HDMR. The solution computed by fine-scale MFEM with the full stochastic space is referred to as the reference solution in this paper. To assess the performance of the various MMsFEM and HDMR combinations, we compute the mean and standard deviation (std) for the quantifies of interest from the two-phase flow model, such as saturation and water-cut. Since the random field in the example does not have strong non-local features, the local MMsFEM (L-MMsFEM) generally gives an accurate approximation [25], and it is used in this example. To choose the most active random dimensions for adaptive HDMR and hybrid HDMR, we choose a threshold constant ζ = 0.9 and use criteria (3.8). This approach produces 31 most active dimensions. The Smolyak sparse-grid collocation method with level 2 is used to tackle the stochastic space. We evaluate the flow model's outputs at the collocation points and use the associated weights to compute the mean and standard deviation of the outputs. Table 5.1 lists the number of deterministic models to be solved for the above four different methods. From the table, we see that MFEM on full random dimensions needs to compute the largest number of deterministic models (fine-scale models), MMsFEM based on hybrid HDMR needs to compute the smallest number of deterministic models (coarse-scale models). This means that MMsFEM based on hybrid HDMR significantly reduces the computational complexity.   Fig. 5.4 shows the mean of water-cut curves, which are all nearly identical. This demonstrates that a good approximation has been achieved. Fig. 5.5 shows the standard deviation of water-cut curves from the four different methods. Here, we see that both L-MMsFEM based on hybrid HDMR and L-MMsFEM based on adaptive HDMR produce a very good approximation for the standard deviation at most time instances. The standard deviation of the water-cut curve from L-MMsFEM based on hybrid HDMR is almost identical to the standard deviation of water-cut curve from L-MMsFEM based on adaptive HDMR.
Next we discuss the relative errors of the domain integrated mean and standard deviation. Let S r /W r , S ms,f /W ms,f , S ms,aH /W ms,aH and S ms,hH /W ms,hH be the saturation/water-cut using MFEM on full random dimensions, MMsFEM on full random dimensions, MMsFEM based on adaptive HDMR and MMsFEM based on hybrid HDMR, respectively. Then we define the relative errors of saturation mean and saturation standard deviation between S r and S ms,f as following, where σ denotes the standard deviation operation. We can similarly define the relative errors E m aH (S) and E std aH (S) between S r and S ms,aH , and the relative errors E m hH (S) and E std hH (S) between S r and S ms,hH . We list these errors on saturation at 0.4 PVI in Table 5.2. From the table, we see that L-MMsFEM based on hybrid HDMR gives the same accuracy of mean as L-MMsFEM based on adaptive HDMR. They are both very close to the error of mean by L-MMsFEM on full random dimensions. The relative error of standard deviation by L-MMsFEM based on hybrid HDMR also has a good agreement with that of L-MMsFEM based on adaptive HDMR. The table implies that the errors only from adaptive HDMR and hybrid HDMR are very small. The main source of errors is from spatial multiscale approximation. relative error of mean relative error of std. L-MMsFEM on full random dim.
1.223540e-002 5.302182e-002 L-MMsFEM based on adaptive HDMR 1.361910e-002 9.389291e-002 L-MMsFEM based on hybrid HDMR 1.361910e-002 9.482789e-002 For water-cut, the relative errors of mean and standard deviation between W r and W ms,f are defined by 4) where [0, 1] is the temporal domain. The relative errors of water-cut, such as E m aH (W ), E std aH (W ), E m hH (W ) and E std hH (W ), are defined in a similar way. We list these relative errors in Table 5.3. The table shows similar behavior for these methods as was observed in the calculation of saturation mean and saturation standard deviation. The performance of L-MMsFEM based on hybrid HDMR is almost the same as the performance of L-MMsFEM based on adaptive HDMR. The number of the most active dimensions may have an important impact on both the computational efficiency and solution accuracy. To this end, we investigate the performance for different threshold constants in (3.8). Specifically, we take ζ = 0.75, 0.8, 0.85, 0.9, 0.95. For these threshold constants, we find that the corresponding number of most active dimensions are dim(Θ J ) = 19, 22, 26, 31, 41.  HDMR increases dramatically. This comparison demonstrates that hybrid HDMR is much more efficient than adaptive HDMR. There are two reasons that adaptive HDMR requires more CPU time: (1) adaptive HDMR needs to solve more deterministic models; (2) to compute the standard deviation (or variance), a large number of stochastic interpolations are involved. and standard deviation (right) for saturation at 0.4 PVI. From the figure, we conclude that increasing the number of the most active dimensions can substantially improve the accuracy of the mean and standard deviation for saturation. We also see that for the saturation approximation the L-MMsFEM based on hybrid HDMR has almost the same accuracy as the L-MMsFEM based on adaptive HDMR. Fig. 5.8 presents the relative errors of water-cut mean (left) and water-cut standard deviation (right) for the different number of the most active dimensions. From the figure, we observe that the reduction in the error with increasing dim(Θ J ) is smaller for water-cut than for mean saturation. Here, the error from L-MMsFEM is dominating the total error, and hence, is overwhelming the impact of improved accuracy in the HDMR itself.

Non-linear transport.
In this subsection, we consider the two-phase flow system (2.5) with the mobilities of water and oil defined by nonlinear functions of saturation, Here µ w /µ o = 0.1, the ratio between viscosity of water and oil. Consequently, the fractional flow function f w (S) of water is given by This results in a non-linear two-phase flow system.
We again consider the random permeability k(x, Θ) = exp a(x, Θ) . Here the covariance function cov[a] of a(x, Θ) is defined in (5.1) with l x = l y = 0.2 and σ 2 = 1. The mean of a(x, Θ) is highly heterogeneous and its map is depicted in Fig. 5.9 (left). It is actually obtained by extracting and rescaling an SPE 10 [8] permeability field (the 45-th layer). We truncate the KLE (5.2) after the first 30 terms to represent the random field a(x, Θ) and so the random field k(x, Θ) is defined in a 30-dimensional random parameter space. Fig. 5.9 shows a realization of the random field a(x, Θ) (right). The stochastic field k(x, Θ) is defined on a 60 × 60 fine grid. We choose a 6 × 6 coarse grid for MMsFEM to compute the velocity. To discretize the temporal variable of the saturation equation, the time step is taken to be 0.01 PVI.  results in this setting, limited global information can improve the accuracy [1,2,25]. In this subsection, we incorporate the hybrid HDMR technique and the adaptive HDMR technique with both local MMsFEM (L-MMsFEM) and global MMsFEM (G-MMsFEM), and compare the performance of these methods. To obtain the most active random dimensions for the truncated HDMR techniques, we choose a threshold constant ζ = 0.9 and use criteria (3.8). This gives rise to 14 most active dimensions in the 30-dim random parameter space. We use the collocation points and weights associated with level 2 Smolyak sparse-grid collocation to compute the mean and standard deviation. truncated HDMR techniques provide good approximations. Fig. 5.11 describes the point-wise standard deviation of saturation at 0.4 PVI for those different methods. Based on the figure, we make two observations: (1) the largest variance occurs around the water front; (2) the hybrid HDMR technique gives almost the same standard deviation map as the adaptive HDMR technique. We note that these observations are consistent with those of Section 5.1. We also compute the water-cut curves for the different methods. To minimize overlap in the visualization, we only present the water-cut curves for four methods: MFEM on full random dim., G-MMsFEM on full random dim., G-MMsFEM based on adaptive HDMR and G-MMsFEM based on hybrid HDMR. Fig. 5.12 shows the mean of the water-cut curves for the four methods. Here, the four curves overlap each other at almost all times. We note that there exists a small fluctuation right after the water breaks through. The reason may be that the value of water-cut changes sharply right after water break-through time. The standard deviation of the water-cut curves for the four methods are illustrated in Fig. 5.13. From the figure, we find that the variance of water-cut rises rapidly at the break-through time, and has a small number of peaks that are likely related to the dominant channelized flow paths in the underlying mean permeability field.
In order to carefully measure the differences caused by L-MMsFEM and G-MMsFEM when integrated with adaptive HDMR and hybrid HDMR, we follow the procedure outlined in Subsection 5.1 and compute the relative errors between the reference solution and the solutions of the various coarse models (MMsFEM with (or without) HDMR techniques). As before, the reference solution is solved by fine-scale MFEM on full random dimensions. The relative errors of saturation and water-cut are defined similarly to (5.3) and (5.4), respectively. Table 5.4 and 5.5 list the relative errors of mean and standard deviation for saturation at 0.4 PVI and water-cut, respectively. Examining this table, we note: (1) limited global information can enhance the accuracy of MMsFEM; (2) the results by adaptive and hybrid HDMR are very close to each other.
Finally we examine the efficiency for the different approaches. To this end, we record the CPU time for the seven different methods: fine-scale MFEM on full random dimensions, L-MMsFEM on full random dimensions, G-MMsFEM on full random  Fig. 5.14 shows the CPU time for the seven different approaches. Here it is apparent that the approaches based on hybrid HDMR need the least CPU time and achieve the best efficiency. The CPU time of hybrid HDMR is only a fraction, about a third, of the CPU time of adaptive HDMR. Moreover, this speedup was shown earlier to increase quickly with the number of most active dimensions. Comparing the performance of the local and global MMsFEM based methods, we find that for each HDMR technique, the CPU time of G-MMsFEMs is only slightly larger than the CPU time of L-MMsFEM. For adaptive HDMR, a significant amount of CPU time is spent in computing variance. However, computation of the variance in hybrid HDMR is straightforward and requires very little CPU time. This is an important advantage of the hybrid HDMR approach.     6. Conclusions. In this paper, we presented a general framework for highdimensional model representations (HDMRs), and proposed a hybrid HDMR technique in combination with a mixed multiscale finite element method (MMsFEM) to simulate two-phase flow through heterogeneous porous media. The hybrid HDMR technique decomposes a high-dimensional stochastic model into a moderate-dimensional stochastic model and a few one-dimensional stochastic models. An optimization criteria was developed that ensures a specified percentage of the total variance is present in the most active dimensions (i.e., the moderate-dimensional stochastic model). In addition, we demonstrated that combining an MMsFEM with hybrid HDMR could significantly reduce the original model's complexity in both the resolution of the physical space and the high-dimensional stochastic space. We also presented MMsFEM based on adaptive HDMR, which has been widely used in stochastic model reduction. Compared with adaptive HDMR, the hybrid HDMR is much more efficient and retains the same (or better) accuracy. To capture strong non-local features in the multiscale models, we have incorporated important global information into the multiscale computation. This can improve the approximation accuracy of the proposed coarse multiscale models. We carefully analyzed the proposed MMsFEM using HDMR techniques and discussed both the computational efficiency and approximation errors. The MMsFEM based on HDMR techniques was applied to two-phase flows in heterogeneous random porous media. Both linear and non-linear saturation dependencies of the mobility were considered. The simulation results confirmed the performance of the proposed approaches.