Flow and transport are the most fundamental phenomena in subsurface porous media associated with various physical processes, e.g., oil and gas flow in petroleum reservoir Terry et al. (2015), CO sequestration Zhang et al. (2015), water pollution dispersion Bear and Cheng (2010), etc. The numerical simulation and analyses of flow and transport in subsurface porous media are highly demanded in practical engineering and mechanism studies. However, the simulation results are always subject to the influence of uncertainties, mainly stemming from the inherent spatial heterogeneity of media properties caused by complex geological processes Boschan and Noetinger (2012). It has been widely recognized that in natural subsurface porous media, most properties, such as permeability, porosity, etc., exhibit an uneven spatial distribution. For example, the hydraulic conductivity can span several orders of magnitude in an aquifer or reservoir. How to quantificationally identify the influence of uncertainties of porous media properties on the flow and transport behaviors in subsurface physical processes has been a research hot spot in recent years.
Therefore, the uncertainty quantification is an essential task in the simulation of practical subsurface flows where porous media properties that unknown or partially known are taken as the input parameters. A possible way to deal with uncertainties of subsurface porous media is to treat porous media properties as random fields, then perform the stochastic simulation on the subsurface flow governing equations with random coefficients to evaluate the quality of interest (QoI). Among commonly-used stochastic simulation methods, e.g. Monte Carlo(MC) method, stochastic finite element method, stochastic collocation method, the MC method demonstrates apparent advantages such as it is a non-intrusive approach that only the realization of coefficients is needed while the original model code remains unchanged, and it is more easily to be implemented. In the standard MC method, the computer-generated (pseudo) random points are used, and in many cases, the computational efficiency is always unsatisfied for large-scale problems. The Quasi-Monte Carlo (QMC) method improve the demerit of MC by using deterministic quasi-random points. These points exhibit lower discrepancy and distribute more uniformly in the probability space. Moreover, to reduce the sample variance and further improve the computational efficiency, the multilevel Monte Carlo (MLMC) method was proposed and developed by HeinrichHeinrich (2001) and Giles Giles (2008). It applies the control variates technique that a series of discretization is adopted with increasing resolution and computes the QoI on each of them, the success of which lies in the effective variance reduction sequentially.
It should be mentioned that in the particular case of subsurface flow with random coefficients, the problem is further aggravated where very detailed geological models are needed (a large number of cells) for an accurate description of the flow. To further alleviate the computational burden connected to the evaluation of random parameter effects on subsurface flow using the MLMC method, in this study, we exploit the similar hierarchies of MLMC and multigrid methods and proposed a full multigrid multilevel (quasi-) Monte Carlo (FMG-MLQMC) approach. In this proposed method, the solution on coarse mesh can be obtained as a byproduct of the full multigrid solution on fine mesh on each mesh level , instead of directly solving the equations on the coarse mesh as the standard MLMC does. The proposed FMG-MLQMC method saves the computation of the . There have been works coupling the multigrid solver with the multilevel framework, see Kumar et al. (2017); Robbe et al. (2018) for example. However, the FMG-MLQMC method we proposed saves the computational cost without modifying the MLMC framework. We exploit the implementation method for upscaling the random coefficient from fine mesh to neighboring coarse mesh. Although in this study we only focus on the simple single-phase subsurface flow with random coefficients, the proposed approach can be applied and extended naturally to multiphase flow and transfer in porous media and any other flow and transport problems associated with uncertainty effect.
The rest of the paper is organized as follows: In Section 2, we give a brief description of the single-phase subsurface flow and then introduce the proposed full multigrid multilevel (quasi-) Monte Carlo method in detail. The methodology on the upscaling method of random coefficients from the fine mesh to the coarse mesh, which preserves the random structure, is presented as well. In Section 3, we verify the effectiveness (a smaller estimator variance and faster convergence rate) of the presented method by comparing with standard MLMC method in two numerical experiments. Finally, in Section 4, we report the concluding remarks of this work along with a brief discussion of future directions.
2.1 Model problem and MLMC method
In this work, we consider the following elliptic problem,
where is the random, spatial-varying coefficient, is the computational domain, is a sample from the probability triple . and are Dirichlet and Neumann boundaries respectively. In single-phase flow context, Eq.1 corresponds the steady-state situation, when and prescribe the pressure and velocity of the fluid at the boundary, then the solution depicts the pressure in the domain .
In this work, we address the random elliptic problem using the multi-level algorithm. Basically, the MLMC method employs a series of control variates, which are often the discretized models with increasing resolution levels. Here, we associate each level with one mesh with given resolution. The approximations of quantity of interest(QoI) on these levels are denoted as , see Figure 1.
We are interested in the approximation on the finest level . The MLMC method not only computes the solution on level itself, but also calculates the solutions on all the preceding meshes. The expectation of such quantity can be expressed by the following telescoping formula,
We can approximate each expectation using the Monte Carlo approach as follows,
Here, we associate level with the -th term in the telescoping formula (2). Notice that on each level , we use the same samples to calculate and . Then and are likely to correlate well, and the variance of will be small, see Eq.4.
where we denote as the variance of the on level .
Also notice that on different levels, independent samples are used so that the variance of the multilevel estimator is the summation of the variance on each level.
where we let . If we write and , then
be an unbiased estimator for,
then the multilevel estimator becomes,
2.2 MLMC Complexity Theory
Let denote a quantity of interest, and denote the corresponding numerical approximation on -th mesh. If we assume that the weak error and the level variance decreases exponentially while the cost per sample on each level increases exponentially, there exist positive constants , and satisfying the following,
where is the cost per sample on level . With the mean square error less than a threshold,
the total computational cost satisfies,
2.3 MLMC Algorithm
This subsection gives a MLMC algorithm initially proposed by M.GilesGiles (2008).
In Algorithm 1, the variances are approximated by the sample variances on the run. The weak error is approximated by Richardson extrapolation .
And here we consider a equal split of the estimator variance and the approximation error, i.e.,
However, it is possible to determine the split factor in an optimal wayCollier et al. (2015).
The optimal number of samples can be obtained by solving a constrained optimization : minimizing the total computational cost subject to the constraint Eq.10.
MLMC algorithm will work under the following three conditions.
The sequence converges. Otherwise, the telescoping equation (2) does not yield a converging result.
and are estimated using the same underlying random sample in equation (3), and are thus well correlated. In this case, the estimator variance is significantly reduced.
The telescoping sum (2) introduces no bias error. Notice that in the telescoping equation, the term for appears twice. However, the two may be evaluated differently. If we denote the -th term in the telescoping equation and by and , respectively, which denote the fine mesh solution and coarse mesh solution on level , then the condition needs to be satisfied in order to introduce no bias error in equation (2). The expectation of fine mesh solution on level should be the same as that of the coarse mesh solution on level .
2.4 MLQMC Algorithm
The QMC method can solve integration problems as well. In contrast to the MC method, the QMC method replaces random points with deterministic points. Figure 2 gives an example of Monte Carlo points, lattice rule, and Sobol’ sequence.
The QMC approximation is given by,
Notice that the points are deterministic. However, deterministic points yield biased estimates. In this work, we use the randomized QMC: random shift and digital scramble for lattice rule and Sobol’ sequence respectively. The interested readers may refer to Dick et al. (2013) for the above mentioned two randomization techniques.
The multigrid method was originally introduced to solve elliptic boundary-value problems efficiently. It has since been developed to solve either linear or non-linear systems. Multigrid methods compute the solution on a sequence of grids. Figure. 3 gives an illustration of the full multigrid scheme.
We observe that, when the full multigrid solver is applied to the MLMC problem, based on the same level hierarchies, the solution on the coarse mesh can be obtained as a byproduct of the multigrid solution on fine mesh on each level . Thus, in our proposed FMG-MLMC method, we have saved the computation for . Also notice that at the red circles in Figure. 3 the solutions are exact, since they are the end point of each V-cycle.
Recall that and are evaluated by the same underlying random sample. In level , we denote and as the coefficients of the fine and coarse models, respectively. Also recall the consistency condition, since we use the same numerical solver for and , we only require that and follow the same distribution.
can be generated by the matrix decomposition method, KL-expansion method or other random field generation methods, see Liu et al. (2019) for example. A way to generate is to coarsen . In order to prevent bias error, should satisfy the same distribution law as . Figure 4 shows a way of coarsening.
In this scheme, the value of the coefficient in the coarse grid is selected to be the corresponding block in the fine grid. We denote the blocks in fine grid and coarse grid by and respectively, then we have the following,
Here we give an algorithm to describe our proposed FMG-MLMC method.
Notice that in our scheme, no correlation is introduced into each levels. Thus the equation (5) still holds true.
3 Numerical Validation
3.1 Problem Statement
Recall the elliptic problem 1 in Section 2, now we focus the physical domain . In this work, we consider cases in two different boundary conditions and quantities of interest, as listed in Table 1. Case I is of the Dirichlet boundary type, with pointwise output quantity. Case II is of the mixed Dirichlet-Neumann boundary condition, whose output is the outflow at the east boundary.
|1||, , ,|
|2||, , ,|
The governing equation (1) is discretized by the finite-volume method on rectangular grids. On level
the degree of freedom is.
(2) Random Fields
We choose the Matérn covariance function
where is the Euclidean distance of two points, is the correlation length, and controls the smoothness of the field.
The parameters of Matérn covariance for four different random fields are given in Table 2.
The random fields are generated using the KL-expansion method. The truncation term is determined when of the variability is captured, meaning that
where the summation of all eigenvalues satisfies the following
Readers of interest can see Ernst et al. (2009) for examples. is the random field region.
(3) QMC Method
In QMC method, the Lattice rule points are generated using the software from Kuo and Nuyens (2016), Sobol’ matrices from Joe and Kuo (2008)
are used to generate Sobol’ sequences. In both cases, 24 randomizations are applied. The confidence intervals are obtained by 10 sets of randomization.
3.2 Numerical results
In this subsection we present the numerical results of the two cases (Tab. 1) with four random field settings (Tab. 2). In each figure, the first and second row corresponds the cases and respectively, while the first column and second column corresponds the cases and . We will first give the results of the first case in the following.
The simulation starts by estimating the asymptotic rates in the assumptions (7). Figs. 5, 6 and 7 plot the variance of the QoI and against level . By comparison, the QMC method reduces in the variance not in the asymptotic variance convergence rate, but in the y-axis offsets.
For the results of MLQMC-Lattice, the variance test (Fig. 6) is presented.
For the results of MLQMC-Sobol’, the variance (Fig. 7) is presented.
The QMC method does not affect the expectation value nor the computational cost. Here we skip the comparison of and , but present the final computational cost of the three methods, given . The results in 4 random fields are presented in Fig. 8.
For case II, we test the QMC convergence rate and then present the computational complexity.
Next the results of MLQMC-Lattice, the convergence test (Fig. 9) variance of level estimator against in each level. The offset between the lines reveals the variance decreases with the levels. The comparison shows Sobol’ sequence’s advantage in decreasing variance. The random parameter settings have small impact on the variances in this case.
Next, the results of MLQMC-Sobol’, the convergence test (Fig. 10)
Finally, we compare the computational cost of the three methods, given . The results in 4 random fields are presented in Fig. 11.
In this work, we combined the MLMC with a full multigrid method. We saved the computation for the coarse grid solution on each level without modifying MLMC hierarchy and introducing correlation between different levels. We applied the consistent coarsening approach such that no bias error was introduced in the telescoping formula (2).
We tested our FMG-MLQMC algorithm on 2-D elliptic PDE with random coefficients for two different types of boundary condition settings and QoIs. The random coefficients were modelled as lognormal random fields with the Matérn covariance function with various parameter settings. We observed that quasi-Monte Carlo approaches have better performance on smoother random fields and problems with more regularity. Also, the comparison of Monte Carlo and quasi-Monte Carlo methods (including Lattice rule and Sobol’ sequence) revealed that QMC outperforms MC due to a smaller estimator variance, and Sobol’ sequence performs slightly better than Lattice rule.
One future work could be the substitution of the geometric multigrid solver with the algebraic multigrid (AMG) solver. In the AMG scheme, the grids are not associated with physical meshes, rather, the grids are fully determined by the matrix entries algebraically.
Another work could be to extend the elliptic model to more sophisticated models, such as two-phase porous flow. The efficient sampling and fast simulation of multi-phase subsurface flow under heterogeneous media could produce practical values.
Further work could extend the multilevel model to a multiscale model, and multiscale meshes rather than geometric meshes would be used. In this case, one level could correspond to one scale, and sampling would be performed on each scale. The literature on multiscale modeling can be found in Jenny et al. (2003); Efendiev and Hou (2009).
The authors gratefully acknowledge the support from the National Natural Science Foundation of China (Nos. 51874262, 51904031) and the Research Funding from King Abdullah University of Science and Technology (KAUST) through the grants BAS/1/1351-01-01.
- Modeling groundwater flow and contaminant transport. Vol. 23, Springer Science & Business Media. Cited by: §1.
- Scale dependence of effective hydraulic conductivity distributions in 3d heterogeneous media: a numerical study. Transport in porous media 94 (1), pp. 101–121. Cited by: §1.
- A continuation multilevel monte carlo algorithm. BIT Numerical Mathematics 55 (2), pp. 399–432. External Links: Cited by: §2.3.
- High-dimensional integration: the quasi-monte carlo way. Acta Numerica 22, pp. 133–288. Cited by: §2.4.
- Multiscale finite element methods: theory and applications. Vol. 4, Springer Science & Business Media. Cited by: §4.
- Efficient solvers for a linear stochastic galerkin mixed formulation of diffusion problems with random data. SIAM Journal on Scientific Computing 31 (2), pp. 1424–1447. Cited by: §3.1.
- Combining sparse grids, multilevel mc and qmc for elliptic pdes with random coefficients. In International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, pp. 265–281. Cited by: §2.4.
- Multilevel quasi-monte carlo path simulation. Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics 8, pp. 165–181. Cited by: §2.4.
- Multilevel monte carlo path simulation. Operations Research 56 (3), pp. 607–617. Cited by: §1, §2.3.
- Multilevel monte carlo methods. In LSSC, Cited by: §1.
- Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. Journal of Computational Physics 187 (1), pp. 47–67. Cited by: §4.
- Constructing sobol sequences with better two-dimensional projections. SIAM Journal on Scientific Computing 30 (5), pp. 2635–2654. Cited by: §3.1.
- A multigrid multilevel monte carlo method using high-order finite-volume scheme for lognormal diffusion problems. International Journal for Uncertainty Quantification 7 (1). Cited by: §1.
- Multilevel quasi-monte carlo methods for lognormal diffusion problems. Mathematics of Computation 86 (308), pp. 2827–2860. Cited by: §2.4.
- Application of quasi-monte carlo methods to elliptic pdes with random diffusion coefficients: a survey of analysis and implementation. Foundations of Computational Mathematics 16 (6), pp. 1631–1696. Cited by: §3.1.
- Advances in gaussian random field generation: a review. Computational Geosciences. External Links: Cited by: §2.5.
- Recycling samples in the multigrid multilevel (quasi-) monte carlo method. arXiv preprint arXiv:1806.05619. Cited by: §1.
- Applied petroleum reservoir engineering. Pearson Education. Cited by: §1.
- Sequentially coupled thmc model for co2 geological sequestration into a 2d heterogeneous saline aquifer. Journal of Natural Gas Science and Engineering 27, pp. 579–615. Cited by: §1.