Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Accurate and efficient uncertainty quantification of subsurface fluid flow via the probabilistic collocation method
(USC Thesis Other)
Accurate and efficient uncertainty quantification of subsurface fluid flow via the probabilistic collocation method
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
ACCURATE AND EFFICIENT UNCERTAINTY QUANTIFICATION OF SUBSURFACE FLUID FLOW VIA THE PROBABILISTIC COLLOCATION METHOD by Heng Li A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) December 2010 Copyright 2010 Heng Li Dedication This dissertation is dedicated to my parents, for their unconditional love and support. ii Acknowledgements I would like to express my deepest gratitude to my advisor, Professor Dongxiao Zhang. I feelfortunate to havehimas myadvisor, who led meon the journeypursuing myPh.D. I have gained invaluable knowledge and experience from his guidance and inspiration that will enlighten my career. Iwouldliketoexpressmyappreciationtothemembersofmydissertationcommittee, ProfessorIrajErshaghiandProfessorRogerGhanem,fortheirtime,expertise,andhelpful comments. IwouldliketothankSchlumbergerandChevronforprovidingmeinternshipopportu- nitieswhereIgainedindustrialexperiences. InparticularatChevron,Ihadopportunities workingonaprojectcloselyrelatedtomyresearchtopic. Iwouldliketothankmymentor Pallav Sarma at Chevron for his help. I would like to thank the professors and staff in the petroleum engineering program at USC, for helping me throughout the academic pursuit. I would also like to thank my fellow graduate students and friends, who have made the journey enjoyable. I am deeply grateful to my parents and brother for their constant support and en- couragement. Last but not least, my special thanks go to my dear wife, for her selfless love and devotion. iii Table of Contents Dedication ii Acknowledgements iii List Of Tables vi List Of Figures vii Abstract xii Chapter 1: Introduction 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Objective of Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Chapter 2: Probabilistic Collocation Method 11 2.1 Polynomial Chaos Expansion . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Probabilistic Collocation Method . . . . . . . . . . . . . . . . . . . . . . . 13 Chapter 3: Gaussian Random Field 15 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Karhunen-Loeve Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Single-Phase Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.1 KL Based Probabilistic Collocation Method . . . . . . . . . . . . . 20 3.3.1.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.1.2 Selection of Collocation Points . . . . . . . . . . . . . . . 22 3.3.1.3 Post-Processing . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.2 Results of PCM and MC . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.3 Comparisons with Other Stochastic Methods . . . . . . . . . . . . 30 3.3.3.1 Polynomial Chaos Expansion Approach . . . . . . . . . . 30 3.3.3.2 KLME Approach. . . . . . . . . . . . . . . . . . . . . . . 35 3.3.4 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3.4.1 Effect of Correlation Length η . . . . . . . . . . . . . . . 37 3.3.4.2 Effect of Large Spatial Variability σ 2 Y . . . . . . . . . . . 39 3.3.4.3 Probability Density Function . . . . . . . . . . . . . . . . 40 iv 3.3.4.4 Illustrative Examples in 2D . . . . . . . . . . . . . . . . . 41 3.3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4 Multi-phase Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4.1 Mathematical Formulations . . . . . . . . . . . . . . . . . . . . . . 49 3.4.2 Case Study: 1D Water Flooding . . . . . . . . . . . . . . . . . . . 52 3.4.2.1 Random Porosity Field . . . . . . . . . . . . . . . . . . . 53 3.4.2.2 Random Permeability Field . . . . . . . . . . . . . . . . . 62 3.4.2.3 Relative Permeability Type . . . . . . . . . . . . . . . . . 71 3.4.3 Case Study: 3D Black Oil Model . . . . . . . . . . . . . . . . . . . 74 3.4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Chapter 4: Non-Gaussian Random Field 86 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.2 Investigation of Approaches for Simulating Non-Gaussian Random Fields 87 4.2.1 Memoryless Transformation Approach . . . . . . . . . . . . . . . . 87 4.2.2 Kernel Principal Component Analysis . . . . . . . . . . . . . . . . 90 4.2.3 Karhunen-Loeve Expansion . . . . . . . . . . . . . . . . . . . . . . 96 4.2.4 Independent Component Analysis . . . . . . . . . . . . . . . . . . 98 4.2.5 Polynomial Chaos Expansion Approach . . . . . . . . . . . . . . . 100 4.3 Modeling Subsurface Flow in Non-Gaussian Random Fields . . . . . . . . 102 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Chapter 5: Conditional Simulation of Random Fields 116 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.2 Conditional Simulation of Gaussian Random Fields . . . . . . . . . . . . . 117 5.3 Conditional Simulation of Non-Gaussian Random Fields . . . . . . . . . . 124 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Chapter 6: Arbitrary Random Variables 133 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.2 Generalized Polynomial Chaos Expansion . . . . . . . . . . . . . . . . . . 135 6.3 Orthogonal Polynomials for Arbitrary Random Variables. . . . . . . . . . 136 6.4 PCM Compared with Experimental Design . . . . . . . . . . . . . . . . . 138 6.4.1 Case Study: Model 1 . . . . . . . . . . . . . . . . . . . . . . . . . 138 6.4.2 Case Study: Model 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 142 6.4.3 Case Study: Model 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 Chapter 7: Conclusions and Recommendations 158 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 7.2 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 References 162 v List Of Tables 6.1 The description of the random parameters for model 1 . . . . . . . . . . . 139 6.2 The description of the random parameters for model 2 . . . . . . . . . . . 146 6.3 Probabilities of the categorical variables for case 2 of model 2 . . . . . . . 149 6.4 The description of random parameters for model 3 . . . . . . . . . . . . . 154 vi List Of Figures 3.1 Series of eigenvalues and their finite sums, for η =4 and η =2 . . . . . . 26 3.2 ComparisonsofthemeanandvarianceofhydraulicheadderivedfromMC, and PCM up to the third order . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 The mean and variance of hydraulic head derived from 4 sets of Monte Carlo simulations, corresponding to 28, 84, 1000, and 10,000 realizations . 29 3.4 Comparisons of head variance derived from PCM and MC, with different spatial variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.5 Comparisons of head variance derived from PCE and MC, with different spatial variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.6 Comparisons of head variance derived from KLME and MC, with different spatial variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.7 Comparisons of head variance derived from PCM, KLME and MC, with different spatial variability for η =2 . . . . . . . . . . . . . . . . . . . . . 42 3.8 Head variance derived from MC and PCM (up to the fourth order) for huge spatial variability, i.e. σ 2 Y =4.0 . . . . . . . . . . . . . . . . . . . . 43 3.9 Comparisons of the probability density function (PDF) of hydraulic head at three positions (i.e., x = 2, 4, and 6) obtained with PCM, KLME, and MC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.10 Comparison of head variance at x 2 =5 in 2D obtained with PCM and MC 45 3.11 Comparisons of mean and variance of hydraulic head along the diagonal line in 2D with a pumping well obtained with PCM and MC . . . . . . . 46 3.12 Series of eigenvalues and their finite sums . . . . . . . . . . . . . . . . . . 55 vii 3.13 Mean and standard deviation of water saturation and pressure obtained from 2nd order PCM and the MC (with 1000 realizations) for the case with standard deviation of log porosity i.e. σ lnϕ = 0.1 and viscosity ratio m=0.5: (a) and (b) for water saturation; (c) and (d) for pressure . . . . 57 3.14 Mean and standard deviation of water saturation and pressure obtained from 2nd order PCM and MC (with 1000 realizations) for the case with standarddeviationoflogporosityi.e. σ lnϕ =0.1andviscosityratiom=2: (a) and (b) for water saturation; (c) and (d) for pressure . . . . . . . . . . 58 3.15 Mean and standard deviation of water saturation and pressure obtained from 2nd order PCM and MC (with 1000 realizations) for the case with standarddeviationoflogporosityσ lnϕ =0.1,correlationlengthη/L=1/5, and viscosity ratio m=2: (a) and (b) for water saturation; (c) and (d) for pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.16 Probability density functions of water saturation S w and pressure at two different locations (x = 150 ft and 350 ft), obtained from 2nd order PCM and MC (with 1000 realizations) for the case with standard deviation of log porosity σ lnϕ = 0.1, correlation length η/L = 2/5, and viscosity ratio m=2: (a) and (c) for water saturation; (b) and (d) for pressure . . . . . 61 3.17 Mean and standard deviation of water saturation and pressure obtained from PCM (2nd and 4th order) and MC (with 1000 realizations) for σ 2 Y = 0.25: (a) and (b) for water saturation; (c) and (d) for pressure . . . . . . 64 3.18 Mean and standard deviation of water saturation and pressure obtained fromPCM(2ndand4thorder)andMC(with1000and10,000realizations) for σ 2 Y =1.0: (a) and (b) for water saturation; (c) and (d) for pressure . . 65 3.19 Mean and standard deviation of water saturation and pressure obtained fromPCM(2ndand4thorder)andMC(with1000and10,000realizations) for σ 2 Y =4.0: (a) and (b) for water saturation; (c) and (d) for pressure . . 66 3.20 Mean and standard deviation of water saturation and pressure obtained from 40 sets of MC (with 210 realizations) for σ 2 Y = 4.0: (a) and (b) for water saturation; (c) and (d) for pressure . . . . . . . . . . . . . . . . . . 69 3.21 Meanandstandarddeviationofwatersaturationandpressure(atx=275 ft) obtained from Monte Carlo simulations versus the number of realiza- tions for σ 2 Y =4.0: (a) and (b) for water saturation; (c) and (d) for pressure 70 3.22 Confidence intervals of the total well oil production at different time for σ 2 Y =1.0, obtained from the 4 th order PCM and MC (with 10,000 realiza- tions) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 viii 3.23 Confidence intervals of the total well oil production at different time for σ 2 Y =4.0, obtained from the 4 th order PCM and MC (with 10,000 realiza- tions) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.24 Mean and standard deviation of water saturation and pressure obtained from 2nd order PCM and MC (with 1000 realizations) for the case with standard deviation of log porosity i.e. σ lnϕ = 0.1, viscosity ratio m = 2 and new relative permeability functions: (a) and (b) for water saturation; (c) and (d) for pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.25 Mean and standard deviation of water saturation and pressure obtained from PCM (2nd and 4th order) and MC (with 1000 and 10,000 realiza- tions) for the case with σ 2 Y = 1.0, viscosity ratio m = 1, and new relative permeability functions: (a) and (b) for water saturation; (c) and (d) for pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.26 Initial oil saturation in the reservoir . . . . . . . . . . . . . . . . . . . . . 79 3.27 Series of eigenvalues in 3 D . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.28 One realization of the horizontal permeability field for σ 2 Y =1.0 . . . . . . 80 3.29 Contour plot of mean and standard deviation of oil saturation obtained from 2nd order PCM and MC (with 1000 realizations): (a) and (b) for PCM; (c) and (d) for MC . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.30 Contour plot of mean and standard deviation of gas saturation obtained from 2nd order PCM and MC (with 1000 realizations): (a) and (b) for PCM; (c) and (d) for MC . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.31 Confidence intervals of reservoir responses for σ 2 Y = 1.0: (a) field cumula- tive oil production (in stock tank barrel STB) normalized by pore volume PV(inreservoirbarrelRB);(b)fieldcumulativegasproduction(inMSCF) over PV (in RB). Blue for PCM, Red for MC. . . . . . . . . . . . . . . . 83 4.1 The PDF of the non-Gaussian random field Y(x) . . . . . . . . . . . . . 93 4.2 The PDF of Y(x) generated with Linear PCA . . . . . . . . . . . . . . . 94 4.3 (a) Eigenvalues for polynomial kernel of order d = 3; (b) finite sums of the eigenvalues for polynomial kernel of order d = 3 . . . . . . . . . . . . . . 95 4.4 Reconstruction of realizations for Gaussian kernel with σ =5 , N =50. . 96 ix 4.5 The marginal PDF of K(x): gamma distribution γ(2,1) with its approxi- mation by the 5th order Hermite polynomial chaos expansion. . . . . . . 104 4.6 Realizations of the random field K(x). . . . . . . . . . . . . . . . . . . . . 105 4.7 Variance of hydraulic head for the case of gamma distribution γ(2,1), ob- tained from PCM and MC for ICA. . . . . . . . . . . . . . . . . . . . . . 107 4.8 (a) PDF of Beta distribution ; (b) head variance for the case of Beta distribution β(4,2), obtained from PCM and MC for ICA. . . . . . . . . . 108 4.9 (a) PDF of Beta distribution ; (b) head variance for the case of Beta distribution β(5,1), obtained from PCM and MC for ICA. . . . . . . . . . 109 4.10 Variance of hydraulic head for the case of gamma distribution γ(2,1), ob- tained with the PCM and MC for PCE, . . . . . . . . . . . . . . . . . . . 110 4.11 (a) PDF of gamma distribution γ(3,2) ; (b) head variance for the case of gamma distribution γ(3,2), obtained from PCM and MC for PCE. . . . . 112 4.12 (a) PDF of distribution (i); (b) head variance from PCM and MC for distribution (i). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.13 (a) PDF of distribution (ii); (b) head variance from PCM and MC for the case of distribution (ii). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.1 Confidence interval of Y(x) . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.2 Variance of hydraulic head, obtained from the PCM and MC . . . . . . . 121 5.3 (a) Confidence interval of Y(x); (b) Variance of hydraulic head, obtained from the PCM and MC for the case of σ 2 Y =2.0. . . . . . . . . . . . . . . 122 5.4 (a) Confidence interval of Y(x); (b) Variance of hydraulic head, obtained from the PCM and MC for the case of σ 2 Y =4.0. . . . . . . . . . . . . . . 123 5.5 The marginal PDF of the hydraulic conductivity K(x) . . . . . . . . . . . 126 5.6 The confidence interval of the hydraulic conductivity K(x) . . . . . . . . 127 5.7 Variance of the hydraulic head after conditioning . . . . . . . . . . . . . . 127 5.8 (a) The mean of the hydraulic conductivity K(x) after conditioning; (b) The standard deviation of K(x) after conditioning. . . . . . . . . . . . . . 129 x 5.9 (a)Themeanofthehydraulicheadafterconditioning,obtainedwithPCM (with 210 simulations); (b) The standard deviation of the hydraulic head after conditioning, obtained with PCM (with 210 simulations). . . . . . . 130 5.10 (a) The mean of the hydraulic head after conditioning, obtained with MC (with 10,000 simulations); (b) The standard deviation of the hydraulic head after conditioning, obtained with MC (with 10,000 simulations). . . 131 6.1 PDFs of cumulative oil production, obtained from PCM and MC for case 1 140 6.2 PDFs of cumulative oil production for case 1 of model 1, obtained from different experimental designs . . . . . . . . . . . . . . . . . . . . . . . . . 141 6.3 The PDF of the multi-modal distribution . . . . . . . . . . . . . . . . . . 143 6.4 The PDF of the multi-modal distribution from sampling . . . . . . . . . . 143 6.5 PDF of cumulative oil production for the case of multi-modal distribution, obtained from PCM and MC . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.6 Comparison of PDF of cumulative oil production for the case of multi- modal distribution, obtained from experimental designs, PCM and MC. For experimental designs: (a) liner; (b) pure quadratic; (c) full quadratic polynomial response surface . . . . . . . . . . . . . . . . . . . . . . . . . . 145 6.7 The configuration of model 2 . . . . . . . . . . . . . . . . . . . . . . . . . 146 6.8 PDFs of cumulative oil production, obtained from PCM and MC . . . . . 147 6.9 PDFs of cumulative oil production, obtained from different experimental designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 6.10 Probability function of transmissibility multiplier 2 and 4. . . . . . . . . . 151 6.11 PDFs of cumulative oil production, obtained from PCM and MC . . . . . 151 6.12 PDFs of cumulative oil production, obtained from experimental designs (with 25 simulations) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 6.13 The configuration of model 3 . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.14 PDFs of cumulative oil production obtained with different experimental designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 6.15 PDFs of cumulative oil production obtained with PCM . . . . . . . . . . 156 xi Abstract Uncertainty quantification of subsurface flow has recently attracted a significant amount of attention. The uncertainty can result from the combination of the formation hetero- geneityandtheincompleteknowledgeofitsproperties. Traditionalflowsimulationstreat the geological formation deterministic, thus resulting in deterministic predictions. How- ever,takinguncertaintyintoaccountnecessitatesastochasticdescriptionoftheformation properties and hence stochastic approaches to flow simulations. This dissertation explores an efficient approach, i.e., the probabilistic collocation method (PCM) for uncertainty quantification of flow in random porous media. In this approach, the dependent random variables are represented by employing the orthogonal polynomial functions (polynomial chaos expansions) as the bases of the random space. Utilizing the collocation technique in the random space directly results in a set of inde- pendent simulations. This independence feature of this stochastic approach allows us to directly employ existing flow simulators. Random fields such as heterogeneous perme- ability (or porosity) fields are parameterized using some dimension reduction techniques such as the Karhunen-Loeve expansion. Applications to both single-phase and multi- phase flows are performed. Other than the typical log-normal probability distribution, other non-Gaussian probability distributions are also considered for characterizing the xii formation properties that are treated as random fields. Conditional simulations of both Gaussian and non-Gaussian random fields when measurements of hydraulic conductivity are available are discussed in this dissertation. This dissertation also considers uncertainty analysis of petroleum reservoir simula- tions, where the uncertain parameters have arbitrary probability distributions. The ex- perimental design methods are widely used for uncertainty analysis in the oil industry. However, the traditional experimental design methods usually have an inherent assump- tionofuniformdistributionsfortherandominputsandtheydonottakeintoaccountthe fullprobabilitydensityfunctions(PDFs)oftheinputrandomparametersconsistently. In this dissertation, orthogonal polynomials for arbitrary PDFs are constructed numerically and the PCM is utilized for uncertainty propagation. Each collocation point captures a different weight of the random variables for a given PDF and thus yields optimal ap- proximation of the responses. Various studied cases reveal that, while the computational efforts are greatly reduced compared to Monte Carlo simulation, the PCM is able to accurately quantify uncertainty of the reservoir performance. Results also reveal that the PCM is more robust, accurate, and efficient than experimental design methods for uncertainty analysis. xiii Chapter 1 Introduction 1.1 Introduction Subsurface fluids such as the groundwater and petroleum hydrocarbon are stored in the interstitial pores of the underground rocks. The geological formations are invisible from the ground surface and the information about the formations is usually acquired through drilling wells through them. The limited measurements make information about the formations inadequate. Further, the geological formations normally exhibit spatial heterogeneity at different scales. Consequently, it is well recognized that there exists significant uncertainty in characterizing the subsurface. The dynamic behaviors of the subsurface fluid flow and transport in porous media are modeled by certain mathematical equations. The model parameters that are used to represent the formation properties are usually considered to be random variables or random fields due to the uncertainty and thus the equations governing the flow and transportbecomestochastic. Hence,uncertaintyquantificationisinvolvedwithstochastic modeling of the fluid flow and transport in porous media. 1 Extensive studies on stochastic modeling of flow and transport in porous media have been studied and many stochastic approaches have been developed in the discipline of hydrogeology in the last two decades [e.g., Dagan, 1989, Gelhar, 1993, Cushman, 1997, Neuman, 1997, Zhang, 2002, Rubin, 2003]. Some approaches have been applied to the discipline of petroleum engineering [e.g., Zhang and Tchelepi, 1999, Zhang et al., 2000]. In most of the studies, the media properties (e.g. the permeability) are assumed to follow log-normal distributions, and most stochastic approaches are developed on the basis of Gaussian random field. The assumption of Gaussian distribution and its broad applications are due to its simplicity and the central limit theorem. A Gaussian random field can be determined merely by its first two statistical moments, e.g. the mean and covariance. In the context of geostatistics, the two-point geostatistics based on the as- sumption of Gaussian random field has gained much popularity in many fields, such as hydrogeology,miningandpetroleumengineering. Thevariogramthatisequivalenttothe covarianceiscalculatedfrommeasurementsandusedtorepresentthespatialcorrelations of the random fields. Based on the assumption of Gaussian distribution, various methods are available to generate realizations of the random fields. However, the nature of the real worldis complex and the properties of the porous me- dia in the subsurface may exhibit other non-Gaussian distributions. The properties may even have discrete patterns and their distributions then behave as discrete distributions. Multiple-pointgeostatisticshasbeendevelopedinordertoaccountforcomplexstructures of the formation such as channels [Strebelle, 2000, Strebelle and Journel, 2001]. Training images of the geological patterns are required as prior knowledge, based on whichrealiza- tions of the geological structures are reconstructed. When these realizations are used as 2 input for modeling flow and transport in porous media, Monte Carlo method is usually considered as a common tool for uncertainty quantification, though it is computationally demanding [Feyen and Caers, 2006]. On the other hand, the computational cost is an important concern in stochastic modeling of flow and transport in porous media. In the literature, few studies have been conducted to efficiently model flow and transport with media properties characterized by non-Gaussian probability distributions. 1.2 Literature Review Stochastic modeling of flow in porous media involves solving the stochastic governing equationsdescribingtheflowbehaviors, whicharemainlystochastic(ordinaryorpartial) differential equations. The Monte Carlo method is the most common and conceptually straightforward technique to solve the stochastic differential equations numerically [e.g., Hassan et al., 1998, Ballio and Guadagnini, 2004]. A random field is represented by an ensemble of equally probable realizations, which can be fed as input into numerical simulators for modeling flow in porous media and quantifying its uncertainty. Monte Carlo method is based on statistical sampling of the random fields and resorts to the law of large numbers. Its accuracy depends on the number of realizations used [e.g., Ballio and Guadagnini, 2004]. Generally a large number of simulations are required to achieve statisticalconvergence,thusthemaindrawbackofMonteCarlomethodisitsrequirement of large computational efforts. Therefore, more efficient stochastic approaches are highly demanded. Despite its high computational cost, Monte Carlo method is usually used to 3 provide benchmark solutions for accessing the accuracy of other stochastic approaches owing to its easy implementation [e.g., Zhang, 2002, Zhang and Lu, 2004]. An alternative to Monte Carlo simulation is the method based on moment equations [e.g. Graham and McLaughlin, 1989, Neuman, 1993, Zhang and Neuman, 1995, Zhang, 1998, 1999, Guadagnini and Neuman, 1999]. It derives a system of deterministic differ- ential equations governing the statistical moments (usually the first two moments) of the random variables, with the method of perturbation or some type of closure approxima- tion. However, thecomputationalcostforthe(conventional)momentequationmethodis still high for large-scale problems. In order to compute the hydraulic head covariance to first-order in terms of the variance of log hydraulic conductivity, one needs to solve a set of deterministic equations on a grid of N nodes for about 2N times: N times for solving the cross-covariance between the hydraulic head and the log hydraulic conductivity, and about N times for the hydraulic head covariance [Zhang, 2002, Lu and Zhang, 2006]. The moment equation approach based on Karhunen-Loeve expansion (KLME) was proposed by Zhang and Lu [2004] for single phase flow in porous media and extended to unsaturated and multiphase flows [Yang et al., 2004, Chen and Zhang, 2005, Chen et al., 2006a], flow in nonstationary random hydraulic conductivity fields [Lu and Zhang, 2007], andsolutetransport[Liuetal.,2007]. Theloghydraulicconductivityisrepresentedwith theKarhunen-Loeve(KL)expansionandtheflowandtransportquantitiessuchasthehy- draulichead,velocity,andconcentrationaredecomposedbytheperturbationexpansions. The KLME leads to uncoupled equations governing the coefficients, which can be used to construct the required statistical moments. With the KLME approach, high order terms of the mean and variance of hydraulic head (or other flow and transport quantities) can 4 be obtained with relatively small computational efforts [Zhang and Lu, 2004] and can be easily implemented with existing simulators [Lu and Zhang, 2006, Liu et al., 2006, 2007]. The KLME approach is generally more efficient than the traditional moment equation approach. Although it is derived with polynomial perturbation expansions, the KLME is found to be applicable to large variability in log hydraulic conductivity for single phase flow. However, its validity range is not yet clear for multi-phase flow. Anotherapproachcalledthepolynomialchaosexpansion(PCE)method,pioneeredby Ghanem and Spanos [1991], has attracted considerable interest for uncertainty quantifi- cationinvariousproblemsofmechanics[GhanemandDham,1998,XiuandKarniadakis, 2002]. This technique includes representing the random variables using polynomial chaos basis and deriving appropriate discretized equations for the expansion coefficients using the Galerkin technique. The PCE method has been applied to uncertainty quantifica- tion of single-phase and multiphase flow in porous media [Ghanem, 1998, Ghanem and Dham, 1998]. The polynomial chaos expansion allows high order approximation of ran- dom variables and possesses the property of fast convergence under certain conditions. However, the deterministic coefficients of the polynomial chaos expansion are governed by a set of coupled equations, which are difficult to solve when the number of coefficients is large (owing to a large random dimensionality or a high polynomial order). Further, the simulators used for deterministic flow simulations cannot be directly employed since the original governing equations are coupled with coefficients of the polynomial chaos expansions. The probabilistic collocation method (PCM) introduced by Tatang et al. [1997] is another efficient stochastic approach. In this method, polynomial chaos expansions are 5 used to represent the random variables of interest. Instead of the Galerkin scheme, a col- locationtechniqueisusedtosolveforthecoefficientsofthepolynomialchaosexpansions, and independent deterministic equations can be derived. The PCM is usually used for parametric uncertainty quantification where the input random variables are independent [Webster et al., 1996, Tatang et al., 1997, Kumar et al., 2005]. However, in stochastic modeling of subsurface flow, random fields may be involved. Sarma et al. [2005] applied the PCM in optimization of the oil reservoir production where permeability is a random field. In the context of uncertainty quantification of flow in porous media, the media prop- erties are usually assumed to follow log-normal distributions. For example, the logarithm of hydraulic conductivity is assumed as a Gaussian random field. Thus most stochastic approaches are developed on the basis of Gaussian random fields. However, stochas- tic modeling of flow in non-Gaussian spatial random fields is engaged with numerical simulation of the non-Gaussian random fields. Simulation of non-Gaussian processes (fields) is mainly based on memoryless trans- form of the standard Gaussian processes. The processes are called translation processes and are restricted to stationary processes. Yamazaki and Shinozuka [1988] proposed an iterative procedure for simulating the translation processes, where the marginal proba- bility distributions and the spectral power density are specified. It is based on a spectral representation method where the target power spectrum function is matched through an iterative algorithm. Ren et al. [1995] proposed a method for conditional simulation of non-Gaussian fields, by utilizing the algorithm for unconditional simulation proposed by Yamazaki and Shinozuka [1988]. In their study, the covariance function instead of the 6 spectral power density of the random fields is specified, and the iterative procedure is thus modified. Grigoriu [1998] presented an algorithm for simulating the non-Gaussian translation processes, by using certain properties of the translation processes. The al- gorithm is based on an analytical relationship between the covariance function of the non-Gaussian process and that of the underlying Gaussian process. Sakamoto and Ghanem [2002] proposed a different approach for generating non- Gaussianandnon-stationaryrandomprocessesbyusingthepolynomialchaosexpansion. In their approach, the non-Gaussian processes are represented through polynomial chaos expansionsofGaussianrandomprocesseswhicharesimulatedthroughtheKLexpansion. Wan and Karniadakis [2009] studied solving elliptic problems with non-Gaussian ran- dom coefficients. They employed the KL expansion directly to represent the correlated non-Gaussianrandomcoefficientsandassumedindependenceamongtherandomvariables in the KL expansion, although they are only mutually uncorrelated. 1.3 Objective of Study Inthisstudy, theprobabilisticcollocationmethod(PCM)isexploredforefficientlyquan- tifying uncertainty of subsurface flow in porous media. The method is first applied in the cases where the formation properties are treated as Gaussian random fields. A Gaus- sian random field is known to have large dimension. An efficient technique, i.e. the Karhunen-Loeve expansion is employed for dimension reduction. Through integrating the Karhunen-Loeve expansion and polynomial chaos expansion, the PCM can be used 7 to efficiently quantify uncertainty of flow in porous media in the cases of Gaussian ran- dom fields. The application is started with the single-phase water flow for deriving the methodology. With the examples of the single-phase flow, the validity of this approach is investigatedforvariousdegreesofspatialvariabilityofthemediaanddifferentcorrelation lengths. This approach is compared with some stochastic approaches for efficiency and accuracy. The examples of multi-phase flow in petroleum reservoirs are also performed. The properties are also considered to be non-Gaussian random fields other than the typical log-normal distribution. The approaches for simulating non-Gaussian random fields are explored and their applicabilities are assessed for the purpose of efficient mod- eling of flow in porous media. When measurements of reservoir properties, such as the hydraulic conductivity, are prescribed, the conditional simulation of flow and transport in porous media is performed for uncertainty reduction. Conditional simulation of both Gaussian and non-Gaussian random fields is discussed in this study. When the uncertain parameters in reservoir simulations are treated as random vari- ables, and their probability distributions are arbitrary, the experimental design method that is widely used for uncertainty quantification may be inappropriate. In this study, an approach based on the PCM is proposed for quantifying uncertainty in the cases of arbitrary random variables. 8 1.4 Dissertation Outline The dissertation is organized as follows: Chapter 1 is the introduction, which describes the motivation of the dissertation and the literature review. Section 2.1 introduces the polynomial chaos expansion, and section 2.2 derives the probabilistic collocation method from the weighed residual method. Chapter 3 derives the KL-based PCM for stochastic modeling of flow associated with Gaussian random fields. Section 3.1 introduces the Karhunen-Loeve(KL) expansion for parameterization of Gaussian random fields. Examples of single-phase flow and multi- phase flow are performed in sections 3.2 and 3.3, respectively. The examples of single- phase flow are used to compare the PCM with other stochastic approaches, such as the polynomial chaos expansion approach, the Karhunen-Loeve expansion based moment equation approach, and the Monte Carlo method. Chapter 4 discusses various approaches for simulating non-Gaussian random fields and modeling flow in porous media characterized by non-Gaussian random fields. The approachessuchasthememorylesstransformationapproach, thekernelprincipalcompo- nent analysis, the Karhunen-Loeve expansion, the independent component analysis, and the polynomial chaos expansion, are described in section 4.2, and their applications to modeling flow in non-Gaussian random fields are presented in section 4.3. Chapter 5 discusses conditional simulation of flow and transport in porous media. Sections 5.1 and 5.2 discuss the conditional simulation of Gaussian and non-Gaussian random fields, respectively. 9 Chapter6 describesan approachfor dealing with arbitraryrandom variables. Section 6.2 introduces the generalized polynomial chaos expansion. Section 6.3 discusses how to generate orthogonal polynomials for arbitrary random variables. Section 6.4 presents the applications of PCM to oil reservoir problems, along with the comparison of PCM with the experimental design methods. Chapter 7 presents the conclusions and recommendations. 10 Chapter 2 Probabilistic Collocation Method 2.1 Polynomial Chaos Expansion Any random variable or random field can be expressed with the polynomial chaos expan- sion. For example, the random field y(x,θ), where x∈ D (physical space) and θ ∈ Θ (a probability space), is expressed as, y(x,θ)=a 0 (x)+ ∞ ∑ i 1 =1 a i 1 (x)Γ 1 (ξ i 1 (θ))+ ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 a i 1 i 2 (x)Γ 2 (ξ i 1 (θ),ξ i 2 (θ)) + ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 i 2 ∑ i 3 =1 a i 1 i 2 i 3 (x)Γ 3 (ξ i 1 (θ),ξ i 2 (θ),ξ i 3 (θ))+..., (2.1) where the coefficients a 0 (x) and a i 1 i 2 ,...,i d (x) are deterministic functions ofx, and Γ d (ξ i 1 , ...,ξ i d ) are orthogonal polynomial chaos of order d with respect to the random variables (ξ i 1 ,...,ξ i d ). For independent standard Gaussian random variables (ξ i 1 ,...,ξ i d ), Γ d (ξ i 1 ,...,ξ i d )=(−1) d e 1 2 T ∂ d ∂ξ i 1 ...∂ξ i d [e − 1 2 T ], (2.2) 11 where is a vector denoting (ξ i 1 ,...,ξ i d ) T . Hermite polynomials form the best orthogonal basis for Gaussian random variables [e.g., Ghanem and Spanos, 1991]. In case of other random distributions, generalized polynomial chaos expansions [Xiu and Karniadakis, 2002] can be used to represent the random field as in equation (2.1). In practice, equation (2.1) is usually truncated by finite terms, and it can be written simply as y(x,θ)= P ∑ j=1 c j (x)Ψ j (), (2.3) where is a vector of dimension N. There is a one-to-one correspondence between the terms in equations (2.1) and (2.3). The total number of terms P is determined by the random dimensionality N and the degree of the polynomial chaos expansion d, P = (N +d)! N!d! . (2.4) The Galerkin scheme is usually used to evaluate the deterministic coefficients in the polynomial chaos expansion approach [e.g., Ghanem, 1998, Xiu and Karniadakis, 2002]. Another approach is the stochastic collocation method (SCM) developed by Mathelin et al. [2005], which has been applied to problems with small number of random variables. If a stochastic process (random field) has to be described with a large number of random variables, the SCM with a tensor product implementation may suffer from the problem of high dimensionality. The probabilistic collocation method (PCM), also called the deterministicequivalentmodelingmethod(DEMM),isintroducedbyTatangetal.[1997] and successfully applied to the uncertainty analysis in some fields [Webster et al., 1996, 12 Kumar et al., 2005]. In their applications, the input parameters are usually independent. In our study, both random fields and independent variables are considered for the inputs. 2.2 Probabilistic Collocation Method Both the Galerkin scheme and the PCM used to determine the coefficients of the polyno- mialchaosexpansionscanbederivedwiththeweightedresidualmethod. Forastochastic differential equation Ly(x,θ)=f(x), (2.5) wherey(x,θ)istheunknownrandomspacefunctionandf(x)isthesourceterm. Theop- eratorLinvolvesdifferentiationsinspaceandcanbenonlinear. Ify(x,θ)isapproximated by the polynomial chaos expansion and the approximation is denoted as ˆ y(x,θ), ˆ y(x,θ)= P ∑ i=1 c i (x)Ψ i ((θ)). (2.6) Define the residual R as R({c i },)=Lˆ y−f. (2.7) The weighted residual method in the random space is expressed as ∫ R({c i },)w j ()p()d =0, (2.8) where w j () is the weighting function, j = 1,...,P, and p() is the joint probability density function of . 13 With the Galerkin scheme, the weighting function is chosen as the basis function of the approximation, w j ()=Ψ j (). (2.9) Substituting equation (2.9) into equation (2.8) and expressing it in the form of ensemble mean yields, ⟨R({c i },)Ψ j ()⟩=0, (2.10) which leads to a coupled system of deterministic equations. A particular example of the Galerkin PCE scheme is given in the following chapter. In the probabilistic collocation method, the weighting function is chosen as the Dirac delta function, w j ()=δ(− j ), (2.11) where j is a particular set selected with certain algorithm out of the random vector . The elements in j are called the collocation points. Then equation (2.8) becomes, R({c i }, j )=0, (2.12) which results in a set of independent equations, evaluated at the given sets of collocation points, j , where j = 1,2,··· ,P. It is seen that P sets of collocation points are needed andthusequation(2.12)hastobesolvedforP timesinordertoobtaintheP coefficients {c i }, where i = 1,2,··· ,P. As indicated in the next chapter, these collocation points maybeselectedwithalgorithmsimilartothatforselectingintegrationpointsinGaussian quadrature. 14 Chapter 3 Gaussian Random Field 3.1 Introduction The hydraulic conductivity or permeability is usually spatially heterogeneous, thus it is considered as a spatial random field. In the literature, the hydraulic conductivity is usually considered to have log-normal distribution. Thus the logarithm of the hydraulic conductivity is a Gaussian random field. The dimension of the Gaussian random field is corresponding to the grid blocks generated in the numerical simulations, and it can be very large in practice. In this chapter, an efficient technique, i.e. the Karhunen-Loeve expansion is utilized to parameterize the Gaussian random field, and the probabilistic collocation method is employed to quantify uncertainty of flow in porous media. Appli- cations to both single- and multi-phase flow are performed. 15 3.2 Karhunen-Loeve Expansion Suppose the hydraulic conductivity K is a spatial random field, and it can be written as a random space function K(x,θ), where x ∈ D and θ ∈ Θ (a probability space). Let Y(x,θ)=lnK(x,θ). One may write Y(x,θ)=Y(x)+Y ′ (x,θ) (3.1) where Y(x) is the mean and Y ′ (x,θ) is the fluctuation. The spatial structure of the random field may be described by the covariance C Y (x,y)=⟨Y ′ (x,θ)Y ′ (y,θ)⟩ (3.2) Since the covariance is bounded, symmetric and positive-definite, it can be decom- posed as [e.g., Ghanem and Spanos, 1991], C Y (x,y)= ∞ ∑ n=1 λ n f n (x)f n (y), (3.3) where λ n and f n (x) are eigenvalues and deterministic eigenfunctions, respectively, and can be solved from the following Fredholm equation: ∫ D C Y (x,y)f(x)dx=λf(y). (3.4) 16 Then the random process Y(x,θ) can be expressed as Y(x,θ)=Y(x)+ ∞ ∑ n=1 √ λ n f n (x)ξ n (θ), (3.5) whereξ n (θ)areorthogonalGaussianrandomvariableswithzeromeanandunitvariance. The expansion in equation (3.5) is called the Karhunen-Loeve (KL) expansion. The KL expansion, which is a spectral expansion, is optimal with mean square convergence when the underlying process is Gaussian [e.g., Ghanem and Spanos, 1991]. Although, in general, the eigenvalue problem of (3.4) has to be solved numerically, there exist analytical or semi-analytical solutions under certain conditions. For a one- dimensional random field with a covariance function C Y (x 1 ,y 1 )=σ 2 Y exp(−|x 1 −y 1 |/η), where σ 2 Y and η are the variance and the correlation length of the process, respectively, the eigenvalues and their corresponding eigenfunctions can be expressed as [e.g., Zhang and Lu, 2004], λ n = 2ησ 2 Y η 2 ω 2 n +1 , (3.6) and f n (x)= 1 √ (η 2 ω 2 n +1)L/2+η [ηω n cos(ω n x)+sin(ω n x)], (3.7) where ω n are positive roots of the characteristic equation (η 2 ω 2 −1)sin(ωL)=2ηωcos(ωL). (3.8) For problems in multi-dimension, if we assume that the covariance function C Y (x,y) is separable, for example C Y (x,y)=σ 2 Y exp(−|x 1 −y 1 |/η 1 −|x 2 −y 2 |/η 2 ) in a rectangular 17 domain D = {(x 1 ,x 2 ) : 0 ≤ x 1 ≤ L 1 ;0 ≤ x 2 ≤ L 2 }, the eigenvalues and eigenfunctions can be obtained by combining those in each dimension. For a non-separate covariance in adomainofarbitraryshape, theeigenvalueproblemof(3.4)hastobesolvednumerically. Furthermore, theKLexpansionisnotlimitedtostationaryrandomfields[LuandZhang, 2007]. From equation (3.5), one can getDσ 2 Y = ∑ ∞ n=1 λ n , whereD is the domain size, which indicates the total variance σ 2 Y is decomposed by an infinite series of eigenvalues λ n . Equation (3.8) has infinite number of positive roots. If the roots ω n are sorted in an increasing order, the related eigenvalues λ n are monotonically decreasing, which allows us to truncate the KL expansion to finite terms. The rate of decay of λ n determines the number of terms that are retained in the Karhunen-Loeve expansion, which determines the random dimensionality of the problem. 3.3 Single-Phase Flow Thesinglephasewaterflowinporousmediaisgovernedthefollowingcontinuityequation and Darcy’s law [e.g., Zhang, 2002]: S s ∂h(x,t) ∂t +∇·q(x,t)=g(x,t), (3.9) q(x,t)=−K(x)∇h(x,t), (3.10) subject to initial and boundary conditions h(x,0)=H 0 (x), x∈D, (3.11) 18 h(x,t)=H(x,t), x∈Γ D , (3.12) q(x,t)·n(x)=Q(x,t), x∈Γ N , (3.13) where h(x,t) is hydraulic head, K(x) is the hydraulic conductivity,q(x,t) is the specific discharge (flux), and g(x,t) is the source (or sink) term. H 0 (x) is the initial head in the domain D, H(x,t) is the prescribed head on Dirichlet boundary segment Γ D , Q(x,t) is the prescribed flux across Neumann boundary segments Γ N , n(x) = (n 1 ,...,n d ) T is an outward unit vector normal to the boundary Γ=Γ D ∪Γ N , and S S is the specific storage. We consider steady-state flow in saturated media satisfying the following equation: ∇·[K(x)∇h(x)]+g(x)=0, (3.14) subject to the boundary conditions h(x)=H(x), x∈Γ D , (3.15) K(x)∇h(x)·n(x)=−Q(x), x∈Γ N . (3.16) In this study, K(x) is treated as a random process (space function), thus equation (3.14) becomes a stochastic differential equation, whose solutions are no longer deter- ministic values but probability distributions or related moments. We aim to efficiently and accurately estimate the statistical properties of the hydraulic head in terms of the statistical moments of log transformed K(x). 19 3.3.1 KL Based Probabilistic Collocation Method 3.3.1.1 Implementation In previous studies [e.g. Webster et al., 1996, Tatang et al., 1997], PCM is used for para- metric sensitivity and uncertainty analysis and is limited to independent scalar random variables. When uncertain parameters are correlated and exist as a random field, the PCM may be combined with the KL expansion. Sarma et al. [2005] made use of the PCM and the KL expansion for optimizing reservoir production in the presence of an uncertain permeability field. In this work we combine the PCM with the KL expansion for analyzing uncertainty propagation for flow in random porous media, where random field and stochastic differential equation have to be considered. Four major steps are involved in implementation of this approach: 1) Representation of the input random field using the KL expansion in terms of a set of independent standard random variables (srv’s); 2) Approximation of the outputs of interest using the polynomial chaos expansion in terms of the set of srv’s; 3)Determinationoftheunknowncoefficientsinthepolynomialchaosexpansionusing the PCM technique; and 4) Evaluation of the statistical properties of the outputs. Substituting the KL expansion of Y(x) in equation (3.5) into the governing equation (3.14) yields ∇· { exp [ Y(x)+ N ∑ n=1 √ λ n f n (x)ξ n (θ) ] ∇h(x) } +g(x)=0. (3.17) 20 Now we take the second-order PCM for example to perform uncertainty analysis in random porous flow. We consider the approximation of the random field h(x) with a second-order polynomial chaos expansion and express it as ˆ h(x), ˆ h(x)=a 0 (x)+ N ∑ i=1 a i (x)ξ i + N ∑ i=1 a ii (x)(ξ 2 i −1)+ N−1 ∑ i=1 N ∑ j>i a ij (x)(ξ i ξ j ), (3.18) where N is the dimensionality of the input random field. This can be written as ˆ h(x)= P ∑ i=1 c i (x)Ψ i (), (3.19) where is a vector of dimension N. Usingtheprobabilisticcollocationmethod,welettheresidualbezeroatsomeselected sets of collocation points, specified as j = (ξ 1,j ,ξ 2,j ,...,ξ N,j ) T , where j = 1,2,...,P. It follows from equation (2.12) that at the jth set of collocation point, ˆ h j (x) satisfies the following equation, ∇· { exp [ Y(x)+ N ∑ n=1 √ λ n f n (x)ξ n,j (θ) ] ∇ ˆ h j (x) } +g(x)=0, (3.20) wherej =1,2,...,P. ItmeansthatweonlyneedtochooseP setsofcollocationpointsfor the random variables = (ξ 1 ,ξ 2 ,...,ξ N ) T , and to solve equation (3.20) P times (one at each set of collocation points j ). It should be noted that these values of hydraulic head ˆ h j (x) are solved independently through equation (3.20), which is identical to the original equation (3.14). That is, the codes for solving deterministic differential equations could be employed directly, like in the Monte Carlo method. 21 With these P sets of collocation points, we can solve equation (3.20) with the finite difference method, or any other methods, to obtain the P sets of the hydraulic head field ˆ h j (x). Here we denote the coefficients in equation (3.19) by c 1 (x), c 2 (x), ..., c P (x) in sequence, and C(x) = [c 1 ,c 2 ,...,c P ] T . The corresponding hydraulic head field at each set of collocation points (like each realization in the Monte Carlo method) are ˆ h 1 (x), ˆ h 2 (x), ..., ˆ h P (x), and h(x) = [ ˆ h 1 , ˆ h 2 ,..., ˆ h P ] T . In PCM, each of these hydraulic head fields is called a “representation” because they do not carry the same weight unlike the realizations in the direct sampling Monte Carlo method. Then we can rewrite (3.19) as, ZC(x)=h(x), (3.21) where Z, of elements Z ji = Ψ i ( j ), is a space-independent matrix of dimension P ×P, consisting of Hermite polynomials evaluated at the selected sets of collocation points. By solving the linear system of equations (3.21), the deterministic coefficients C(x) could be obtained readily. Note that the sets of collocation points should be so selected that the matrix Z satisfies the condition of rank(Z)=P. 3.3.1.2 Selection of Collocation Points The performance of PCM strongly depends on the choice of the collocation points. One particular scheme is, in analogy to Gaussian quadrature, to select the collocation points atagivenorderofpolynomialsfromtherootsofthenexthigherorderorthogonalpolyno- mialforeachuncertainparameter[e.g.,Websteretal.,1996,Tatangetal.,1997,Villadsen and Michelsen, 1978]. If the polynomial order is d, then the number of collocation points 22 availableis(d+1) N , whichisalwayslargerthanthenumberofcollocationpointsneeded. For example, the third-order Hermite polynomial is H 3 (ξ)=ξ 3 −3ξ, and the three roots ofthispolynomialare0,− √ 3, √ 3, respectively. Forthesecond-orderPCM,theP setsof collocationpoints(ξ 1,1 ,ξ 2,1 ,...,ξ N,1 ),(ξ 1,2 ,ξ 2,2 ,...,ξ N,2 ),...,(ξ 1,P ,ξ 2,P ,...,ξ N,P )arechosen from all the combinations of the three roots. Here we demonstrate how to choose the P setsofcollocationpoints. Werankthesethreerootsinorderofdecreasingprobability,(0, − √ 3, √ 3), since 0 has the highest probability for the standard Gaussian random variable with zero mean and unit variance. The first set of collocation point, i.e., (0,0,...,0) will alwayscontaintherootswiththehighestprobabilityforeachrandomvariable. Theother sets are selected by keeping as many of the variables of high probability as possible. In order to satisfy the condition of rank(Z) = P in (3.21), each new set j+1 of proposed collocation points must possess the property that the (j +1) th row of Z is linearly inde- pendent with the previous j rows in Z. Otherwise, the proposed set is rejected and the next set with equal or lower probability is tested. The last step is to make sure that all of itscolumnsarealsolinearlyindependentforthepurposeofsatisfyingtheconditionoffull rankforthematrixZ. WhenN ordislarge,theselectionprocesscanbecomputationally demanding. However, for a given type of polynomial chaos of known N and d, the P sets of collocation points are independent of physical problems. As for integration points in Gaussianquadrature, theselectionofcollocationpointsonlyhastobedoneonceandcan hencebetabulatedbeforehand. Analternative,whichcircumventstherequirementoffull rank for the matrix, is to estimate the set of coefficients in the polynomial chaos expan- sion via a regression based method called stochastic response surface method [Isukapalli et al., 1998]. In the regression method, a set of sample points is chosen and the model 23 outputs are equated to the estimates of the approximation at these points. However, the number of sample points must be more than the number of unknown coefficients to be estimated, resulting in a high computational effort. 3.3.1.3 Post-Processing Once we obtain the coefficients of equation (3.18), we could easily evaluate the statistical quantities, such as various moments and the probability density function, of the outputs based on expansion (3.18) or (3.19) by certain sampling methods such as Monte Carlo. Sincetheexpansionis reducedtoapolynomialform anditdoesnotinvolvesolving equa- tions, the evaluation is computationally efficient. Alternatively, one may directly derive the statistical moments of h(x) from equation (3.19), such as the mean and variance, ⟨h(x)⟩=c 1 (x), (3.22) σ 2 h = P ∑ j=2 c j (x) 2 ⟨Ψ 2 j ⟩. (3.23) Higher order moments can be obtained similarly. 3.3.2 Results of PCM and MC This and the following sections show some examples to illustrate the KL based PCM and examine its validity and applicability to flow in porous media. We first consider steady state flow in a one-dimensional domain of length L = 10 [L] (where [L] denotes any consistent length unit) and assume the forcing term to be zero. The boundary conditions are prescribed heads at the two ends, H 0 =7 [L] and H L =5 [L]. The correlation length 24 is given as η = 4[L]. The mean of the log hydraulic conductivity is given as ⟨Y⟩ = 0.0. Two cases with different spatial variability (i.e. σ 2 Y =1.0 or 2.0) are performed. Theeigenvalueandeigenfunctionλ n andf n (x), n=1,2,...,canbedeterminedanalyt- ically by solving equations (3.6) and (3.7). The eigenvalues are monotonically decreasing as illustrated in Figure 3.1(a), for cases with different correlation lengths (η = 2 and 4). Figure 3.1(b) shows the sum of eigenvalues as a function of number of terms included. Owing to the rapid decay, only the first 6 terms are retained in the KL expansion for η = 4. That is, the random dimensionality of {ξ n } is N = 6. One can see from Figure 3.1(b) that about 90% energy of the process is preserved in this case. For the second- order PCM, the total number of collocation points is P = 28. More precisely, the PCM involves 28 sets of collocation points, Inordertotestthesensitivityofthenumberofcollocationpoints,weconsiderahigher order PCM. As we know, based on the polynomial chaos expansion, the total number of collocation points is P = (N+d)! N! d! , where N and d are random dimension and polynomial order, respectively. It means that P grows fast with the random dimensionality and the order of the expansion. If we consider the third order polynomial chaos expansion, and the number of retained terms in K-L expansion is also 6, then the total number of collocation points grows to 84. In this study, the finite difference method is used to solve equation (3.20) and 151 physicalnodesarechosen. Figure3.2(a)and(b)showthemeanandvarianceofhydraulic headrespectively,whereresultsobtainedfromthesecond-orderandthird-orderPCMand MonteCarlomethodareallshown. NotethatinourMonteCarlomethod, weimplement the K-L expansion to generate the random field of the log hydraulic conductivity, based 25 1 5 10 15 20 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 n λ n /L η = 4.0 η = 2.0 1 10 20 30 40 50 60 0.4 0.5 0.6 0.7 0.8 0.9 1 n Σλ n /L η = 4.0 η = 2.0 (a) Eigenvalues (b) Sum of eigenvalues Figure 3.1: Series of eigenvalues and their finite sums, for η =4 and η =2 26 onequation(3.5)with60terms, whichisfoundtobeadequate. Wecanseethatboththe second-orderandthird-orderPCMsolutionsareclosetothosefromMonteCarlomethod. A comparison on computational effort and solution accuracy can be made between PCM and Monte Carlo method. The computational effort in solving equation (3.20) for each set of collocation points in the PCM is almost the same as solving equation (3.17) for each realization in the Monte Carlo method. We only need to solve equation (3.20) for 28 times in the second-order PCM and 84 times in the third-order PCM, much less than the number of realizations needed in Monte Carlo method. Figure 3.3 illustrates 4 setsofMonteCarlosimulations,correspondingto28, 84, 1000and10,000realizations. It is seen that Monte Carlo simulations with 28 or 84 realizations fail to obtain statistically accurate results and that even 1000 realizations are not enough compared to the Monte Carlo results with 10,000 realizations, which are believed to have converged. Obviously, PCM is far more efficient than Monte Carlo method. Both the Monte Carlo method and PCM involve sampling. The difference is that in the Monte Carlo method realizations aregeneratedrandomly, whereas inthePCM astructuralexpression (i.e. the polynomial chaos expansion) of the output random field is generated at first and then the collocation techniqueisadoptedaccordingtotheGaussianquadraturerule. Owingtothefastconver- gence of the polynomial chaos expansion and the high accuracy of Gaussian quadrature, PCM can yield quite good results with much less samplings than Monte Carlo method. To further test the validity of PCM, we have examined two more cases of different spatial variability, while keeping other conditions the same as in the first case. Owing to the particular boundary conditions in our examples, the mean head obtained from different approaches are very close to each other. We thus focus our discussion only on 27 0 1 2 3 4 5 6 7 8 9 10 5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 x <h> PCM, 2 nd order PCM, 3 rd order MC (10,000) η = 4.0, σ Y 2 = 1.0 (a) The mean of hydraulic head 0 1 2 3 4 5 6 7 8 9 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 x σ h 2 η = 4.0, σ Y 2 = 1.0 PCM, 2 nd order PCM, 3 rd order MC (10,000) (b) The variance of hydraulic head Figure 3.2: Comparisons of the mean and variance of hydraulic head derived from MC, and PCM up to the third order 28 0 1 2 3 4 5 6 7 8 9 10 5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 x <h> MC (28) MC (84) MC (1000) MC (10,000) η = 4.0, σ Y 2 = 1.0 (a) The mean of hydraulic head 0 1 2 3 4 5 6 7 8 9 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 x σ h 2 MC (28) MC (84) MC (1000) MC (10,000) η = 4.0, σ Y 2 = 1.0 (b) The variance of hydraulic head Figure 3.3: The mean and variance of hydraulic head derived from 4 sets of Monte Carlo simulations, corresponding to 28, 84, 1000, and 10,000 realizations 29 the head variance in the following sections. In Figure 3.4, we present the comparisons of the head variance derived from the second-order and third-order PCM and Monte Carlo method for σ 2 Y = 0.3 and 2. For small spatial variability, i.e., σ 2 Y = 0.3, the results from thesecond-orderPCMandthird-orderPCMarepracticallythesameandbotharealmost identical to the Monte Carlo result. As the spatial variability becomes large, i.e., σ 2 Y =2, the second-order and third-order PCM agree with the Monte Carlo result fairly well in spiteofsomedeviationsamongthethree. Thisobservationisencouragingasthevariance ofσ 2 Y =2,beingequivalenttothecoefficientofvariationofhydraulicconductivityCv K = σ K /⟨K⟩=253%, represents a large variability in hydraulic conductivity. We will discuss the case of huge spatial variability in section 3.3.4. 3.3.3 Comparisons with Other Stochastic Methods 3.3.3.1 Polynomial Chaos Expansion Approach Inthissection, wecomparePCMwiththetraditionalpolynomialchaosexpansion(PCE) approach. The traditional PCE approach utilizes the Galerkin scheme to determine the coefficientsofthePCE[GhanemandSpanos,1991,Ghanem,1998]. Forthesameproblem as in the previous section, we also use the KL expansion to represent the random field of the log hydraulic conductivity in the PCE approach. The governing equation is (3.17) and the PCE of hydraulic head is expressed as equation (3.19). 30 0 1 2 3 4 5 6 7 8 9 10 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 x σ h 2 PCM, 2 nd order PCM, 3 rd order MC (10,000) η = 4.0, σ Y 2 = 0.3 (a) 0 1 2 3 4 5 6 7 8 9 10 0 0.05 0.1 0.15 0.2 0.25 x σ h 2 η = 4.0, σ Y 2 = 2.0 PCM, 2 nd order PCM, 3 rd order MC (10,000) (b) Figure 3.4: Comparisons of head variance derived from PCM and MC, with different spatial variability 31 Fortheparticularone-dimensionalexampleonhand,assumingY(x)=0,substituting equation (3.19) into equation (3.17) and simplifying yields P ∑ j=1 d 2 c j (x) dx 2 Ψ j + P ∑ j=1 ( N ∑ i=1 √ λ i f ′ i (x)ξ i ) dc j (x) dx Ψ j =0. (3.24) WiththeGalerkintechnique,multiplyingequation(3.24)byeachoftheΨ k ,k =1,2,3,...,P, and taking the inner product, one has the following equation: d 2 c k (x) dx 2 ⟨Ψ 2 k ⟩+ P ∑ j=1 N ∑ i=1 √ λ i f ′ i (x) dc j (x) dx ⟨ξ i Ψ j Ψ k ⟩=0, (3.25) where the orthogonality of Ψ j , i.e., ⟨Ψ j Ψ k ⟩ = ⟨Ψ 2 k ⟩δ jk , is utilized. The finite difference method is also used to discretize this equation in physical space. Equation (3.25) results in a coupled system of P equations. Note that ⟨ξ i Ψ j Ψ k ⟩ should be calculated ahead of time, according to each term of ξ i , Ψ j and Ψ k . After solving the coupled system of equations for the deterministic coefficients c i (x), the statistical properties, such as the mean and variance, of the hydraulic head can be evaluated with equations (3.22)–(3.23). The post-processing is the same as in PCM. The head variances derived from the second- and third-order polynomial chaos ex- pansion approach, denoted as PCE, are presented in Figure 3.5, where the results from Monte Carlo simulations are added as the benchmarks. Three cases with different spatial variabilities are shown. It is seen that when the spatial variability is small, i.e. σ 2 Y =0.3, the results from both the second- and third-order PCE are in good agreement with those 32 from Monte Carlo method. As the spatial variability grows, the results from the second- order PCE have some deviations from those from Monte Carlo method. And the results from the third-order PCE are closer to those from Monte Carlo method. Since both the PCM and PCE use the same polynomial chaos expansion, their orders have the same meaning and the difference between them is how the coefficients of the expansion are evaluated. PCE utilizes the Galerkin scheme involved with the integration on each of the basis functions, which is accurate up to the given order of polynomials (subject to discretization errors). However, the PCM utilizes the optimal collocation techniqueaccordingtotheGaussianquadrature,withwhichtheaccuracyofthesolutions can be equal to or even higher than the given order of polynomials [e.g., Villadsen and Michelsen, 1978]. The superior performance of the second-order PCM relative to the second-orderPCEmaystemfromtheoptimalcollocationschemeutilizedintheGaussian quadrature. Withthesamenumberofterms,P,inthesamepolynomialchaosexpansion, PCM only needs to run the deterministic differential equation independently forP times; in PCE the P terms of coefficients are governed by P coupled equations, which lead to a large computational effort, especially for large-scale problems. And computing all the terms of ⟨ξ i Ψ j Ψ k ⟩ in PCE takes of order N ×P 2 operations, where N and P are the retained number of terms in the KL expansion and the number of terms in polynomial chaos expansion, respectively. In addition, since the PCM approach is non-intrusive in that the resulting equations are the same as the original equations, it can be easily implemented with existing codes and can be naturally parallelized. This advantage of PCM is even more pronounced for nonlinear problems. 33 0 1 2 3 4 5 6 7 8 9 10 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 x σ h 2 η = 4.0, σ Y 2 = 0.3 PCE, 2 nd order PCE, 3 rd order MC (10,000) (a) 0 1 2 3 4 5 6 7 8 9 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 x σ h 2 η = 4.0, σ Y 2 = 1.0 PCE, 2 nd order PCE, 3 rd order MC (10,000) (b) 0 1 2 3 4 5 6 7 8 9 10 0 0.05 0.1 0.15 0.2 0.25 x σ h 2 η = 4.0, σ Y 2 = 2.0 PCE, 2 nd order PCE, 3 rd order MC (10,000) (c) Figure 3.5: Comparisons of head variance derived from PCE and MC, with different spatial variability 34 3.3.3.2 KLME Approach Anotherefficientapproach,themoment-equationapproachbasedontheKarhunen-Loeve expansion (KLME) is recently developed by [Zhang and Lu, 2004]. Here we compare the KLME approach with the PCM. In the KLME approach, the log hydraulic conductivity Y(x) is expressed by the KL expansion and the hydraulic head h(x,t) is represented by the perturbation expansion, h(x) = h (0) +h (1) +h (2) +···. In this series, the order of each term is with respect to σ Y , the standard deviation of Y(x). And each term h (m) (x) (m≥1)isexpandedintermsoftheproductofmorthogonalGaussianrandomvariables, i.e., h (m) (x)= ∞ ∑ i 1 ,i 2 ,...,im=1 m ∏ j=1 ξ i j h (m) i 1 ,i 2 ,...,im (x), (3.26) where h (m) i 1 ,i 2 ,...,im are deterministic coefficients to be determined. Deriving from equations (3.9)–(3.13),onecanfinallyobtainthefollowingequationswithboundaryconditionsthat govern the coefficients [Zhang and Lu, 2004]: ∇ 2 h (0) (x)+∇⟨Y(x)⟩·∇h (0) (x)=− g(x) K G (x) , (3.27) h (0) (x)=H(x), x∈Γ D , (3.28) ∇h (0) (x)·n(x)=−Q(x)/K G (x), x∈Γ N , (3.29) and for m≥1, ∇ 2 h (m) i 1 ,i 2 ,...,im (x)+∇⟨Y(x)⟩·∇h (m) i 1 ,i 2 ,...,im (x)= (−1) m+1 g(x) m!K G (x) m ∏ j=1 f i j (x) 35 − 1 m ∑ P i 1 ;:::;im ∇f i 1 (x)·∇h (m−1) i 2 ,...im (x), (3.30) h (m) i 1 ,i 2 ,...,im (x)=0, x∈Γ D , (3.31) ∇h (m) i 1 ,i 2 ,...,im (x)·n(x)= (−1) m+1 Q(x) m!K G (x) m ∏ j=1 f i j (x), x∈Γ N , (3.32) where K G (x) = exp(⟨Y(x)⟩) and f n (x) is used to represent ˜ f n (x) = √ λ n f n (x) for sim- plicity. It is seen that the equations for different order terms of the hydraulic head are recursive and the equations with the same order are independent. To obtain the coeffi- cientsh (m) i 1 ,i 2 ,...,im (x),wherei j =1,...,n,thenumberoftimesrequiredtosolveindependent differentialequationsisS m =n(n+1)...(n+m−1)/m!. Oncethesecoefficientsaresolved, onecandirectlycalculatethemean,variance,andthecovarianceofhydraulichead[Zhang and Lu, 2004]. We solve the same problem as in the previous sections using the KLME approach and present the results in Figure 3.6. We also examine three cases with different spatial variabilities σ 2 Y and add the Monte Carlo results as benchmarks. Up to second-order (in σ 2 Y )correctionforheadvariance,oneneedstosolveforthecoefficientsh (1) i ,h (2) ij ,andh (3) ijk , the numbers of terms included, which are chosen as 6, 4 and 4, respectively, i.e., index i in h (1) i running up to 6 and each index in h (2) ij running up to 4, and so on. Owing to the symmetry of the coefficients such as h 2 ij , the computational times are reduced. The total number of times to solve the deterministic equations for h (1) i , h (2) ij , and h (3) ijk is 6 + 10 + 20 = 36. From Figure 3.6, one can see that the first-order (in σ 2 Y ) approximations of head vari- ance have some deviations from the Monte Carlo results. When σ 2 Y is small or moderate, 36 i.e., σ 2 Y = 0.3 or 1.0, the second-order head variances from the KLME approach are in good agreement with the Monte Carlo results. However, when the spatial variability becomes large, i.e., σ 2 Y = 2.0, the second-order head variance from the KLME approach deviates from the Monte Carlo result. We try to increase the number of terms included in h (i) (x), i = 1,2,..., consequently bringing about a larger computational effort, but do not find much improvement. It has been found by Zhang and Lu [2004] that for similar cases the approximation can be improved with third-order corrections for σ 2 Y =2. It should be noted that the order in the KLME is different from that in the PCM (or PCE).TheorderinKLMEcomesfromtheperturbationexpansionandiswithrespectto σ Y , whereas that in PCM or PCE is with respect to the degree of the polynomial chaos expansion. If all terms in the perturbation expansion are collected, one can find that the KLME has the same construction of expansion as the polynomial chaos expansion with the PCM or PCE. In the KLME, the number of terms included in approximating h (i) (x) can be distinctly selected for each order i [Zhang and Lu, 2004]. Generally, the number decreases as the increase of order i in h (i) (x). While in the PCM or PCE, for different order polynomial chaos expansions the same number of terms in the KL expansion is selected as the random dimensionality. 3.3.4 Discussions 3.3.4.1 Effect of Correlation Length η AsshowninFigure3.1,therateofdecayintheeigenvaluesisdependentonthecorrelation length η relative to the domain length L. To further test the effect of correlation length on the recently developed approaches, i.e., the PCM and the KLME, we perform 3 cases 37 0 1 2 3 4 5 6 7 8 9 10 0 0.01 0.02 0.03 0.04 0.05 0.06 x σ h 2 η = 4.0, σ Y 2 = 0.3 KLME, 1 st order KLME, 2 nd order MC (10,000) (a) 0 1 2 3 4 5 6 7 8 9 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 x σ h 2 η = 4.0, σ Y 2 = 1.0 KLME, 1 st order KLME, 2 nd order MC (10,000) (b) 0 1 2 3 4 5 6 7 8 9 10 0 0.06 0.12 0.18 0.24 0.3 0.36 x σ h 2 η = 4.0, σ Y 2 = 2.0 KLME, 1 st order KLME, 2 nd order MC (10,000) (c) Figure 3.6: Comparisons of head variance derived from KLME and MC, with different spatial variability 38 for η = 2, with different spatial variabilities, i.e., σ 2 Y = 0.3, 1.0 and 2.0. We also include the Monte Carlo simulation results for comparison. The head variance obtained from those approaches are presented in Figure 3.7. Note that the number of physical nodes is also chosen as 151 and the number of realizations in Monte Carlo simulations is 10,000. A smaller number of physical nodes or realizations in Monte Carlo simulations could lead to unstable results. The retained number in the KL expansion for PCM is also chosen as 6, the same as that for η = 4. And for the second-order KLME, the numbers of terms included in h (1) i , h (2) ij , and h (3) ijk are also chosen as 6, 4, and 4, respectively. From Figure 3.7, we observe that for small or moderate spatial variability, i.e., σ 2 Y = 0.3 or 1, the results obtained from both the second-order PCM and the second-order KLME agree well with those from Monte Carlo method. When σ 2 Y is large, i.e., σ 2 Y = 2, the second-order PCM result is still close to the Monte Carlo result while the second-order KLME result exhibits some deviation. 3.3.4.2 Effect of Large Spatial Variability σ 2 Y As illustrated in the previous sections, for moderate spatial variability, the second-order or third-order PCM (or PCE) and the second-order KLME can obtain quite accurate results. However, for large variability there are some deviations. In this section, we examine some cases with a huge spatial variability, i.e, σ 2 Y = 4.0, corresponding to the coefficient of variation of hydraulic conductivity Cv K =732% [Zhang, 2002]. We explore the higher order PCM solutions for the purpose of improving the accuracy. For the cases of σ 2 Y = 4.0 and with different correlation lengths η, we perform the PCM from the second order up to the fourth order. The retained number in the KL expansion is 6, 39 thus the total number of collocation points for the fourth-order PCM reaches 210. Figure 3.8 (a) and (b) show the head variances for η = 4 and η = 2, respectively, derived from the second-order to fourth-order PCM and Monte Carlo method with 10,000 realizations. From Figure 3.8(a) and (b), it is seen that the head variances from both the second- order and fourth-order PCM are quite close to those from Monte Carlo simulations, for both η = 4 and η = 2. However, the third-order solutions deviate significantly from their Monte Carlo counterparts. The inferior performance of the third-order PCM is related to the fact that its collocation points are selected from the roots of fourth order Hermitepolynomial, whichdonotcontaintheoriginpointthatcorrespondstotheregion of highest probability for a standard Gaussian random variable. This indicates that an evenorderPCMmayperformbetterthanthenexthigheroddorderPCMinthepresence of huge variability in hydraulic conductivity when such a scheme is invoked for selecting collocation points. 3.3.4.3 Probability Density Function Once the dependent random field such as the hydraulic head h(x) is approximated by eitherthepolynomialchaosexpansion(3.18)inthePCM(orthePCE)ortheperturbation polynomial expansion (3.26) in the KLME, the probability density function (PDF) of the randomfieldcanbesimulatedwithcertainsamplingmethods. Takingthefirstcaseofη = 4andσ 2 Y =1.0forexample,wecalculatethePDFofh(x)atthreepositions(i.e.,x=2,4, and6), withthePCM,KLME,andMonteCarlomethod. Thecomparisonsarepresented inFigure3.9. Onecanseethattheresultsobtainedfromboththesecond-orderPCMand the second-order KLME are close to those from Monte Carlo simulations. Note that one 40 needs to solve the differential equations for 10,000 times with the Monte Carlo method. However, with the PCM or the KLME approach, one only needs to solve the differential equations for 28 or 36 times in order to construct the polynomial expansions, based on which 10,000 samplings are then drawn to approximate the PDF. Since the dependent random field is explicitly expressed in a polynomial form, the computational time is significantlyreducedcomparedtosolvingoriginaldifferentialequationswithMonteCarlo method. 3.3.4.4 Illustrative Examples in 2D In this section, we consider two cases in a two-dimensional domain of saturated het- erogeneous medium, which is a square of size L 1 = L 2 = 10 [L], uniformly discretized into 40 × 40 square elements. The non-flow conditions are prescribed at two lateral boundaries. The hydraulic head is prescribed at the left and right boundaries as 10.5 [L] and 10.0 [L], respectively. The mean of the log hydraulic conductivity is given as ⟨Y⟩ = 0.0. Assume the covariance function of the log hydraulic conductivity is C Y (x,y) = σ 2 Y exp(−|x 1 −y 1 |/η 1 −|x 2 −y 2 |/η 2 ), where η 1 = η 2 = 4.0 and σ 2 Y = 1.0. For the first case, there is no source term, i.e. g(x) = 0. For the second case, there is a pumping well extracting water at the flow rate of 1 [L 3 /T] at the center of the domain. For the first case without well, Figure 3.10 shows the head variance along the section of x 2 =5. Forthesecondcasewithwell,Figure3.11(a)and(b)showthemeanandvariance of hydraulic head along the diagonal line, respectively. In both cases, we found it suffice to retain 7 modes (N = 7) in KL expansion for the two-dimensional log transformed 41 0 1 2 3 4 5 6 7 8 9 10 0 0.01 0.02 0.03 0.04 0.05 0.06 x σ h 2 η = 2.0, σ Y 2 = 0.3 PCM, 2 nd order KLME, 2 nd order MC (10,000) (a) 0 1 2 3 4 5 6 7 8 9 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x σ h 2 η = 2.0, σ Y 2 = 1.0 PCM, 2 nd order KLME, 2 nd order MC (10,000) (b) 0 1 2 3 4 5 6 7 8 9 10 0 0.05 0.1 0.15 0.2 0.25 0.3 x σ h 2 η = 2.0, σ Y 2 = 2.0 PCM, 2 nd order KLME, 2 nd order MC (10,000) (c) Figure 3.7: Comparisons of head variance derived from PCM, KLME and MC, with different spatial variability for η =2 42 0 1 2 3 4 5 6 7 8 9 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 x σ h 2 η = 4.0, σ Y 2 = 4.0 PCM, 2 nd order PCM, 3 rd order PCM, 4 th order MC (10,000) (a) 0 1 2 3 4 5 6 7 8 9 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x σ h 2 η = 2.0, σ Y 2 = 4.0 PCM, 2 nd order PCM, 3 rd order PCM, 4 th order MC (10,000) (b) Figure 3.8: Head variance derived from MC and PCM (up to the fourth order) for huge spatial variability, i.e. σ 2 Y =4.0 43 4.5 5 5.5 6 6.5 7 7.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 h(2) PDF η = 4.0, σ Y 2 = 1.0 PCM, 2 nd order KLME, 2 nd order MC (10,000) (a) 4 4.5 5 5.5 6 6.5 7 7.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 h(4) PDF η = 4.0, σ Y 2 = 1.0 PCM, 2 nd order KLME, 2 nd order MC (10,000) (b) 4.5 5 5.5 6 6.5 7 7.5 8 0 0.2 0.4 0.6 0.8 1 1.2 h(6) PDF η = 4.0, σ Y 2 = 1.0 PCM, 2 nd order KLME, 2 nd order MC (10,000) (c) Figure 3.9: Comparisons of the probability density function (PDF) of hydraulic head at three positions (i.e., x=2, 4, and 6) obtained with PCM, KLME, and MC 44 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 x 10 −3 x 1 σ h 2 2 nd PCM MC 10,000 η = 4.0, σ Y 2 = 1.0 Figure 3.10: Comparison of head variance at x 2 =5 in 2D obtained with PCM and MC hydraulic field, thus the total number of terms in the second-order polynomial chaos ex- pansion is 36. Results from the second-order PCM with 36 simulations agree quite well with those from the Monte Carlo method with 10,000 simulation runs. It is seen that the computational efficiency and accuracy observed in one dimension also carries to two dimensions. 3.3.5 Summary The probabilistic collocation method (PCM) combined with the Karhunen-Loeve (KL) expansion is proposed for uncertainty analysis of flow in random porous media. The underlying (input) Gaussian random field is represented by the KL expansion and the dependent random field such as the hydraulic head is expressed by the polynomial chaos expansion. ThePCMisusedtodeterminethecoefficientsinthepolynomialchaosexpan- sionbyevaluatingthehydraulicheadatdifferentsetsofcollocationpoints. Thisapproach 45 0 1 2 3 4 5 6 7 8 9 10 9.6 9.7 9.8 9.9 10 10.1 10.2 10.3 10.4 10.5 10.6 x 1 <h> 2 nd PCM MC 10,000 η = 4.0, σ Y 2 = 1.0 (a) Mean hydraulic head 0 1 2 3 4 5 6 7 8 9 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 x 1 σ h 2 2 nd PCM MC 10,000 η = 4.0, σ Y 2 = 1.0 (b) Variance of hydraulic head Figure 3.11: Comparisons of mean and variance of hydraulic head along the diagonal line in 2D with a pumping well obtained with PCM and MC 46 is appealing because it results in independent deterministic differential equations, which, similartotheMonteCarlomethod,canbeimplementedwithexistingcodes. ThePCMis applied to several cases of flow in random porous media in one and two dimensions, with differentspatialvariabilitiesandcorrelationlengths. ItiscomparedwiththeMonteCarlo method, the traditional polynomial chaos expansion (PCE) approach based on Galerkin scheme, and the moment-equation approach based on Karhunen-Loeve expansion. The following conclusions can be addressed: 1. With a relatively small computational effort, the probabilistic collocation method (PCM) based on Karhunen-Loeve expansion is feasible for accurately quantifying uncertainty associated with flow in random porous media, where the random field and stochastic differential equation have to be considered. 2. While both the PCM and the Galerkin PCE approach rely on polynomial chaos ex- pansion to approximate the hydraulic head, the difference depends on how the coefficients of the expansion are evaluated. The traditional PCE approach uses a Galerkin scheme and leads to a set of coupled equations whereas the PCM uses a collocationapproachandresultsinindependentdeterministicdifferentialequations. The PCM is thus generally more efficient than the Galerkin PCE approach. 3. SimilartotheMonteCarlomethod,thePCMcanbeeasilyimplementedwithexisting codes and naturally parallelized. Both the direct sampling Monte Carlo method and the PCM involve sampling. The Monte Carlo method resorts to the law of large numbers by generating a large number of equally probable realizations and is thus computationally demanding; the PCM builds upon an optimal approximation 47 evaluatedatselectedsetsofcollocationpointsandhencereducesthecomputational efforts significantly. 4. The moment-equation approach based on Karhunen-Loeve expansion (KLME) first decomposes the hydraulic head into terms of different order of σ Y with the per- turbation expansion and then expresses each term with a polynomial expansion. If all terms in the perturbation expansion are collected, the KLME has the same construction as the polynomial chaos expansion, which is used by the PCM and the PCE approach. The difference is how the approximations are made in all these methods. In the KLME, the equations for different order terms of the hydraulic headarerecursive,i.e.,theequationsforthehigherordertermsdependonthelower orderterms; theequationsatthesameorderareindependent. Thenumberofterms retained in the KL expansion of the log hydraulic conductivity can be different for different order terms of hydraulic head. 5. The PCM, the Galerkin PCE approach, and the KLME are all fairly accurate at least forthespatialvariabilityσ 2 Y aslargeas2.0(correspondingtothecoefficientofvari- ation of hydraulic conductivity as 253%). And the PCM and the KLME are more efficient than the Galerkin PCE of coupled equations. With small computational efforts, the probability density functions (PDF) of the dependent random fields can be obtained with the PCM and the KLME. In general, for moderate spatial variability the second-order PCM (or KLME) is enough to obtain quite accurate results while for larger spatial variability, a higher order PCM (or KLME) may be 48 needed. It should be noted that PCM is not restricted to the flow type and config- uration considered in this study. Next section discusses the applications of PCM in multi-phase flow in porous media. 3.4 Multi-phase Flow 3.4.1 Mathematical Formulations The three-phase black oil model is expressed by the following continuity equations [Chen et al., 2006b], ∂ ∂t ( ϕS i B i ) −∇· [ k(x)k ri µ i B i (∇p i −ρ i ℘∇z) ] =q i , i=w,o (3.33) ∂ ∂t [ ϕ ( S g B g + S o R so B o )] − ∇· [ k(x)k ro R so µ o B o (∇p o −ρ o ℘∇z)+ k(x)k rg µ g B g (∇p g −ρ g ℘∇z) ] =q g (3.34) where x and t denote the position and time, respectively; w,o, and g denote the three phases, i.e., water, oil, and gas, respectively; ϕ is the porosity of the media; B i and ρ i (i = w,o,g) are the formation volume factor and density of phase i, respectively; R so is the solubility of gas in oil; q i is the source or sink term; S i and p i are the saturation and pressure of fluid i, respectively; k(x) is the absolute (intrinsic) permeability; k ri is the relative permeability of fluid i, which is a function of S i ; µ i is the viscosity of fluid i; z is the depth; and ℘ is the gravitational acceleration. 49 Equations (3.33) and (3.34) are coupled with S w +S o +S g =1 (3.35) p cwo (S w )=p o −p w (3.36) p cgo (S g )=p g −p o (3.37) where p cwo which is a function of S w , is the capillary pressure between water and oil, and p cgo which is a function of S g , is the capillary pressure between gas and oil. Based on the governing equations (3.33)-(3.37) subject to certain initial and boundary conditions, one can solve for the fluid saturation and pressure. In this study, the properties of the heterogeneousporousmedia,suchasthepermeabilityandporosityaretreatedasrandom functions, thus the governing equations become stochastic partial differential equations, whosesolutionsarenolongerdeterministicvaluesbutprobabilitydistributionsorrelated statistical moments. The permeability or porosity is considered as a random field which has log-normal distribution. The KL expansion can be used to represent the Gaussian random field with known covariance function. The fluid saturation and pressure are dependent on medium properties, such as the permeability and porosity. While the covariance of the dependent randomprocessesareyettobefound,theKLexpansioncannotbeusedtorepresenttheir random structures. Instead, the polynomial chaos expansion can be used to effectively 50 express the dependent random fields. For example, we express the output random fields S i (x,t) and p i (x,t) with the polynomial chaos expansions, S i (x,t,θ)=a 0 (x,t)+ ∞ ∑ i 1 =1 a i 1 (x,t)Γ 1 (ξ i 1 (θ))+ ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 a i 1 i 2 (x,t)Γ 2 (ξ i 1 (θ),ξ i 2 (θ)) + ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 i 2 ∑ i 3 =1 a i 1 i 2 i 3 (x,t)Γ 3 (ξ i 1 (θ),ξ i 2 (θ),ξ i 3 (θ))+..., (3.38) p i (x,t,θ)=b 0 (x,t)+ ∞ ∑ i 1 =1 b i 1 (x,t)Γ 1 (ξ i 1 (θ))+ ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 b i 1 i 2 (x,t)Γ 2 (ξ i 1 (θ),ξ i 2 (θ)) + ∞ ∑ i 1 =1 i 1 ∑ i 2 =1 i 2 ∑ i 3 =1 b i 1 i 2 i 3 (x,t)Γ 3 (ξ i 1 (θ),ξ i 2 (θ),ξ i 3 (θ))+..., (3.39) where the coefficientsa 0 (x,t),···, a i 1 i 2 ,...,i d (x,t) andb 0 (x,t),···, b i 1 i 2 ,...i d (x,t) are deter- ministicfunctionsofx,andΓ d (ξ i 1 ,...,ξ i d )areorthogonalpolynomialchaosoforderdwith respect to the random variables (ξ i 1 ,...,ξ i d ). For independent standard Gaussian random variables (ξ i 1 ,...,ξ i d ), Γ d (ξ i 1 ,...,ξ i d ) are the multi-dimensional Hermite Polynomials of degree d expressed as Γ d (ξ i 1 ,...,ξ i d )=(−1) d e 1 2 T ∂ d ∂ξ i 1 ...∂ξ i d [e − 1 2 T ], (3.40) where is a vector denoting (ξ i 1 ,...,ξ i d ) T . Equations (3.38) and (3.39) can be expressed as S i (x,t)= P ∑ j=1 c j (x,t)Ψ j () (3.41) 51 p i (x,t)= P ∑ j=1 d j (x,t)Ψ j () (3.42) where c j (x,t) and d j (x,t) are deterministic coefficients and Ψ j () are polynomials of random vector . There are P =(N+d)!/(N!d!) terms in each of the polynomials chaos expansions, where N is the random dimensionality and d is the degree of polynomial chaos expansion. The P terms of coefficients can be determined through the PCM . The PCM is demonstrated for multi-phase flow with the commercial Eclipse black oil simulator. Different scenarios of one-dimensional, (oil and water) two-phase flows are used to analyze the impact of various factors, such as type of variability, variance, correlation length, viscosity ratio, and relative permeability type. On the basis of the 9 th SPE comparative project, a three-dimensional three-phase reservoir model with a moderate grid size is designed for reservoir performance forecast. 3.4.2 Case Study: 1D Water Flooding At first a one-dimensional (water and oil) system is considered. The reservoir depth is 4000ft,andthelengthLis1000ft,uniformlydividedby40gridblocks. Assumeconstant water injection rate 100 STB/day for the injection well at the inlet, and constant bottom hole pressure 3000 psi for the production well at the outlet. The capillary pressure is neglected for simplicity. The following Corey type relative permeability functions are considered: k rw =(S ∗ ) n , (3.43) k ro =(1−S ∗ ) m , (3.44) 52 S ∗ = S w −S wc 1−S or −S wc , (3.45) where S wc and S or are the connate water saturation and residual oil saturation, respec- tively. Because the PCM makes a direct use of existing simulators, it is not restrictive to particular types of relative permeability curves or tabulations. The following illustrative examples employ the linear Corey type relative permeability curves, i.e., n = m = 1 in Equations (3.43) and (3.44), unless specified otherwise. Both the connate water satura- tion and residual oil saturation are equal to zero, S wc =S or =0. 3.4.2.1 Random Porosity Field We first assume the porosity ϕ(x) to be a random field while keeping the permeability deterministic. Let β(x)=lnϕ(x). The covariance of β(x) satisfies the exponential form, C β (x 1 ,x 2 )=σ 2 β exp(−|x 1 −x 2 |/η), (3.46) where σ 2 β and η are the variance and the correlation length of the process, respectively. β(x) is represented by the Karhunen-Loeve (KL) expansion, β(x)=⟨β(x)⟩+ N ∑ n=1 √ λ n f n (x)ξ n , (3.47) 53 Substituting equation (3.47) into (3.33) yields, ∂ ∂t exp [ ⟨β(x)⟩+ ∑ N n=1 √ λ n f n (x)ξ n ] S i B i − ∇· { k(x)k ri µ i B i (∇p i −ρ i ℘∇z) } =q i , i=w,o (3.48) With the PCM, we only need to choose certain sets of collocation points for random variables j =(ξ 1 ,ξ 2 ,...ξ N ) T j , wherej =1,2,··· ,P. Foreachset j , thewatersaturation and pressure fields are solved using the reservoir simulator, in a manner similar to Monte Carlo simulations. The coefficients c j and d j in (3.41)-(3.42) are then solved from the P sets of saturation and pressure fields, respectively. Statistical properties of the output random fields can then be evaluated with these coefficients using equations (3.22)-(3.23). Here we assume ⟨β(x)⟩ as -0.9613, corresponding to the mean porosity about 0.4, σ β as 0.1, and η/L as 2/5. The eigenvalue and eigenfunction λ n and f n (x), n=1,2,..., can be determined by solving equations (3.6) and (3.7). The eigenvalues are monotonically decreasing as illustrated in Figure 3.12(a), for cases with different correlation lengths (η/L = 1/5 and 2/5). Figure 3.12(b) shows the sum of eigenvalues as a function of number of terms included. Owing to the rapid decay, only the first 6 terms are retained intheKLexpansionforη/L=2/5. Thatis, therandomdimensionalityof{ξ n }isN =6. OnecanseefromFigure3.12(b)thatabout90%energyoftheprocessispreservedinthis case. For the second-order PCM, the total number of collocation sets is P = 28. More precisely, the PCM involves 28 sets of collocation points, each of which requires one run of the simulator. 54 1 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 n λ n /L η/L = 2/5 η/L = 1/5 (a) 1 2 4 6 8 10 12 14 16 18 20 0.4 0.5 0.6 0.7 0.8 0.9 1 n Σλ n /L η/L = 2/5 η/L = 1/5 (b) Figure 3.12: Series of eigenvalues and their finite sums 55 Different cases with different viscosity ratios of oil to water, defined as m = µ o /µ w , are performed. Figures 3.13 and 3.14 depict the mean values of the water saturation and pressure and the associated standard deviations at 100 days for different viscosity ratios (m = 0.5 and 2). The results obtained from the second-order PCM and Monte Carlo (MC) simulations with 1000 realizations are presented together for comparison. Excellentagreementsareachievedinallthecases,bothforwatersaturationandpressure. However, the computational cost of the PCM is much less than the MC. The PCM just needs to perform 28 deterministic simulations while the MC requires 1000 simulations. It is observed from the figures that the standard deviations of pressure are small compared to their mean values. This indicates that the random heterogeneity of porosity has a small effect on the pressure variations for this particular setup of porosity variability and boundary conditions. And the viscosity ratio affects distributions of the mean and standard variation of the water saturation: the spreading of the mean saturation is larger (indicating larger transition zone) when oil is more viscous than water, and the standard deviation profile of water saturation has a higher peak and is more compact when oil is less viscous. This is consistent with the finding of [Zhang and Tchelepi, 1999]. As shown in Figure 3.12, the rate of decay in the eigenvalues is dependent on the correlation length η relative to the domain length L. When the correlation length is smaller, the convergence of the eigenvalues in the KL expansion becomes slower. To further test the effect of correlation length on the PCM, we perform another case with smaller correlation length (η/L = 1/5), keeping other conditions the same as in the previouscaseofviscosityratiom=2. Inthiscase,tokeepthenumberoftermsretainedin the KL expansion as 6, about 80% of the total energy of the random process is preserved, 56 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> m = 0.5 PCM MC (a) 0 200 400 600 800 1000 0 0.05 0.1 0.15 x, ft σ Sw m = 0.5 PCM MC (b) 0 200 400 600 800 1000 3000 3100 3200 3300 3400 x, ft <Pressure>, psi m = 0.5 PCM MC (c) 0 200 400 600 800 1000 0 1 2 3 4 x, ft σ Pressure , psi m=0.5 PCM MC (d) Figure 3.13: Mean and standard deviation of water saturation and pressure obtained from 2nd order PCM and the MC (with 1000 realizations) for the case with standard deviationoflogporosityi.e. σ lnϕ =0.1andviscosityratiom=0.5: (a)and(b)forwater saturation; (c) and (d) for pressure 57 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> m = 2 PCM MC (a) 0 200 400 600 800 1000 0 0.01 0.02 0.03 0.04 0.05 x, ft σ Sw m = 2 PCM MC (b) 0 200 400 600 800 1000 3000 3200 3400 3600 3800 4000 4200 x, ft <Pressure>, psi m = 2 PCM MC (c) 0 200 400 600 800 1000 0 2 4 6 8 10 x, ft σ Pressure , psi m = 2 PCM MC (d) Figure3.14: Meanandstandarddeviationofwatersaturationandpressureobtainedfrom 2nd order PCM and MC (with 1000 realizations) for the case with standard deviation of log porosity i.e. σ lnϕ = 0.1 and viscosity ratio m = 2: (a) and (b) for water saturation; (c) and (d) for pressure 58 as seen from Figure 3.12 (b). The mean and standard deviation of water saturation (and pressure) at 100 days obtained with the second-order order PCM and the MC are presentedinFigure3.15. ItcanbeseenthatthePCMresultsarestillinperfectagreement with the MC results. Once the dependent random field such as the water saturation is approximated by the polynomial chaos expansion (3.41), the probability density functions (PDFs) of the random field can be simulated with certain sampling methods. By using the case with η/L = 2/5 and viscosity ratio m = 2 as an example, we calculate the PDFs of water saturation and pressure at several locations, with the PCM and direct sampling MC simulations. Figures 3.16 (a) and (c) show the PDFs of water saturation at 100 days at two locations (x = 150 ft and 350 ft), respectively. And the PDFs of the pressure at the two locations are shown in Figures 3.16 (b) and (d), respectively. It can be seen that the PCM results agree with the MC results very well. With the direct sampling Monte Carlo simulations, 1000 realizations are solved with the simulator and then the results are statistically evaluated to obtain the PDFs. However, with the PCM, one only needstorun28deterministicsimulationsinordertoconstructthepolynomialexpansions, based on which 10,000 samplings are then drawn to approximate the PDFs. Since the dependent random field is explicitly expressed in a polynomial form, the computational time is significantly reduced compared to running the reservoir simulator with the Monte Carlo method. 59 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> m = 2, η/L = 1/5 PCM MC (a) 0 200 400 600 800 1000 0 0.01 0.02 0.03 0.04 0.05 x, ft σ Sw m=2, η/L = 1/5 PCM MC (b) 0 200 400 600 800 1000 3000 3200 3400 3600 3800 4000 4200 x, ft <Pressure>, psi m = 2, η/L = 1/5 PCM MC (c) 0 200 400 600 800 1000 0 2 4 6 8 x, ft σ Pressure , psi m = 2, η/L = 1/5 PCM MC (d) Figure3.15: Meanandstandarddeviationofwatersaturationandpressureobtainedfrom 2nd order PCM and MC (with 1000 realizations) for the case with standard deviation of log porosity σ lnϕ =0.1, correlation length η/L=1/5, and viscosity ratio m=2: (a) and (b) for water saturation; (c) and (d) for pressure 60 0.2 0.3 0.4 0.5 0.6 0.7 0 2 4 6 8 10 Sw (at x= 150 ft) PDF PCM MC (a) 3880 3900 3920 3940 3960 0 0.01 0.02 0.03 0.04 0.05 0.06 Pressure (at x= 150 ft), psi PDF PCM MC (b) 0 0.05 0.1 0.15 0.2 0 5 10 15 20 25 Sw (at x= 350 ft) PDF PCM MC (c) 3730 3735 3740 3745 3750 0 0.05 0.1 0.15 0.2 Pressure (at x= 350 ft), psi PDF PCM MC (d) Figure 3.16: Probability density functions of water saturation S w and pressure at two different locations (x = 150 ft and 350 ft), obtained from 2nd order PCM and MC (with 1000 realizations) for the case with standard deviation of log porosity σ lnϕ = 0.1, correlation lengthη/L=2/5, and viscosity ratiom=2: (a) and (c) for water saturation; (b) and (d) for pressure 61 3.4.2.2 Random Permeability Field Next we consider the permeability k(x) to be a random field while keeping the porosity deterministic. The covariance of Y(x)=lnk(x) satisfies the exponential form, C Y (x 1 ,x 2 )=σ 2 Y exp(−|x 1 −x 2 |/η) (3.49) where σ 2 Y and η are the variance of Y(x) and correlation length of the random field, respectively. Substituting the KL expansion for Y(x) into (3.33) yields, ∂ ∂t ( ϕS i B i ) − ∇· exp [ ⟨Y(x)⟩+ ∑ N n=1 √ λ n f n (x)ξ n ] k ri µ i B i (∇p i −ρ i ℘∇z) =q i , i=w,o (3.50) The PCM can be implemented as described in the previous section. By choosing P sets of collocation points for random variables j , we can obtain the P number of water saturation and pressure fields with the reservoir simulator and then compute the coefficients c j and d j in (3.41)-(3.42), from which statistical properties of the output random fields are evaluated. Assume ⟨Y(x)⟩ = 5.0 [+ln md], corresponding to the mean permeability of about 148 md. Consider the same reservoir conditions as in the previous case of correlation lengthη/L=2/5andviscosityratiom=1. Differentspatialvariabilitiesofpermeability 62 (i.e., σ 2 Y = 0.25,1.0, and 4.0) are explored. The mean (and standard deviation) of water saturation and pressure at 100 days, obtained from the second-order and fourth-order PCM and Monte Carlo simulations, are presented in Figures 3.17 - 3.19. The number of terms retained in the KL expansion is 6, thus the second-order PCM needs 28 sets of collationpointsandthefourth-orderPCMrequires210setsofcollocationpoints. Forthe cases of σ 2 Y = 0.25, both the second-order and fourth-order PCM agree quite well with Monte Carlo simulations of 1000 realizations, both for the water saturation and pressure. When the spatial variability of the permeability becomes larger, i.e., in the cases of σ 2 Y =1.0, and4.0, sodoestheuncertaintyinwatersaturationandpressure. MonteCarlo simulations with up to 10,000 realizations are performed for comparison. It is seen from Figures3.18-3.19thatthestandarddeviationofwatersaturationaswellasthemeanand standarddeviationofpressurefromtheMCwith1000realizationsaredifferentfromthose with 10,000 realizations. The convergence of MC results shall be investigated in more details. In general, it is easier to accurately estimate the mean and standard deviation of saturation than those of the pressure as the log permeability variability increases. At σ 2 Y = 1.0, corresponding to the permeability coefficient of variation being 131%, there is a good agreement between the second-order PCM and the MC approach for the mean and standard deviation of saturation while there exist large deviations for the pressure statistics. It is found that the fourth-order PCM gives better agreements with the MC of 10,000 realizations than does the second-order PCM. Even at σ 2 Y =4.0, corresponding to the permeability coefficient of variation being 732%, the fourth-order PCM results in reasonable approximations for the saturation and pressure statistics. 63 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> σ Y 2 = 0.25 2 nd order PCM 4 th order PCM MC 1000 (a) 0 200 400 600 800 1000 0 0.5 1 1.5 x 10 −3 x, ft σ Sw σ Y 2 = 0.25 2 nd oder PCM 4 th order PCM MC 1000 (b) 0 200 400 600 800 1000 3000 3200 3400 3600 3800 x, ft <Pressure>, psi σ Y 2 = 0.25 2 nd order PCM 4 th order PCM MC 1000 (c) 0 200 400 600 800 1000 0 50 100 150 200 250 x, ft σ Pressure , psi σ Y 2 = 0.25 2 nd oder PCM 4 th order PCM MC 1000 (d) Figure 3.17: Mean and standard deviation of water saturation and pressure obtained from PCM (2nd and 4th order) and MC (with 1000 realizations) for σ 2 Y = 0.25: (a) and (b) for water saturation; (c) and (d) for pressure 64 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (a) 0 200 400 600 800 1000 0 1 2 3 4 5 x 10 −3 x, ft σ Sw σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (b) 0 200 400 600 800 1000 3000 3200 3400 3600 3800 4000 x, ft <Pressure>, psi σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (c) 0 200 400 600 800 1000 0 100 200 300 400 500 600 700 800 x, ft σ Pressure , psi σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (d) Figure 3.18: Mean and standard deviation of water saturation and pressure obtained fromPCM(2ndand4thorder)andMC(with1000and10,000realizations)forσ 2 Y =1.0: (a) and (b) for water saturation; (c) and (d) for pressure 65 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> σ Y 2 = 4.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (a) 0 200 400 600 800 1000 0 0.01 0.02 0.03 0.04 0.05 x, ft σ Sw σ Y 2 = 4.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (b) 0 200 400 600 800 1000 3000 4000 5000 6000 7000 8000 x, ft <Pressure>, psi σ Y 2 = 4.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (c) 0 200 400 600 800 1000 0 2000 4000 6000 8000 10000 x, ft σ Pressure , psi σ Y 2 = 4.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (d) Figure 3.19: Mean and standard deviation of water saturation and pressure obtained fromPCM(2ndand4thorder)andMC(with1000and10,000realizations)forσ 2 Y =4.0: (a) and (b) for water saturation; (c) and (d) for pressure 66 It is seen that the low order PCM approximations require a smaller number of sim- ulations but give reasonable results compared to the high-resolution MC method. One apparent advantage of the MC is flexibility in that it directly utilizes the existing simu- lators and it can be performed with any number of simulations. In practice, only a small numberofrealizationsaresolvedowingtothelargecomputationaleffortrequiredforeach of them. The salient question is whether the MC approach gives consistent and accurate resultsifasmallnumberofsimulationsareperformed. Forthecaseofσ 2 Y =4.0,weevalu- ated the mean and standard deviation of water saturation and pressure using 40 different setsofMonteCarlosimulations,eachwith210realizations. Therefore,thecomputational effort for running each of the MC sets is approximately equivalent to that for the fourth- order PCM because the number of simulations is the same in the two approaches. Figure 3.20 shows that the standard deviation of saturation, the mean and standard deviation of pressure obtained from the Monte Carlo simulations have large spreads. That is, the MC results are highly sensitive to the particular choice of realizations when the number of realizations is small and before statistical convergence is achieved. However, as seen in Figure 3.19, the fourth-order PCM (with the same number of simulation runs as the MC of 210 realizations) yields a single set of results, which are closer to the MC results with 10,000 realizations. The PCM is consistent and robust because the results are fixed for a given polynomial order and a given dimensionality in the KL expansion, unlike the MC in which the results depend on the choice of realizations (which in turn are dependent on the random seed and how the realizations are generated). In order to further investigate the statistical convergence and robustness of the MC results, we plot the mean and standard deviation of water saturation and pressure at a 67 specified location (x = 275 ft) obtained from Monte Carlo simulations versus the number ofrealizationsinFigure3.21. Itshowsthatafterusing1000realizations, themeanvalues start to converge while the standard deviations have fluctuations even with 5,000 realiza- tions. With the direct sampling Monte Carlo simulations, resorting to the law of large numbers, normally a large number of realizations are needed to achieve the statistically accurate results. Both the Monte Carlo method and the PCM involve sampling. The difference is that in the direct sampling Monte Carlo method the realizations are equally probable whereas in the PCM a structural expression (i.e. the polynomial chaos expan- sion) of the output random field is generated at first and then the collocation technique is adopted according to the Gaussian quadrature rule. Owing to the fast convergence of the polynomial chaos expansion and the high accuracy of Gaussian quadrature, the PCM can yield quite robust results with much less samplings than the direct sampling Monte Carlo method. It should be noted that the statistical properties of other flow-related quantities of interest can also be quantified with the PCM. For example, the total (cumulative) oil production from the well can be treated as a random function of time and represented with the polynomial chaos expansion. When the expansion coefficients are obtained with the P sets of simulations, the mean and standard deviation of the total oil production can be determined, from which the confidence interval can be estimated. For the cases of σ 2 Y = 1.0 and 4.0, Figures 3.22 and 3.23, respectively, exhibit the confidence interval of total well oil production as a function of time, obtained from the fourth-order PCM and the MC with 10,000 realizations. The confidence intervals are constructed by the mean plus and minus one standard deviation. The PCM and MC results agree well with 68 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> σ Y 2 = 4.0 40 sets of MC (210) (a) 0 200 400 600 800 1000 0 0.01 0.02 0.03 0.04 0.05 x, ft σ Sw σ Y 2 = 4.0 40 sets of MC (210) (b) 0 200 400 600 800 1000 3000 4000 5000 6000 7000 8000 x, ft <Pressure>, psi σ Y 2 = 4.0 40 sets of MC (210) (c) 0 200 400 600 800 1000 0 2000 4000 6000 8000 10000 x, ft σ Pressure , psi σ Y 2 = 4.0 40 sets of MC (210) (d) Figure 3.20: Mean and standard deviation of water saturation and pressure obtained from 40 sets of MC (with 210 realizations) forσ 2 Y =4.0: (a) and (b) for water saturation; (c) and (d) for pressure 69 1 2000 4000 6000 8000 10000 0.47 0.48 0.49 0.5 0.51 0.52 0.53 Number of realizations <Sw> σ Y 2 = 4.0 (a) 1 2000 4000 6000 8000 10000 0 0.01 0.02 0.03 0.04 0.05 Number of realizations σ Sw σ Y 2 = 4.0 (b) 1 2000 4000 6000 8000 10000 4000 6000 8000 10000 12000 Number of realizations <Pressure>, psi σ Y 2 = 4.0 (c) 1 2000 4000 6000 8000 10000 0 1000 2000 3000 4000 5000 6000 7000 Number of realizations σ Pressure , psi σ Y 2 = 4.0 (d) Figure 3.21: Mean and standard deviation of water saturation and pressure (at x = 275 ft)obtainedfromMonteCarlosimulationsversusthenumberofrealizationsforσ 2 Y =4.0: (a) and (b) for water saturation; (c) and (d) for pressure 70 each other and can be used as the measures of the production uncertainty. Figure 3.22 showsthat whenσ 2 Y =1.0the uncertaintyof the totaloil productionis quite small, while Figure 3.23 indicates a larger uncertainty of the total oil production when σ 2 Y =4.0. 3.4.2.3 Relative Permeability Type To further investigate the effect of relative permeability functions of the fluids on the statistical quantities of the two-phase flow, two more cases with different types of relative permeability functions are performed in this section. Let n=m=2 in Equations (3.43) and (3.44). Consequently, the following Corey type relative permeability functions are obtained with both connate water and residual oil saturation equal to zero: k rw =S 2 w , (3.51) k ro =(1−S w ) 2 . (3.52) With the above relative permeability functions, two cases with random log porosity and permeability, respectively, are conducted. The first one has the same conditions as the previous case shown in Figure 3.14 with the standard deviation of log porosity σ lnϕ = 0.1, correlation length η/L = 2/5, and viscosity ratio m = 2. Figure 3.24 depicts the mean and standard deviation of water saturation and pressure obtained with the second-order PCM and the MC with 1000 realizations. It can be seen that the results are still closely matched between the PCM and MC. It is also observed that the shape of the water saturation profile changes with different relative permeability functions, and that the standard deviation profile of water saturation appears to be more compact and have 71 0 20 40 60 80 100 120 0 2000 4000 6000 8000 10000 12000 t, days Total Well Oil Production, STB σ Y 2 = 1.0 4 th order PCM (a) 0 20 40 60 80 100 120 0 2000 4000 6000 8000 10000 12000 t, days Total Well Oil Production, STB σ Y 2 = 1.0 MC 10,000 (b) Figure 3.22: Confidence intervals of the total well oil production at different time for σ 2 Y =1.0, obtained from the 4 th order PCM and MC (with 10,000 realizations) 72 0 20 40 60 80 100 120 0 2000 4000 6000 8000 10000 12000 t, days Total Well Oil Production, STB σ Y 2 = 4.0 4 th order PCM (a) 0 20 40 60 80 100 120 0 2000 4000 6000 8000 10000 12000 t, days Total Well Oil Production, STB σ Y 2 = 4.0 MC 10,000 (b) Figure 3.23: Confidence intervals of the total well oil production at different time for σ 2 Y =4.0, obtained from the 4 th order PCM and MC (with 10,000 realizations) 73 a higher peak in the water front zone in this case with the quadratic Corey type relative permeability functions compared to the previous case of Figure 3.14 with linear Corey type relative permeability functions. The profile of pressure standard deviation shows a similar effect. TheothercasehasthesameconditionsasthepreviouscaseshowninFigure3.18with the variance of log permeability σ 2 Y = 1.0, correlation length η/L = 2/5, and viscosity ratiom=1. Themeanandstandarddeviationofwatersaturationandpressureobtained from the PCM (with both the second and fourth order) and the MC (with both 1000 and 10,000 realizations) are shown in Figure 3.25. It is found that the standard deviation of pressureissimilartothatinFigure3.18withthelinearrelativepermeabilityfunctionsand thatthepeakvalueofthestandarddeviationofwatersaturationisonemagnitudehigher thanthatinFigure3.18. Thesecond-orderPCMyieldsreasonableresultscomparedtothe high-resolutionMCresultswith10,000realizationswhilethefourth-orderPCMenhances the accuracy. The two cases performed in this section indicate that the type of relative permeabilityfunctionsdoesaffectthetwo-phaseflowresponsesandhencetheirstatistical quantities. 3.4.3 Case Study: 3D Black Oil Model In this section, we apply the PCM to a three dimensional (3D) three-phase model with moderate grid size, modified from the 9 th SPE comparative project. The 9 th SPE project wasoriginallyorganizedtoinvestigatethecomplicationsforblack-oilreservoirsimulation brought about by a high degree of heterogeneity in a geostatistically-based permeability 74 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> m = 2, σ ln ϕ = 0.1 2 nd order PCM MC 1000 (a) 0 200 400 600 800 1000 0 0.05 0.1 0.15 x, ft σ Sw m=2, σ ln ϕ = 0.1 2 nd order PCM MC 1000 (b) 0 200 400 600 800 1000 3000 3200 3400 3600 3800 4000 4200 x, ft <Pressure>, psi m = 2, σ ln ϕ = 0.1 2 nd order PCM MC 1000 (c) 0 200 400 600 800 1000 0 2 4 6 8 10 12 x, ft σ Pressure , psi m = 2, σ ln ϕ = 0.1 2 nd order PCM MC 1000 (d) Figure3.24: Meanandstandarddeviationofwatersaturationandpressureobtainedfrom 2nd order PCM and MC (with 1000 realizations) for the case with standard deviation of logporosityi.e. σ lnϕ =0.1,viscosityratiom=2andnewrelativepermeabilityfunctions: (a) and (b) for water saturation; (c) and (d) for pressure 75 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 x, ft <Sw> σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (a) 0 200 400 600 800 1000 0 0.005 0.01 0.015 x, ft σ Sw σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (b) 0 200 400 600 800 1000 3000 3200 3400 3600 3800 4000 4200 x, ft <Pressure>, psi σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (c) 0 200 400 600 800 1000 0 200 400 600 800 1000 x, ft σ Pressure , psi σ Y 2 = 1.0 2 nd order PCM 4 th order PCM MC 1000 MC 10,000 (d) Figure 3.25: Mean and standard deviation of water saturation and pressure obtained from PCM (2nd and 4th order) and MC (with 1000 and 10,000 realizations) for the case with σ 2 Y = 1.0, viscosity ratio m = 1, and new relative permeability functions: (a) and (b) for water saturation; (c) and (d) for pressure 76 field [Killough, 1995]. The problem involves water flooding with natural water encroach- mentfromanaquiferatthebottomofthedippingreservoir. Themodelhasdimensionof 7200 ft × 7500 ft × 360 ft, which is represented by a 24×25×15 grid. Cell (1,1,1) is at a depth of 9000 ft subsea and the remaining cells dip in the x direction at an angle of 10 degrees. Thereisnodipinthey direction. Thereisonewaterinjectionwellinthecorner grids (i = 24, j = 25), and 25 somehow randomly distributed producers in the reservoir. The initial oil phase pressure is 3600 psia at a depth of 9035 ft subsea, which is the saturation pressure of oil. The oil-water contact is 9950 ft subsea. And there is no initial free gas in the reservoir. Figure 3.26 depicts the configuration of the reservoir with one injection well (I1) and 25 production wells (P2-P26). The initial oil saturation is shown in Figure 3.26. More detailed descriptions of the reservoir can be found in the literature [Killough, 1995]. The 9 th SPE project just considered one realization of heterogeneous permeability field and performed deterministic simulation. In this study, we adopt the same reservoir and fluid descriptions, however impose a random permeability field to the reservoir, and perform stochastic simulations to predict uncertainty in the reservoir. In our simulations, the bottom hole pressure of the injection well is controlled as 4000 psi, and all production wells are controlled as 1000 psi. The vertical permeability is assumed 1% of the horizontal permeability. We assume the log horizontal permeability Y(x) = lnk(x), is a Gaussian random field with the separate exponential covariance function: C Y (x,y)=σ 2 Y exp(−|x 1 −y 1 |/η 1 −|x 2 −y 2 |/η 2 −|x 3 −y 3 |/η 3 ) (3.53) 77 where σ 2 Y is the variance of log permeability, and η 1 , η 2 , η 3 are the correlation lengths in x, y and z direction, respectively. Here we assume the mean of the log horizontal permeability ⟨Y⟩ = 4 [+ln md], vari- ance σ 2 Y =1.0, and the correlation lengths are 1/4 of the domain sizes in each direction. Y(x) is represented by the Karhunen-Loeve (KL) expansion (3.5). Owing to the rapid decayoftheeigenvaluesasshowninFigure3.27,onlythefirst20termsareretainedinthe KL expansion for PCM. Thus the total number of collocation points for the second-order PCM is 231. The formulation for the PCM in 3 D can be given similar to (3.50) as in 1 D, and the implementation is the same. The PCM is performed with the commercial reservoir simulator Eclipse and compared to Monte Carlo (MC) simulations. One realization of the horizontal permeability field is shown in Figure 3.28, indicating the high heterogeneity of the reservoir. Figure 3.29 shows the contour plots of the mean (and standard deviation) of oil saturation at the end of simulation (600 days) obtained from the PCM and MC. And Figure 3.30 is for gas saturation. It reveals that the statis- tics of both oil and gas saturation achieve good agreement between the PCM and MC. The spatial distributions of the statistical properties (such as the mean and standard deviation) can be regarded as important visual measures of the uncertainty rising from the randomness of the heterogeneous permeability. The PCM is more efficient than the MC, since the PCM is involved with 231 independent simulations and the MC requires 1000. Besidestheflow-relatedquantitiessuchasfluidsaturation, theoutputvariablescould be any other reservoir responses generated from the reservoir simulation. For instance, if thefieldoilandgasproductionarethequantitiesofinterest, wecanobtainthestatistical 78 Figure 3.26: Initial oil saturation in the reservoir moments (e.g., mean and variance) of them, based on which the confidence intervals can then be constructed. Figure 3.31(a)-(b) exhibit the confidence intervals (using mean plus or minus the standard deviation) of the cumulative field oil and gas production, respectively, both normalized by the reservoir pore volume. It can be seen that the results obtained from the PCM agree well with those from the MC. It indicates that, for such a three phase model with high nonlinearity and spatial variability, the PCM can quantify uncertainty of the reservoir performances efficiently, with a smaller number of simulations compared to the MC simulations. 79 1 5 10 15 20 25 30 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 n λ n Figure 3.27: Series of eigenvalues in 3 D x y z permx: 4 20 40 60 80 100 200 300 800 md Figure 3.28: One realization of the horizontal permeability field for σ 2 Y =1.0 80 x y z <So>: 0.15 0.25 0.35 0.45 0.55 0.65 0.7 (a) x y z STD_So: 0.001 0.01 0.02 0.03 0.04 0.05 0.06 (b) x y z <So>: 0.15 0.25 0.35 0.45 0.55 0.65 (c) x y z STD_So: 0.001 0.01 0.02 0.03 0.04 0.05 0.06 (d) Figure3.29: Contourplotofmeanandstandarddeviationofoilsaturationobtainedfrom 2nd order PCM and MC (with 1000 realizations): (a) and (b) for PCM; (c) and (d) for MC 81 x y z <Sg>: 0 0.04 0.08 0.12 0.16 (a) x y z STD_Sg: 0 0.01 0.02 0.03 0.04 (b) x y z <Sg>: 0 0.04 0.08 0.12 0.16 (c) x y z STD_Sg: 0 0.01 0.02 0.03 0.04 (d) Figure 3.30: Contour plot of mean and standard deviation of gas saturation obtained from 2nd order PCM and MC (with 1000 realizations): (a) and (b) for PCM; (c) and (d) for MC 82 0 100 200 300 400 500 600 700 0 0.05 0.1 t, days Field CUM. OIL PROD. /PV, STB/RB (a) 0 100 200 300 400 500 600 700 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 t, days Field CUM. Gas PROD./PV, MSCF/RB (b) Figure 3.31: Confidence intervals of reservoir responses forσ 2 Y =1.0: (a) field cumulative oil production (in stock tank barrel STB) normalized by pore volume PV (in reservoir barrel RB); (b) field cumulative gas production (in MSCF) over PV (in RB). Blue for PCM, Red for MC. 83 3.4.4 Summary The PCM is applied to modeling multiphase flow in random heterogeneous media. The underlyingrandomfieldssuchaspermeabilityandporosityarerepresentedbytheKarhunen- Loeve expansion, and the dependent (output) random fields (e.g., fluid saturation and pressure) or variables (e.g., hydrocarbon production) are expressed by polynomial chaos expansions. The PCM is utilized to determine the coefficients in the polynomial chaos expansion and evaluate the statistical properties of the random outputs, by running the deterministic reservoir simulation models independently. The commercial Eclipse reser- voir simulator is incorporated for uncertainty quantification in a straightforward manner. Applications of the PCM to different scenarios of multi-phase flow are performed, in one and three physical dimensions. The one dimensional examples are used to illustrate the approach and study its sensitivity to various factors such as type of variability, variance, correlation length, viscosity ratio, and relative permeability type. A three-dimensional, three-phase reservoir model is modified from the 9 th SPE comparative project for the purpose of reservoir performance forecast. The statistical properties of the flow-related quantities of interest, such as the fluid saturation, pressure, and the cumulative oil (and gas) production, are predicted and compared with those from Monte Carlo simulations. The conclusions are summarized as follows: 1. The PCM is found to be a feasible approach for accurately and efficiently quanti- fying uncertainty associated with reservoir simulations at the field scale. With a smaller computational effort, the PCM can achieve a good agreement with the Monte Carlo results from a large number of realizations. 84 2. The second-order PCM has been found to be consistently efficient in different cases with moderate correlation lengths of the underlying random porosity field and in caseswithmildspatialvariabilityoftherandompermeabilityfield. Whenthevarianceof random log permeability becomes large (e.g, with the coefficient of variation being larger than 253%), a high order of polynomial chaos may be needed to achieve higher accuracy. 3. The statistical moments (e.g., mean and standard deviation) and probability den- sity functions of water saturation and pressure are important measures of uncertainty rising from the randomness of the heterogeneous porosity or permeability. Confidence intervals of the cumulative oil (and gas) production can be constructed to show the un- certainty in production. The viscosity ratio of oil to water and the relative permeability type affect distributions of the mean and standard variation of the water saturation: the spreading of the mean water saturation is larger (indicating larger transition zone) when oil is more viscous than water, and the standard deviation profile of water saturation has ahigherpeakandismorecompactwhenoilislessviscous; withthequadraticCoreytype relative permeability functions in contrast to the linear relative permeability functions, the spreading of mean water saturation is smaller, while the standard deviation profile of water saturation has a higher peak and is more compact. 85 Chapter 4 Non-Gaussian Random Field 4.1 Introduction Simulation of non-Gaussian random fields has gained more and more interest in many engineering problems, especially in the disciplines of engineering mechanics and struc- tural engineering since the last two decades. Non-Gaussian random fields exist in many engineering problems. For example, the Young’s modulus cannot be negative thus its distribution cannot be treated as Gaussian. Many phenomena in the real world exhibit non-Gaussian distributions, such as the elevation of the ocean waves and the acceleration of the seismic ground motions. Permeability of the geological formations is non-Gaussian and it is usually assumed to have log-normal distribution. In this chapter, we consider othernon-Gaussiandistributionsforthepermeability(orhydraulicconductivity). Unlike the Gaussian random field whose joint probability density functions are fully character- ized by the first- and second-order statistics, characterization of non-Gaussian random field requires higher order statistics that are usually difficult to obtain from limited data. 86 Therefore, the marginal probability distribution and the covariance function are often utilized to characterize the non-Gaussian random field in practice. 4.2 Investigation of Approaches for Simulating Non- Gaussian Random Fields 4.2.1 Memoryless Transformation Approach First the memoryless transformation approach is explored for simulating a stationary non-Gaussian random field. A Gaussian field is generated at first, and it is then trans- formed to the non-Gaussian field through a mapping between the marginal probability distributions of the non-Gaussian field and the Gaussian distribution. Let g(x) and n(x) be the stationary Gaussian and non-Gaussian field, respectively. They have the following relationship, n(x)=F −1 (Φ(g(x))) (4.1) whereΦistheGaussiandistributionfunctionandF isthemarginaldistributionfunction of the n(x). The output based on the nonlinear memoryless transformation is called translation process [Grigoriu, 1998]. Yamazaki and Shinozuka [1988] proposed an iterative algorithm for simulating the translation processes with the aid of spectral representation of Gaussian random fields. The marginal CDF and spectral density function of the non-Gaussian random field n(x), i.e., S T nn (κ) are given, where κ is the wave number. A Gaussian field g(x) is sought such that the non-Gaussian field can be generated from the relationship (4.1). 87 At the first step, the spectral density function of g(x), S 1 gg (κ) is set to be S T nn (κ). One realization of the Gaussian field g(x) prescribed with the spectral density function S 1 gg (κ) is generated using the spectral representation method. Then it is transformed to a sample non-Gaussian field n(x) through equation (4.1). However, the spectral density function from the sample field n(x) is generally different from the target, because of the nonlinearity of the transformation (4.1). An iterative procedure is thus introduced to update the spectral density function of g(x), S gg (κ): S i+1 gg (κ)= S i gg (κ) S i nn (κ) S T nn (κ) (4.2) whereS i gg (κ) andS i+1 gg (κ) are the spectral density functions of the Gaussian random field g(x) at thei th and (i+1) th iteration step, respectively, andS i nn (κ) is the spectral density function of the sample random field n(x) at the i th step. S i nn (κ) is calculated from the sample filed n(x) generated at the ith step: S i nn (κ)= 1 2πL ∫ L 0 n(x)exp(−iκx)dx (4.3) Equation (4.2) is used to update until the spectral density function of n(x) coincides with the target S T nn (κ). This iterative algorithm has been applied to generate sample functions of non-Gaussian random fields with slightly skewed marginal probability func- tions, with sufficient accuracy. Deodatis and Micaletti [2001] addressed the difficulties of thealgorithmwhenthemarginalprobabilityfunctionsofthenon-Gaussianrandomfields 88 are highly skewed. The iterative equation (4.2) was modified as follows [e.g., Deodatis and Micaletti, 2001] S i+1 gg (κ)=S i gg (κ) [ S T nn (κ) S i nn (κ) ] β (4.4) where the exponent β is introduced to correct the spectral density functions. Bocchini and Deodatis [2008] critically reviewed the algorithms originating from the Yamazaki and Shinozuka algorithm, and made some improvements on this class of algorithms in terms of accuracy and computational efficiency. However, this class of algorithms is based on the spectral density functions and the spectral representation method is used for simulating Gaussian random fields. It is developed for the purpose of generating realizations of non-Gaussian random fields. If efficiency is taken into consideration for uncertainty quantification, this class of algorithms may not be a good choice since it is limited to Monte Carlo simulations. BasedontheYamazakiandShinozukaalgorithm,Renetal.[1995]proposedamethod for conditional simulation of non-Gaussian fields. In their study, the covariance function instead of the spectral power density of the random fields is specified. The algorithm for generating non-Gaussian random fields is described as follows. Suppose the covariance function of the non-Gaussian field n(x) is given as ρ n (∆x)=E[(n(x 1 )−n(x))(n(x 2 )−n(x))] (4.5) where ∆x = x 2 −x 1 . The covariance function ρ g (∆x) of the underlying Gaussian field g(x)isobtainedbyaniterativeprocedure. Theρ g (∆x)issettobeρ n (∆x)atinitialstep, 89 and realizations of the Gaussian field g(x) can be generated by any existing techniques, then realizations of the non-Gaussian field n(x) can be obtained through equation (4.1). The sampling covariance function of the generated non-Gaussian field is compared to the targetcovariancefunction. Iftheydonotmatch,ρ g (∆x)isupdatedthroughthefollowing equation, ρ g (∆x)=ρ g (∆x)[1+k(ρ (T) n (∆x)−ρ (S) n (∆x))] (4.6) where k is a factor analogous to the relaxation factor in iteration methods of solution of linear equations [Ren et al., 1995]. The iteration procedure proceeds until the target and simulating covariance function of n(x) match within an acceptable error. In order to calculate the sampling covariance ρ (S) n , a large number of realizations of n(x) and thus of g(x) are required at each iteration step. The accuracy of generating the Gaussian random field g(x) will influence that of ρ (S) n and the stability of the iteration. One potential problem of the method is that the positive definite condition of the simulated covariance function of the Gaussian field may violate during the iteration. 4.2.2 Kernel Principal Component Analysis The principal component analysis (PCA) is usually used for simulating Gaussian random field, where only the first and second order statistical moments, i.e., the mean and co- variance are used. PCA is a linear representation of the random field. To account for the nonlinear relationship, the kernel principal component analysis (KPCA) was developed and applied in signal and image processing [Scholkopf et al., 1998]. Sarma et al. [2008] 90 applied the KPCA for parameterization of the non-Gaussian random field in the context of optimization of the oil production. In their study, the KPCA was used to represent the non-Gaussian permeability field, which is estimated in an inverse problem. However, the properties of the non-Gaussian field represented by the KPCA are not verified in the statistical manner. In this study, we will investigate the statistical properties of the random field generated with KPCA. In KPCA, data are first mapped into some feature space via a nonlinear function φ , and then the linear PCA is performed on the mapped data in the feature space. Assume a set of Nr centered observations x k (with zero mean), k = 1,,N r , x k ∈ R N . A kernel matrix is defined as a dot product of two feature expansions, k(x i ,x j )=⟨φ(x i ),φ(x j )⟩ (4.7) The kernel matrix K needs to be normalized to centered matrix, ˜ K =K−JK−KJ +JKJ (4.8) where J is the matrix with each element equal to 1/Nr and Nr is number of realizations. Perform eigenvalue decomposition for the centered matrix, N r λ ˜ K = ˜ Ka (4.9) where λ is the eigenvalue, and a is the eigenvector. 91 ThePCA(orKLexpansion)isperformedinthefeaturespacetogeneraterealizations. A realization in the original data space needs to be inverted from the generated feature expansion, y = φ −1 (Y). This is called the pre-image problem. However, the inverse of the nonlinear projection φ −1 is usually not explicit. A minimization problem is solved to seek for an optimal approximation of the pre-image, minρ(y)=||φ(y)−Y|| 2 (4.10) The critical part of the kernel PCA is the pre-image problem. Pre-image is involved with reconstructing the realization from the feature space. Since it is an inverse process, thepre-imagemaynotexist. Evenifitexists,itmaynotbeunique. Scholkopfetal.[1998] proposed a fix-point iteration scheme to solve the pre-image problem for one particular typeofkernelfunction,i.e. theGaussian(orRBF)kernel,wherethediagonalisconstant. However,adrawbackoftheiterationschemeistheinstability. Rathietal.[2006]proposed an alternative to reconstruct a unique approximate pre-image. The algorithm does not require iteration, while it is modified from the iteration scheme. Applications of KPCA mainly focused on the data extraction (denoising) problems. Our aim is to generate realizations based on a low random dimensionality. Sarma et al. [2008] utilized the polynomial kernels to re-parameterize non-Gaussian fields. The fix-point iteration algorithm for solving pre-images was utilized and modified. We consider some numerical examples to test the applicability of KPCA for simulating non-Gaussian fields. Y(x) is a bimodal random field, whose distribution is shown in Figure 4.1. Suppose we have 1000 realizations of Y(x). It can be seen in Figure 4.2 that 92 −3 −2 −1 0 1 2 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 PDF Y Figure 4.1: The PDF of the non-Gaussian random field Y(x) the reconstruction of the input random field by linear PCA is Gaussian. It indicates that the linear PCA is not applicable for the non-Gaussian random field. WeexaminetheKPCAusingthefix-pointiterationalgorithm. Thepolynomialkernel is used at first. Polynomial kernel function is expressed as, k(x i ,x j )=(x i ·x j ) d (4.11) The polynomial kernel may have a problem that the eigenvalues are very large as the exponent d increases. This may lead to some problems for the implementation of the pre-image problems. The eigenvalues of d = 3 are shown in Figure 4.3(a), and their finite sums are shown in Figure 4.3(b). It can be seen from Figure 4.3(b) that a very large number of terms of eigenvalues need to be retained to preserve enough level of energy. The number of retained eigenvalues determines the random dimensionality 93 −4 −3 −2 −1 0 1 2 3 4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Y PDF Figure 4.2: The PDF of Y(x) generated with Linear PCA of the problem which is related to the computational efficiency for further uncertainty quantification. If the number is too large, the corresponding computational cost for uncertainty quantification will be excessive. The Gaussian kernel is expressed as, k(x i ,x j )=exp ( − ||x i −x j || 2 2σ 2 ) (4.12) For this kernel, we observe that σ needs to be selected properly. If σ is too small, the kernel matrix is singular, and if σ is too large, the kernel matrix will be an identity matrix. How to choose σ is still not clear. Take the previous random field Y(x) as an example. The number of nodes in the one-dimensional domain is 100. The number of realizations of Y(x) is 1000. If σ = 5 in the Gaussian kernel, and the truncated number in the KL expansion is chosen as N = 50, we observe that some of the iterations do 94 1 200 400 600 800 1000 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 n λ n (a) 1 100 200 300 400 500 600 700 800 900 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 n Σλ n /Σλ N (b) Figure 4.3: (a) Eigenvalues for polynomial kernel of order d = 3; (b) finite sums of the eigenvalues for polynomial kernel of order d = 3 95 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Y PDF Figure 4.4: Reconstruction of realizations for Gaussian kernel with σ =5 , N =50. not converge during the fix-point iteration process. The results that fail to converge are eliminated and the remaining realizations are kept. Figure 4.4 shows the PDF of the generated realizations, which behaves like a Gaussian distribution. 4.2.3 Karhunen-Loeve Expansion Insteadoftransformingthenon-GaussianrandomfieldstoGaussianrandomfields,Phoon et al. [2002] proposed an approach to directly simulate the non-Gaussian random fields using the Karhunen-Loeve (KL) expansion. The KL expansion is well known for repre- senting a Gaussian random field through a set of independent Gaussian random variables (See Section 3.2). A random field Y(x) can be expressed as Y(x,θ)=Y(x)+ ∞ ∑ n=1 √ λ n f n (x)ξ n (θ) (4.13) 96 ξ n (θ) can be expressed as ξ n = 1 √ λ n ∫ D (Y(x)−⟨Y(x)⟩)f n (x)dx (4.14) which satisfies, ⟨ξ n (θ))⟩=0 (4.15) ⟨ξ i (θ))ξ j (θ))⟩=δ ij (4.16) Therandomvariablesξ n (θ)intheKLexpansionarenotGaussianifisanon-Gaussian random field, while satisfying equations (4.15) and (4.16). Phoon et al. [2002] proposed an approach to generate realizations of a non-Gaussian random field through an iterative procedure for updating the KL random variables. The algorithm is designed to match the specified covariance function and the marginal cumulative distribution function F. Sample functions of the non-Gaussian random field are generated with the truncated KL expansion, Y(x,θ)=Y(x)+ N ∑ n=1 √ λ n f n (x)ξ n (θ) (4.17) At the first step, ξ n are initially set as identically F-distributed random variables satisfying equations (4.15) and (4.16). Based on the sample functions, the sampling covariance and marginal CDFs can be calculated. Satisfyingequations(4.15)and(4.16), thesamplingcovariancefunctionwillapproach the target covariance function if the sampling number N s is sufficiently large. However, the sampling CDFs calculated from the initially generated sample functions are generally different from the target CDFs. The random variables ξ n need to be modified in order 97 to match the sampling and target CDFs. An iterative equation is established to update the random variables ξ n , ξ k+1 n = 1 √ λ n ∫ D ( Y k (x)− ⟨ Y k (x) ⟩) f n (x)dx (4.18) The iteration is carried out until the sampling CDF matches the target. The non- GaussianKLexpansioncanbeusedtogeneratenon-Gaussianrandomfieldswithdifferent types of marginal PDFs and covariance functions. However, it is limited to Monte Carlo simulation if it is used for the uncertainty quantification. 4.2.4 Independent Component Analysis The independent component analysis (ICA) has recently attracted growing interest and has been widely used in signal processing, neural network research and other disciplines. Itcanbeusedforextractingtheindependentcomponentsfromthemixedsignals. Khhalil andSarkar[2008]addressedthepossibilityofrepresentingnon-GaussianfieldsusingICA. However, the statistics of the generated non-Gaussian fields were not examined by the authors and ICA was not used for uncertainty quantification. A random field Ξ(x) can be represented by the KL expansion Ξ(x)=⟨Ξ(x)⟩+ ∞ ∑ n=1 √ λ n f n (x)ζ n (4.19) where⟨·⟩denotesthemean,λ n andf n (x)arethedeterministiceigenvaluesandeigenfunc- tions of the covariance of Ξ(x), respectively. Owing to the fast decay of the eigenvalues, 98 the expansion can be truncated to finite terms, say, M, which is referred to as the proper random dimension. The uncorrelated random variables ζ n can be formulated as ζ n = 1 √ λ n ∫ D (Ξ(x)−⟨Ξ(x)⟩)f n (x)dx (4.20) The uncorrelated random variables are independent when Ξ(x) is a Gaussian random field. However, for a non-Gaussian random field, the uncorrelated variables ζ n are not necessarily independent and their joint distributions are yet to be found. Presuming mutualindependenceoftherandomvariablescansimplifytheproblemforthegeneration of random fields and the application of many stochastic models. On the other hand, simply assuming the random variables in the KL expansion to be independent can result in loss of information about the joint statistics among these variables. The ICA can be employed to transform the uncorrelated random variables to statis- tically independent ones. The components {ζ n } M n=1 of a random vector are assumed to be linear mixtures of the independent components {ζ n } M n=1 denoted by . The mixing model is represented by ζ =Aη (4.21) where A is the unknown mixing matrix that maps between ζ and η . And then the independent components can be obtained through η =A −1 ζ (4.22) 99 where A −1 is the inverse of A. There are several algorithms for implementing ICA, and the popular FastICA algorithm proposed in [Hyvarinen and Oja, 2000] is adopted in our study. Eachoftheindependentcomponentsη n ,(n=1,,M)canbemappedtoastandard Gaussian random variable {η n } M n=1 via η n =G n (ξ n )=F −1 n (Φ(ξ n )) (4.23) where Φ is the cumulative distribution function (CDF) of the standard Gaussian random variableandF n istheCDFofη n ,whichcanbeobtainedfromtherealizationsofη n . With thetwomappingfunctionsAandG n , theuncorrelatedrandomvector isrepresentedby the independent, standard Gaussian vector via ζ =ζ(η(ξ)), where ξ denotes {ξ n } M n=1 . For a non-Gaussian random field, the KL expansion (4.19) can now be written as Ξ(x)=⟨Ξ(x)⟩+ M ∑ n=1 √ λ n f n (x)ζ n (η(ξ)) (4.24) With the independent, standard Gaussian random vector , the KL expansion can be readily used for generating non-Gaussian random fields and conveniently integrated with stochastic models for uncertainty propagation. It should be noted that the proper random dimensionality (i.e., M in equation (4.24) strongly depends on the stochastic method employed for uncertainty quantification. 4.2.5 Polynomial Chaos Expansion Approach Sakamoto and Ghanem [2002] proposed an approach for simulating non-Gaussian fields basedonthepolynomialchaosexpansion. Thisapproachisapplicableforbothstationary 100 and non-stationary non-Gaussian fields. A random field u(x) can be represented by the polynomial chaos expansion, u(x)= P ∑ i=0 U i (x)Γ i (x) =U 0 (x)+U 1 (x)γ(x)+U 2 (x)(γ(x) 2 −1)+U 3 (x)(γ(x) 3 −3γ(x))+... (4.25) where U i (x) are deterministic coefficients, Γ i (x) are the orthogonal Hermite polynomials with respective to the standard Gaussian random field γ(x), and P is the highest order of the polynomial chaos expansion. Equation (4.25) provides an expression relating the non-Gaussian field to the standard Gaussian field. The coefficients in the equation can be evaluated by the expression, U i = ⟨uΓ i ⟩ ⟨ Γ 2 i ⟩ (4.26) The covariance function of the non-Gaussian field u(x) can be expressed as, ρ u (x 1 ,x 2 )= P ∑ j=1 P ∑ i=1 U i (x 1 )U j (x 2 )⟨Γ i (x 1 )Γ j (x 2 )⟩ (4.27) We have the following expressions, ⟨Γ i (x 1 )Γ j (x 2 )⟩= i!⟨γ(x 1 )γ(x 2 )⟩ i ,i=j 0,i̸=j (4.28) 101 So equation (4.27) becomes, ρ u (x 1 ,x 2 )= P ∑ i=1 U i (x 1 )U i (x 2 )i!⟨γ(x 1 )γ(x 2 )⟩ i (4.29) Equation (4.29) provides the relationship between the covariance of the non-Gaussian field and that of the standard Gaussian random field. The computation of the covariance function of the Gaussian field from equation (4.29) is more stable than the iterative pro- cedureinthememorylesstransformationapproachintroducedforthetranslationprocess. 4.3 Modeling Subsurface Flow in Non-Gaussian Random Fields We have reviewed several approaches for simulating non-Gaussian random fields in the previous sections. However, our main purpose is to perform uncertainty quantification, usually involved with solving stochastic differential equations. In this section, we will perform some applications for uncertainty quantification associated with non-Gaussian random fields. As discussed before, both the ICA and KPCA are lack of mathematical justifica- tion for generating realizations of non-Gaussian random fields. Other approaches have been proven for generating realizations of non-Gaussian random fields with honoring the corresponding statistics. The difference among different approaches for generating ran- dom fields in terms of the computational efficiency can be almost ignored, compared to the computational efforts for solving the governing equations in the further step. The 102 memoryless transformation approaches that originate from the Yamazaki and Shinozuka algorithm make use of the spectral representation of Gaussian random fields. The ap- proach employing the KL expansion for directly generating realizations of non-Gaussian random fields is based on iteration of the non-Gaussian KL random variables. The ICA andPCEcanbeusedforparameterizationofnon-Gaussianrandomfieldsthroughasetof independent Gaussian random variables, which enables them to be integrated into some efficient stochastic approaches for uncertainty quantification, such as the stochastic finite element method and the probabilistic collocation method. The probabilistic collocation method (PCM) is employed in our study. Weconsiderthesteady-stateflowinporousmediaasanexample. Theone-dimensional domain length is L=10. The boundary conditions are prescribed hydraulic head at the two ends, H 0 = 7, and H L = 5. The hydraulic conductivity K(x) is a considered as a non-Gaussian random field, and the statistical moments of the hydraulichead h(x) are to be evaluated. Marginal PDFs of K(x) are prescribed. The covariance function of K(x) is assumed to satisfy the exponential form, C K (x 1 ,x 2 )=σ 2 K exp ( − |x 1 −x 2 | τ ) (4.30) where x 1 and x 2 denote two locations, and τ is the correlation length, here equal to 2. Different types of the marginal PDF for the hydraulic conductivity are considered. In the first case, a skewed non-Gaussian distribution, i.e., the gamma distribution γ(2,1) showninFigure4.5isconsideredtobethemarginalPDFofK(x). AsshowninFigure4.5, the marginal PDF of K(x) can be approximated accurately with the 5th order Hermite 103 0 2 4 6 8 10 12 14 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 K PDF 5th order HPC Gamma Figure 4.5: The marginal PDF of K(x): gamma distribution γ(2,1) with its approxima- tion by the 5th order Hermite polynomial chaos expansion. polynomial chaos expansion. We first generate 1000 realizations of K(x), using the PCE approach. 20 realizations of the non-Gaussian field K(x) are shown in Figure 4.6. Note that the realizations can be from any available tool for generating non-Gaussian random fields. These realizations of the non-Gaussian random field K(x) can be used for Monte Carlo simulation, whose solutions will be used as benchmark for comparing with other approaches. We employ ICA to represent the non-Gaussian random field at first. Once the non- Gaussian field is parameterized by ICA, the PCM can be employed for efficient uncer- tainty quantification. The major steps involved in the proposed approach for uncertainty quantification are summarized as follows: 104 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 x K Figure 4.6: Realizations of the random field K(x). 1) representing the (input) non-Gaussian field with the KL expansion in terms of a set of uncorrelated random variables {ζ i } N i=1 , and obtaining realizations of {ζ i } N i=1 from those of the underlying non-Gaussian field; 2) obtaining two mapping functions: a) the mixing matrix A using the ICA tech- nique,formappingfromtheindependentcomponents{η i } N i=1 totheuncorrelatedrandom variables {ζ i } N i=1 , and b) the function G i for mapping between each η i and a standard Gaussian random variable ξ i ; 3)selectingcollocationpointsfor andrepresentingtheinputrandomfieldaccording to equation (4.24); 4) propagating uncertainty in the underlying physical models through approximation of the dependent (output) random variables (or fields) with Hermite PCE in terms of independent, standard Gaussian random vector ; 105 5) determining the unknown PCE coefficients and evaluating the statistics of outputs with the PCM. The 1000 realizations of input non-Gaussian random field K(x) are used for imple- menting the KL expansion in ICA. N = 6 terms are retained in the KL expansion to perform PCM. The 4th order PCM that requires 210 simulations is performed. Two typesofPCMareexplored, oneemployingICAforextractingnon-Gaussiancomponents, andtheothersimplyassumingtherandomvariablesintheKLexpansiontobeGaussian. MonteCarlo(MC)simulationsusingtheoriginal1000realizationsofK(x)areconducted for comparison. Another MC simulation with 10,000 realizations is also run for examin- ing statistical convergence of the solution. Since the mean values of hydraulic head h(x) have similar results, we just compare the variance of hydraulic head (σ 2 h ) in the whole domain. Figure4.7showstheheadvarianceobtainedfromPCM(withandwithoutICA) and MC (with 1000 and 10,000 realizations). The deviation of the 1000 MC from the 10,000 MC results indicates that even 1000 realizations are not enough in this case to achieveconvergence. Asithasbeencheckedtogiveaconvergentvariance,the10,000MC is used as the benchmark. Note that the computational cost for each PCM simulation is roughly the same as each of the MC simulations. Figure 4.7 also shows that although with much less computational effort, the PCM (with ICA) with 210 simulations gives results identical to those from the 10,000 MC. And the PCM (without ICA) results have some deviations. With the same correlation coefficient for the hydraulic conductivity field, other types of marginal PDFs are considered. Figures 4.8 (a) and 4.9(a) show the PDFs for the beta distributions β(4,2) and β(5,1), respectively. From the simulation results for the 106 0 2 4 6 8 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 x σ h 2 PCM 210 (w/ ICA) MC 1000 MC 10,000 PCM 210 (w/o ICA) Figure4.7: Varianceofhydraulicheadforthecaseofgammadistributionγ(2,1),obtained from PCM and MC for ICA. hydraulic head as shown in Figures 4.8(b) and 4.9(b), the conclusions from the previous case of gamma distribution can also be drawn for the two beta distributions. As discussed before, the PCE approach can be used to represent a non-Gaussian ran- dom field through a set of independent Gaussian random variables. This also enables the use of the PCM for efficient propagation of uncertainty through the flow models. The PCM combined with the KL expansion has been proposed for efficient stochastic simu- lation in the case of Gaussian random field. For a non-Gaussian random field, once the covariance function of the underlying Gaussian field in the PCE approach is calculated, the KL expansion can be used to represent the underlying Gaussian field and the PCM can be incorporated for uncertainty propagation in the flow simulation. Consider the same example where the marginal PDF of K(x) is the gamma distri- bution γ(2,1). We perform PCM and compare it with the MC simulation. Figure 4.10 107 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 K PDF (a) 0 2 4 6 8 10 0 0.005 0.01 0.015 0.02 0.025 x σ h 2 PCM 210 (w/ ICA) MC 1000 MC 10,000 PCM 210 (w/o ICA) (b) Figure 4.8: (a) PDF of Beta distribution ; (b) head variance for the case of Beta distri- bution β(4,2), obtained from PCM and MC for ICA. 108 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 K PDF (a) 0 2 4 6 8 10 0 0.002 0.004 0.006 0.008 0.01 x σ h 2 PCM 210 (w/ICA) MC 10,000 PCM 210 (w/o ICA) (b) Figure 4.9: (a) PDF of Beta distribution ; (b) head variance for the case of Beta distri- bution β(5,1), obtained from PCM and MC for ICA. 109 0 2 4 6 8 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 x σ h 2 PCM 28 PCM 210 MC 1000 MC 10,000 Figure 4.10: Variance of hydraulic head for the case of gamma distribution γ(2,1), ob- tained with the PCM and MC for PCE, shows the variance of hydraulic head obtained with the PCM and MC. As shown in Fig- ure 4.10, the results from MC 1000 and MC 10,000 have some difference, which indicates that 1000 realizations are not enough to achieve the statistical convergence for MC. The 2nd order and 4th order PCM are performed. While the results from the 2nd order PCM (with 28 simulations) have larger deviation from those of MC 10,000, those from the 4th order PCM (with 210 simulations) are close to them. Another gamma distributionγ(3,2) is considered for the marginal PDF ofK(x). The covariance of K(x) has the exponential form with correlation length η = 4. Figure 4.11 (a)showsthemarginalPDFofK(x),whichcanbeapproximatedaccuratelywiththe5th order Hermite polynomial chaos expansion. PCM and MC are performed for uncertainty 110 quantification. As seen in Figure 4.11 (b), the results obtained from the PCM 210 are close to those from the MC 10,000. The above examples demonstrate that the PCM combined with either the ICA or the PCE seems to produce good results compared to the Monte Carlo simulation, with much smaller computational efforts. The range of values for the beta distributions is from 0 to 1, thus the two beta distributions in the previous examples have small variability. We next investigate the effect of larger variability of the non-Gaussian random field. We consider two more non-Gaussian PDFs for the hydraulic conductivity. The following two distributions (named distributions (i) and (ii), respectively) have similar shapes with the beta distributions, β(4,2) and β(5,1), respectively. However, their ranges of values are both from 0 to 10. The PDFs of the two distributions are shown in Figures 4.12 (a) and 4.13(a), respectively. Both the ICA and PCE are used to represent the non-Gaussian field K(x) and they are combined with PCM for uncertainty quantification. The Monte Carlo method with 10,000 simulations is also performed for comparison. The results for distributions (i) and (ii) are presented in Figures 4.12 (b) and 4.13 (b), respectively. For both cases, the PCM with ICA has significant deviations from the MC results, and the PCM with PCE yields results close to the true solutions. We suspect that the ICA is questionable for preserving the statistics of the non-Gaussian random fields especially when they have large variability. The linear relationship imposed in the ICA may not be valid for an arbitrary non-Gaussian field and thus may cause loss of statistics. 111 0 5 10 15 20 25 30 35 40 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 K PDF (a) 0 2 4 6 8 10 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 x σ h 2 PCM 210 MC 10,000 (b) Figure 4.11: (a) PDF of gamma distribution γ(3,2) ; (b) head variance for the case of gamma distribution γ(3,2), obtained from PCM and MC for PCE. 112 (a) 0 2 4 6 8 10 0 0.005 0.01 0.015 0.02 0.025 0.03 x σ h 2 MC 10,000 PCM(ICA) 210 PCM (PCE) 210 (b) Figure 4.12: (a) PDF of distribution (i); (b) head variance from PCM and MC for distribution (i). 113 (a) 0 2 4 6 8 10 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 x σ h 2 MC 10,000 PCM (ICA) 210 PCM (PCE) 210 (b) Figure 4.13: (a) PDF of distribution (ii); (b) head variance from PCM and MC for the case of distribution (ii). 114 4.4 Summary The memoryless transformation methods originating from the Yamazaki and Shinozuka algorithm are widely used for generating non-Gaussian fields. However, they are only ap- plicableforstationaryGaussianrandomfields,andlimitedtotheMonteCarlosimulation foruncertaintyquantification. TheKLexpansioncanbeusedfordirectsimulationofnon- Gaussian random fields by iteratively updating the random variables in the expansion. This approach is limited to Monte Carlo simulation if it is used for uncertainty quantifi- cation. KPCA and ICA are two techniques involved with non-Gaussian distributions. It is questionable for both techniques to represent non-Gaussian random fields by honoring the statistics. The PCE approach is able to simulate both stationary and non-stationary non-Gaussianrandomfields. ThePCEapproachprovidesawaytoparameterizethenon- Gaussian field through a set of independent Gaussian random variables, which makes it possible to integrate it into the PCM for efficient uncertainty quantification. In this ap- proach, the KL expansion is employed for representing an underlying Gaussian random field with dimension reduction. Numerical results reveal that, with a much smaller num- ber of simulations than the Monte Carlo method, the PCM is able to accurately quantify uncertainty for flow associated with non-Gaussian hydraulic conductivity fields. 115 Chapter 5 Conditional Simulation of Random Fields 5.1 Introduction The previous chapters discussed the stochastic simulation without conditioning on mea- sured data. Conditional simulation of the flow and transport in heterogeneous porous media has been widely studied [e.g., Dagan, 1982, Rubin, 1991, Lu et al., 2002, Lu and Zhang, 2004]. When some measurements of the hydraulic conductivity are prescribed in the process of stochastic simulation, it is called the conditional simulation. The condi- tioning on the measured hydraulic conductivity can reduce the overall uncertainty of the hydraulic conductivity field, and that of responses of the flow and transport. In the liter- ature of conditional simulation of flow and transport in porous media, the log hydraulic conductivity is usually a Gaussian random field, and the Monte Carlo simulation and the moment equation approach are the two most common approaches. In this chapter, we discuss the conditional simulation of both Gaussian and non-Gaussian random fields. And the PCM is proposed for efficient uncertainty quantification. 116 5.2 Conditional Simulation of Gaussian Random Fields The conditional simulation of Gaussian random fields is based on the kriging technique [e.g., Ren et al., 1995]. Suppose the Gaussian random field Y(x). The measurements of Y(x) at some locations x i (i=1,2,...N) are known. The objective is to simulate Y(x) at other locations without measurements. Let Y k (x) denote the kriging estimate of Y(x). The simulated field of Y(x) is expressed as follows, Y(x)=Y k (x)+ϵ(x) (5.1) where ϵ(x) is the estimation error at location x. The kriging estimate Y k (x) can be expressed as follows, Y k (x)= N ∑ i=1 λ i (x)Y(x i ) (5.2) where λ i (x) are the Lagrangian multipliers. The kriging technique tries to find the best linear unbiased estimate. The following σ 2 is to be minimized, σ 2 =E{[Y(x)−Y k (x)] 2 } (5.3) Substitution of equation (5.2) into (5.3) yields, σ 2 =C(x,x)−2 N ∑ i=1 λ i (x)C(x i ,x)+ N ∑ i=1 N ∑ j=1 λ i (x)λ j (x)C(x i ,x j ) (5.4) 117 where C(x i ,x j ) is the covariance function of Y(x) between x i and x j . The minimization of σ 2 yields, N ∑ j=1 λ j (x)C(x i ,x j )=C(x i ,x), i=1,2,...,N. (5.5) The kriging weights λ i (x) can be obtained from equation (5.5). And the statistics of ϵ(x) can be derived as follows, E[ϵ(x r )]=0 (5.6) E[ϵ(x r )ϵ(x s )]=C(x r ,x s )− N ∑ i=1 λ i (x s )C(x r ,x i ). (5.7) where x r and x s are any two unrecorded locations. The simulation of Y(x) at any unmeasured locations can be accomplished by cal- culating Y k (x) and simulating ϵ(x). The kriging estimate Y k (x) can be obtained from equation (5.2). Since a Gaussian random field is considered, ϵ(x) is a Gaussian random field whose statistics are shown in equations (5.6) and (5.7). The Karhunen-Loeve (KL) expansion can be used to simulate ϵ(x), ϵ(x)=ϵ(x)+ M ∑ n=1 √ λ n f n (x)ξ n (5.8) where M is the number of terms truncated in the KL expansion. Through the KL expansion, the Y(x) can be eventually represented through a set of independent Gaussian random variables {ξ n }, and the PCM can then be employed for uncertainty quantification. The single-phase flow in porous media considered in section 3.3.2 is used as an ex- ample. The flow configuration is the same as considered in section 3.3.2. The length of 118 the one dimensional domain is L = 10. The fixed hydraulic heads are prescribed at two boundaries, H 0 =7, and H L =10. The only difference from the previous example is that somemeasurementsofhydraulicconductivityaretakenatseverallocationsinthecurrent example. The measurements are taken at three locations x = 0, 3, and 8, respectively. The log hydraulic conductivity, Y(x)=lnK(x), is a Gaussian random field. Y(x) has its covariance with the exponential form, C Y (x 1 ,x 2 )=σ 2 Y exp( −|x 1 −x 2 | η ), where the variance σ 2 Y =1.0, and the correlation length η =4. Using the kriging technique introduced above, the confidence interval of the log hy- draulic conductivity Y(x) can be constructed through the kriging estimate and its the standard deviation. Figure 5.1 shows the confidence interval of the log hydraulic conduc- tivity Y(x). As shown in 5.1, the spatial variability of the random field Y(x) is reduced after conditioning to the measurements. Both the PCM and the MC are performed for uncertainty quantification of the flow responses. Figure 5.2 shows the variance of the hy- draulic head obtained from the PCM and MC. The Monte Carlo method requires 10,000 simulations to achieve converged resulsts, while the PCM only requires 210 simulations. As can be seen in Figure 5.2, the PCM can obtain almost the same solutions as the MC, and it is obviously more efficient than the MC. The effect of the spatial variability is investigated. Two more cases of larger variance of the log hydraulic conductivity σ 2 Y are considered. The σ 2 Y is increased to 2.0, while other conditions are the same as in the previous example. Figure 5.3 (a) shows the confidence interval of the log hydraulic conductivity for σ 2 Y = 2.0, and Figure 5.3 (b) shows the variance of hydraulic head obtained from PCM and MC. The PCM (with 210 simulations) yields results close to the MC (with 10,000 simulations), while there are 119 0 2 4 6 8 10 −2 −1.5 −1 −0.5 0 0.5 1 x Y(x) Figure 5.1: Confidence interval of Y(x) slight deviations. Another case of σ 2 Y = 4.0 is also performed. Figure 5.4 (a) and (b) show the confidence interval of the log hydraulic conductivity and variance of hydraulic head, respectively. For the case of σ 2 Y = 4.0, where the coefficient of variation of the hydraulic conductivity is very large, i.e., Cv K = 732%, the PCM (with 210 simulations) can obtain results close to the MC (with 10,000 simulations), although the deviation also becomes larger. By comparing the examples performed in this section with the previous examplesinsection3.3.2whereconditioningisnotconsidered,onecanseethat,theoverall uncertainty of the hydraulic head is reduced after conditioning to the measurements of hydraulic conductivity. 120 0 2 4 6 8 10 0 0.01 0.02 0.03 0.04 0.05 0.06 x var(h) PCM 210 MC 1000 MC 10,000 Figure 5.2: Variance of hydraulic head, obtained from the PCM and MC 121 0 2 4 6 8 10 −2 −1.5 −1 −0.5 0 0.5 1 1.5 x Y(x) (a) 0 2 4 6 8 10 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 x σ 2 h PCM 210 MC 10,000 (b) Figure 5.3: (a) Confidence interval of Y(x); (b) Variance of hydraulic head, obtained from the PCM and MC for the case of σ 2 Y =2.0. 122 0 2 4 6 8 10 −3 −2 −1 0 1 2 3 x Y(x) (a) 0 2 4 6 8 10 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 x var(h) PCM 210 MC 10,000 (b) Figure 5.4: (a) Confidence interval of Y(x); (b) Variance of hydraulic head, obtained from the PCM and MC for the case of σ 2 Y =4.0. 123 5.3 Conditional Simulation of Non-Gaussian Random Fields For conditional simulation of a non-Gaussian random field, Ren et al. [1995] utilized the technique for conditional simulation of Gaussian random field, together with the method for non-conditional simulation of non-Gaussian random field. In their work, the non- conditional simulation of the non-Gaussian random field K(x) is based on the mapping between the non-Gaussian random field to an Gaussian random field, K(x)=F −1 {Φ[g(x)]} (5.9) where g(x) is the underlying Gaussian random field, and F and Φ are the cumulative distribution functions of the K(x) and g(x), respectively. However, this approach is limited to the Monte Carlo simulation. In this study, we employ the Hermite polynomial chaos expansion for simulation of non-Gaussian random field, as introduced in section 4.2.5. K(x) is expressed with the Hermite polynomial chaos expansion, K(x)= P ∑ i=0 U i (x)Γ i (x) =U 0 (x)+U 1 (x)g(x)+U 2 (x)(g(x) 2 −1)+U 3 (x)(g(x) 3 −3g(x))+... (5.10) where U i (x) are deterministic coefficients, Γ i (x) are the orthogonal Hermite polynomials with respective to the standard Gaussian random field g(x), and P is the highest order 124 of the polynomial chaos expansion. The coefficients can be evaluated by the following equation, U i = ⟨KΓ i ⟩ ⟨ Γ 2 i ⟩ (5.11) The covariance function of K(x) is given as ρ K (x 1 ,x 2 ), and it can be expressed by the covariance of g(x) denoted by ρ g (x 1 ,x 2 ), ρ K (x 1 ,x 2 )= P ∑ i=1 U i (x 1 )U i (x 2 )i!ρ g (x 1 ,x 2 ) i (5.12) ρ g (x 1 ,x 2 ) then can be solved from equation (5.12). Onceρ g (x 1 ,x 2 ) is known, and the conditionalvaluesofK(x)aretransferredtog(x)basedontheequation(5.9),conditional simulation of the underlying Gaussian random field g(x) can be accomplished through the approaches introduced in the previous section. Either the PCM or the MC can be employed for uncertainty quantification. The flow scenario considered in the previous section is used for application. The flow configuration is the same except that the hydraulic conductivity is a non-Gaussian randomfield. Supposethegammadistributionγ(2,1)showninFigure5.5isthemarginal PDF of the hydraulic conductivity K(x). The covariance of K(x) is given as C K (x 1 ,x 2 ) = σ 2 K exp( −|x 1 −x 2 | η ) where σ 2 K is the variance of K(x), and the correlation length η = 4. There are 5 measurements taken at x = 0, 1, 5, 7, and 10. With the approach introduced above, the confidence interval of the hydraulic conductivity can be constructed, as shown in Figure 5.6. Both the PCM and MC are performed for uncertainty quantification. Figure 5.7 shows the variance of hydraulic head after conditioning. It can be seen that the PCM (with 210 simulations) 125 0 5 10 15 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 K PDF Figure 5.5: The marginal PDF of the hydraulic conductivity K(x) cangetresultsclosetotheMC(with10,000simulations),whiletherearesomedeviations. Another case of groundwater flow in a two-dimensional domain is performed. The domainisasquareofsizeL 1 =L 2 =10,uniformlydiscretizedinto20×20squareelements. The non-flow conditions are prescribed at two lateral boundaries. The hydraulic head is prescribed at the left and right boundaries as 10.5 and 10.0, respectively. Suppose the gamma distribution γ(2,1) shown in Figure 5.5 is the marginal PDF of the hydraulic conductivity K(x). Assume the covariance function of the hydraulic conductivity has the form C K (x,y) = σ 2 K exp(−|x 1 −y 1 |/η 1 −|x 2 −y 2 |/η 2 ), where σ 2 K is the variance of K(x), and the correlation lengths are η 1 = η 2 = 5.0. There are 9 conditional points regularly distributed in the square domain. The mean and standard deviation of the hydraulic conductivity after conditional simulation are shown in Figure 5.8 (a) and (b), respectively. The overall uncertainty of the hydraulic conductivity field is reduced. Both 126 0 2 4 6 8 10 1 2 3 4 5 6 7 8 9 10 11 x K(x) Figure 5.6: The confidence interval of the hydraulic conductivity K(x) 0 2 4 6 8 10 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 x σ h 2 PCM (210) MC (10,000) Figure 5.7: Variance of the hydraulic head after conditioning 127 thePCMandMCareperformedforuncertaintyquantificationofthehydraulichead. The meanandthestandarddeviationofthehydraulicheadobtainedwiththePCM(with210 simulations) are shown in Figure 5.9 (a) and (b), respectively. And those obtained from the MC (with 10,000 simulations) are presented in Figure 5.10 (a) and (b), respectively. It can be seen that, both the mean and the standard deviation of the PCM are in good agreement with those from the MC, although the central strip of high standard deviation obtained from the MC is slightly wider than that from the PCM. 5.4 Summary This chapter discusses the conditional simulations of flow in heterogeneous porous media when measurements of hydraulic conductivity are available. When the hydraulic conduc- tivityischaracterizedasalog-normalrandomfield, themethodisbasedonthealgorithm ofthekrigingtechniqueforconditionalsimulationofGaussianrandomfields. TheKLex- pansionisintegratedintothealgorithmforrepresentingtheunderlyingGaussianrandom field, and then the PCM is proposed for uncertainty quantification. Different degrees of the variance of log hydraulic conductivity are investigated. When the hydraulic conduc- tivity is characterized as a non-Gaussian random field other than the typical log-normal field, thepolynomialchaosexpansionisusedtorepresentthenon-Gaussianrandomfield, and the technique for conditional simulation of Gaussian random fields is incorporated. ThePCMandMCareperformedandcomparedforconditionalsimulationsofflow,where thehydraulicconductivityischaracterizedwithboththeGaussianandnon-Gaussianran- dom fields. Results show that the PCM is computationally more efficient than the MC. 128 x y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 2 3 4 5 6 7 8 9 10 11 (a) x y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 (b) Figure 5.8: (a) The mean of the hydraulic conductivity K(x) after conditioning; (b) The standard deviation of K(x) after conditioning. 129 x y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 10 10.05 10.1 10.15 10.2 10.25 10.3 10.35 10.4 10.45 10.5 (a) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 (b) Figure 5.9: (a) The mean of the hydraulic head after conditioning, obtained with PCM (with 210 simulations); (b) The standard deviation of the hydraulic head after condition- ing, obtained with PCM (with 210 simulations). 130 x y 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 10 10.05 10.1 10.15 10.2 10.25 10.3 10.35 10.4 10.45 10.5 (a) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 (b) Figure 5.10: (a) The mean of the hydraulic head after conditioning, obtained with MC (with 10,000 simulations); (b) The standard deviation of the hydraulic head after condi- tioning, obtained with MC (with 10,000 simulations). 131 The conditioning on the measurements of the hydraulic conductivity reduces the overall uncertainty of the hydraulic conductivity field and the hydraulic head. 132 Chapter 6 Arbitrary Random Variables 6.1 Introduction Inthepreviouschapters,thepropertiesofthegeologicalformations,suchasthehydraulic conductivity(orpermeability)andporosity,areconsideredasrandomfields. However,for many real-world applications of uncertainty assessment in the oil industry, the uncertain parameters associated with the reservoir simulation are usually considered as random variables. Those parameters could be the permeability and porosity of the reservoir, the connate water saturation, water-oil contact, gas-oil contact, and so on. When treated as random variables, they can have different types of probability distributions. In order to perform accurate uncertainty analysis, a large number of simulations are often required, as in the well-known Monte Carlo method. However, it is usually prohibitive to do so sinceevenasinglesimulationofpracticallarge-scalesimulationmodelsmaybequitetime consuming. The experimental design method associated with different response surface method- ologiesiswidelyusedintheoilindustryforassessinguncertaintiesinreservoirproduction 133 andeconomicappraisal[e.g.,VanElketal.,2000,Venkataraman,2000,Friedmannetal., 2003, Yeten et al., 2005, Tipping et al., 2008]. Experimental design provides a proxy of reservoir models and it has been used to speed up the history matching process [Amudo et al., 2008, Schaaf et al., 2008]. The experimental design method is more efficient than the direct Monte Carlo method. Most experimental design methods use the two-level design (with low and high values), or three-level design (including the middle value). Multilevel designs increase the number of simulations significantly and some mixed de- signs are more efficient [Kalla and White, 2007]. However, a key disadvantage of the experimental design approach is that it does not take into account the full probability distributions of the parameters consistently while creating the response surface, that is, thefullPDFsarenotusedforsamplinganddesign,butonlyusedduringpost-processing. Further, because all samples are equally weighted for response surface generation, there is an inherent assumption that the distributions of these parameters are uniform. As a result, the experimental design method may not be appropriate when parameter distri- butions are arbitrary, which is common in real-world applications. However, it is still broadly used in the industry and these limitations are usually disregarded. Recently, Lawal [2009] discussed some limitations of experimental design methods for modeling uncertainty in the oil industry, such as the inconsistency and non-uniqueness of proxy models, and the impact of input distributions on the output, etc. In this study, we apply the probabilistic collocation method (PCM) for uncertainty quantification, and compare it with the experimental design methods. We propose an approach to deal with arbitrary input PDFs with the PCM. Orthogonal polynomials for arbitrary PDFs are constructed numerically, and then the probabilistic collocation 134 method is utilized for uncertainty propagation. As a result, PCM can be applied effi- ciently for any arbitrary numerical or analytical distribution of the input parameters. The approach is also extended to discrete distributions. 6.2 Generalized Polynomial Chaos Expansion As introduced before, the original polynomial chaos expansion is based on the Hermite polynomials in terms of Gaussian random variables. The generalized polynomials chaos expansion,basedondifferentpolynomialsintheAskeyschemefordifferenttypesofinput random variables, has been developed by Xiu and Karniadakis [2002]. Not restricted to theGaussiandistribution,theinputrandomvariablescanbeanyoftheclassicprobability distributions in the Askey scheme. Different input distributions determine the types of orthogonal polynomials used. For example, Legendre polynomial chaos expansion is used for the uniform distribution, Laguerre polynomial chaos expansion for the gamma distribution,Jacobipolynomialchaosexpansionforthebetadistribution,etc. Foragiven input distribution, the specified polynomial chaos expansion yields better convergence rate than other expansions. In reservoir models, the probability distributions of input random variables could be arbitrary. For such distributions, the associated orthogonal polynomial expansions need to be constructed. For a random variable y of interest, the polynomial chaos expansion can be expressed in the following general form: y = ∑ k≤d a k Γ k () (6.1) 135 where d is the highest order of the expansion, a k are the deterministic coefficients and Γ k () are the k th order multi-dimensional polynomials of the random variables = (ξ 1 ,ξ 2 ,···,ξ N ) T . In this study, we construct the orthogonal polynomial chaos expansion numerically with respect to arbitrary probability density functions. 6.3 Orthogonal Polynomials for Arbitrary Random Variables In this section, we briefly describe the construction of orthogonal polynomials for arbi- trary random variables, including discrete and discretized continuous random variables. The orthogonality is usually with respect to a positive weight function on an interval or some discrete positive weights on specified points. The orthogonal polynomials can be constructed numerically [Gautschi, 1994]. For a weight function ω(t), the inner product of polynomials u(t) and v(t) is defined as, (u,v)= ∫ R u(t)v(t)ω(t)dt (6.2) The i th order orthogonal polynomials {γ i (t)} are defined as follows, γ i (t)=t i + lower-order terms, i=0,1,2,... (6.3) (γ i ,γ j )=0, if i̸=j (6.4) 136 They have the distinct feature, that they satisfy the three-term recurrence relations, γ k+1 (t)=(t−α k )γ k (t)−β k γ k−1 (t), (6.5) γ 0 (t)=1, γ −1 (t)=0, (6.6) StieltjesprocedureortheLanczosalgorithmisusuallyusedtocomputetherecurrence coefficients α i and β i . The Stieltjes procedure uses the fact that the coefficients can be expressed as follows, α i = (tγ i ,γ i ) (γ i ,γ i ) , i=0,1,2,... (6.7) β 0 =(γ 0 ,γ 0 ) (6.8) β i = (γ i ,γ i ) (γ i−1 ,γ i−1 ) , i=1,2,... (6.9) In the case of discrete weights, the inner product (6.2) has the discretized form: (u,v)= M ∑ k−1 w k u(x k )v(x k ) (6.10) where x k are M distinct points and w k are the corresponding weights. Let = (ξ 1 ,ξ 2 ,· · ·,ξ N ) T be a random vector consisting of the random variables. Once the one-dimensional orthogonal polynomials γ k i (ξ i ) are constructed, the multi- dimensional polynomial expansion Γ k (ξ) can be obtained by tensor product of the one- dimensional polynomials, Γ k ()=Π N i=1 γ k i (ξ i ) (6.11) 137 where k = ∑ N i=1 k i . 6.4 PCM Compared with Experimental Design Based on the constructed orthogonal polynomials, the PCM can be employed for un- certainty quantification. In this section, the probabilistic collocation method is applied to three reservoir models, and the results are compared to experimental design methods and Monte Carlo simulations. Through the examples, application of PCM to continuous, discrete and mixed distributions is demonstrated. 6.4.1 Case Study: Model 1 The first reservoir model is a synthetic model used for comparing different experimental design methods and response surfaces [Yeten et al., 2005]. It has 5×10×3 = 150 active cells, with a bottom drive aquifer and a gas cap. There are 4 vertical producers and 1 horizontal water injector. The simulation period is 11 years. There are 7 uncertain parametersinthismodel, asshowninTable6.1. Theparametersareassumedasuniform random variables. We compare the PDFs of the cumulative field oil production obtained from the PCM and different experimental design methods. Figure 6.1 shows the PDF of the cumulative oil production at the end of simulation time, obtained with the 2nd order PCM (with 36 simulations). The Monte Carlo simula- tions (MC) with 1250 simulations are performed for comparison. As shown in Figure 6.1, the PCM with only 36 simulations agrees with the MC with 1250 simulations quite well. Next, we perform experimental designs for this model, and compare them to MC. Three- level design is used. There are several different experimental design methods, such as the 138 Parameter ID Parameter description Min Max SWC Connate water sat., fraction 0.15 0.25 WOC Water oil contact depth, ft 8650 8750 GOC Gas oil contact depth, ft 7650 7750 RCOMP Rock compressibility, 1/psi 8×10 −7 8×10 −6 PORO Global porosity, fraction 0.15 0.25 XPERM Global k x , md 50 550 ANISMLT Global k x /k z multiplier 0.01 1 Table 6.1: The description of the random parameters for model 1 D-optimal design, space filling design, exhaustive design. Six types of response surfaces, i.e., the linear polynomial, pure quadratic polynomial, full quadratic polynomial, neural network,splines,andkrigingresponsesurfacesaregeneratedforeachdesignmethod. The PDFs are evaluated by performing Monte Carlo sampling on the constructed response surfaces. In thisstudy, thenumberforthe MonteCarlosampling forallthe experimental design is chosen as 10,000, which has been proven to be sufficient to achieve statistical convergence. Figure 6.2 (a) shows the PDFs of the cumulative field oil production using the D-optimal design with 36 simulations, the same number used for the PCM. All the response surfaces are presented together. As can be seen in Figure 6.2 (a), the pure quadratic polynomial response surface agrees with the MC very well, the splines and kriging are close to the MC, and other response surfaces have larger deviations. Figure 6.2 (b) shows the PDFs of the cumulative field oil production using the space filling de- sign (with 36 simulations) with various response surfaces. Figure 6.2 (b) indicates that the full quadratic polynomial response surface for the design has extraordinarily large discrepancy, though other response surfaces are close to the MC. Figure 6.2 (c) shows the results for the exhaustive design that uses all the design points. 2187 simulations are 139 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 6 0 0.5 1 1.5 2 2.5 x 10 −6 Cum. Oil Prod., STB PDF PCM 36 MC 1250 Figure 6.1: PDFs of cumulative oil production, obtained from PCM and MC for case 1 requiredfortheexhaustivedesign. AsshowninFigure6.2(c), allotherresponsesurfaces are close to the MC results, except for the linear polynomial response surface. We learn from this case that, for a particular problem, different experimental design methods or different response surfaces for a particular design method can yield different results. Itindicatesthattheexperimentaldesignmethodsarenotrobust. Althoughsome combinationsofadesignmethodandaresponsesurfacemaygivegoodresults,othersmay have bad solutions. In contrast, the PCM gives robust results, since the PCM is based on a particular polynomial chaos expansion according to the probability distributions of the random parameters. In this case, the Legendre polynomial chaos expansion is used for the uniform random parameters. 140 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 6 0 0.5 1 1.5 2 2.5 x 10 −6 Cum. Oil Prod, STB PDF linear Pure qaud Full quad Neural Net Spline Kriging PCM MCS (a) D-optmial design (36 simulations) −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 6 0 0.5 1 1.5 2 2.5 x 10 −6 Cum. Oil Prod., STB PDF Linear Pure quad Full quad Neural net Spline Kriging PCM MCS (b) Space-filling design (36 simulations) 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 6 0 0.5 1 1.5 2 2.5 x 10 −6 Cum. Oil Prod. STB PDF Linear Pure quad Full quad Neural net Spline Kriging PCM MCS (c) Exhaustive design (2187 simulations) Figure 6.2: PDFs of cumulative oil production for case 1 of model 1, obtained from different experimental designs 141 For the same reservoir model, we consider more complicated probability distributions for the input random parameters. Let the random inputs have the multi-modal distri- bution shown in Figure 6.3. The PDF can be also evaluated from sampling, as shown in Figure 6.4. Figure 6.5 shows the comparison of PDFs of PCM (using both explicit PDF and sampling PDF) and MC. The results of PCM using explicit PDF and sampling PDF are almost the same, and close to those of MC. Figure 6.6(a)-(c) show the compar- isonsofPDFobtainedfromdifferentexperimentaldesignschemeswithdifferentresponse surfaces, i.e., linear, pure quadratic and full quadratic polynomials, respectively. Both linearandpurequadraticresponsesurfacesdeviatefromMCresults, evenfortheexhaus- tive design with 2817 simulations. For the full quadratic polynomial response surfaces, D-optimal design with 125 simulations and exhaustive design are close to PCM with 36 simulationsandMCwith2000simulations,whereasboththeD-optimalandSpaceFilling design with 36 simulations run have large deviations. 6.4.2 Case Study: Model 2 The second reservoir model is modified from a field model used for comparing different experimental design methods and response surfaces [Yeten et al. 2005]. It is a sector model of a real compartmentalized oil reservoir, with 1 water injector and 2 producers. The model has 27×46×36 = 44,712 cells, of which 19,922 are active. The reservoir configuration is shown in Figure 6.7. There are a number of faults in the system, the transmissibility multipliers of two of which are uncertain. Each compartment is modeled as different equilibrium regions. Reservoir rock properties are populated using 3 facies distributions. Total simulated time period is 10 years. Five reservoir parameters are 142 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 PDF Figure 6.3: The PDF of the multi-modal distribution 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 1.2 1.4 1.6 PDF Figure 6.4: The PDF of the multi-modal distribution from sampling 143 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 6 0 0.5 1 1.5 2 2.5 x 10 −6 Cum. oil prod., STB PDF PCM 36 (Explicit PDF) PCM 36 (Sampling PDF) MC 2000 Figure 6.5: PDF of cumulative oil production for the case of multi-modal distribution, obtained from PCM and MC 144 (a) (b) (c) Figure 6.6: Comparison of PDF of cumulative oil production for the case of multi-modal distribution, obtained from experimental designs, PCM and MC. For experimental de- signs: (a) liner; (b) pure quadratic; (c) full quadratic polynomial response surface 145 Parameter ID Parameter description Min Max PORV1 Pore vol. multiplier for equilibrium region 1 1 10 PORV2 Pore vol. multiplier for equilibrium region 2 and 6 0.5 1.5 TRANSMULT2 Transmissibility multiplier for fault 2 0.0 10.0 TRANSMULT4 Transmissibility multiplier for fault 4 0.0 10.0 ANISMLT kx/kz multiplier for facies 3 0.01 5.0 Table 6.2: The description of the random parameters for model 2 Figure 6.7: The configuration of model 2 considered uncertain, as shown in Table 6.2. Three cases are considered for the same model: case 1 for all parameters with the same random distributions, and case 2 and 3 for parameters mixed with continuous and discrete distributions. In the first case, all 5 parameters are assumed to have uniform distributions, with the ranges shown in Table 6.2. The aim of the study is to quantify uncertainty in the cumulative field oil production. In the following cases, the numbers of MC simulations have been shown to be sufficient, unless specified otherwise. 146 (a) The 2nd order (21 simulations) and 4th order (126 simulations) PCM 3 3.5 4 4.5 x 10 6 0 1 2 3 4 5 6 x 10 −6 Cum. Oil Prod., STB PDF PCM 31 (6th pure) MC 1422 (b) The 6th pure order PCM (31 simulations) Figure 6.8: PDFs of cumulative oil production, obtained from PCM and MC 147 Figure 6.8 shows the PDFs of the cumulative field oil production at the end of simu- lation time, obtained from PCM (with different orders) and MC (with 1422 simulations). As shown in Figure 6.8 (a), the 2nd order PCM requires 21 simulations and shows some deviation from MC, and the 4th order PCM that requires 126 simulations agrees with MC quite well. The number of total simulations in PCM increases fast as the order of polynomial chaos expansion increases. To increase efficiency, we may try to keep some higher order terms, and eliminate the crossing terms in the polynomial chaos expansion. For example, if we increase the order up to 6 but remove cross terms, only 31 simulations are required. The PCM associated with this reduced order polynomial chaos expansion is referred to as a pure PCM. Figure 6.8 (b) shows that the pure PCM matches MC very well. Notethatitiscurrentlynotclearhowtobestremovetermsfromthefullpolynomial expansion such that accuracy is least affected and maximum efficiency is gained. This will be investigated in future work. Figure 6.9 (a) shows the PDFs of the cumulative field oil production using the D- optimal design with 21 simulations. Results for all response surfaces based on the same design are shown together. As seen from Figure 6.9 (a), different response surfaces give different results, and none of them matches MC. When the number of simulations for the D-optimal design is increased to 100, Figure 6.9 (b) shows that the PDFs perform better for some response surfaces than the D-optimal design with 21 simulations, although they still deviate from MC results and the linear polynomial response surface is far away from MC. Figure 6.9 (c) shows the results for the space-filling design with 31 simulations, and Figure 6.9 (d) is for the exhaustive design that uses the complete set of design points 148 Parameters \ Probability 0.3 0.4 0.3 TRANSMULT2 (transmissibility multiplier for fault 2) 0 5.0 10.0 TRANSMULT4 (transmissibility multiplier for fault 4) 0 5.0 10.0 Table 6.3: Probabilities of the categorical variables for case 2 of model 2 (243) and requires 243 simulations. None of the experimental designs agrees with MC results. Incase2,thetransmissibilitymultipliersforfault2and4areassumedtobecategorical variables with three possible discrete values whose values and probabilities are shown in Table 6.3. Other parameters are uniform and their ranges are shown in Table 6.2. Figure6.10showsthediscreteprobabilityfunctionofthetransmissibilitymultipliers. For implementingPCMonmixedrandomdistributions, polynomialchaosexpansionsneedto beconstructedonthebasisofthespecificdistributionofeachrandomvariable. ThePCM with different order are performed and compared with MC (with 1323 simulations). The PDFsofcumulativeoilproductionobtainedwithPCMandMCareshowninFigure6.11. As seen from Figure 6.11, the 2nd order PCM with 21 simulations has some deviation from MC, but the 6th order (pure) PCM with only 25 simulations agree with MC quite well. Two different experimental designs are performed, i.e., the popular D-optimal and spacefillingdesign. 25reservoirsimulationsareconductedforbothdesigns, anddifferent response surfaces are used to analyze the cumulative oil production. The results are shown in Figure 6.12. For each design, different response surfaces exhibit different results and deviate from MC results, except that the neural network response surface for the space filling design yields better results. 149 3 3.5 4 4.5 5 x 10 6 0 1 2 3 4 5 6 7 8 x 10 −6 Cum. Oil Prod., STB PDF linear Pure quad Full quad Neural net Spline Kriging MC 1422 (a) D-optimal design (21 simulations) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 x 10 6 0 1 2 3 4 5 6 x 10 −6 Cum. Oil Prod, STB PDF linear Pure quad Full quad Neural net Spline Kriging MC 1422 (b) D-optimal design (100 simulations) 3 3.5 4 4.5 5 x 10 6 0 1 2 3 4 5 6 x 10 −6 Cum. Oil Prod., STB PDF linear Pure quad Full quad Neural net Spline Kriging MC 1422 (c) Space-filling design (31 simulation) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 x 10 6 0 1 2 3 4 5 6 x 10 −6 Cum. Oil Prod., STB PDF linear Pure quad Full quad Neural net Spline Kriging MC 1422 (d) Exhaustive design (243 simulations) Figure 6.9: PDFs of cumulative oil production, obtained from different experimental designs 150 Figure 6.10: Probability function of transmissibility multiplier 2 and 4. 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 x 10 6 0 1 2 3 4 5 6 x 10 −6 Cum. Oil Prod., STB PDF PCM 21 (2nd full) PCM 25 (6th pure) MC 1323 Figure 6.11: PDFs of cumulative oil production, obtained from PCM and MC 151 3 3.5 4 4.5 x 10 6 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 −6 Cum. Oil Prod., STB PDF Linear Pure Quad Full Quad Neural Net Splines Kriging MC 1323 (a) D-optimal design (25 simulations) 3 3.5 4 4.5 5 5.5 x 10 6 0 1 2 3 4 5 x 10 −6 Cum. Oil Prod., STB PDF Linear Pure Quad Full Quad Neural Net Splines Kriging MC 1323 (b) Space filling design (25 simulations) Figure 6.12: PDFs of cumulative oil production, obtained from experimental designs (with 25 simulations) 152 6.4.3 Case Study: Model 3 The third reservoir model is a large field model with grid size 96×73×150=1,051,200. The reservoir configuration is shown in Figure 6.13. The simulation time period is 4 years. There are 5 random parameters as shown in Table . All parameters are assumed to have uniform distributions. The quantity of interest is the PDF of the cumulative field oil production. MC is still performed as benchmark for comparison. 863 simulations are performed for MC and converged results are obtained. Note that the MC simulation is time consuming for such a large model. Each simulation takes about 40 min, thus MC with 863 simulations requires about 575 hours of CPU time. Different experimental designsusingvariousresponsesurfacesareperformedandtheresultsareshowninFigure 6.14. As can be seen from Figure 6.14 (a), all the response surfaces for the Folded Plackett-Burman (FPB) design yield results far away from the MC results. The FPB designrequires16simulations, andFigure6.14(a)indicatesthatthe16designpointsare not enough to capture the model uncertainty. We increase the number of design points to 30, and use two different design methods, i.e. the D-optimal design and the space filling design, whose results are shown in Figures 6.14 (b) and (c), respectively. Both the D-optimal and space filling design have deviations for different response surfaces, while the space filling design is closer to the MC results. Figure 6.14 (d) is for the exhaustive design which requires 243 simulations. It is noted that the exhaustive design that uses the full combination of design points does not guarantee good results, although it usually requiresalargenumberofsimulations. As showninFigure6.14(d), none oftheresponse surfaces for the exhaustive design gives satisfactory results compared to the MC solution. 153 Parameter ID Parameter description Min Max SWC Connate water saturation 0.08 0.18 RCOMP Rock compressibility, 1/psi 1.1×10 −6 1.1×10 −5 XPERMMULT k x multiplier 0.5 3 YPERMMULT k y multiplier 0.5 3 ZPERMMULT k z multiplier 0.5 3 Table 6.4: The description of random parameters for model 3 Figure 6.13: The configuration of model 3 Figure 6.15 shows the comparison of PCM and MC. The 6th order pure PCM with 31 simulations can get results quite close to the MC solution. 6.5 Summary A stochastic uncertainty quantification approach to deal with arbitrary input PDFs is proposed with the aid of probabilistic collocation method (PCM) and the construction of numerical orthogonal polynomials. Further, PCM and experimental design methods 154 5 6 7 8 9 10 11 12 13 x 10 7 0 1 2 3 4 5 6 7 x 10 −8 Cum. Oil Prod., STB PDF Linear Pure quad Neural net Splines Kriging MC 863 (a) FPB design (16) 5 6 7 8 9 10 11 12 13 x 10 7 0 1 2 3 4 5 6 7 x 10 −8 Cum. Oil Prod., STB PDF Linear Pure quad Full quad Neural net Splines Kriging MC 863 (b) D-optimal design (30) 5 6 7 8 9 10 11 12 13 14 x 10 7 0 1 2 3 4 5 6 7 x 10 −8 Cum. Oil Prod., STB PDF Linear Pure quad Full quad Neural net Splines Kriging MC 863 (c) Space filing design (30) 5 6 7 8 9 10 11 12 13 x 10 7 0 1 2 3 4 5 6 7 x 10 −8 Cum. Oil Prod., STB PDF Linear Pure quad Full quad Neural net Splines Kriging MC 863 (d) Exhaustive design (243) Figure 6.14: PDFs of cumulative oil production obtained with different experimental designs 155 5 6 7 8 9 10 11 12 13 x 10 7 0 1 2 3 4 5 6 7 x 10 −8 Cum. Oil Prod., STB PDF PCM 31 (6th pure) PCM 41 (8th pure) MC 863 Figure 6.15: PDFs of cumulative oil production obtained with PCM are compared for a number of cases. The applicability, accuracy, and efficiency of the methods are explored. From the cases studied, the following conclusions can be drawn: 1. Mostexperimentaldesignmethodsutilizetheminimum,maximumplusthemiddle valueoftheunderlyingparameterstodesignexperimentsandconstructresponsesurfaces based on the experiments. A key disadvantage of the approach is that it does not take intoaccountthefullprobabilitydistributionoftheparametersconsistentlywhilecreating the response surface, and there is an inherent assumption that the distribution of these parameters is uniform. As a result, experimental design method may be inappropriate when parameter distributions are arbitrary. 2. For a particular reservoir model, different design methods give different results. And even for a particular experimental design method, different response surfaces can 156 yield different probability density distributions for the model output. Some of the re- sponse surfaces even have significant discrepancy. Additionally, the design methods or response surfaces suitable for one problem may not be proper for other problems. How to choose the best design method or response surfaces is out of the scope of this study. 3. PCM is more efficient than experimental design methods, due to the fact that the orthogonal polynomial expansions (polynomial chaos expansions) capture the nature of the input random variables and generate accurate ”polynomial chaos proxies”. 4. For a particular type of PDF of the input parameters, the associated polyno- mial chaos expansions are used to achieve best convergence. For arbitrary PDFs, the orthogonal polynomial expansions need to be constructed numerically. 5. Foragivencomputationaleffort,thePCMgivesmorerobustresults. UnlikeMonte Carlo simulations which depend on the random samples or the experimental design in which the design points rely on the specific design strategy, the collocation points in the PCM are associated with the given PDFs of the input parameters. Overall, the results demonstrate that PCM is more robust, accurate, and efficient than experimental design for uncertainty analysis. 157 Chapter 7 Conclusions and Recommendations 7.1 Conclusions The detailed conclusions for each chapter are presented at the end of each chapter. They are summarized as follows: Anefficientapproachisproposedforstochasticmodelingofflowinspatiallyheteroge- neousporousmedia. TheGaussianrandomfieldsareefficientlyparameterizedintermsof some independent random variables through the Karhunen-Loeve (KL) expansion. The probabilistic collocation method (PCM) combining with the KL expansion provides an efficient approach for solving stochastic differential equations where random fields or ran- dom processes exist. The approach is non-intrusive to the flow models such that existing codesorsimulatorscanbeemployed. Theexamplesofbothsingle-phaseandmulti-phase flowshowthat,witharelativelysmallnumberofsimulations,theapproachcanaccurately quantify uncertainty of flow in heterogeneous porous media characterized by log-normal random fields. 158 The reservoir properties are also considered as non-Gaussian random fields. Several approaches for simulating non-Gaussian random fields are investigated and their appli- cabilities for uncertainty quantification are studied. Numerical examples of single phase flow demonstrate that the polynomial chaos expansion can be used to simulate the non- Gaussianrandomfields,anditcanbeintegratedwiththeprobabilisticcollocationmethod for efficient uncertainty quantification. The integration of the KL expansion and the PCM into the conditional simulation of both Gaussian and non-Gaussian random fields yields computational efficiency in the process of uncertainty quantification. The conditioning on hydraulic conductivity re- duces the overall uncertainty of the hydraulic conductivity field and the uncertainties in prediction of the flow responses. For reservoir simulations where the uncertain parameters in the reservoir models are consideredasrandomvariableswitharbitraryprobabilitydistributions, theexperimental designmethodscommonlyusedforuniformrandomvariablesareinappropriateforuncer- tainty quantification. The probabilistic collocation method based on the the generalized polynomialschaosexpansionisappropriateforarbitraryrandomvariables. Theexamples with different reservoir models demonstrate that the PCM is more robust, accurate, and efficient than the experimental design methods. 159 7.2 Recommendations Some recommendations are presented as follows: The total number (P) of the collocation points in the PCM is essential since it repre- sents the efficiency of the method. It is equal to the number of terms in the polynomial chaos expansion. It is shown that the total number depends on the the number of input randomparameters(orthenumberofrandomdimensionalityoftheinputrandomfields), N, and the order of the polynomial chaos expansion, d. The accuracy of the solutions generally increases with the increase of N and d. However, the number P and the com- putational effort would also increase rapidly. For a specific problem, the sufficient N or d may be examined numerically by comparing with the next-level PCM. Such a proce- dure could increase more computational efforts. For future research, a posterior error estimator may be developed to guide the choice of the proper N and d. The hydraulic conductivity considered in this study is a stationary non-Gaussian random field. Simulation of non-stationary non-Gaussian fields and the subsequent un- certainty quantification will be a new topic. The marginal PDFs of the non-Gaussian random fields considered in this study are skewed continuous functions. When the PDFs are discrete or multi-modal, e.g., for dis- crete patterns of rock types, multiple-point geostatistics provides a choice for generating realizations of the random field, however the uncertainty quantification is restricted to Monte Carlo method, which is computationally expensive. Parameterization of the ran- dom fields with discrete distributions as well as the efficient uncertainty quantification will be a challenge. 160 This study focuses on the forward modeling of uncertainty quantification. When the measurements of flow responses, such as the hydraulic head of groundwater and the production data from oil wells are also available, estimation of the formation properties such as the hydraulic conductivity or permeability is another topic. That is the inverse modeling (named the history matching in petroleum engineering). The computational efficiency of the PCM could accelerate the process if it is integrated to the algorithms of inverse modeling. The applications of subsurface flow and transport are of great breadth. The efficiency of the PCM allows it to be applied to various field-scale problems. It can play an vi- tal role in uncertainty quantification of the real-world engineering problems, such as the management of groundwater, remediation of underground contamination transport, ge- ological sequestration of carbon dioxide, and exploration of oil and gas as well as other non-conventionalfossilenergy. Riskassessmentwiththeaidofuncertaintyquantification of the subsurface flow and transport would be an indispensable step in the activities of underground exploration, production, and management. 161 References C. Amudo, T. Graf, N.R. Harris, R. Dandekar, F. Ben Amor, and R.S. May. Exper- imental design and response surface models as a basis for stochastic history match-a Nigerdeltaexperience. IntheInternationalPetroleumTechnologyConference,Kuala Lumpur, Malaysia, 3-5 December 2008. F. Ballio and A. Guadagnini. Convergence assessment of numerical Monte Carlo simulations in groundwater hydrology. Water Resour. Res., 40:W04603, doi:10.1029/2003WR002876, 2004. P. Bocchini and G. Deodatis. Critical review and latest development of a class of simulation algorithms for strongly non-Gaussian random fields. Probabilistic Engi- neering Mechanics, 23:393–407, 2008. M.ChenandD.Zhang. Stochasticanalysisoftwophaseflowinheterogeneousmedia by combining Karhunen-Loeve expansion and perturbation method. Water Resour. Res., 41:W01006, doi:10.1029/2004WR003412, 2005. M. Chen, A. A. Keller, D. Zhang, Z. Lu, and G. A. Zyvoloski. Stochastic analysis of transient two-phase flow in heterogeneous porous media. Water Resour. Res., 42: W03425, doi:10.1029/2005WR004257, 2006a. Z. Chen, G. Huan, and Y. Ma. Computational Methods for Multiphase Flows in Porous Media. SIAM, 2006b. J. H. Cushman. The Physics of Fluids in Hierarchical Porous Media: Angstroms to Miles. Kluwer Academic, Norwell, MA., 1997. G.Dagan. Stochasticmodelingofgroundwaterflowbyunconditionalandconditional probabilities: 1. conditional simulation and the direct problem. Water Resour Res, 18(4):813–33, 1982. G. Dagan. Flow and Transport in Porous Formations. Springer, New York, 1989. G. Deodatis and R.C. Micaletti. Simulation of highly skewed non-Gaussian stochas- tic processes. J. Engng. Mech, 127(12):1284–95, 2001. L. Feyen and J. Caers. Quantifying geological uncertainty for flow and transport modeling in multi-modal heterogeneous formations. Advances in Water Resources, 29:912–929, 2006. 162 F. Friedmann, A. Chawathe, and D. K. Larue. Assessing uncertainty in channelized reservoirs using experimental designs. SPE Reservoir Evaluation and Engineering, 2003. W. Gautschi. Algorithm 726: ORTHPOL- A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Transactions on Mathematical software, 20(1):21–62, 1994. L. W. Gelhar. Stochastic Subsurface Hydrology. Prentice-Hall, Englewood Cliffs, NJ., 1993. R. Ghanem. Scales of fluctuation and the propagation of uncertainty in random porous media. Water Resour. Res., 34(9):2123–2136, 1998. R. Ghanem and S. Dham. Stochastic finite element analysis for multiphase flow in heterogeneous porous media. Tranp. Porous Media, 32:239, 1998. R. Ghanem and P. Spanos. Stochastic Finite Element: A Spectral Approach. Springer, New York, 1991. W. D. Graham and D. McLaughlin. Stochastic analysis of nonstationary subsurface solute transport: 1. unconditional moments. Water Resour. Res., 25(2):215–232, 1989. M.Grigoriu. Simulationofstationarynon-Gaussiantranslationprocesses. J. Engng. Mech., 124(2):121–6, 1998. A. Guadagnini and S. P. Neuman. Nonlocal and localized analysis of conditional mean steady state flow in bounded, randomly nonuniform domains: 1. theory and computational approach. Water Resour. Res., 35(10):2999–3018, 1999. A. E. Hassan, J. H. Cushman, and J. W. Delleur. A Monte Carlo assessment of Eulerian flow and transport perturbation models. Water Resour. Res., 34(5):1143– 1163, 1998. A. Hyvarinen and E. Oja. Independent component analysis: algorithms and appli- cations. Neural Networks, 13:411–430, 2000. S. S. Isukapalli, A. Roy, and P. G. Georgopoulos. Stochastic response surface meth- ods (SRSMs) for uncertainty propagation: Application to environmental and biolog- ical systems. Risk Analysis, 18(3):351–363, 1998. S. Kalla and C.D. White. Efficient design of reservoir simulation studies for devel- opment and optimization, SPE-95456-PA. SPE Res Eval Eng, 10(6):629–637, 2007. M. Khhalil and A. Sarkar. Independent component analysis for uncertainty repre- sentationofstochasticsystems. In49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Schaumburg, USA, 2008. J.E.Killough. NinthSPEcomparativesolutionproject: areexaminationofblack-oil simulation. SPE 29110, 1995. 163 Y. S. Kumar, J. Li, C. Talarico, and J. Wang. A probabilistic collocation method based statistical gate delay model considering process variations and multiple input switching. InProceedings of Design, Automation and Test in Europe,pages770–775, 7-11, March 2005. K.A. Lawal. Modeling subsurface uncertainty with experimental design: some ar- guments of non-conformists, SPE-128350. In the 33rd Annual SPE International Technical Conference and Exhibition, Abuja, Nigeria, 3-5, August 2009. G. Liu, D. Zhang, and Z. Lu. Stochastic uncertainty analysis for unconfined flow systems. Water Resour. Res., 42:W09412, doi:10.1029/2005WR004766, 2006. G. Liu, Z. Lu, and D. Zhang. Stochastic uncertainty analysis for solute transport in randomly heterogeneous media using a Karhunen-Loeve based moment equation approach. Water Resour. Res., 43:W07427, doi:10.1029/2006WR005193, 2007. Z. Lu and D. Zhang. Conditional simulations of flow in randomly heterogeneous porous media using a KL-based moment-equation approach. Advances in Water Resources, 27:859–874, 2004. Z. Lu and D. Zhang. Accurate, effective quantification of uncertainty for flow in heterogeneous reservoirs using the KLME approach. SPE J., 12:239–247, 2006. Z. Lu and D. Zhang. Stochastic simulations for flow in nonstationary randomly heterogeneous porous media using a KL-based moment-equation approach. SIAM Multiscale Modeling and Simulation, 6(1):228–245, 2007. Z. Lu, SP. Neuman, A. Guadagnini, and DM. Tartakovsky. Conditional moment analysis of steady-state unsaturated flow in bounded, randomly heterogeneous soils. Water Resour Res, 38(4):10.1029/2001WR000278., 2002. L. Mathelin, M. Y. Hussaini, and T. A. Zang. Stochastic approaches to uncertainty quantification in CFD simulations. Numerical Algorithms, 38:209–236, 2005. S. P. Neuman. Eulerian-Lagrangian theory of transport in spacetime nonstationary velocity fields: Exact nonlocal formalism by conditional moments and weak approx- imations. Water Resour. Res., 29:633–645, 1993. S. P. Neuman. Stochastic approach to subsurface flow and transport: A view to the future. In G. Dagan and S. P. Neuman, editors, Subsurface Flow and Transport: A Stochastic Approach, pages 231–241. Cambridge Univ. Press, New York, 1997. K. K. Phoon, S.P. Huang, and S.T. Quek. Simulation of second-order processes using Karhunen-Loeve expansion. Computers and Structures, 80:1049–1060, 2002. Y. Rathi, S. Dambreville, and A. Tannenbaum. Statistical shape analysis using kernel PCA. In SPIE, Electronic Imaging, 2006. Y. J. Ren, I. Elishakoff, and M. Shinozuka. Conditional simulation of non-Gaussian random fields for earthquake monitoring systems. Chaos, Solitons and Fractals, 5 (1):91–101, 1995. 164 Y. Rubin. Prediction of tracer plume migration in disordered porous media by the method of conditional probabilities. Water Resour Res, 27(6):1291–308, 1991. Y. Rubin. Applied Stochastic Hydrogeology. Oxford University Press, New York, 2003. S. Sakamoto and R. G. Ghanem. Polynomials chaos decomposition for simulation of non-Gaussian non-stationary stochastic processes. ASCE J. Engrg. Mech., 128(2): 190–200, 2002. P. Sarma, L. J. Durlofsky, and K. Aziz. Efficient closed-loop production optimiza- tion under uncertainty, SPE-94241. In the SPE Europec/EAGE Annual Conference, Madrid, Spain, 13-16, June, 2005. P. Sarma, L. J. Durlofsky, and K. Aziz. Kernel principal component analysis for efficient, differentiable parameterization of multipoint geostatistics. Math. Geosci., 40:3–32, 2008. T.Schaaf,B.Coureaud,andN.Labat. Assistedhistorymatchingtoolsandbayesian framework to get probabilistic production forecasts. In the Europec/EAGE Confer- ence and Exhibition, Rome, Italy, 9-12, June 2008. B.S.Scholkopf,A.J.Mika,G.Smola,K.Ratsch,andR.Muller. KernelPCApattern reconstruction via approximate pre-images. In the 8th international conference on artificial neural networks, Skovde, Sweden, 1998. S.B. Strebelle. Conditional simulation of complex geological structures using multiple-point statisticss. Mathematical Geology, 34:1–22, 2000. S.B. Strebelle and A.G. Journel. Reservoir modeling using multiple-point statis- tics, SPE-71324. In SPE Annual Technical Conference and Exhibition held in New Orleans, Louisiana, 30 September-3 October, 2001. M. A. Tatang, W. Pan, R. G. Prinn, and G. J. McRae. An efficient method for parametric uncertainty analysis of numerical geophysical models. J. Geophy. Res., pages 21925–21931, 1997. E. D. Tipping, Deimbacher X. F., D. Kovyazin, A. Medvedev, and R. Valencia. Experimental design reduces risks of a Siberian development with subsurface and financial uncertainties. In the SPE Russian Oil and Gas Technical Conference and Exhibition, 28-30, Octobor 2008. J. F. Van Elk, L. Guerrera, K. Vijayan, and R. Gupta. Improved uncertainty man- agement in field development through the application of the experimental design method to the multiple realizations approach, SPE-64462. In the SPE Annual Tech- nical Conference and Exhibition, Dallas, Texas, USA, 1-4, Octobor 2000. R. Venkataraman. Application of the method of experimental design to quantify uncertainty in production profiles, SPE-59422. In the SPE Asia Pacific Conference on Integrated Modelling for Asset Management,Yokohoma,Japan,25-26,April2000. 165 J. Villadsen and M. L. Michelsen. Solution of Differential Equation Models by Poly- nomial Approximation. Prentice-Hall, 1978. X. Wan and G.E. Karniadakis. Solving elliptic problems with non-Gaussian spatially-dependent random coefficients. Comput. Methods Appl. Mech. Engrg., 198: 1985–1995, 2009. M. Webster, M. A. Tatang, and G. J. Mcrae. Application of the probabilistic col- location method for an uncertainty analysis of a simple ocean model. MIT Joint Program on the Science and Policy of Global Change Report Series No. 4, Mas- sachusetts Institute of Technology, 1996. D. Xiu and G. E. Karniadakis. Modeling uncertainty in steady state diffusion prob- lems via generalized polynomial chaos. Comput. Methods Appl. Mech. Eng., 191: 4927–4948, 2002. F.YamazakiandM.Shinozuka. Digitalgenerationofnon-Gaussianstochasticfields. J. Engng. Mech., 114(7):1183–1197, 1988. J. Yang, D. Zhang, and Z. Lu. Stochastic analysis of saturated-unsaturated flow in heterogeneous media by combining Karhunen-Loeve expansion and perturbation method. J. hydrol., 294:18–38, 2004. B. Yeten, A. Castellini, B. Guyaguler, and W. H. Chen. A comparison study on experimental design and response surface methodologies, SPE-93347. In the 2005 SPE Reservoir Simulation Symposium, Houston, Texas, 31 January -2 February 2005. D. Zhang. Numerical solutions to statistical moment equations of groundwater flow in nonstationary, bounded heterogeneous media. Water Resour. Res., 34:529–538, 1998. D. Zhang. Nonstationary stochastic analysis of transient unsaturated flow in ran- domly heterogeneous media. Water Resour. Res., 35:1127–1141, 1999. D.Zhang. StochasticMethodsforFlowinPorousMedia: CopingWithUncertainties. Academic Press, San Diego, CA, 2002. D. Zhang and Z. Lu. An efficient, high-order perturbation approach for flow in random porous media via Karhunen-Loeve and polynomial expansions. J. Comput. Phys., 194:773–794, 2004. D. Zhang and S. P. Neuman. Eulerian-Lagrangian analysis of transport conditioned on hydraulic data, 1. analytical-numerical approach. Water Resour. Res., 31:39–51, 1995. D. Zhang and H. Tchelepi. Stochastic analysis of immiscible two-phase flow in het- erogeneous media. SPE Journal, 4(4):380–388, 1999. D. Zhang, L. Li, and H. A. Tchelepi. Stochastic formulation for uncertainty analysis of two-phase flow in heterogeneous reservoirs. SPE Journal, 5(1):60–70, 2000. 166
Abstract (if available)
Abstract
Uncertainty quantification of subsurface flow has recently attracted a significant amount of attention. The uncertainty can result from the combination of the formation heterogeneity and the incomplete knowledge of its properties. Traditional flow simulations treat the geological formation deterministic, thus resulting in deterministic predictions. However, taking uncertainty into account necessitates a stochastic description of the formation properties and hence stochastic approaches to flow simulations.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Inverse modeling and uncertainty quantification of nonlinear flow in porous media models
PDF
Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems
PDF
Stochastic oilfield optimization under uncertainty in future development plans
PDF
Deep learning architectures for characterization and forecasting of fluid flow in subsurface systems
PDF
Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing
PDF
An extended finite element method based modeling of hydraulic fracturing
PDF
The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media
PDF
Efficient control optimization in subsurface flow systems with machine learning surrogate models
PDF
A method for characterizing reservoir heterogeneity using continuous measurement of tracer data
PDF
Uncertainty management for complex systems of systems: applications to the future smart grid
PDF
Subsurface model calibration for complex facies models
PDF
Pore-scale characterization and fluid-flow simulation: application of the Lattice Boltzmann method to conventional and unconventional rocks
PDF
Sparse feature learning for dynamic subsurface imaging
PDF
Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation
PDF
History matching for fractured reservoirs with Hough-transform-based parameterization
PDF
Assessing induced seismicity rate increase based on deterministic and probabilistic modeling
PDF
Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow
PDF
Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation
PDF
Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection
PDF
Design optimization under uncertainty for rotor blades of horizontal axis wind turbines
Asset Metadata
Creator
Li, Heng
(author)
Core Title
Accurate and efficient uncertainty quantification of subsurface fluid flow via the probabilistic collocation method
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Petroleum Engineering
Publication Date
10/12/2010
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
OAI-PMH Harvest,porous media,probabilistic collocation method,random field,uncertainty quantification
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Zhang, Dongxiao (
committee chair
), Ershaghi, Iraj (
committee member
), Ghanem, Roger Georges (
committee member
)
Creator Email
hengli@usc.edu,liheng2005@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m3504
Unique identifier
UC1476670
Identifier
etd-Li-4164 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-418996 (legacy record id),usctheses-m3504 (legacy record id)
Legacy Identifier
etd-Li-4164.pdf
Dmrecord
418996
Document Type
Dissertation
Rights
Li, Heng
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
porous media
probabilistic collocation method
random field
uncertainty quantification