University of Southern California Dissertations and Theses

Risk assessment, intrinsic interpolation and computationally efficient models for systems under uncertainty

(USC Thesis Other)

Risk assessment, intrinsic interpolation and computationally efficient models for systems under uncertainty

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content RISK ASSESSMENT, INTRINSIC INTERPOLATION AND COMPUTATIONALLY EFFICIENT MODELS FOR SYSTEMS UNDER UNCERTAINTY by Charanraj Thimmisetty 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 (CIVIL AND ENVIRONMENTAL ENGINEERING) May 2016 Copyright 2016 Charanraj Thimmisetty To family and friends. ii Acknowledgments Foremost, I would like to express my sincere gratitude to my advisor Prof. Roger Ghanem for the time he spent with me, his honest recommendations and criticism from time to time, which inspired me to excel in my work. I am grateful to him for his encouragement and positive attitude towards my work and showing me the right path throughout my PhD. I would like to thank him for providing me many opportunities to present my work on dierent platforms which increased my condence and leadership skills. I would like to thank Prof. Sami Masri, Prof. Fred Aminzadeh, Prof. Fei Sha, Prof. Vincent Lee, Dr. Ketan Savla and Dr. Felipe De Barros for serving in my qualifying, and screening committees and providing valuable suggestions for my research. I would like to thank my collaborators, Dr. Arman Khodabakhshnejad, Dr. Nima Jabbari and Prof. Fred Aminzadeh for sharing their knowledge and helping me to build reservoir models. I would like to thank Dr. Kelly Rose, Dr. Jennifer Bauer and Dr. Corinne Disenhof of the National Energy Technology Laboratory (NETL) for their valu- able inputs and sharing their knowledge on the Bureau of Ocean Energy Management (BOEM) database. I would like to thank Nicolas Leseur of the Saudi Aramco company for providing me the permeability and well-log data and giving me the insights into the data. I would like to thank Dr. Venkat Aitharaju, John N. Owens and Selina Zhao of the General Motors corporation for their valuable suggestions and experimental data of resin transport modeling (RTM). Also, I would like to thank Praveen Pasupuleti, Dr. Xiaoshi iii Jin and Dr. Hang Yu of ESI corporation for helping me to build the deterministic models of RTM and teaching me how to use ESI software PAM-RTM. I gratefully acknowledge the nancial support received from the US Department of Energy and National Science Foundation during this work. I would like to thank the Viterbi school of engineering for oering me Provost fellowship. Also, I would like to thank Sonny Astani Department of Civil and Environmental Engineering for oering me teaching assistantship. I would like to extend my sincere thanks to all my colleagues and friends who en- couraged directly or indirectly during this work. I would like to thank Dr. Ramakrishna Tippireddy and Dr. Evangelia Kalligiannaki for their guidance in the earlier phase of my PhD. My sincere thanks to my former and current fellow group members Dr. Vahid Keshavarzzadeh, Dr. Bedrick Sousdik, Dr. Hadi Meidani, Dr. Hamed Haddadzadegan, Iman Yadegaran, Ruda Zhang, Nastaran Bassamzadeh, Dr. Loujaine Mehrez and Pana- giotis Tsilis for sharing their knowledge. I would like to thank my friends Dr. Prasanth Koganti and Dr. Deval Modi for proofreading my documents and their encouragement throughout my PhD. Also, I would like to thank my friends Viraj Shirodkar, Jayasankar Malladi and Nageswara Rao Janjanam for their support and encouragement. I would like to thank Prof. C. S. Manohar, who introduced me to the eld of un- certainty quantication, guiding my research work at the Indian Institute of Science. I am grateful to all my teachers who taught and inspired me throughout my student life. Last but not the least, I would like to thank my family: my father Anjaneyulu, my mother Koteswaramma, sisters Sunitha & Gowthami and my uncle Siva Nageswara Rao for their unconditional love and sacrices without which this work would not have been completed. iv Contents Dedication ii Acknowledgments iii List of Tables viii List of Figures ix Abstract xii Chapter 1: Introduction 1 Chapter 2: Multiscale Stochastic Representation in High-Dimensional Data using Gaussian Processes and Polynomial Chaos Expansion 5 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Motivation and challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Summary of the Methodology . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.1 Reservoir QoI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.1.1 PCE of reservoir QoI . . . . . . . . . . . . . . . . . . . . 13 2.3.2 Wellbore QoI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.2.1 Well in ow performance relationship (IPR) curve . . . . 15 2.3.2.2 Well Tubing Performance Relationship (TPR) curve . . . 16 2.3.2.3 PCE of wellbore QoI . . . . . . . . . . . . . . . . . . . . 17 2.3.3 Reservoir reduced order model . . . . . . . . . . . . . . . . . . . . 17 2.3.4 Wellbore reduced model . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.1 PCE of QoI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4.2 Reservoir ROM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4.2.1 RROM with known parameters at the target site . . . . . 26 2.4.2.2 RROM with unknown parameters at the target site . . . 28 2.4.3 Wellbore ROM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4.3.1 Normal scenario . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.3.2 Damaged scenario with known damage location . . . . . 32 2.4.3.3 Damaged scenario with unknown damage location . . . . 39 2.5 Data analytics on the BOEM database . . . . . . . . . . . . . . . . . . . . 39 v 2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Chapter 3: High-Dimensional Intrinsic Interpolation using Gaussian Pro- cess Regression 43 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2.1 Gaussian process regression . . . . . . . . . . . . . . . . . . . . . . 46 3.2.2 Diusion maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.3 Intrinsic interpolation using Gaussian process regression . . . . . . 50 3.3 Case study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1 Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.2 Diusion map of the dataset . . . . . . . . . . . . . . . . . . . . . 52 3.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Chapter 4: An Ecient Basis Adaptation in Homogeneous Chaos Space for Stochastic Initial-boundary Value Problem with Specic Quantities of Interest 61 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.3 Stochastic dimensions of the problem . . . . . . . . . . . . . . . . . . . . . 64 4.4 Polynomial chaos decomposition . . . . . . . . . . . . . . . . . . . . . . . 65 4.5 Basis adaptation to QoI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.5.1 Gaussian adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.6 Ecient computation of the isometry for initial boundary value problems for specic quantities of interest . . . . . . . . . . . . . . . . . . . . . . . . 68 4.7 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.7.1 Basis adaptation using proposed method . . . . . . . . . . . . . . . 72 4.7.2 Computational eciency of the proposed method . . . . . . . . . . 76 4.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Chapter 5: Convergence Enhancement of Basis Adaptation in Homoge- neous Chaos Spaces 78 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.2 Modeling the RTM process . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2.1 Flow in the porous media . . . . . . . . . . . . . . . . . . . . . . . 80 5.2.2 Thermal phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.2.3 The kinetics of the resin polymerization . . . . . . . . . . . . . . . 82 5.3 Polynomial chaos decomposition . . . . . . . . . . . . . . . . . . . . . . . 82 5.3.1 Basis adaptation to QoI . . . . . . . . . . . . . . . . . . . . . . . . 83 5.3.1.1 Gaussian adaptation . . . . . . . . . . . . . . . . . . . . . 84 5.4 Convergence enhancement of basis adaptation . . . . . . . . . . . . . . . . 85 5.5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.5.1 Basis adaptation of the modied Rosenbrock function . . . . . . . 87 5.5.2 Uncertainty propagation in the RTM process . . . . . . . . . . . . 89 vi 5.5.2.1 Uncertainty propagation by treating permeability as a random process . . . . . . . . . . . . . . . . . . . . . . . 90 5.5.2.2 Uncertainty propagation in 18 dimensions . . . . . . . . . 93 5.5.2.3 Uncertainty propagation in 13 dimensions . . . . . . . . . 94 5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Bibliography 97 Appendix A 102 vii List of Tables 2.1 Few attributes in BOEM database . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Random input parameters used for the reservoir model. . . . . . . . . . . 23 2.3 Mean square error for dierent models for training and initial test reservoirs. 27 2.4 Statistics of the ow rate on day one (QoI) from the actual model, pre- dictions using GPR with known & unknown parameters and preliminary estimate models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Flow values for dierent breach scenarios. . . . . . . . . . . . . . . . . . . 36 2.6 Equivalent normal scenarios for failure scenarios 1 and 2. . . . . . . . . . 36 2.7 Equivalent normal scenarios for failure scenario 3. . . . . . . . . . . . . . . 36 2.8 Relation between breach location and the ow rate for P out = 1; 844: . . . 37 2.9 Relation between breach location and the ow rate for P out = 2; 144: . . . 38 2.10 Parameters used for the data analytics . . . . . . . . . . . . . . . . . . . . 40 2.11 Number of reservoirs in each cluster of the BOEM data . . . . . . . . . . 40 3.1 Mean square error (MSE) for dierent correlation functions and regression models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.1 Relative permeability table for oil and water. . . . . . . . . . . . . . . . . 71 4.2 Pressure volume temperature (PVT) of the dead oil. . . . . . . . . . . . . 71 4.3 Computational cost associated with PCE, Gaussian adaptation and the proposed method, here c t is computational time to solve the PDE up to time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.1 Computational cost associated with dierent uncertainty propagation meth- ods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.2 Mean and coecient of variations of the random variables considered in the example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3 Distributions of the random variables considered in the example . . . . . 94 viii List of Figures 2.1 Subgroups in BOEM Database . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Spatial location of boreholes over the Gulf of Mexico; Source of the data [1] 10 2.3 An example of IPR and TPR curves used for the simulations. . . . . . . . 16 2.4 Spatial location of boreholes over the Gulf of Mexico and dierent groups in the datasets: red circles, yellow triangles, and green stars respectively indicate nal test set, initial test and training reservoirs. . . . . . . . . . . 21 2.5 Distribution of the centered (a) log permeability (b) log porosity for few of the elds under study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Probability distribution function (pdf) of oil production rate for few train- ing reservoirs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.7 Predictions of the rst PCE coecient using GPR with squared exponen- tial autocorrelation at nal test reservoirs. . . . . . . . . . . . . . . . . . . 27 2.8 QoI: Flow rate on day one in bbl (a) Solid green line indicates pdf of QoI at prediction location obtained from available data at that location; red dashed line indicates predicted pdf of QoI using Gaussian Process regression over parameter space; (b) red dashed and dotted line indicates samples of the pdf for QoI, each sample pdf re ects a realization from the Gaussian Process regression over parameter space; The black dashed line indicates most likely pdf at the prediction location obtained as an average over Gaussian Process regression. . . . . . . . . . . . . . . . . . . . . . . . 28 2.9 QoI: Flow rate on day one in bbl (a) Solid green line indicates pdf of QoI at prediction location obtained from available data at that location; Dashed red line indicates predicted pdf of QoI using Gaussian Process regression over interpolated parameter space; (b) Dashed and dotted lines indicate samples of pdf of QoI, each sample pdf re ects a realization from the Gaussian Process regression over the interpolated parameter space; Each color in dotted line indicates dierent realization of the parameters; 29 2.10 Probability Density Function for various quantities of interest at sea oor; undamaged pipe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.11 Normal scenario: Geometry of the pipe system. . . . . . . . . . . . . . . . 32 2.12 Schematic of ow in drillpipe system. . . . . . . . . . . . . . . . . . . . . . 33 2.13 Failure scenario 1: Geometry of the pipe system. . . . . . . . . . . . . . . 34 2.14 Failure scenario 2: Geometry of the pipe system. . . . . . . . . . . . . . . 34 2.15 Failure scenario 3-a: Geometry of the pipe system. . . . . . . . . . . . . . 35 2.16 Failure scenario 3-b: Geometry of the pipe system. . . . . . . . . . . . . . 35 ix 2.17 Probability Density Function of oil ow rate. Pipe diameter of 7" is con- sidered nominal with other pipe diameters construed as deviations from normal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.18 Plot of 36D BOEM data embedding into the rst two diusion map co- ordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1 Spatial locations of the wells. . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.2 (a) Eigenvalue decay of the diusion kernel (b) vs L() curve . . . . . . . 53 3.3 Scatter plot of the rst three dominant directions in (a) Diusion space (b) Principle component space. . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Histogram of the measured permeability in (a) Original space (b) Trans- formed space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5 Blind test results for the well 1 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.6 Blind test results for the well 2 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.7 Blind test results for the well 23 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.8 Blind test results for the well 29 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.9 Blind test results for the well 106 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.10 Blind test results for the well 107 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.11 Blind test results for the well 116 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.12 Blind test results for the well 236 using (a) Proposed method (b) Gaussian process regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1 (a) Spatial discretization of the reservoir (b) A realization of the perme- ability obtained from the chosen covariance kernel. . . . . . . . . . . . . . 70 4.2 A realization of the (a) Pressure (b) Oil saturation at the end of one year obtained from solving the numerical reservoir model with chosen covariance kernel for the permeability. . . . . . . . . . . . . . . . . . . . . 72 4.3 Probability density function of cumulative oil production in MSTB for year one, ve, ten, fteen and twenty obtained using the proposed method. 73 4.4 Evaluation of the 1 over the (a) years one to ten (b) years eleven to twenty 74 4.5 First nine PCE coecients of the permeability with respect to (a) original Gaussian variable (b) one dimensional Gaussian variable 1 . . . . . . . 75 5.1 Probability distribution function of the modied Rosenbrock function for dierent dimensional adaptations obtained using isometry with initialA 0 based on the identity matrix. . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.2 Probability distribution function of the modied Rosenbrock function, for dierent dimensional adaptations obtained using the proposed method. . 89 x 5.3 Initial A 0 (top left) and isometry A (top right) based on the identity matrix. InitialA 0 (bottom left) and isometryA (bottom right) from the proposed method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.4 RTM setup used for the demonstration . . . . . . . . . . . . . . . . . . . . 90 5.5 Probability distribution function of lltime from 1D, 2D, 3D adaptations: obtained using the proposed method (left), and with initialA 0 based on the identity matrix (right). . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.6 Initial A 0 (top left) and isometry A (top right) based on the identity matrix. InitialA 0 (bottom left) and isometryA (bottom right) from the proposed method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.7 Histogram of the lltime obtained by sampling using: Monte Carlo sim- ulations (top left), Gaussian part of PCE (top right), PCE from the 1D adaptation (bottom left), PCE from the 2D adaptation (bottom right). . 94 5.8 Probability density function of the ll time obtained using basis adapta- tions in 1D and 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 xi Abstract This dissertation focuses on developing probabilistic models for predicting the behavior of physical systems in the presence of uncertainty. Novel contributions are made to distinct aspects of the problem, namely, 1) risk assessment and risk mitigation models for the ultra-deepsea drilling in the Gulf of Mexico, 2) regression models for intrinsic variables constrained onto a manifold embedded in an ambient space, and, 3) ecient isometry computation for the stochastic basis adaptation method to detect lower dimensional stochastic space adapted to quantities of interest. In particular, the rst part of the thesis addresses an important component of risk mitigation for ultra-deepsea drilling in the Gulf of Mexico, namely the probabilistic char- acterization of uid uxes at the sea oor from future drilling operations. In the process, a stochastic representation of functions is dened on a high-dimensional space condi- tional on their marginal statistics over a denumerable subset of this space, together with their global correlation structure. The representation leverages a particular structure of the functional dependence of interest which exhibits scale separation. Specically, we construct a polynomial chaos representation for the scalar quantity of interest (QoI) whose coecients are themselves random. The intrinsic randomness of the polynomial chaos expansion (PCE) re ects local uncertainty and captures dependence on a subset (say S 1 ) of the parameters. While randomness in the PCE coecients captures a global structure of the uncertainty and dependence on the remaining parameters (say S 2 ) in xii the high-dimensional space (let S =S 1 US 2 ). This construction is demonstrated by pre- dicting wellbore signatures in the Gulf of Mexico (GoM) where a 120-dimensional table is populated at several thousand wellbore locations throughout the GoM. Physics-based models of multiphase ow in porous media are used to calculate the PCE representations at the sites where data is available. In this context, random parameters describing the subsurface dene the parameter set (S 1 ). A Gaussian process model is then developed for each coecient in these representations, construed as a function onS 2 . The combined probabilistic representation permits the delineation of separate stochastic in uences on predictions of interest. In the second part, a regression procedure is proposed for intrinsic variables con- strained onto a manifold embedded in an ambient space. The procedure is meant to sharpen high-dimensional interpolation by introducing constraints delineated from within the data being interpolated. Proposed method augments manifold learning procedures with a Gaussian process regression (GPR). This method, rst identies, using diusion maps, a low dimensional manifold embedded in an ambient high dimensional space asso- ciated with the data. The diusion distance associated with this construction is used to dene proximity between points on the manifold. This distance is then used to compute the correlation structure of a Gaussian process that describes the dependence of quanti- ties of interest in the high dimensional ambient space. The proposed method can be used to predict any quantity on the arbitrary high-dimensional dataset; a specic example is considered to characterize the permeability of rock formations from petrophysics data. In the nal part, two ecient isometry computation methods for stochastic basis adaptation procedure to detect lower dimensional stochastic space adapted to quantities of interest are proposed. Specically, the rst method is useful to the stochastic initial- boundary value problem with specic quantities of interest. It uses the solution from the Gaussian part of the QoI at time t T to compute the isometry of the QoI at time T . The proposed method is used to compute a basis adapted for cumulative ow rate in a petroleum reservoir with random permeability. The second method uses the xiii sensitivity obtained from the Gaussian part of the solution to enhance the convergence properties of the adaptation. This method is used to develop computationally ecient models for uncertainty propagation and calibration of various stages in a typical carbon ber composite manufacturing processes such as lling, curing and distortion. xiv Chapter 1 Introduction Although deterministic models satisfactorily simulate the real world processes for a large number of applications, they are often not adequate to explain observed behavior in physical systems. For instance, sometimes the underlying physics is unknown, such as in weather predictions. Other times, the models are so complex that we have to resort to simpler models due to limitations in computational resources, for example, many nonlinear phenomena are described using linear models. Sometimes, we cannot accurately measure the model parameters, such as micro-scale permeability in a petroleum reservoir. In such instances, probabilistic models allow us to capture the discrepancy between the associated model and physical behavior of the system. This dissertation focuses on developing probabilistic models for predicting the behavior of physical systems in the presence of uncertainty. Novel contributions are made to distinct aspects of the problem, namely, 1) risk assessment and risk mitigation models for the ultra-deepsea drilling in the Gulf of Mexico, 2) regression models for intrinsic variables constrained onto a manifold embedded in an ambient space, and, 3) ecient isometry computation for the stochastic basis adaptation method to detect lower dimensional stochastic space adapted to quantities of interest. Multiscale stochastic representation in high-dimensional data using Gaussian processes This work is motivated by uncertainty associated with predictions of wellbore signatures due to anthropogenic events in hydrocarbon production such as the one encountered in Deepwater Horizon oil spill. During this event attempts were made to control and 1 cease the blowout. The scarcity of data, as well as uncertainty in leakage and seepage mechanisms, made associated predictions challenging. Several investigations have been conducted to model the blowout and estimate the damage. It took almost two months to get reliable wellbore signatures which were used to estimate the degree of contamination in the Gulf of Mexico (GoM). In collaboration with the National Energy Technology Laboratory (NETL) and colleagues from the petroleum engineering department here at USC, an important component of risk mitigation for ultra-deepsea drilling in the GoM is addressed, namely the probabilistic characterization of uid uxes at the sea oor from future drilling operations [2]. In the process, a stochastic representation of functions is dened on a high-dimensional space conditional on their marginal statistics over a denu- merable subset of this space, together with their global correlation structure. The repre- sentation leverages a particular structure of the functional dependence of interest which exhibits scale separation. Specically, a polynomial chaos representation for the scalar quantity of interest (QoI) whose coecients are themselves random is constructed. The intrinsic randomness of the polynomial chaos expansion (PCE) re ects local uncertainty and captures dependence on a subset (say S 1 ) of the parameters. While randomness in the PCE coecients captures a global structure of the uncertainty and dependence on the remaining parameters (say S 2 ) in the high-dimensional space (let S = S 1 US 2 ). This construction is demonstrated by predicting wellbore signatures in the GoM where a 120-dimensional table is populated at several thousand wellbore locations throughout the GoM. Physics-based models of multiphase ow in porous media are used to calculate the PCE representations at the sites where data is available. In this context, random parameters describing the subsurface dene the parameter set (S 1 ). A Gaussian pro- cess model is then developed for each coecient in these representations, construed as a function on S 2 . The combined probabilistic representation permits the delineation of separate stochastic in uences on the prediction of the QoI. 2 High-Dimensional Intrinsic Interpolation using Gaussian Process Re- gression This research work is motivated by challenges faced in the characterization of the per- meability of rock formations in a petroleum reservoir, due to microscale uctuations, scarcity, and uncertainty associated with the data collection. For this work, a regres- sion procedure is proposed for intrinsic variables constrained onto a manifold embedded in an ambient space. The procedure is meant to sharpen high-dimensional interpola- tion by introducing constraints delineated from within the data being interpolated. Pro- posed method augments manifold learning procedures with a Gaussian process regression (GPR). This method, rst identies, using diusion maps, a low dimensional manifold embedded in an ambient high dimensional space associated with the data. The diusion distance associated with this construction is used to dene proximity between points on the manifold. This distance is then used to compute the correlation structure of a Gaussian process that describes the dependence of quantities of interest in the high dimensional ambient space. An ecient basis adaptation in homogeneous chaos space for stochastic initial-boundary value problem with specic quantities of interest This research work is motivated by the randomness of the material properties from one point of the domain to the other, leading to high stochastic dimension (N). Even with the advent of powerful computers, numerically solving the partial dierential equations of corresponding physics is computationally expensive due to huge mesh sizes. For this problem, an ecient isometry computation method for the stochastic basis adaptation procedure to detect lower dimensional (nN) stochastic space adapted to quantities of interest is proposed. Especially, this method is useful for the stochastic initial-boundary value problem with specic quantities of interest. It uses the solution from the Gaussian part of the QoI at time t T to compute the isometry of the QoI at time T . The 3 proposed method is used to compute a basis, adapted for cumulative ow rate in a petroleum reservoir with random permeability. An ecient model for uncertainty propagation and calibration of the carbon ber composite manufacturing processes This work is motivated by the diculties associated with uncertainty propagation and calibration of the various stages of the composite manufacturing processes such as ll- ing, curing, and distortion. Randomness in the material properties, computational cost associated with the solution of corresponding physics and dependence of the uncertainty from one stage to the other, makes this problem challenging. For this work, a method was proposed for ecient isometry computation of the stochastic basis adaptation pro- cedure to detect lower dimensional stochastic space adapted to the quantities of interest. Specically, the proposed method uses the sensitivity obtained from the Gaussian part of the solution to enhance the convergence properties of the adaptation. 4 Chapter 2 Multiscale Stochastic Representation in High-Dimensional Data using Gaussian Processes and Polynomial Chaos Expansion This chapter addresses an important component of risk mitigation for ultra-deepsea drilling in the Gulf of Mexico, namely the probabilistic characterization of uid uxes at the sea oor from future drilling operations. In the process, a stochastic representation of functions dened on high-dimensional space conditional on their marginal statistics over a denumerable subset of this space, together with their global correlation structure is de- veloped. The representation leverages a particular structure of the functional dependence of interest which exhibits scale separation. Specically, a polynomial chaos representa- tion for the scalar quantity of interest (QoI) whose coecients are themselves random is constructed. The intrinsic randomness of the polynomial chaos expansion (PCE) re- ects local uncertainty and captures dependence on a subset (say S 1 ) of the parameters. While randomness in the PCE coecients captures a global structure of the uncertainty and dependence on the remaining parameters (say S 2 ) in the high-dimensional space (let S = S 1 US 2 ). This construction is demonstrated by predicting wellbore signatures in the Gulf of Mexico (GoM) where a 120-dimensional table is populated at several 5 thousand wellbore locations throughout the GoM. Physics-based models of multiphase ow in porous media are used to calculate the PCE representations at the sites where data is available. In this context, random parameters describing the subsurface dene the parameter set (S 1 ). A Gaussian process model is then developed for each coecient in these representations, construed as a function on S 2 . The combined probabilistic rep- resentation permits the delineation of separate stochastic in uences on predictions of interest. 2.1 Introduction Interpolation of knowledge is a ubiquitous challenge in many problems of science and engineering. The form it takes depends on which knowledge is being interpolated and which data is available. In many cases of interest, only statistical knowledge is available. In some other cases, knowledge is predicated on rst principles, embodied in equations of mathematical physics. Recently, an explosion of data acquisition capabilities has enabled the rational fusion of data-driven and model-driven inferences. This chapter addresses the problem of conditional regression of a stochastic pro- cess q(x;y;) based on the available generator functions in the form of PCE of QoI, q i (x i ;y i ; i ) at i = 1; 2 ;n locations. Herex2R 3 refers to spatial locations,y and are, respectively, deterministic and stochastic inputs to the underlying physical process. Our interest is to estimate generator functions of q s (x s ; y s ; s ) at a location "x s " with partial or no information ony s and marginal densities at n spatial locations. One approach to this problem is to pose this question as a conditional density regres- sion problem, which is to estimatef qs (q s jfq i g n i=1 ). Heref qs refers to marginal probability density function ofq at locations, andq i refers to samples ofq at the locationi, obtained from the generator function. For this approach, methods for simulating conditional sta- tionary non-Gaussian process has been proposed in the literature [3, 4]. These methods compute a non-linear transformation of the underlying Gaussian process using iterative 6 procedures to match the spectral density function of the process or its correlation func- tion. Recently, Sakamoto and Ghanem proposed a method for simulating non-stationary and non-Gaussian process using polynomial chaos expansion (PCE) and Karhunen-Lo eve representation [5]. More recently, density regression techniques using Dirichlet process are proposed [6, 7]. In this chapter, stochastic processes of interest describe wellbore signatures in the Gulf of Mexico (GoM) and exhibit dierent character of statistical uctuations at distinct spatial scales. Knowledge of these processes is gained by collection and interpretation of various datasets such as sonar imagery, digital seismic-re ections, and well-logs. The uncertainty of this knowledge will be reduced over the time, as more data is collected and assimilated into the database. Also, these wellbore signatures are governed by micro (reservoir) scale, meso (eld) scale and macro (GoM) scale processes. Our objective is to develop a method for conditional regression of these stochastic processes with two main features. First, it should leverage the multiscale dependence of QoI on available knowledge. Such a hierarchical model will enhance the knowledge on dierent scales of uncertainty, the value in reducing these uncertainties and the worth of associated information. Second, it should use all the current assimilated knowledge, this allows reducing the uncertainty of the QoI as more data is available to the decision maker. Section 2.2 introduces the motivation and challenges associated with the problem, Section 2.3 summary of the methodology and Section 2.4 numerical demonstration of the proposed method in the GoM. 2.2 Motivation and challenges An increase in the world hydrocarbon consumption pushed the mankind to explore new sources of hydrocarbons such as deep and ultra-deep oshore regions in the GoM. The 7 NETL/BOEM Geospatial database Geological attributes Example : Spa- tial references Macro scale attributes Example : Initial pressure Average micro scale attributes Example : Permeability Standard nomen- clature attributes Example : API Number Production attributes Example : GOR Figure 2.1: Subgroups in BOEM Database recent Deepwater Horizon oil spill highlighted the uncertainty associated with predic- tions of wellbore signatures due to anthropogenic events in hydrocarbon production in these regions. During this event attempts were made to control and cease the blowout. The scarcity of data, as well as uncertainty in leakage and seepage mechanisms, made associated predictions challenging. Several investigations have been conducted to model the blowout and estimate the damage. It took almost two months to get reliable wellbore signatures which were used to estimate the degree of contamination in the GoM [8]. Wellbore sand data is released annually with a 3-year data moratorium. A repository of wellbore data from the GoM is curated by the Bureau of Ocean and Energy Manage- ment (BOEM) and can be queried online (http://www.boem.gov). This data pertains to reservoirs that are currently in production and no longer in production. Figure 2.2 shows the spatial distribution of one of such available databases. Table 2.1 shows few attribute names (entire list is shown in Appendix A) in the dataset. Attributes in the data set are divided into geological, macro scale, average microscale, standard nomenclature and production attributes as shown in Figure 2.1. This data is complex and features dierent types of uncertainty at the micro (reser- voir) scale, meso (eld) scale and macro (GoM) scale. Though, more data is available to the private operators, publicly accessible data is very limited. In the current study, we restrict ourselves to BOEM database, which adds more uncertainty to the QoI. The interplay between data entries, spatial locations, subsurface physical and geophysical 8 processes is quite complex. Attempting to predict all aspects of oil & gas production throughout the GoM is a task that, we believe, cannot be carried out with any level of credibility. This chapter focuses attention its to specic QoI because of the complexities and uncertainties associated with the data. This will simultaneously help to develop reduced order models (ROM), which are better and more eciently adapted to specic objectives. Constructed ROM has the hierarchical structure, that enhances our understanding of the role of dierent scales of uncertainty play in reducing predictability, as well as the value in reduction of each of these uncertainties and the worth of associated information. This ROM can be updated with currently assimilated knowledge, which reduces the uncertainty of QoI as more data is available to the decision maker. It can rapidly estimate wellbore signatures or QoI based on the currently available database of wellbore data over the GoM. It also serves as a rapid risk assessment tool in case of anthropogenic events such as Deepwater Horizon oil spill, anomaly detection tool by doing hypothesis testing with normal conditions and a stochastic prior model for the sites with limited data. The objective here is to develop a ROM to estimate a QoI at sites where partial or no data has been collected. As a rst step, the QoI is computed at sites where the data is available by evaluating a reservoir simulation code that takes a subset of the data from BOEM with predominant local in uence, such as permeability. Typically, measurements are taken from multiple wellbores within a reservoir and a model for permeability as a stochastic process is constructed using geostatistical methods. Although, permeability data with higher resolution is available to private operators throughout the Gulf of Mexico, publicly available data is limited to the BOEM dataset which provides a single permeability value averaged along the length of the wellbore. Sparsity and uncertainty levels associated with the measured data makes this task of regression quite challenging. Here, two ROMs to transport the evidence available through the BOEM databases to credible uxes at the sea oor are developed. These ROMs pertain, respectively, to ow behavior within a reservoir, and ow through the wellbore. The reservoir ROM 9 Figure 2.2: Spatial location of boreholes over the Gulf of Mexico; Source of the data [1] (RROM) provides a stochastic representation of ows within the reservoir at a target location where BOEM data is not available. It essentially aims to interpolate the BOEM database to that target location. A challenge in achieving this pertains to the fact that the entries in the BOEM database perceive space and time dierently and cannot nec- essarily be interpolated with uniformly credible accuracy. We thus rely on models that, in a rst step, transform some of the BOEM data (such as permeability) into derived quantities (such as ow within reservoir), taking into account the physical processes in the neighborhood of each BOEM wellbore. We then interpolate those derived quantities, across the GoM, conditional on other parameters from the BOEM dataset. To account for uncertainties stemming both from missing data or the missing physics, both the lo- cal models and the GoM-wide interpolation are inherently probabilistic with statistics that re ect available information. The RROM is developed by synthesizing statistical realizations from stochastic reservoir simulation (using polynomial chaos expansion with a reservoir simulator) models and multi-dimensional Gaussian processes. The second 10 Attribute name Attribute name Attribute name Attribute name Sand sequence number Chronozone Total area Oil API gravity Sand name Play type Total volume Initial gas formation volume factor Assessed Proved or Unproved Oil average net thickness Initial oil formation volume factor Sand discovery date Sand type Oil total area Recoverable oil per acre-foot Sand discovery year Water depth Oil total volume Recoverable gas per acre-foot Sand discovery date high Proved recoverable oil Gas average net thickness Oil in place Sand discovery year high Proved recoverable gas Gas total area Gas in place Discovery well Proved recoverable BOE Gas total volume Oil recovery factor Field name Cumulative oil produced Drive type Oil reservoir's recoverable oil Reservoir type Porosity Water saturation Permeability Initial pressure Initial temperature Sand pressure gradient Sand temperature gradient Initial solution gas-oil ratio Yield Proportion oil Gas-oil ratio Table 2.1: Few attributes in BOEM database ROM, the Wellbore ROM (WROM) characterizes the multiphase ow at the outlet of a drillpipe. Conditions within the reservoir at the location of the pipe are rst obtained by running the RROM. The WROM accounts for possible loss of integrity along the wellbore by including provisions for a random pipe diameter. The range for the diameter is obtained by running several detailed ow analysis corresponding to dierent owpaths within the wellbore. For each such analysis, an equivalent \healthy" pipe was computed with a corresponding pipe diameter. The input for the WROM are few outputs from the RROM and also details about the pipe geometry such as nominal drillpipe diameter and length. 11 2.3 Summary of the Methodology In the present setting, the interplay between data entries, spatial locations, subsurface physical and geophysical processes is quite complex. Attempting to predict all aspects of oil & gas production throughout the GoM is a task that, we believe, cannot be carried out with any level of credibility. To mitigate the challenge presented by this complexity and associated uncertainties, we focus attention on specic QoI. This will simultaneously help us develop ROMs that are better and more eciently adapted to specic objectives. 2.3.1 Reservoir QoI In designing the QoI, it is noticed that spatially resolved data, at the reservoir scale is not available. Thus characterizing the heterogeneity within the reservoir would introduce a signicant uncertainty into the prediction process, signicantly reducing the credibility of the predictions. Assuming a homogeneous medium within each reservoir, while highly unrealistic, only carries implications for long times at which the medium away from the well begins to feel the pressure from the well. For early times, the assumption of homogeneity is suciently justied, as it needs parameters only in the immediate vicinity of the wellbore. QoI for the reservoir ROM (RROM), therefore, was designed to re ect the conditions in the reservoir very soon after the drillpipe intercepts the reservoir. While this may seem to be a signicant restriction on our predictive capabilities, we note that it is indeed in these stages of reservoir exploration that issues of uid stability and well- integrity are most signicant. Furthermore, as the time progresses, uncertainty about local conditions decreases with direct access to local data. For these reasons, a QoI that consists of cumulative productions of gas and oil, during early stages of the production is used. These are obtained by solving uid ow (Darcy's) equations along with conservation laws. Since BOEM data have insucient data, to 12 simulate detailed compositional reservoir models, following block-oil partial dierential equations [9] are solved numerically to obtain QoI, @ @t Ws B w S w +r Ws B w ~ u w q W = 0 (2.1) @ @t Os B o S o +r Os B o ~ u o q O = 0 (2.2) @ @t R so Gs B o S o + Gs B g S g +r R so Gs B o u o + Gs B g u g q G = 0 (2.3) Here denotes the porosity;S ; ;p ;B ; s ;q are the saturation, viscosity, pressure, formation volume factor, density and volumetric rate of each phase (water, gas, and oil); R so is the gas solubility; u w ; u g ; u o are the Darcy velocities obtained by solving following the Darcy's law for each phase u = k r (rp rz) =o;w;g; (2.4) where ; k r is the mass density and relative permeabilities of each phase, is the magnitude of the gravitational acceleration, and z is the depth. Corresponding volumetric ow rates and cumulative productions are obtained using Q = ZZ A u dA =o;w;g; (2.5) C = Z t Q dt =o;w;g: (2.6) Here Q o ; Q w ; Q g are volumetric ow rates and C o ; C w ; C g are cumulative productions of oil, gas and water; A is the cross-sectional area, and t is the time. From here onwards, C o ; C w ; C g are reservoir QoI and symbol "q" is used to represent these quantities. 2.3.1.1 PCE of reservoir QoI Uncertainty associated with permeability, porosity and some other parameters renders the QoI at each site random. PCE of QoI are constructed to represent these random 13 quantities, using numerical reservoir simulations with Eclipse [10], a commercial reservoir simulator. The PCE represents the QoI, q, as a multidimensional polynomial in the input random variables [11, 12] with deterministic coecients, p th order PCE of q can be written as q p (x;y;) = X j q j (x;y) j (): (2.7) Here () are the orthonormal Hermite polynomials up to order p; is the vector of independent standard Gaussian random variables that represent the uncertain input pa- rameters to the numerical reservoir model; x2 R 3 refers to spatial locations, y is the vector of deterministic inputs to the numerical reservoir model; q j (x;y) are determinis- tic coecients of the expansion computed using non-intrusive Galarikin projection and quadrature rule, q j (x;y) = <q p (x;y;) ()> < () ()> = nq P k=1 q p k (x;y; k ) j ( k )w k nq P k=1 j ( k ) j ( k )w k ; (2.8) q p k (x;y; k ) is the QoI evaluated at k'th quadrature point k , and nq is the number of quadrature points. 2.3.2 Wellbore QoI The QoIs for wellbore ROMs (WROM) are sea oor uxes obtained using well in ow performance relationship (IPR) and tubing performance relationship (TPR) curves. In addition to sea oor uxes, it is important to be able to characterize loss of integrity in the wellbore following an accident. So, QoI can be specied either as ow rates at sea oor conditioned on a specic pipe geometry (in which case we are certain of the pipe condition), or ow rates at sea oor under conditions of uncertain pipe diameter. The range of pipe diameters investigated is obtained by computing equivalent "healthy" pipes corresponding to a statistical distribution of breaches along the length of the casing. 14 In this manner, following an accident, it may be possible to determine the most likely owpath associated with a given level of observed uxes at the sea oor. 2.3.2.1 Well in ow performance relationship (IPR) curve IPR curve is used in production engineering to obtain ideas on the well performance by plotting the ow pressure (aka bottom-hole pressure) versus the ow rate. It is also used to model the uid ow from the reservoir, through the formation, and into the well. Here, we use the linear well productivity index (PI) model and the Vogel model [13] for a pressure range above and below the bubble point. The PI model for IPR curve is given by, Q L =J L (p ws p wf ): (2.9) Here, Q L is stock-tank oil rate, p ws is well static (or reservoir) pressure and p wf is well owing (or bottom hole) pressure. The Vogel model for IPR curve is given by, Q =Q max 1 (1C) p wf p ws C p wf p ws 2 ! : (2.10) Here, Q is liquid ow rate, Q max is the absolute open hole ow potential (the ow rate when the bottom hole pressure is 0), an output from the reservoir simulation, one of the reservoirs QoI; p ws is the well static (or reservoir) pressure, wf is well owing (or bottom hole) pressure andC is the Vogel coecient (0.8). The corrected Vogel equation for bottom hole pressure smaller than bubble point pressure is as follows, QQ bp =Q max 1 (1C) p wf p ws C p wf p ws 2 ! ; (2.11) where Q bp is the ow rate at the bubble point pressure. 15 Figure 2.3: An example of IPR and TPR curves used for the simulations. 2.3.2.2 Well Tubing Performance Relationship (TPR) curve TPR curve is calculated using Hagedorn and Brown correlation model [14], this cor- relation is used for pressure loss and hold up. Using this model, pressure gradient is calculated as follows: 144 p h = g g c [ L H L + g (1H L )] 1 + f( v SL + v SG ) 2 2gD + ( v SL + v SG ) 2 2gh : (2.12) Here p is the pressure (psi), h is the depth (ft); g is the gravitational acceleration, ft/sec 2 , L is the integrated average density of liquid at owing condition (lb/ft 3 ); g is the integrated average density of gas at owing condition (lb/ft 3 ); H L is the liquid- holdup factor; f is the Darcy-Weisbach friction factor; v SL is the average velocity of supercial liquid at owing condition (ft/sec); v SG is the average velocity of supercial gas at owing condition (ft/sec), andD is the pipe diameter (ft). The intersection point of the IPR and the TPR curves as shown in Figure 2.3 gives the well owrate q w . 16 2.3.2.3 PCE of wellbore QoI Uncertainty associated withQ max and the diameter (representation of the failure scenario will be detailed in Section 2.4.3) renders the wellbore QoI W q at each site random. The p th order PCE of W q can be written as W p q (x;y;) = X j W j (x;y) j (): (2.13) Here () are the orthonormal Hermite polynomials up to order p; is the vector of independent standard Gaussian random variables; x2R 3 refers to spatial locations,y is the vector of deterministic parameters to IPR and TPR curves; W j (x;y) are deter- ministic coecients of the expansion, computed using non-intrusive Galarikin projection and quadrature rule, W j (x;y) = <W p (x;y;) ()> < () ()> = nq P k=1 W p k (x;y; k ) j ( k )w k nq P k=1 j ( k ) j ( k )w k ; (2.14) W p k (x;y; k ) and W q 's are the wellbore owrates at quadrature point k computed by constructing the IPR and TPR curves, and nq is the number of quadrature points. 2.3.3 Reservoir reduced order model Reservoir reduced order model is obtained using PCE in conjunction with Gaussian Process Regression (GPR) [15, 16, 17, 18] over a parameter space dened by the content of the database. To compute the RROM, rst PCE of q is computed at all the sites where BOEM data is available. Denoting these number of sites by R, we thus have, q(r; r ) = X j q j (r) j ( r ); r = 1; 2; 3;:::R; (2.15) 17 where r are independentR d -valued Gaussian random variables (d is equal to the number of uncertain inputs to the numerical reservoir model), and are the multidimensional Hermite polynomials. The zeroth order polynomial is equal to 1, signifying that the zeroth order coecient is the average of q. We will assume in the sequel that the ne scale uncertainty, described by the random variables r are identical across all reservoirs. At site "s", knowingfq j (s)g, then we can estimate of q(s) as q(s; r ) = X j q j (s) j ( r ) : (2.16) To get fq j (s)g, an interpolatory model is constructed over the parameter set S 2 (macroscale attributes) using GPR. Interpolation with GPR makes PCE coecients random at all locations where no data was collected, including spatial locations where only partial data is available. This leads to a stochastic expansion (PCE) with random coecients. If some of the data is not available at a spatial site of interest, this data is itself interpolated using GPR over the spatial domain (R 3 ), thus adding another level of uncertainty in the model leading to a stochastic expansion with three levels of un- certainty. Since S 2 is not observed at the site "s," its content is interpolated over GoM with respect to the parameters in the set S 1 (Geological attributes) and with respect to geospatial references (referred to as set D 2 ). This interpolation task is denoted symboli- cally by, ^ S 2 (s;S 1 ;D 2 ) =I g (s;S 1 jS 1 ;D 2 ); (2.17) whereI g is an interpolation operator, which in our case is the Gaussian process regres- sion. This leads to a PCE representation of entries in ^ S 2 , ^ S 2 (s; g ) = X j S 2j (s) j ( g ) ; (2.18) where g are new Gaussian random variables that re ect the fundamental uncertainty that controls variability at the scale of the Gulf. These random variables are assumed 18 to be independent of the set r . The PCE coecients q j (s) are thus interpolated as a function of the interpolated ^ S 2 and the PCE coecients obtained at all the observed reservoir sites q j (r). This is expressed as, ^ q j (^ s 2 (s; g );q j (r)) =K g (q j (s)j^ s 2 (s; g );q j (r)) (2.19) whereK g is an interpolation operator, which in this case is again given as Gaussian process regression over the parameter space. This can be described in the following PCE, ^ q j (D s 2 (s; g ); f ) = X i q ij (D s 2 (s; g )) i ( f ) (2.20) in which a new set of f , "fundamental" random variables are introduced. Putting all the representations together leads to ^ q (s; r ; f ; g ) = X i;j q ij X k s 2 k (s) k ( g ) ! i ( r ) j ( f ) : (2.21) where the parentheses on the right-hand side refers to functional dependence which can further be unfolded yielding the following expression, ^ q (s; r ; f ; g ) = X i;j;k q ijk k ( g ) i ( r ) j ( f ) : (2.22) This model has three levels of uncertainty, which arise from three dierent physical scales. The rst level of uncertainty is in the reservoir model because of the lack of detailed information on the microscale features such as permeability. This uncertainty can be reduced by having detailed information of microscale features, i.e. by having more wellbores inside the reservoir. The second level of uncertainty is due to the uncertainty of local geological features such as initial pressure, which can be reduced by having an understanding of the geology and the variability of the parameter set D 2 in local scale, i.e. inside the eld where the reservoir "s" is located. The third level of uncertainty is 19 attributed to geological uctuations across the GoM, which can be reduced by detailed information on the geology of the entire GoM. If we know S 2 at site "s," the above model reduces to ^ q (s; r ; g ) = X i;j q ij (s) i ( r ) j ( f ): (2.23) Note that in case if the BOEM database is updated, PCE in the Equation 2.15 can be used as a surrogate model (assuming updated knowledge reduced the uncertainty in the input parameters) for the reservoir model and PCE will be updated to q(r; r ) = X j q j (r) j ( r ); r = 1; 2; 3;:::R (2.24) by updating the PCE coecients from q j to q j , then, the updated PCE will be used to build the RROM. 2.3.4 Wellbore reduced model Once RROMq (s; r ; f ; g ) computed WROM at location "s" is obtained using PCE as W q (s;y; r ; f ; g ; w ) = X j W j (s;y; r ; f ; g ) j ( w ): (2.25) Here j ( w ) are the orthonormal Hermite polynomials; w is the vector of indepen- dent standard Gaussian random variables that represent the uncertainty in the wellbore model; y is the vector of deterministic parameters to IPR and TPR curves at loca- tion "s"; W j (s;y; r ; f ; g ) are random coecients of the expansion, computed using non-intrusive Galarikin projection and quadrature rule W j (s;y; r ; f ; g ) = <W (s;y; r ; f ; g ) ( w )> < ( w ) ( w )> = nq P k=1 W k (s;y; r ; f ; g ) j ( w k )w k nq P k=1 j ( w k ) j ( w k )w k ; (2.26) 20 W k (s;y; r ; f ; g ) is the wellbore owrates at quadrature point w k computed by con- structing the IPR and TPR curves, and nq is the number of quadrature points. 2.4 Numerical simulations Figure 2.4: Spatial location of boreholes over the Gulf of Mexico and dierent groups in the datasets: red circles, yellow triangles, and green stars respectively indicate nal test set, initial test and training reservoirs. To show the capability of the developed procedure, simulations are carried at few elds in the GoM. These elds are divided into training, initial testing and nal test- ing reservoirs, which consists of 131, 16 and 22 reservoirs respectively. Figure 2.4 shows the spatial location of boreholes over the GoM. Here red circles, yellow triangles, and green stars respectively indicate nal test set, initial test and training reservoirs. Training reservoirs fall in Gunnison, Popeye, Bullwinkle, Neptune (At), Allegheny, Nansen, Pom- pano, Troika, Genesis, Brutus, Jolliet, Mars-Ursa, and Front Runner elds. Initial test reservoirs fall in Devils Tower, and Prince elds which are used for the model selection of GPR and nal test reservoirs fall in Auger, and Thunder Horse elds. These reservoirs demonstrate the predictive value of our ROM methodologies. More specically, we will 21 demonstrate using the given data from various wells located in the training and initial test reservoirs, how we can anticipate ows at some arbitrary well, "s" in the test reser- voirs. The credibility of our predictions will depend on the type of information available within training sites as well as on the information available for well "s" in testing sites. Information about "s" consists of spatial coordinates, but could also include additional geological, geophysical and subset of parameters of S 2 . 2.4.1 PCE of QoI Numerical reservoir simulations are carried out using Eclipse [10], an inclusive numerical simulator to model dierent scenarios of uid ow in the oil and gas reservoirs. Since BOEM data do not yield sucient information to allow for numerical reservoir sim- ulations, the following additional assumptions are made to enable numerical reservoir simulation, Each reservoir is assumed to be bounded by a cube with dimensions specied by the area of the reservoir, and the average thickness (both available in the BOEM dataset), which is consistent with our stated QoI. Permeability inside the eld is considered as a random variable. Using the samples of permeability in the eld, a pdf of this random variable was obtained. Permeabil- ity in each reservoir is assumed to follow this distribution but shifted towards the mean of that reservoir. The pdf of centered permeability for few elds generated by this procedure is shown in Figure 2.5(a). Similar procedure is followed to generate a distribution of the porosity. The pdf of centered porosity for few elds generated by this procedure is shown in Fig- ure 2.5(b). Permeability along the z direction is considered to be 0.1 to 0.9 times of perme- ability along the x direction. 22 No. Random variable Distribution 1 Permeability along the X and Y directions Generated by the measurements in the eld 2 Permeability along the Z direction Permeability along the X direction * U(0.1, 0.9) 3 Porosity Generated by the measurements in the eld 4 Bottom hole pressure U(0.43 psi/ft, 1 psi/ft) 5 Gas saturation rate U(0.6017, 1.027) 6 Irreducible water saturation U(0.1, 0.4) Table 2.2: Random input parameters used for the reservoir model. Bottom hole pressure is assumed to be a uniform random variable with lower and upper limits as hydro-static (1psi/ft) and litho-static (0.433 psi/ft) pressures based on [19]. PVT data was obtained from the correlations developed in [20] with random gas saturations varying from 0.6017 to 1.0270. Relative permeability coecients are calculated using Wyllie and Gardner corre- lations with irreducible water saturation (S wc ) varying from 0.1 to 0.4. These assumptions lead to the reservoir model with random input parameters listed in Table 2.2. Using these assumptions, PCE of dimension 6 and order 4 is created for QoI. Probability distribution functions (pdf) of QoI for some of the training reservoirs are shown in Figure 2.6. 23 −0.5 0 0.5 0 2 4 Centered Log Porosity Density −0.5 0 0.5 0 5 10 Centered Log Porosity Density −0.5 0 0.5 0 2 4 Centered Log Porosity Density −0.5 0 0.5 0 2 4 6 Centered Log Porosity Density −0.5 0 0.5 0 2 4 Centered Log Porosity Density −0.5 0 0.5 0 2 4 6 Centered Log Porosity Density −0.2 0 0.2 0 5 10 Centered Log Porosity Density −1 0 1 0 1 2 3 Centered Log Porosity Density −0.5 0 0.5 0 5 10 Centered Log Porosity Density (a) −0.5 0 0.5 0 2 4 Centered Log Porosity Density −0.5 0 0.5 0 5 10 Centered Log Porosity Density −0.5 0 0.5 0 2 4 Centered Log Porosity Density −0.5 0 0.5 0 2 4 6 Centered Log Porosity Density −0.5 0 0.5 0 2 4 Centered Log Porosity Density −0.5 0 0.5 0 2 4 6 Centered Log Porosity Density −0.2 0 0.2 0 5 10 Centered Log Porosity Density −1 0 1 0 1 2 3 Centered Log Porosity Density −0.5 0 0.5 0 5 10 Centered Log Porosity Density (b) Figure 2.5: Distribution of the centered (a) log permeability (b) log porosity for few of the elds under study. 24 0 1 2 3 4 5 6 7 x 10 4 0 1 2 x 10 −4 QoI: Oil production rate on day one in bbl/day pdf of QoI Figure 2.6: Probability distribution function (pdf) of oil production rate for few training reservoirs. 25 2.4.2 Reservoir ROM Assuming only partial information of attributes regarding well "s" at the testing reser- voirs is available, QoI at this site is predicted by using PCE coecients from the training reservoirs. GPR of PCE coecients is performed over the subset of S 2 parameter space since these parameters are used in numerical reservoir simulations. Here, we used subset formed by porosity, water saturation, initial pressure, thickness of the reservoir, oil API and area of the reservoir. If these parameters are known at well "s", they are directly used for GPR, otherwise, they are computed using GPR by interpolating over the spatial reference and used for the interpolation of the PCE coecients. 2.4.2.1 RROM with known parameters at the target site Gaussian Process regression of PCE coecients is carried over the parameters porosity, water saturation, initial pressure, thickness of the reservoir, oil API and area of the reservoir using Python toolkit Scikit-learn [21]. Model selection for GPR is done by training the models on the training reservoirs and evaluating the performance of the each model on the initial test reservoirs. For example, mean square error (MSE) of the rst polynomial coecient using dierent autocorrelation functions for training and initial test reservoirs are shown in Table 2.3. The squared exponential autocorrelation function is chosen based on MSE metric, the performance of this model on the nal test reservoirs is shown in Figure 2.7. This gure shows that for most of the nal test reservoirs, predicted rst PCE coecient fell within the 95% condence interval calculated from the chosen model. Using the similar procedure, prediction of other PCE coecients is carried out. GP regression gives distribution for each PCE coecient at prediction site, each realization from these distributions gives one pdf for QoI. For one of the wells "s" in Auger eld, the realizations of these pdf are shown in Figure 2.8 (b). The red dashed and dotted line indicates samples of the pdf for QoI, each sample pdf re ects a realization from the GP regression. All the samples from these realizations are collected and single 26 Figure 2.7: Predictions of the rst PCE coecient using GPR with squared exponential autocorrelation at nal test reservoirs. No Kernel MSE on training reservoirs MSE on initial test reservoirs 1 Exponential 9.0330172e-21 6.107499e+06 2 Squared Exponential 2.533343e-17 4.041412e+06 3 Linear 1.974906e-16 5.097380e+06 4 Cubic 2.9689354e-21 5.244330e+06 Table 2.3: Mean square error for dierent models for training and initial test reservoirs. pdf is obtained from these samples which serves as a pdf of QoI at the prediction site. The black dashed line here indicates most likely pdf at the prediction location obtained as an average over GP regression. In Figure 2.8 (a) the solid green line indicates pdf of QoI at the prediction location obtained from available data (actual) at that location. Here, the dashed red line indicates predicted pdf of QoI using GP regression over the parameter space. 27 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 0 0.2 0.4 0.6 0.8 1 1.2 x 10 −4 QoI: Oil flow rate on day one in bbl/day Probabilty density of QoI Predicted pdf of QoI Actual pdf of QoI (a) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x 10 −4 QoI: Oil flow rate on day one in bbl/day Probabilty density of QoI Relizations of predicted pdf’s QoI Predicted most probable pdf of QoI (b) Figure 2.8: QoI: Flow rate on day one in bbl (a) Solid green line indicates pdf of QoI at prediction location obtained from available data at that location; red dashed line indicates predicted pdf of QoI using Gaussian Process regression over parameter space; (b) red dashed and dotted line indicates samples of the pdf for QoI, each sample pdf re ects a realization from the Gaussian Process regression over parameter space; The black dashed line indicates most likely pdf at the prediction location obtained as an average over Gaussian Process regression. 2.4.2.2 RROM with unknown parameters at the target site Since parameters porosity, water saturation, initial pressure, the thickness of the reser- voir, oil API and the area of the reservoir at the well prediction locations are not available, they are obtained by doing GPR over the spatial reference. This makes these parame- ters at this location random, for each realization of these parameters, regression of PCE coecients are carried over the parameters similar to Section 2.4.2.1. For one of the wells "s", realizations of these pdf of QoI are shown in Figure 2.9 (b) the dashed and dotted lines indicate samples of the pdf of QoI, each sample pdf re ects a realization from the Gaussian Process regression over the interpolated parameter space. Here, each color in dotted line re ects pdf of QoI, corresponding to dierent realizations of the parameters. All the samples from these realizations are collected and single pdf is obtained from these samples which serves as a pdf of QoI at the site "s". In Figure 2.9 28 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 0 0.5 1 1.5 x 10 −4 QoI: Oil flow rate on day one in bbl/day Probabilty density of QoI Actual pdf of QoI Predicted pdf of QoI (a) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 4 0 0.5 1 1.5 2 2.5 3 3.5 x 10 −4 QoI: Oil flow rate on day one in bbl/day Probabilty density of QoI (b) Figure 2.9: QoI: Flow rate on day one in bbl (a) Solid green line indicates pdf of QoI at prediction location obtained from available data at that location; Dashed red line indi- cates predicted pdf of QoI using Gaussian Process regression over interpolated parameter space; (b) Dashed and dotted lines indicate samples of pdf of QoI, each sample pdf re- ects a realization from the Gaussian Process regression over the interpolated parameter space; Each color in dotted line indicates dierent realization of the parameters; (a) solid green line indicates pdf of QoI at the prediction location obtained from available data (actual model) at that location. The dashed red line indicates predicted pdf of QoI using GP regression over the parameter space. Table 2.4 shows few statistics from the actual model, predictions using GPR with known & unknown parameters and preliminary estimate models (obtained by collecting all the realizations from training reservoirs). Table 2.4 and Figure 2.8 & Figure 2.9 show that our proposed methodology is able to provide reliable probabilistic estimates of QoI at the target sites. 29 Model Mean of QoI (in bbl) Standard deviation of QoI Actual model 9680 9724 Prediction model with known parameters 10566 11029 Prediction model with unknown parameters 6214 7620 Preliminary estimate model 3898 17753 Table 2.4: Statistics of the ow rate on day one (QoI) from the actual model, predictions using GPR with known & unknown parameters and preliminary estimate models. 2.4.3 Wellbore ROM Two reduced order models are pursued for the well ROM (WROM), each of which is adapted to a particular QoI. The rst WROM aims to estimate the ow at the sea oor that is consistent with conditions within the reservoir estimated using the RROM de- scribed above. The Hagerdon-Brown correlation model is adopted for this WROM. The operating conditions for this model are at the intersection of an IPR (In ow performance relation) and a TPR (Tubing performance relation). This condition is evaluated for a statistical ensemble describing the output of the RROM at the target site. Wellbore geometry (pipe length and diameter) also needs to be specied in this calculation. We provide two means for specifying these quantities. First, these quantities may be known from current operating conditions, in which case they are input to the executable. Al- ternatively, these conditions may yet be unknown. In this case, we also allow for the eventuality of a breach occurring somewhere along the length of the wellbore, therefore aecting the owpath within it. A statistical model is thus adopted for the well diameter (under conditions of known length). PIPESIM [22] is used to evaluate the ow rates through the pipe, with a breach occurring at some uniformly random height along the pipe. An equivalent "healthy" pipe is then estimated through an optimization task. A second WROM consists of essentially the same building blocks as the rst WROM. Its use is essentially dierent. Specically, the WROM allows the evaluation of PDF for multiphase ow at the sea oor for dierent conditions of pipe integrity. A test of 30 4500 5000 5500 6000 6500 7000 7500 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x 10 −3 QoI: Flow rate of oil in bbl/Day Probabilty density of QoI Figure 2.10: Probability Density Function for various quantities of interest at sea oor; undamaged pipe. hypothesis is then conducted for observed ows against each of these PDFs. This permits the rapid identication of breach location that is most consistent with observed ows. 2.4.3.1 Normal scenario This scenario corresponds to the normal operating conditions and variability in sea oor predictions stems from variability in predicting reservoir conditions, the ow through the pipe being deterministic. Figure 2.11 shows the snapshot of the corresponding model in PIPESIM, and Figure 2.12 shows the schematic of ow in drillpipe system. Figure 2.10 shows pdf of the cumulative oil ow rate after day one at the sea oor level. 31 Figure 2.11: Normal scenario: Geometry of the pipe system. 2.4.3.2 Damaged scenario with known damage location In this scenario, it is assumed that there is a leakage in the wellbore with known damage location. Here, owrate corresponds to two outlet pressures and deterministic reservoir (mean) conditions are assumed to demonstrate the eect of the damage location on the sea oor uxes. In case one, the outlet pressure is equal to the pressure of water column above the well head and the second one is the same water pressure plus a pressure drop value (assumed to be 300 psi), (P O ) 1 = P w = 4:289 0:43 = 1; 844psi (P O ) 2 = P w + P = 1; 844 + 300 = 2; 144psi: It should be noted that the drilling pipe diameter is assumed to be 5 inches, and the pipe system is thought to have two sets of casings, one with diameter of 12 inches and the other of 10 inches. A schematic, in Figure (2.12), drawn for scenario 3-b illustrates the physics of the problem. The damage scenarios are as described below: Scenario 1 (lower limit): Annular ow between the 10.5" casing and 12" liner. 32 Figure 2.12: Schematic of ow in drillpipe system. Scenario 2 (upper limit): Tubing ow through the drill pipe and annular ow be- tween the drill pipe and casing (drill pipe is modeled as a 5" tubing inside the 10" casing). Scenario 3 (intermediate case): This scenario is divided into two sub-scenarios with dierent breach locations: (a) Breach located at 100 ft below the drill pipe (dominant ow is the annular one), uid ows up in the area between 10.5" casing and 12" liner up to the breach point, afterwards it goes into the 10.5" casing and goes up (both in the casing and drill pipe). (b) Breach located at 100 ft above the perforations (100 ft above the bottom), the ow path is similar to case (a). In terms of behavior, case (a) shows a close relation to scenario 1, whereas case (b) behaves similarly to scenario 2. The pipe congurations for damaged scenarios are modeled with PIPESIM, as shown in the Figures (2.13), (2.14), (2.15) and (2.16), respectively. 33 Figure 2.13: Failure scenario 1: Geometry of the pipe system. Figure 2.14: Failure scenario 2: Geometry of the pipe system. The cumulative oil ow rate at the end of the day one, for each scenario is computed by doing the nodal analysis with PIPESIM and the results are shown in Table 2.5. Finding diameter for normal operation equivalent to the damaged scenarios Here, rst the damage location or the scenario is assumed and the ux at the sea oor is computed using the PIPESIM. Then, IPR and TPR curves are calculated for dier- ent diameters assuming no leakage in the system. Using this procedure, an equivalent 34 Figure 2.15: Failure scenario 3-a: Geometry of the pipe system. Figure 2.16: Failure scenario 3-b: Geometry of the pipe system. diameter has chosen such that resulting ow rate is the same as the damaged scenario. Table 2.6 and 2.7 shows the eective diameters for the scenarios 1, 2 and 3 used in the previous section. Similarly, assuming the ratio of the breach location to the pipe length is uniform random variableU[0; 1]. The equivalent diameters for the normal scenario are calculated for the P out 1844 psi and 2144 psi. These results are shown in Table 2.8 and 2.9. 35 P 0 = 1; 844psi Q(stb=day) Scenario 1 3661 Scenario 2 5208 Scenario 3-a 3737 Scenario 3-b 4864 P 0 = 2; 144psi Q(stb=day) Scenario 1 3255 Scenario 2 4530 Scenario 3-a 3221 Scenario 3-b 4256 Table 2.5: Flow values for dierent breach scenarios. Equivalent Normal Scenario Failure Scenario Flow Rate (stb/d) Diameter (inch) Flow Rate (stb/d) SC1 (annular) 3,661 2.06 3,650 SC2 (all inside the casing) 5,208 5.48 5,208 Table 2.6: Equivalent normal scenarios for failure scenarios 1 and 2. Scenario 3: annular+tubing Equivalent Normal Scenario Breach point MD (ft) Flow Rate (stb/d) Diameter (inch) Flow Rate (stb/d) 10624 4864 3.08 4869 9624 4587 2.69 4591 8624 4356 2.48 4359 7624 4158 2.34 4162 6624 3987 2.23 3980 5624 3836 2.15 3832 4889 3737 2.10 3733 Table 2.7: Equivalent normal scenarios for failure scenario 3. 36 P out = 1; 844 psi Breach Loca- tion / Pipe length Breach Depth Flow Path Flow Rate (stb/d) Equivalent Diameter (inch) 0 11,624 Tubing 5,208 10.5 0.1 10,891 Mixed 4,948 3.14 0.15 10,524 Mixed 4,834 2.96 0.2 10,157 Mixed 4,728 2.83 0.25 9,790 Mixed 4,629 2.73 0.3 9,424 Mixed 4,537 2.65 0.35 9,057 Mixed 4,451 2.58 0.4 8,690 Mixed 4,370 2.52 0.45 8,323 Mixed 4,293 2.46 0.5 7,957 Mixed 4,221 2.41 0.55 7,590 Mixed 4,152 2.36 0.6 7,223 Mixed 4,087 2.32 0.65 6,856 Mixed 4,025 2.28 0.7 6,490 Mixed 3,966 2.25 0.75 6,123 Mixed 3,910 2.21 0.8 5,756 Mixed 3,855 2.18 0.85 5,389 Mixed 3,804 2.16 0.9 5,023 Mixed 3,754 2.13 0.95 4,656 Mixed 3,724 2.11 1 No Breach Annulus 3,661 2.10 Table 2.8: Relation between breach location and the ow rate for P out = 1; 844: 37 P out = 2; 144 psi Breach Loca- tion / Pipe length Breach Depth Flow Path Flow Rate (stb/d) Equivalent Diameter (inch) 0 11,624 Tubing 4,530 10.5 0.05 11,257 Mixed 4,422 3.73 0.1 10,891 Mixed 4,324 3.28 0.15 10,524 Mixed 4,233 3.03 0.2 10,157 Mixed 4,146 2.86 0.25 9,790 Mixed 4,066 2.73 0.3 9,424 Mixed 3,990 2.64 0.35 9,057 Mixed 3,920 2.56 0.4 8,690 Mixed 3,853 2.49 0.45 8,323 Mixed 3,789 2.43 0.5 7,957 Mixed 3,729 2.38 0.55 7,590 Mixed 3,672 2.33 0.6 7,223 Mixed 3,618 2.29 0.65 6,856 Mixed 3,566 2.25 0.7 6,490 Mixed 3,516 2.22 0.75 6,123 Mixed 3,469 2.19 0.8 5,756 Mixed 3,423 2.16 0.85 5,389 Mixed 3,379 2.13 0.9 5,023 Mixed 3,336 2.11 0.95 4,656 Mixed 3,308 2.09 1 No Breach Annulus 3,255 2.06 Table 2.9: Relation between breach location and the ow rate for P out = 2; 144: 38 2.4.3.3 Damaged scenario with unknown damage location As described in the previous section, breaches at dierent locations along the pipe restrict the ow and are thus transformed into equivalent pipe diameters. To better understand the eect of changes in pipe diameters on ow conditions at the sea oor, we carried out the parametric study shown in Figure (2.17). The gure depicts the cumulative oil ow rate at the sea oor with dierent pipe diameters. It is clear that changes in pipe diameter have a clear impact on the centrality and the spread of the associated probability density function. This observation implies that changes in well-integrity condition can be readily detected through a standard statistical test of hypothesis. Figure 2.17: Probability Density Function of oil ow rate. Pipe diameter of 7" is consid- ered nominal with other pipe diameters construed as deviations from normal. 2.5 Data analytics on the BOEM database BOEM 2009 database consists of 120 parameters as listed in Appendix A collected from 13056 wellbores. Data analytics techniques were used to understand the underlying struc- ture in the database. Here, we used diusion maps on 36-dimensional parameter space 39 Sub- group Attributes in subgroup 1 Spatial references X, Spatial references Y, RKB Elevation, Bore hole total measured depth, True vertical depth, Surface North/South distance, Surface East/West distance, Bottom North/South distance, Bottom East/West distance, Depth of water, Surface longitude, Bottom longitude, Surface latitude, Bottom latitude. 2 Discovered BOE, Subsea Depth, Total average net thickness, Total area, Total volume, Oil average net thickness, Oil total area, Oil total volume, Gas average net thickness, Gas total area, Gas total volume, Porosity, Water saturation, Initial pressure, Initial temperature, Sand pressure gradient, Sand temperature gradient, Initial solution gas-oil ratio, Gas specic gravity, Oil API gravity, Initial gas formation volume factor, Initial oil formation volume factor Table 2.10: Parameters used for the data analytics Cluster Production type of reservoirs Number of reservoirs in cluster 1 Only oil 1940 2 Both oil & gas 2474 3 Only gas 8642 Total=13056 Table 2.11: Number of reservoirs in each cluster of the BOEM data (listed in Table 2.10) using a Gaussian kernel with = 1 to identify the underlying database. Based on the eigenvalue decay diusion kernel, dimension is reduced to 15. To visualize manifold, diusion coordinates corresponding to rst two dimensions are plotted in Figure 2.18. This plot shows that there are three clusters present in the data. Reservoirs corresponding to each cluster are identied based on simple linear separation. Careful observation of data in each cluster lead to the conclusion that these clusters correspond to reservoirs producing gas, oil and oil & gas. Table 2.11 shows the number of reservoirs in each cluster. A sample of 65 reservoirs was used to test the classica- tion accuracy, which was 100% for this set since there is a clear separation between the clusters. 40 −4 −3 −2 −1 0 1 2 −6 −4 −2 0 2 4 Φ 1 (x) Φ 2 (x) Only oil reservoirs Both oil and gas reservoirs Only gas reservoirs All 65 reservoirs Figure 2.18: Plot of 36D BOEM data embedding into the rst two diusion map coor- dinates. 41 2.6 Conclusions A procedure has been developed for conditional density regression which produced a reduced order model useful for interpolating wellbore signatures over the Gulf of Mex- ico. These are useful for rapid risk assessment in case of anthropogenic events such as the Deepwater Horizon oil spill. The construction uses polynomial chaos expansion in conjunction with a Gaussian process regression and a database of wellbore data over the GoM. The nal model has three levels of uncertainties which are arising from three dierent physical scales. Despite these uncertainties, the proposed model is able to pro- vide reliable probabilistic estimates of wellbore signatures at unobserved sites. This is demonstrated by carrying out numerical simulations reservoirs in a few elds in the Gulf of Mexico. 42 Chapter 3 High-Dimensional Intrinsic Interpolation using Gaussian Process Regression This chapter introduces a regression procedure for intrinsic variables constrained onto a manifold embedded in an ambient space. The procedure is meant to sharpen high- dimensional interpolation by introducing constraints delineated from within the data being interpolated. Our method augments manifold learning procedures with a Gaussian process regression (GPR). Proposed method rst identies, using diusion maps, a low dimensional manifold embedded in an ambient high dimensional space associated with the data. It relies on the diusion distance associated with this construction to dene proximity between points on the manifold. This distance is then used to compute the correlation structure of a Gaussian process that describes the dependence of quantities of interest in the high dimensional ambient space. The proposed method can be used to predict any quantity on the arbitrary high-dimensional dataset; a specic example is considered to characterize the permeability of rock formations from petrophysics data. 3.1 Introduction As natural resources become scarcer, their recovery becomes more crucially dependent on illuminating the subsurface. While seismic data acquisition, processing, and interpreta- tion are well-established [23], the value of the information it gleans about the subsurface 43 continues to improve as a function of the sensor and computational technologies. New modalities of subsurface exploration, however, solicit coupled interactions, thus compli- cating the task of signal deconvolution and interpretation. The interdependence between these dierent observables carry the signature of multi-scale spatial and temporal scales. Recent data analytic techniques permit such analysis by rst considering an ambient Euclidean space containing the data, and then discovering a low dimensional manifold embedded in this space. Several methods have been developed to detect the low dimen- sional manifold [24], these methods have grown in rigor and sophistication, and they are applied to a broad range of problems from across science and engineering, including image analysis, shape analysis and medical imaging [25, 26, 27]. The present work extends diusion on manifold algorithms to the analysis and in- terpolation of data constrained on the manifold. The proposed method can be used to predict any quantity on the arbitrary high-dimensional dataset; a specic example is considered to characterize the permeability of rock formations from petrophysics data. Point-by-point permeability over the reservoir are relevant several purposes. For example, to develop completion strategies and input to the geocellular model and ow simulation models. Permeability here is dened as being the matrix single-phase absolute perme- ability (as opposed to fracture permeability or multi-phase eective permeability). For a static geomodel, matrix permeability is an input to saturation height models which enables the petrophysicist to create 3D models of hydrocarbon saturation. It is also used to create rock type classes which relate to ow capacity. Coupled with uid saturation, both properties are usually used as a proxy (heuristics) to inll drilling (lateral variations of the sweet spots) [28, 29, 30]. For dynamic ow simulation models, matrix permeability is input to Navier-Stokes equation. It is needed to make the link between total ow rates (barrels/day) and draw- down pressure ( P, dierence between bottom hole pressure and reservoir pressure). In general, estimates of eective permeabilities are required rather than the relative per- meabilities as they are only used by simulation engineers. The permeability will also 44 control critical quantities such as water arrival/breakthrough time [31, 32]. Finally, the lateral distribution of permeability plays an important role in dening the volume of the reservoir swept by water injection and, therefore, impacts recovery factors [33, 34]. Measurements of the permeability are obtained at very few locations using a variety of the techniques such as well-log evaluation, core measurements, and well-testing [35, 36]. Classical approaches to predict the permeability at unmeasured locations use GPR [37] also known as Kriging [15, 16, 17, 18] and the other variants of the Kriging. Even though GPR is the best linear unbiased predictor, as the dimension of the predictors increases, the uncertainty associated with predictions also increases leading to uninformative pre- dictions. Several soft computing methods such as neural network, fuzzy logic, support vector machines [38, 39, 23, 40, 41, 42] were used in the literature to improve the predic- tions. However, unlike GPR, these methods will not provide a full probabilistic predic- tions and an estimate of uncertainty in the prediction at each location. These estimates are necessary to characterize the uncertainty associated with the predictions, and also to build the prior models for the uncertainty calibration based on the measurements of other physical quantities. The proposed method overcomes the problem with the dimen- sionality in GPR, by detecting the lower dimensional embedding associated with the data using manifold detection techniques. Then, GPR is carried out in low dimensional space with intrinsic metric. The paper starts by brie y reviewing Gaussian process regression, and Diusion maps before demonstrating how they are integrated for the characterization of permeability at unexplored site conditional on available eld data. 3.2 Formulation Let ( ;F; P ) be the probability triple, q(x;!) be a random process. Suppose, we have a set of observationsD =fq(x i ;!) atx i g k i=1 , of the quantity of interest q(x;!). Here, x2R n denotesn attributes which may include spatial location and other geological data. 45 Our interest, is to predict the value of q atx =x a R n -coordinate based on available dataD. Classical approaches for predictions of q(x;!) considered in this work such as permeability, uses, GPR with covariates by assuming the q(x;!) as a weak stationary Gaussian stochastic process. In this case, the aim is to nd conditional distributions of q at locationx =x given the dataD, that is specied by p q (q jD). If the underlying random eld is not Gaussian, a common choice is to use a transformed Gaussian random eld model [43]. 3.2.1 Gaussian process regression Gaussian process regression assumes following functional relation between q(x;!) and x, q(x;!) =f(x) +z(x;!): (3.1) Here,f(x) is regression model or trend function inx, andz(x;!) is a zero-mean random process, with zero mean and covariance between any two pointsx;y dened by C(x;y) =Efz(x;:)z(y;:)g (3.2) When the processz(x;!) is second-order stationary, a correlation function ofz(x),R(d), can be introduced in the form 2 R(d(x;y)) =C(x;y) (3.3) where 2 is the variance of the process and d(x;y) is a distance function between two points x and y2 R n . A common choice, for the correlation model, is associated with the Mat ern class [44] with correlation function, R(d) = 1 ()2 1 p 2 d ! K p 2 d ! : (3.4) 46 Where is the gamma function, and are non-negative parameters of the covari- ance and K is the modied Bessel function of the second kind. Parameters and , respectively, control the smoothness of the random process and its decay with respect to distance. The correlation parameter , the variance 2 , and the regression parameter are obtained following a maximum likelihood argument. In the present development, we express the corresponding likelihood function, L(; 2 ;), in terms of its logarithm, in the form lnL(; 2 ;) = 1 2 [n ln( 2 ) + ln(detR D ) + (qf D ) T )R 1 D (qf D ) T )= 2 ]: (3.5) Here, R D is the n x n matrix of correlations of the n design points with R D ij = R(d(x i ;x j )) 1 i n; 1 j n, q =fq(x i ;!)g n i=1 and f D =ff(x i )g n i=1 . Once the parameters ; 2 ; are computed, q atx is estimated as, ^ q(x ;!) = ^ f(x ) +r T (x )R 1 D (qf D ): (3.6) Here,r T (x ) = [R(d(x 1 ;x )); ;R(d(x n ;x ))] is the vector of correlations between the design points andx. It can be shown that ^ q(x ;!) obtained using GPR is the best linear unbiased predictor (BLUP) of q(x ;!) [16]. 3.2.2 Diusion maps A central idea of the present work is to construct a GPR on manifoldM that is embedded within a high dimensional ambient Euclidean space. The manifoldM is identied using manifold learning techniques and specically diusion maps. A m dimensional manifold M is a topological space, in which the neighborhood of each point is homeomorphic to R m . There are a plethora of methods available for manifold detection [24] in the arbitrary dataset. Here, we use diusion maps to identify the underlying manifold and diusion 47 distance as the metric. Diusion maps have the property of preserving the intrinsic (i.e. local) geometry of the data. For parameter = 1 of the diusion maps, it gives optimal embedding, as it approximates the Laplace-Beltrami operator, a generalization of the Laplace operator, operating on surfaces and manifolds. Also, the diusion distance dened using diusion maps is more robust than the geodesic distance to the noise and small perturbations in the data [45]. Diusion maps [46] is based on dening the Markov random walk on the data. By doing the random walk on the data, manifold of the corresponding data is obtained. Identication of low dimensional embedding using anisotropic diusion maps algorithm involves following steps: 1. Based on n points, x (1) ;x (2) ;x (n) , construct the adjacency matrix (W ) of the graph. The components of W are the weights along the edges connecting corre- sponding nodes. For example, the weight of the edge connecting two pointsx i and x j is determined by kernel k(x i ;x j ), W ij =k(x i ;x j ) = exp kx i x j k 2 : (3.7) Herek:k is the Euclidean norm and is the kernel parameter that controls the diusion process or the step of size of the random walk. A relatively small and large values of lead to very little diusion and most diusion process, respectively. In the random walk point of view, choosing small makes random walker wanders around the starting point and large makes random walker moves outside the manifold. So, the value of the should be in between these two extremes, one of the schemes to pick is as proposed by Singer et al. [47]. This scheme computes a parameter L() as a function of dened as L() = n X j=1 n X j=1 W ij (): (3.8) 48 Then will be chosen such that the logarithmic plot between and L() is linear. 2. Renormalize the W to get anisotropic kernel W : w ij = w ij q i q j with q = n X j=1 k ji : (3.9) The choice of the parameter is related to the data point density on the innites- imal transition of the diusion [48]. For example, the cases = 0, = 0:5 and = 1 corresponds to the Markov chain approximates of graph Laplacian normal- ization, Fokker-Planck diusion and Laplace{Beltrami operator, respectively. For = 1, we will recover Riemannian geometry of the dataset for any arbitrary sam- pling distribution. Since, petrophysical data comes from arbitrary distribution, we choose = 1 for our problem. 3. Compute the diusion operator P using, P =D 1 W and D ii = n X j=1 W ji ; (3.10) where the entires of the matrix P represents the transition kernel of the Markov chain. In other words, it represents the transition probability of going from state i to j in one time step. 4. Compute Eigen decomposition diusion operator after t times step P t P t i;j = X l t l l (x i ) l (x j ); (3.11) here, 1 = 0 > 1 2 are sequence of the Eigenvalues,f l (x i )g andf l (x j )g are sequence of biorthogonal right and left eigenvectors, respectively. These eigen- vectors give the low dimensionally embedding. Since, the graph is fully connected, 49 the largest eigenvalue is trivial (viz. 1 = 1), hence the eigenvector 1 is ignored. Then, the low dimensional embedding after one time step is given by, (x) = ( 2 2 (x); 3 3 (x);:::; m+1 m+1 (x)): (3.12) Attention of the eigenvalues dictates the dimension m of the embedding. Diusion distance between any two points on the manifold is evaluated using the Eu- clidean distance between diusion coordinates . D (x i ;x j ) 2 = X l 2 l ( l (x i ) l (x j )) 2 : (3.13) Denition 1. The diusion map (x) :x!R m with is dened by (x) := 0 B B B B B B B @ 2 2 (x) 3 3 (x) . . . m+1 m+1 (x) 1 C C C C C C C A (3.14) 3.2.3 Intrinsic interpolation using Gaussian process regression Once the low dimensional embeddingM is identied using the diusion map , corre- sponding diusion coordinates of all the sample pointsf i (x i )g k i=1 are obtained and the functional relation between q( (x);!) and (x) is dened as q( (x);!) =f( (x)) +z( (x);!): (3.15) 50 Here, f( (x)) is regression model or trend function in diusion coordinates (x), and z( (x); !) is a zero-mean random process, with zero mean and covariance between any two diusion coordinates (x); (y) is dened as, C ( (x); (y)) =Efz( (x);:);z( (y);:)g: (3.16) When the process is second-order stationary, a correlation function of z( (x)),R (d ), can be introduced in the form 2 R (d ( (x); (y)) =C ( (x); (y)) (3.17) Here, 2 is the variance of the process and d ( (x); (y) is the intrinsic distance between two diusion coordinates (x) and (y). The correlation parameters , the variance 2 and regression parameters are obtained following maximizing the like- lihood argument. In the present development, we express the corresponding likelihood function, L ( ; 2 ; ), in terms of its logarithm, in the form lnL ( ; 2 ; ) = 1 2 [nln( 2 )+ln(detR D )+(qf D ) T )R 1 D (qf D ) T )= 2 ]: (3.18) Here,R D is then xn matrix of correlationsR( (x i ); (x j )) of then design points 1 i n; 1 j n, q =fq i ( (x);!)g n i=1 and f D =ff( (x i )g n i=1 . Once the parameters , 2 ; are identied using the optimization, q at (x ) is given by, ^ q( (x );!) = ^ f( (x )) +r T ( (x ))R D ( (x )) 1 (qf D ( (x ))): (3.19) Here,r T ( (x ))) = [R ( (x 1 ); (x )); ;R ( (x n ); (x ))] is the vector of correla- tions between the design points and (x ). Theq( (x );!) obtained using the proposed method is the best linear unbiased predictor for the intrinsic random variable on the manifold. 51 3.3 Case study The proposed regression method is validated by predicting the permeability using the core log data, by leave-one-out cross-validation test or blind test. Specically, permeabil- ity is predicted at one well, by training the models at other wells, and compared with the measured data. 3.3.1 Dataset 0 20 40 60 80 Spatial coordinate 1 −2 0 2 4 6 8 10 12 Spatial coordinate 2 1 2 106 107 116 23 236 29 Figure 3.1: Spatial locations of the wells. The core log dataset considered in the example consists of eight wells numbered as 1, 2, 23, 29, 106, 107, 116, and 236. Figure 3.1 shows the spatial locations of these wells. The attributes in this log are 3D spatial coordinates, porosity from logs, acoustic slowness attribute, radioactivity and rock density. Since the attributes have dierent scales and dimensions, each component in the dataset is normalized using x 0 = x min(x) max(x) min(x) : (3.20) 3.3.2 Diusion map of the dataset Low dimensional embedding of the dataset is obtained using diusion maps. The param- eter of the kernel chosen as 1.125 based on the plot of vsL() is shown in Figure 3.2 52 2 4 6 8 Index, n 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Eigenvalue, λ n (a) 10 -10 10 0 10 10 ǫ 10 4 10 5 10 6 L(ǫ) ǫ vs L( ǫ) ǫ=1.125 (b) Figure 3.2: (a) Eigenvalue decay of the diusion kernel (b) vs L() curve . (b). Figure 3.2 (a) shows the eigenvalue decay of the diusion kernel, based on the attenuation of the eigenvalues, the dimension of the embedding is restricted to six. Figure 3.3 (a) & (b) depict the scatter plot of rst three dominant directions in diu- sion space and principle component (PC) space, for all the wells. These plots reveal that clusters obtained from the diusion maps are dierent from the clusters in the PC space. The spread of the data points in each cluster or individual well in the diusion space is minimum compared to the spread in the PC space. Even after data normalization, the cluster arrangement in the principle component space is similar to the spatial locations as shown in Figure 3.1, but in diusion space, it is dierent. 53 Diffusion coordinate 1 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 Diffusion coordinate 2 −0.4 −0.2 0.0 0.2 0.4 0.6 Diffusion coordinate 3 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 1 2 106 107 116 23 236 29 (a) PC coordinate 1 −1.5 −1.0 −0.5 0.0 0.5 1.0 PC coordinate 2 −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 PC coordinate 3 −1.0 −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 1 2 106 107 116 23 236 29 (b) Figure 3.3: Scatter plot of the rst three dominant directions in (a) Diusion space (b) Principle component space. 3.3.3 Results Permeability is a positive quantity and predictions using GPR and proposed method can produce negative values, hence, the prediction is carried out on the log permeability. Figure 3.4 shows the histogram of the measured permeability and log measured perme- ability. The following procedure is followed to validate the proposed method using the blind tests: Data points from all the wells except the well of interest are grouped together to form training set and data points from the well of interest is treated as nal testing test. Hundred random samples are drawn from the training set and treated as an initial testing set, and remaining samples are taken as an initial training set. 54 0 500 1000 1500 2000 2500 Permeability 0 200 400 600 800 1000 1200 Frequency (a) −10 −8 −6 −4 −2 0 2 4 6 8 Log permeability 0 100 200 300 400 500 600 700 Frequency (b) Figure 3.4: Histogram of the measured permeability in (a) Original space (b) Transformed space. No Correlation function Regression model MSE for GPR MSE for proposed method 1 Exponential Constant 2.88e+27 2.65e+33 2 Squared Exponential Constant 1.12e+37 1.95e+89 3 Exponential Linear 4.27e+33 3.16e+34 4 Squared Exponential Linear 6.02e+37 9.86e+99 Table 3.1: Mean square error (MSE) for dierent correlation functions and regression models. Correlation function and regression model are chosen based on mean square error (MSE) of the prediction on the initial test set, with the models trained using the initial training set. Once the correlation function and the regression model are chosen, the nal regres- sion model is obtained by computing the parameters of the correlation function and regression model with the training set. The nal regression model is used to predict the permeability at the well of interest. The above procedure is repeated for all the wells. Table 3.1 shows the MSE of various correlation functions and regression models for the initial training set. Similar behavior is observed for all initial training sets, hence, 55 0 1 2 3 4 5 6 Log measured permeability 0 1 2 3 4 5 6 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 Log measured permeability 0 1 2 3 4 5 6 7 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.5: Blind test results for the well 1 using (a) Proposed method (b) Gaussian process regression. constant regression function and non-isotropic exponential correlation function are used. Parameters of the regression model and correlation function are calculated by doing several one-dimensional optimizations instead of running one single multi-dimensional optimization as proposed in [49]. Figure 3.5 to 3.12 depict the blind test results for the wells 1, 2, 23, 29, 106, 108, 116, 136 using GPR and proposed approach. Here, the red circular dots are scatter plots of the predicted mean versus measured permeability and dotted red line is the plot of measured versus measured permeability. The starting and ending of the blue line at each location shows the rst and third quantiles of the predictions. These results show that the proposed approach gave a considerable improvement in the prediction of the permeability except for the wells 106 and 107 because these wells are located very far from other wells in diusion space as shown in Figure 3.3. 56 0 1 2 3 4 5 6 Log measured permeability −1 0 1 2 3 4 5 6 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 Log measured permeability −3 −2 −1 0 1 2 3 4 5 6 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.6: Blind test results for the well 2 using (a) Proposed method (b) Gaussian process regression. 0 1 2 3 4 5 6 7 Log measured permeability −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 7 Log measured permeability −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.7: Blind test results for the well 23 using (a) Proposed method (b) Gaussian process regression. 0 1 2 3 4 5 6 7 Log measured permeability −1 0 1 2 3 4 5 6 7 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 7 Log measured permeability −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.8: Blind test results for the well 29 using (a) Proposed method (b) Gaussian process regression. 57 0 1 2 3 4 5 6 7 8 Log measured permeability −10 −8 −6 −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 7 8 Log measured permeability −10 −8 −6 −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.9: Blind test results for the well 106 using (a) Proposed method (b) Gaussian process regression. 0 1 2 3 4 5 6 7 8 Log measured permeability −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 7 8 Log measured permeability −6 −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.10: Blind test results for the well 107 using (a) Proposed method (b) Gaussian process regression. 0 1 2 3 4 5 6 7 8 Log measured permeability 0 1 2 3 4 5 6 7 8 9 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 7 8 Log measured permeability 0 1 2 3 4 5 6 7 8 9 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.11: Blind test results for the well 116 using (a) Proposed method (b) Gaussian process regression. 58 0 1 2 3 4 5 6 7 8 Log measured permeability −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (a) 0 1 2 3 4 5 6 7 8 Log measured permeability −4 −2 0 2 4 6 8 Log predicted permeability Mean prediction Quantiles of prediction Measurement (b) Figure 3.12: Blind test results for the well 236 using (a) Proposed method (b) Gaussian process regression. 59 3.4 Conclusions We introduced a regression procedure for intrinsic variables constrained onto a mani- fold embedded in an ambient space. The proposed method augments manifold learning procedures with a Gaussian process regression (GPR). The superiority of the proposed method over GPR was shown by considering a case study, to characterize the permeabil- ity of rock formations from petrophysics data. The proposed method gave considerable improvements in the predictions over the GPR. 60 Chapter 4 An Ecient Basis Adaptation in Homogeneous Chaos Space for Stochastic Initial-boundary Value Problem with Specic Quantities of Interest This chapter introduces an ecient isometry computation method for the stochastic basis adaptation procedure to detect a low dimensional subspace adapted to quantities of interest. Especially, this method is useful for the stochastic initial-boundary value problem with specic quantities of interest. It uses the solution from the Gaussian part of the QoI at timetT to compute the isometry of the QoI at the timeT . Computational eciency of the proposed method is demonstrated to detect a low dimensional subspace adapted to the cumulative ow rate in a petroleum reservoir with a random permeability eld. 4.1 Introduction Recent advances in science have enabled us to include new physical processes to emu- late the real world processes more accurately. Nevertheless, the epistemic uncertainty associated with the coecients of the corresponding partial dierential equations (PDE) 61 renders the discrepancy between observations and models. One way to understand this discrepancy is by using the design of computer experiments (DACE) [50] such as cen- tral composite design [51], Box-Behnken design [52], and full and fractional factorial design [53]. Another way is to use stochastic sampling methods such as Monte Carlo, and Latin hypercube sampling (LHS). These two methods allow us to perform: un- certainty quantication to characterize the uncertainty associated with corresponding coecients; inverse analysis to correct for bias in the simulations by updating the co- ecients; data assimilation to combine data from multiple simulations and physics to enhance the model; systems design to identify the sensitivity of the inputs to produce optimal system performance. Even with the advent of the computers, the numerical solution to the PDE is compu- tationally expensive because of the mesh size and step size of the spatial and temporal discretization. Also, the epistemic randomness and costs associated with the data col- lection render the coecients vary from one point domain to another, leading to very high dimensional N stochastic eld representation of the coecients. Recently, meth- ods of white-noise calculus [54, 12, 55] have been adapted to reduce the computational burden associated with the DACE and sampling. These methods also become computa- tionally prohibitive as the dimension increases. More recently, a method was developed by Tippireddy et al. [56] to nd a low-dimensional subspace adapted to the quantity of interest. This paper addresses the computational burden associated with the Gaussian adap- tation procedure developed by Tippireddy [56] by introducing an ecient isometry com- putation to detect lower dimensional (nN) subspace adapted to QoI. Especially, this method is useful for the stochastic initial-boundary value problem with specic quanti- ties of interest. The central idea of the proposed method is to use the solution from the Gaussian part of the QoI at time tT to compute the isometry of the QoI at the time T . 62 4.2 Formulation Let ( ;F;P ) be a complete probability space, constructed by sample space ,-algebra of eventsF 2 , and probability measureP :F! [0; 1]. Consider multi-phase block-oil ow model [9] with random permeability @ @t Ws B w S w +r Ws B w ~ u w q W = 0 (4.1) @ @t Os B o S o +r Os B o ~ u o q O = 0 (4.2) @ @t R so Gs B o S o + Gs B g S g +r R so Gs B o u o + Gs B g u g q G = 0 (4.3) Here denotes the porosity;S ; ;p ;B ; s ;q are the saturation, viscosity, pressure, formation volume factor, density and volumetric rate of each phase (water, gas, and oil);R so is the gas solubility;u w ; u g ; u o are the Darcy velocities obtained by solving the Darcy's law for each phase u = k r (x;!) (rp rz) =o;w;g; x2 ; !2 ; (4.4) where ;k r is the mass density and relative permeabilities of each phase, is the magnitude of the gravitational acceleration, and z is the depth. The quantities of the interest (QoI) represented byq is the scalar valued functions of the time such as volumetric ow rates, and cumulative production, which are obtained using Q = ZZ E u dE; C = Z t Q dt; (4.5) here E is the cross-sectional area of the wellbore. 63 4.3 Stochastic dimensions of the problem In many real world scenarios, microscale uctuations lead to random permeability val- ues from one point domain to another, and their uncertainty should be described as a random process with given covariance structure. For simplicity, letk, the absolute perme- ability tensor be the diagonal tensor andk 2;2 (x) andk 3;3 (x) are deterministic functions of thek 1;1 (x). For example,k 1;1 (x) =k 2;2 (x) = 10k 3;3 (x) andY (x;!) = ln[k 1;1 (x)] be a Gaussian random process. Since the covariance functionC Y (x;y) =<Y (x;!)Y (y;!)> is bounded, symmetric and positive denite, it allows the following decomposition [57] C Y (x;y) = 1 X i=1 i f i (x)f i (y); (4.6) where are the eigenvalues, f(x) and f i (y) are eigenfunctions, which are deterministic and orthogonal functions. Z D f i (x)f j (y)dx = ij ; i;j 1 (4.7) Using the Karhunen{Lo eve expansion, the random process Y (x;!) can be expressed in terms of f i (x) as Y (x;!) =<Y (x;!)> + 1 X i=1 i (!) p i f i (x); (4.8) where (!)'s are orthonormal Gaussian random variables with zero mean, i.e., < i (!) >= 0 and < i (!) j (!) >= ij . The eigenvalues , and eigenfunctions f n (x) are obtained by solving following integral equation either analytically or numerically, Z D C Y (x;y)f(x)dx =f(y): (4.9) 64 The attenuation of the eigenvalues allows truncation of the innite sum in the Equa- tion 4.8 up to N terms, Y (x;!) =<Y (x;!)> + N X i=1 i (!) p i f i (x); (4.10) here N is the stochastic dimension. 4.4 Polynomial chaos decomposition Consider q :R N !R such that R R N (q(x)) 2 e 1 2 x T x <1, then q2L 2 ( ) in general, the second order moment of QoI, q2L 2 ( ), represented in a Polynomial Chaos decomposi- tion of the form [55] q p () = X jjp q (); (4.11) Here () are the Hermite polynomials, = ( 1 ;:::; N ) is the multi-index,jj = D P i=1 i , p is the order of the polynomials, lim p!1 q p ()!q() and = ( 1 ;::: D ). The normalization Hermite polynomials are given by, () =H = p (); kH k 2 L 2 ( )=! (4.12) where! = N Q i=1 i . The coecients in the chaos decomposition of q can be evaluated as an integral which is then approximated through quadrature as, Z R N q() ()e 1 2 T d = r X i=1 q( r )w r ; (4.13) where q( r ) is the solution to PDE Equation 4.1 at the quadrature point r and w r is the weight of the quadrature point. The computational burden of carrying out the above 65 high-dimensional quadrature increases with stochastic dimension N andjj and quickly becomes prohibitive. 4.5 Basis adaptation to QoI To mitigate the computational burden basis adaptation in homogeneous chaos space was proposed by Tippireddy et al. [56]. This method identies a low dimensional subspace adapted to the quantity of interest q using isometry and the projection. Let A be the isometry onR N and let be dened as, =A; AA T =I: (4.14) Also, let q A be the image of q under this mapping, q A () =q(), since is the basis in H, homogeneous Chaos of order N gives A () = (); (4.15) and the expansions of the q and q A are respectively given by, q() = X 2Ip q (); q A () = X 2Ip q A () (4.16) leads to q A = X 2Ip q ( ; A ): (4.17) Here is the multi-index and ( ; A ) is the inner product dened as the expectation of the product relative to the Gaussian measure. 66 Consider a subspace V l of L 2 ( ) spanned byf A ;2Ig, whereII p ,I p is the set of d-dimensional multi-indices of degree less than or equal to p. The projection of q A on V l denoted by q A;l can be expressed in the form q A;l () = X 2I X 2Ip q ( ; A ) (): (4.18) The projection ofq A onV l alternatively can be expressed with respect tof ()g leads to q A;l () = X 2I q 0 (); (4.19) leading to q 0 = X 2I X 2Ip q ( ; A )( A ; ) (4.20) 4.5.1 Gaussian adaptation For the problem of interest, the viable way to compute the isometry A is using the Gaussian adaptation procedure [56]. In this procedure, Gaussian part of the solution q is computed using lower order quadrature, let that Gaussian part be G 1 = X 2I 1 q = N X i=1 q g i i ; (4.21) then the rst row of the isometry A is taken as A 1;: =fq e 1 ;:::;q e N g, where q e i = q g i =jq g j, here q g is the vector (q g 1 ;:::;q g N ). Let the remaining rows of the isometry A be constructed through a Gram-Schmidt procedure, then =A and the rst dominant direction 1 = N P i=1 q e i i . In this case q A () is given by q A () =q A 0 +q A e 1 1 + X 1<jjp q A () = q A 0 +q A e 1 1 + X 1<jjp 2E 1 q A () + X 1<jjp = 2E 1 q A () (4.22) 67 and the q A;l on the subspace V l is given by q A;l () =q A 0 + X 2E 1 q A (): (4.23) The error in the above representation is given by q A ()q A;l () = X 1<jjp = 2E 1 q A (): (4.24) 4.6 Ecient computation of the isometry for initial bound- ary value problems for specic quantities of interest Even with the advent of the computers, the numerical solution to Equation 4.1 is compu- tationally expensive because of the mesh and step sizes of the spatial and temporal dis- cretization. Even though the Gaussian adaptation procedure reduced the computational burden, there is still a small burden in computing the Gaussian part of the solution. To lower this burden, we introduce an ecient isometry computation method for the basis adaptation procedure to detect lower dimensional (n N) subspace adapted to quantities of interest. Especially, this method is useful for the stochastic initial-boundary value problem with specic quantities of interest. The idea here is to use the solution from the Gaussian part of the QoI at time tT to compute the isometry of the QoI at the time T . Let the Gaussian part of the solutionq at timetT be computed using lower order quadrature, by solving Equation 4.1 up to the time stept, and let that Gaussian part be G 1 (t) = X 2I 1 q (t) = N X i=1 q g i (t) i ; (4.25) then the rst row of the isometry at timet,A(t) is taken asA 1;: =fq e 1 ;:::;q e N g, where q e i =q g i (t) =jq g(t) j, hereq g(t) is the vector (q g 1 (t);:::;q g N (t)). The remaining rows of the 68 isometry A(t) are constructed through a Gram-Schmidt procedure, then (t) =A(t) and the rst dominant direction 1 (t) = N P i=1 q e i (t) i . In this case q A (t)() is given by q A(t) ((t)) =q A(t) 0 +q A(t) e 1 1 (t) + X 1<jjp q A(t) ((t)) = q A(t) 0 +q A(t) e 1 1 (t) + X 1<jjp 2E 1 q A(t) ((t)) + X 1<jjp = 2E 1 q A(t) ((t)) (4.26) and the q A(t);l on the subspace V l is given by q A(t);l ((t)) =q A(t) 0 + X 2E 1 q A(t) ((t)): (4.27) The error in the above representation is given by q A(t) ((t))q A(t);l () = X 1<jjp = 2E 1 q A(t) ((t)): (4.28) Note here that using all the dimensions in the (t) will lead to exact representation of q A(t);l . 4.7 Numerical example A two-phase ow problem in a petroleum reservoir is considered to demonstrate the computational eciency of the proposed method. The quantity of interest considered here is the cumulative ow rate at the end of twenty years. The corresponding partial dierential equations are numerically solved using Eclipse [10], a commercial reservoir simulator. The size of the physical domain considered is 3600 m x 3600 m x 15 m. The domain is discretized with the grid size of 30 m x 30 m x 15m leading to total 14400 69 (a) (b) Figure 4.1: (a) Spatial discretization of the reservoir (b) A realization of the permeability obtained from the chosen covariance kernel. grid points, as depicted in Figure 4.1 (a). A time step of one day is considered for the analysis leading to a total of 7200 time steps. Rock pressure is considered to be 6000 psi, the rock compressibility 1 10 6 psi 1 , initial pressure of the reservoir set at 6000 psi, the initial water saturation is taken to be 0.30 and oil specic gravity assumed to be 0.83. Assumed relative permeability coecients for the water and oil at dierent saturation rates are shown in Table 4.1. Pressure-volume-temperature (PVT) data of the dead oil used for the simulations are given in Table 4.2. The reservoir consists of one injection well on the top left corner and three production wells at the remaining corners. Permeability is assumed to be log-normal process with exponential correlation func- tion, isotropic length scale (correlation length) of 150 m amplitude scale (standard de- viation) of 2.25 and mean 3 mD. A Karhunen{Lo eve expansion is used to discretize the Gaussian process associated with the permeability, with 100 terms retained in the ex- pansion (90% contribution from the corresponding eigenmodes). Figure 4.1(b) depicts a realization obtained from the corresponding process. A realization of the pressure and oil 70 Water saturation Water relative permeability Oil relative permeability when only oil and water are present Water-oil capillary pressure (psi) 0.200 0.000 1.000 20.000 0.250 0.007 0.840 17.564 0.300 0.028 0.694 15.255 0.350 0.063 0.563 13.074 0.400 0.111 0.444 11.021 0.450 0.174 0.340 9.095 0.500 0.250 0.250 7.296 0.550 0.340 0.174 5.626 0.600 0.444 0.111 4.083 0.650 0.563 0.063 2.667 0.700 0.694 0.028 1.379 0.750 0.840 0.007 0.218 0.800 1.000 0.000 -0.814 Table 4.1: Relative permeability table for oil and water. Oil phase pressure (psia) Oil formation volume factor (rb/stb) Oil viscosity (cP) oil and water are present 300 1.05 2.85 800 1.02 2.99 8000 1.01 3 Table 4.2: Pressure volume temperature (PVT) of the dead oil. saturation at the end of year one, obtained from solving the numerical reservoir model are presented in Figure 4.2. 71 (a) (b) Figure 4.2: A realization of the (a) Pressure (b) Oil saturation at the end of one year obtained from solving the numerical reservoir model with chosen covariance kernel for the permeability. 4.7.1 Basis adaptation using proposed method First row of the isometry 1 (t) is computed based on the Gaussian part of the QoI at t by constructing the multiple PCE of order one up to the year one. The rst subplot of Figure 4.4 (a) depicts a variation of the 1 (t) within the each year. This Figure shows that 1 (t) remained unchanged after few months, hence, 1 (t) at the end of the year one is taken as the rst row of the isometry A(t). The remaining rows of the isometry are computed using Gram-Schmidt procedure. Using this A(t), PCE of QoI at the end of year one, ve, ten, fteen and twenty are computed using 1D and 2D adaptation with order two and three using sparse grid cubature in 1D and 2D. Figure 4.3 depicts the probability density function of the QoI at the end of year one, ve, ten, fteen and twenty obtained by sampling PCE constructed from the proposed method. To verify the convergence of 1 (t), multiple PCE of order one are constructed up to the T , twenty years. The Figure 4.4 (a) & (b) show the evaluation of the 1 over the time. This plot shows that 1 remained constant from year two to twenty. Even though, 72 45 50 55 60 65 70 75 80 85 Cumulative production (QoI) 0.00 0.02 0.04 0.06 0.08 0.10 0.12 pdf of QoI Year one 150 200 250 300 350 400 Cumulative production (QoI) 0.000 0.005 0.010 0.015 0.020 pdf of QoI Year five 350400450500550600650700750800 Cumulative production (QoI) 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 pdf of QoI Year ten 500 600 700 800 900 1000 1100 Cumulative production (QoI) 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 pdf of QoI Year fifteen 700 800 900 10001100120013001400 Cumulative production (QoI) 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 pdf of QoI Year twenty 1D Adaptation, order=2 2D Adaptation, order=2 1D Adaptation, order=3 2D Adaptation, order=3 Figure 4.3: Probability density function of cumulative oil production in MSTB for year one, ve, ten, fteen and twenty obtained using the proposed method. there is a change in the 1 at one of the dimensions in the last year, basis obtained 1 (t) at the end year one still provided the ecient isometry. Figure 4.5 (a) depicts the rst nine polynomial chaos coecients of permeability associated with an expansion of the lognormal permeability process. Note the global features in these coecients. Figure 4.5 (b) shows the same quantities for the problem adapted to the QoI specied by cumulative ow rate at the end of the year one. Note the localization around each well in the polynomial chaos coecients induced by the quantity of interest. 73 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year one 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year two 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year three 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year four 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year five 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year six 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year seven 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year eight 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year nine 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year ten Month:1 Month:2 Month:3 Month:4 Month:5 Month:6 Month:7 Month:8 Month:9 Month:10 Month:11 Month:12 (a) 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year eleven 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year twelve 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year thirteen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year fourteen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year fifteen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year sixteen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year seventeen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year eighteen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year nineteen 0 20 40 60 80 100 Dimension (i) −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2 η 1(i) Year twenty Month:1 Month:2 Month:3 Month:4 Month:5 Month:6 Month:7 Month:8 Month:9 Month:10 Month:11 Month:12 (b) Figure 4.4: Evaluation of the 1 over the (a) years one to ten (b) years eleven to twenty 74 (a) (b) Figure 4.5: First nine PCE coecients of the permeability with respect to (a) original Gaussian variable (b) one dimensional Gaussian variable 1 75 4.7.2 Computational eciency of the proposed method Method Computational time PCE with sparse grid 20601*20c t =412020c t Gaussian adaptation (201+8+21)*20c t =4600c t Proposed method 201c t +(8+21)*20c t =760c t Table 4.3: Computational cost associated with PCE, Gaussian adaptation and the pro- posed method, here c t is computational time to solve the PDE up to time t. Letc t sec be the computational time required to obtain the QoI at the timet =T=20, PCE with the sparse grid cubature requires 20601 (spare-grid level two) numerical eval- uations of the PDE up to the time T , thus requires the total computational time of 412020c t sec. Gaussian adaptation requires 201 (sparse grid level one) computations of the PDE up to the timeT to compute Gaussian part and 29 (8, 21 for 1D, 2D) evaluations to verify the convergence of the adaptation, thus requires total 4600 c t sec of compu- tational time. The proposed method requires 201 (sparse grid level one) evaluations of the PDE up to the time t to compute Gaussian part and 29 (8, 21 for 1D, 2D adapta- tions) evaluations of PDE up to timeT to verify the convergence of the adaptation, thus requires total 760c t sec of computational time. In conclusion, using the Gaussian adap- tation computational requirement of the model has reduced by 98.88% and the proposed approach has further increased the eciency by 83.02%. 76 4.8 Conclusions This chapter introduced an ecient isometry computation method for the stochastic basis adaptation procedure to detect lower dimensional (n N) subspace adapted to quantities of interest. Especially, this method is useful for the stochastic initial-boundary value problem with specic quantities of interest. It uses the solution from the Gaussian part of the QoI at time tT to compute the isometry of the QoI at the time T . The eciency of this method is demonstrated to detect low dimensional subspace adapted to the cumulative ow rate in a petroleum reservoir with a random permeability eld. The Gaussian adaptation computational requirement of the model has reduced by 98.88% and the proposed approach has further increased the eciency by 83.02%. 77 Chapter 5 Convergence Enhancement of Basis Adaptation in Homogeneous Chaos Spaces Basis adaptation in homogeneous chaos spaces identies the low-dimensional subspaces associated with the quantities of interest (QoI). Sometimes, the dimension of this sub- space is greater than one and other times, even if there is one-dimensional subspace does exist, it may be computationally expensive to nd the corresponding isometry. For these scenarios, a computationally ecient method is proposed to compute the sub- space. Specically, the proposed method uses the sensitivity obtained from the Gaussian part of the solution adaptation to enhance the convergence properties of the adaptation. Convergence eciency of this method is demonstrated on a multidimensional modied Rosenbrock function and also in the context of uncertainty propagation in composite manufacturing processes. Randomness in the material properties and computational cost associated with the solution of corresponding physics and makes this uncertainty prop- agation challenging. 5.1 Introduction The increase in the world energy consumption with depletions in the natural resources is one of the problems faced by the developed and developing countries in the world. 78 Although, there is an active research to increase the renewable sources of energy contri- bution, their usage projections by the year 2040 is limited to 10% [58]. The consumption of petroleum products can be decreased by reducing the motor gasoline consumption of on-road light-duty vehicles (LDVs) by increasing their average eciency. It has been estimated that 10% reduction in the weight of the LDVs leads to 6.5% increase in their fuel eciency. At present, there is an active usage of carbon ber (CF) composites for the Class-A vehicle body panels. Currently, researchers are focusing to expand the use of CF composites to replace the structural parts of the Class-A vehicles. Composite man- ufacturing involves many physical processes such as ber draping, resin injection, and curing. Numerical modeling of these processes involves solving many partial dierential equations (PDE) involving micro and macro scale parameters. Epistemic uncertainty associated with these parameters render the PDEs to be stochastic. The uncertainty propagation and calibration using these stochastic PDEs allow us to increase the safety and eciency of the Class-A vehicles. This work focuses on the uncertainty propagation of the resin transfer molding (RTM). Numerical modeling of the RTM process is helpful to detect air bubbles, dry spots and to spot the resin rich areas and also zones of the high porosity [59, 60, 61]. This modeling also helps to nd the optimal positions of injection ports and vents. Also, it is economical to design the mold based on numerical simulations before constructing the actual mold. Randomness in the material properties and uncertainties associated with manufacturing process render the quantities of the interest (QoI) random. Since, the numerical models for the RTM involve solving coupled PDEs with multi-physics, associ- ated computational cost is very high. Uncertainty propagation using traditional sampling methods such as Monte Carlo sampling, Latin hyper cube sampling [62], and expansion methods [54, 12, 55] such as polynomial chaos expansion (PCE) becomes prohibitive as the dimensionality of the random variables increases. To mitigate this, a basis adaptation in homogeneous chaos spaces to identify the low- dimensional subspaces associated with the quantities of interest (QoI) is proposed by 79 Tippireddy at el. [56]. Sometimes, the dimension of this subspace is greater than one and other times, even if there is one-dimensional subspace does exist, it may be compu- tationally expensive to nd the corresponding isometry. For these scenarios, a computa- tionally ecient method is proposed to compute the subspace. Specically, the proposed method uses the sensitivity obtained from the Gaussian part of the solution adaptation to enhance the convergence properties of the adaptation. Convergence eciency of this method is demonstrated on a multidimensional modied Rosenbrock function, also in the context of uncertainty propagation in composite manufacturing processes 5.2 Modeling the RTM process Numerical simulation of the RTM process involves modeling three physical processes: ow in the porous media, thermal phenomena and the kinetics of the resin polymerization. 5.2.1 Flow in the porous media In the RTM process, the resin is injected and allowed to ow through the ber rein- forcement by creating the pressure dierence. Here the ber reinforcement acts like a porous media. This ow is governed by the permeability of the porous media, viscosity of the resin and pressure dierence. The corresponding ow phenomena is modeled by the Darcy's law which is given by V = K rP; in the domain : (5.1) Here K is permeability tensor of the porous media, is the resin viscosity, V is the vector of the Darcy's velocities, P is the pressure andr is the gradient operator. For the resin mass balance, the Darcy's velocities have to satisfy the following equation, r:V = 0: (5.2) 80 Combining the Equations 5.1 and 5.2 leads to r: K rP = 0: (5.3) The boundary conditions required to solve the above equation can be specied as imposed pressure condition or Dirichlet conditionp =p(x), or imposed ow rate condition or the Neumann conditionV:~ n =Q(x) on the boundary d . Here~ n is the unit normal vector and Q is the ow rate. 5.2.2 Thermal phenomena In the RTM process, the brous reinforcement is kept inside the cavity of the mold, and the resin is injected to ll the mold and then the part is progressively solidied. Thermal phenomena have a strong in uence on this process during the lling and curing stages. Since the viscosity of the resin is a function of the temperature, it in uences the ow during the mold lling. Also, the reactivity of the polymerization reaction is a function of the temperature, hence it in uences the curing process. Three thermal phenomenon can cause the changes in the temperature: rst, the heat conduction phenomena between the bers and the resin, second, convection transport between the cavity and the resin and nally, heat production due to the exothermic chemical reaction of resin polymerization. These phenomena are modeled by the following equation, C p @T @t + r c pr V:rT =r:(k:rT ) r h d dt ; (5.4) where is the density, C p is the specic heat, T is the temperature, t is the time, ~ V is the vector of the Darcy's velocities, k is the heat conduction coecient sensor, h is the total enthalpy of the polymerization of the resin and is the degree of the cure. Here, subscript r indicates the resin. The above equation is solved with the prescribed temperature T =T 0 or the heat ux @T @n =H or the heat convection @T @n =h(T 1 T ). 81 5.2.3 The kinetics of the resin polymerization During the curing phase of the RTM process, the injected liquid resin will be solidied, and the corresponding kinetics of the resin polymerization is modeled using Kamal- Sourour [63] model. The Kamal-Sourour equation with n components is given by, = n X i=1 C i i ; (5.5) d i dt =K i (T ): m i i :(1 i ) p i : (5.6) Here, for each component i, i is the degree of the cure, d i dt is the cure rate, m i and p i are the constants describing the sensitivity of autocatalytic reaction. The K i 's are dened by the Arrhenius law K i =A i exp(E i =RT ); (5.7) hereA i is the prefactor,E i is the activation energy,R is the universal gas constant, and T is the temperature. 5.3 Polynomial chaos decomposition Most of the parameters in the Equation 5.1 to Equation 5.6 are microscale parameters, epistemic uncertainty associated with these parameters renders the coecients in these equations random. Consider ( ;F;P ) be a probability triplet, with the sample space , and the -algebra of eventsF 2 , and a probability measure P :F! [0; 1]. Let q : R N ! R such that R R N (q(x)) 2 e 1 2 x T x <1, be our quantity of the interest (QoI). Since q2 L 2 ( ) q2 L 2 ( ), represented in a Polynomial Chaos decomposition of the form [55] q p () = X jjp q (); (5.8) 82 Here () are the Hermite polynomials, = ( 1 ;:::; N ) is the multi-index,jj = D P i=1 i , p is the order of the polynomials, lim p!1 q p ()!q() and = ( 1 ;::: D ). The normalization Hermite polynomials are given by, () =H = p (); kH k 2 L 2 ( )=! (5.9) where! = N Q i=1 i . The coecients in the chaos decomposition of q can be computed with an integral which is then approximated using the quadrature as, Z R N q() ()e 1 2 T d = r X i=1 q( r )w r ; (5.10) where q( r ) is the solution to Equation 5.1 to Equation 5.6 at the quadrature point r andw r is the weight of the quadrature point. The computational burden associated with the above quadrature increases with the number of stochastic parameters N and with thejj and quickly becomes prohibitive. 5.3.1 Basis adaptation to QoI Recently, to reduce the computational burden, basis adaptation in homogeneous chaos space was proposed by Tippireddy et al. [56]. This method identies a low dimensional subspace adapted to the quantity of interest q using isometry and the projection. LetA be the isometry onR N and let be dened as, =A; AA T =I: (5.11) Also, let q A be the image of q under this mapping, q A () =q(), since is the basis in H, homogeneous Chaos of order N, gives A () = (); (5.12) 83 and the expansions of the q and q A are respectively given by, q() = X 2Ip q (); q A () = X 2Ip q A () (5.13) leads to q A = X 2Ip q ( ; A ): (5.14) Here is the multi-index and ( ; A ) is the inner product dened as the expectation of the product relative to the Gaussian measure. Consider a subspace V l of L 2 ( ) spanned byf A ;2Ig, whereII p ,I p is the set of d-dimensional multi-indices of degree less than or equal to p. The projection of q A on V l , denoted by q A;l can be expressed in the form q A;l () = X 2I X 2Ip q ( ; A ) (): (5.15) The projection ofq A onV l alternatively can be expressed with respect tof ()g, leads to q A;l () = X 2I q 0 (); (5.16) leading to q 0 = X 2I X 2Ip q ( ; A )( A ; ) (5.17) 5.3.1.1 Gaussian adaptation A computationally viable way to nd the isometryA for the problem of interest is using the Gaussian adaptation procedure introduced by Tippireddy et al. [56]. In this proce- dure, Gaussian part of the solution q is computed by using the lower order quadrature rules, let that Gaussian part be 84 G 1 = X 2I 1 q = N X i=1 q g i i ; (5.18) then the rst row of the isometry A is taken as A 1;: =fq e 1 ;:::;q e N g, where q e i = q g i =jq g j, here q g is the vector (q g 1 ;:::;q g N ). Let the remaining rows of the isometry A be constructed through a Gram-Schmidt procedure, then =A and the rst dominant direction 1 = N P i=1 q e i i . In this case q A () is given by q A () =q A 0 +q A e 1 1 + X 1<jjp q A () = q A 0 +q A e 1 1 + X 1<jjp 2E 1 q A () + X 1<jjp = 2E 1 q A () (5.19) and the q A;l on the subspace V l is given by q A;l () =q A 0 + X 2E 1 q A (): (5.20) The error in the above representation is given by q A ()q A;l () = X 1<jjp = 2E 1 q A (): (5.21) 5.4 Convergence enhancement of basis adaptation Basis adaptation method described above nds the low-dimensional subspaces associated with the quantities of interest (QoI). Sometimes, the dimension of this subspace is greater than one and other times, even if there is one-dimensional subspace does exist, it may be computationally expensive to nd the corresponding isometry. In these scenarios, 85 the PCE is associated with higher dimensional subspace in Equation 5.20 obtained by including the other dimensions, q A;l () =q A 0 + X 2( n S i=1 E i ) q A (); (5.22) heren is the dimension of the associated subspace. The error in the above representation is given by q A ()q A;l () = X 1<jjp = 2( n S i=1 E i ) q A (): (5.23) The rst, second, . . .n th 's are obtained based on the isometryA constructed using A 1;: =fq e 1 ;:::;q e N g as the rst row and arbitrary values for the remaining rows, then using Gram-Schmidt procedure, complete isometry ofA is obtained. One of the choices for the initial starting point ofA is based on the identity matrix, that is A 0 ij = 8 > > > > > > < > > > > > > : q e j if i = 1; 0 if 1<i< N &i6=j; 1 if 1<i< N &i =j: (5.24) The choice of A will aect the dimension of the subspace V l , thus aecting the con- vergence of the basis adaptation procedure. Here, we introduce a new choice of A 0 based on the sensitivity obtained from the Gaussian part of the solution. The set c 1 =fq e 1 ;:::;q e N g computed from the Gaussian part gives an indication of the sen- sitivity q on each dimension. Let c 2 be the vector obtained by sorting absolute values of c 1 in descending order and c 3 be the corresponding indices that is [c 2 ;c 3 ] = sort(jc 1 j; descending): (5.25) 86 Based on c 2 and c 3 , new starting point ofA is dened as A 0 ij = 8 > > > > > > < > > > > > > : q e j if i = 1; 0 if 1<i< N &j6=c 3 (j); c 2 (j) if 1<i< N &j =c 3 (j): (5.26) TheA 0 ij in the Equation 5.26 allows faster convergence of the basis adaptation procedure since the associated subspace includes the sensitivity information. 5.5 Numerical examples To show the numerical eciency of the basis adaptation method and also convergence enhancement of the basis adaptation using the proposed method, multiple numerical examples will be considered here. First, a multidimensional modied Rosenbrock function is used to demonstrate the convergence. Then, uncertainty propagation in the RTM process with 70, 18 and 13 uncertain parameters are considered. 5.5.1 Basis adaptation of the modied Rosenbrock function As a rst numerical example, a modied Rosenbrock function, as shown in the Equa- tion 5.27, is used to demonstrate the convergence of the adaptation using the proposed method. The rst row of the isometry is constructed from the Gaussian adaptation. Fig- ure 5.1 shows the probability distribution function of the modied Rosenbrock function, for dierent dimensional adaptations obtained using isometry with initial A 0 based on the identity matrix. This gure reveals that an 8D adaptation is required to represent the modied Rosenbrock function with this method. Figure 5.2 shows the probabil- ity distribution function of the modied Rosenbrock function for dierent dimensional adaptations obtained using the proposed method, by constructing the isometry using the sensitivity obtained from the Gaussian part. This gure reveals that a 6D adaptation 87 −600 −500 −400 −300 −200 −100 0 100 q 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 pdf of q D=1 D=2 D=3 D=4 D=5 D=6 D=7 D=8 D=9 D=10 Figure 5.1: Probability distribution function of the modied Rosenbrock function for dierent dimensional adaptations obtained using isometry with initial A 0 based on the identity matrix. can represent the modied Rosenbrock function, showing the faster convergence of the proposed method. q(x) = 9 X i=1 (2i + 1)(x i+1 x 2 i ) 2 + (1x i ) 2 where x = [x 1 ;:::;x 10 ]2R 10 (5.27) The top left, and the top right images in Figure 5.3 show the A 0 and isometry A based on the identity matrix. The bottom left, and the bottom right images in Figure 5.3 show theA 0 and isometryA obtained from the proposed method, by constructing the isometry using the sensitivity calculated from the Gaussian part. The isometry computed from the proposed method is non-sparse in nature by allowing the contribution from the sensitive parameters in low dimensional subspace. 88 −600 −500 −400 −300 −200 −100 0 100 q 0.000 0.005 0.010 0.015 0.020 0.025 0.030 pdf of q D=1 D=2 D=3 D=4 D=5 D=6 D=7 D=8 D=9 D=10 Figure 5.2: Probability distribution function of the modied Rosenbrock function, for dierent dimensional adaptations obtained using the proposed method. Figure 5.3: InitialA 0 (top left) and isometryA (top right) based on the identity matrix. InitialA 0 (bottom left) and isometryA (bottom right) from the proposed method. 5.5.2 Uncertainty propagation in the RTM process Figure 5.4 depicts the RTM setup used in the following examples. Using this setup, a square plaque of size 44.5 cm x 44.5 cm x 1.8 mm is manufactured. A quasi-isotropic laminate with 0 0 =45 0 = 45 0 =90 0 = 90 0 =45 0 =45 0 =0 0 carbon ber lay-up is used here. The locations of the inlet and outlet are also shown in Figure 5.4. The resin is pumped from the inlet with a gauge pressure of 55 psi, and the outlet is kept at the atmospheric 89 Figure 5.4: RTM setup used for the demonstration pressure creating a 40 psi pressure dierence. The quantity of interest considered here is the lltime, which is dened as the time required for the resin to reach from the inlet to outlet. The corresponding physics are modeled in the PAM-RTM [64], a commercial composite manufacturing simulator. 5.5.2.1 Uncertainty propagation by treating permeability as a random pro- cess For this example, 0 0 layer permeability is assumed to be a log-normal random process of exponential correlation function, with isotropic length scale (correlation length) 4.45 cm, amplitude scale (standard deviation) 2.89x10 12 m 2 and the mean 2.89x10 11 m 2 . A Karhunen{Lo eve expansion is used to discretize the Gaussian process associated with 90 48 50 52 54 56 58 60 62 Fill-time 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 pdf of the fill-time 1D 2D 3D 53 54 55 56 57 58 59 Fill-time 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 pdf of the fill-time 1D 2D 3D Figure 5.5: Probability distribution function of lltime from 1D, 2D, 3D adaptations: obtained using the proposed method (left), and with initial A 0 based on the identity matrix (right). the permeability, with 70 terms retained in the expansion (90% contribution from the corresponding eigenmodes). Figure 5.5 shows the probability distribution function of lltime from 1D, 2D and 3D adaptations obtained using: proposed method (left) and isometry with initialA 0 based on the identity matrix (right). This gure reveals that a 2D adaptation can represent the QoI with the proposed method, showing the faster convergence of the proposed method. The top left, and the top right images in Figure 5.6 show the initial A 0 and nal isometryA based on the identity matrix. The bottom left, and the bottom right images in Figure 5.6 show the A 0 and nal isometry A obtained using the proposed method, by constructing the isometry using the sensitivity obtained from the Gaussian part. The isometry computed from the proposed method is non-sparse in nature by allowing the contribution from the sensitive parameters in low dimensional subspace. Table 5.1 shows the cost associated with dierent uncertainty propagation methods and the computa- tional eciency of the PCE with adapted basis. 91 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 Figure 5.6: InitialA 0 (top left) and isometryA (top right) based on the identity matrix. InitialA 0 (bottom left) and isometryA (bottom right) from the proposed method. Method Number of samples Comments 1 MCMC Orders of thousands Computationally expensive 2 PCE with Quadrature n d [2 70 ] Map between input and output, curse of dimensionality 3 PCE with sparse grid <n d [10221 for level 2] Map between input and output, curse of dimensionality 4 PCE with adapted basis Cost of Gaussian part computation + cost of (1D +2D+3D) adaptations [208] Map between input and output, very ecient Table 5.1: Computational cost associated with dierent uncertainty propagation methods 92 5.5.2.2 Uncertainty propagation in 18 dimensions For this example, uncertainty propagation in 18 dimensions is considered. Table 5.2 shows the mean and coecient of variations of the random variables considered in this example. Gaussian part of the solution is computed using sparse grid level one. Using this, basis adaptation in 1D and 2D is carried out. Figure 5.7 shows the histogram of the lltime obtained by sampling using: Monte Carlo simulations (top left), Gaussian part of PCE (top right), PCE of order 2 from the 1D adaptation (bottom left), and the 2D adaptation (bottom right). No. Random Variable Deterministic values COV 1 Permeability along direction one 4.235E-10 m 2 0.75 2 Permeability along direction two 4.235E-10 m 2 0.75 3 Permeability of the racetrack 3.000E-08 m 2 0.5 4 Permeability of the channel 3.320E-06 m 2 0.5 5 Density of the resin 1100 kg=m 3 0.005 6 Specic heat of the resin 1205 J=kg-k 0.01 7 Thermal conductivity of the resin 0.11 W=m-K 0.01 8 Density of the fabric 1700 kg=m 3 0.005 9 Specic heat of the fabric 700 0.01 10 Thermal conductivity of the fabric (k 1) 2 W=m-K 0.01 11 Thermal conductivity of the fabric (k 2) 0.5 W=m-K 0.01 12 Specic heat of the laminate 700 J=kg-k 0.01 13 Thermal conductivity of the laminate 0.1 W=m-K 0.01 14 Thickness of each layer 2.25E-04 m 0.01 15 Porosity of each layer 0.52 0.01 16 Density of the channel 1000 kg=m 3 0.01 17 Specic heat of the channel 700 J=kg-k 0.01 18 Thermal conductivity of the channel 0.02 W=m-K 0.01 Table 5.2: Mean and coecient of variations of the random variables considered in the example 93 20 40 60 80 100 120 Fill time in seconds 0 20 40 60 80 100 Frequency 0 10 20 30 40 50 60 70 80 90 Fill time in seconds 0 5000 10000 15000 20000 25000 30000 Frequency 0 20 40 60 80 100 120 140 Fill time in seconds 0 10000 20000 30000 40000 50000 Frequency 0 20 40 60 80 100 120 140 160 Fill time in seconds 0 10000 20000 30000 40000 50000 60000 Frequency Figure 5.7: Histogram of the lltime obtained by sampling using: Monte Carlo simulations (top left), Gaussian part of PCE (top right), PCE from the 1D adaptation (bottom left), PCE from the 2D adaptation (bottom right). 5.5.2.3 Uncertainty propagation in 13 dimensions For this example, uncertainty propagation is carried out in 13 dimensions. The random variables considered are 0 0 layer permeability in the principle direction one, two, three, the and permeability of the channel and porosity. Also, the plaque is divided into 22.25 cm x 22.25 cm zones on the top, and bottom and the temperature in these zones is assumed to be a random variable. No. Random Variable Distribution 1 0 0 layer permeability (in m 2 ) in principle direction one N (2:89x10 11 ; 2:89x10 12 ) 2 0 0 layer permeability (in m 2 ) in principle direction two N (2:89x10 11 ; 2:89x10 12 ) 3 0 0 layer permeability(in m 2 ) in principle direction three N (2:89x10 11 ; 2:89x10 12 ) 4 Channel permeability (in m 2 ) N (3x10 8 ; 3x10 9 ) 5 Porosity N (0:5; 0:01) 6-13 Temperatures (in 0 K) zone 1 to 8 N (403:1; 1) Table 5.3: Distributions of the random variables considered in the example 94 40 50 60 70 80 90 100 110 120 Filltime (sec) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Density 1D Adaptation 2D Adaptation Figure 5.8: Probability density function of the ll time obtained using basis adaptations in 1D and 2D. Table 5.3 shows the distributions of the random variables considered in this example. The quantity of interest considered here is the lltime. Gaussian part of the solution is computed using sparse grid level one. Using this, basis adaptation in 1D and 2D is computed, Figure 5.8 depicts the probability density function of the ll time obtained from these adaptations. This gure shows the 1D adaptation can represent the QoI eectively. 95 5.6 Conclusions A computationally ecient method to compute the low dimensional subspace using the basis adaptation method is proposed. Specically, the proposed method uses the sen- sitivity obtained from the Gaussian part of the solution to enhance the convergence properties of the adaptation. The convergence eciency of this method is demonstrated on a multidimensional modied Rosenbrock function and also in the context of uncer- tainty propagation in the resin transfer molding process. 96 Bibliography [1] G. Jim, R. Kelly, B. Jennifer, D. Corinne, J. Chris, N. Jake, R. Chris, S. Lawrence, , and V. Kaylyn, \Integration of spatial data to support risk and impact assessments for deep and ultra-deepwater hydrocarbon activities in the gulf of mexico," EPAct Technical Report Series,Morgantown, WV; p 36., vol. NETL-TRS-4-2012, 2012. [2] C. Thimmisetty, A. Khodabakhshnejad, N. Jabbari, F. Aminzadeh, R. Ghanem, K. Rose, J. Bauer, and C. Disenhof, \Multiscale stochastic representation in high- dimensional data using gaussian processes with implicit diusion metrics," in Dy- namic Data-Driven Environmental Systems Science, pp. 157{166, Springer, 2015. [3] M. Shinozuka and G. Deodatis, \Simulation of stochastic processes by spectral rep- resentation," Applied Mechanics Reviews, vol. 44, no. 4, pp. 191{204, 1991. [4] F. Yamazaki and M. Shinozuka, \Digital generation of non-gaussian stochastic elds," Journal of Engineering Mechanics, vol. 114, no. 7, pp. 1183{1197, 1988. [5] S. Sakamoto and R. Ghanem, \Simulation of multi-dimensional non-gaussian non- stationary random elds," Probabilistic Engineering Mechanics, vol. 17, no. 2, pp. 167{176, 2002. [6] P. M uller, A. Erkanli, and M. West, \Bayesian curve tting using multivariate nor- mal mixtures," Biometrika, vol. 83, no. 1, pp. 67{79, 1996. [7] D. B. Dunson, N. Pillai, and J.-H. Park, \Bayesian density regression," Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 69, no. 2, pp. 163{183, 2007. [8] M. K. McNutt, R. Camilli, T. J. Crone, G. D. Guthrie, P. A. Hsieh, T. B. Ryerson, O. Savas, and F. Shaer, \Review of ow rate estimates of the deepwater horizon oil spill," Proceedings of the National Academy of Sciences, vol. 109, no. 50, pp. 20260{ 20267, 2012. [9] Z. Chen, G. Huan, and Y. Ma, Computational methods for multiphase ows in porous media, vol. 2. Siam, 2006. [10] Schlumberger, \Eclipse-industry reference reservoir simulator," 97 [11] R. G. Ghanem and P. D. Spanos, Stochastic nite elements: A spectral approach, vol. 41. Springer, 1991. [12] D. Xiu and G. E. Karniadakis, \The Wiener{Askey polynomial chaos for stochas- tic dierential equations," SIAM Journal on Scientic Computing, vol. 24, no. 2, pp. 619{644, 2002. [13] J. Vogel, \In ow performance relationships for solution-gas drive wells," Journal of Petroleum Technology, vol. 20, no. 1, pp. 83{92, 1968. [14] A. Hagedorn and B. Kermit, \Experimental study of pressure gradients occurring during continuous two-phase ow in small-diameter vertical conduits," Journal of Petroleum Technology, vol. 17, no. 4, pp. 475{484, 1965. [15] N. Cressie, \The origins of kriging," Mathematical Geology, vol. 22, no. 3, pp. 239{ 252, 1990. [16] E. Isaaks and R. Srivastava, \Applied geostatistics," London: Oxford University, 2011. [17] G. Matheron, \Principles of geostatistics," Economic geology, vol. 58, no. 8, pp. 1246{1266, 1963. [18] D. Krige, \A statistical approach to some basic mine valuation problems on the witwatersrand," Jnl. C'hem. Met. and Min. Soc. S. Afr, 1951. [19] L. A. Burke, S. A. Kinney, R. F. Dubiel, and J. K. Pitman, \Regional map of the 0.70 psi/ft pressure gradient and development of the regional geopressure-gradient model for the onshore and oshore gulf of mexico basin, usa," 2012. [20] B. Dindoruk, P. G. Christman, et al., \Pvt properties and viscosity correlations for gulf of mexico oils," SPE Reservoir Evaluation & Engineering, vol. 7, no. 06, pp. 427{437, 2004. [21] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blon- del, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Courna- peau, M. Brucher, M. Perrot, and E. Duchesnay, \Scikit-learn: Machine learning in Python," Journal of Machine Learning Research, vol. 12, pp. 2825{2830, 2011. [22] Schlumberger, \Pipesim user manual, version 2011.1," 2011. [23] M. Nikravesh and F. Aminzadeh, \Past, present and future intelligent reservoir characterization trends," Journal of Petroleum Science and Engineering, vol. 31, no. 2, pp. 67{79, 2001. [24] L. J. van der Maaten, E. O. Postma, and H. J. van den Herik, \Dimensionality reduction: A comparative review," Journal of Machine Learning Research, vol. 10, no. 1-41, pp. 66{71, 2009. 98 [25] H. S. Seung and D. D. Lee, \The manifold ways of perception," Science, vol. 290, no. 5500, pp. 2268{2269, 2000. [26] A. Srivastava, S. H. Joshi, W. Mio, and X. Liu, \Statistical shape analysis: Clus- tering, learning, and testing," Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 27, no. 4, pp. 590{602, 2005. [27] B. Fischl, A. Liu, and A. M. Dale, \Automated manifold surgery: constructing geo- metrically accurate and topologically correct models of the human cerebral cortex," Medical Imaging, IEEE Transactions on, vol. 20, no. 1, pp. 70{80, 2001. [28] J. S. Gomes, M. T. Ribeiro, C. J. Strohmenger, S. Naghban, M. Z. Kalam, et al., \Carbonate reservoir rock typing-the link between geology and scal," in Abu Dhabi International Petroleum Exhibition and Conference, Society of Petroleum Engineers, 2008. [29] N. M. Gomaa, S. R. Azer, D. E. Ouzzane, O. Saif, M. O. Okuyiga, D. F. Allen, D. Rose, R. Ramamoorthy, E. Bize, et al., \Case study of permeability, vug quan- tication, and rock typing in a complex carbonate," in SPE Annual Technical Con- ference and Exhibition, Society of Petroleum Engineers, 2006. [30] S. Masalmeh and X. Jing, \Carbonate scal: characterisation of carbonate rock types for determination of saturation functions and residual oil saturations," in SCA-08 presented at the SCA 2004 conference, Abu Dhabi, 2004. [31] A. C. Del Rey, M. Peres, A. Marcello, A. B. Barreto Jr, S. R. de Almeida, C. L. Sombra, et al., \Direct permeability estimation using production log," in SPWLA 50th Annual Logging Symposium, Society of Petrophysicists and Well-Log Analysts, 2009. [32] D. Moore et al., \Impact of super permeability on completion and production strate- gies," in Middle East Oil Show, Society of Petroleum Engineers, 1989. [33] B. Rashid, A. Muggeridge, A.-L. Bal, G. J. Williams, et al., \Quantifying the impact of permeability heterogeneity on secondary-recovery performance," SPE Journal, vol. 17, no. 02, pp. 455{468, 2012. [34] C. Uguru, H. Obiuwevbi, J. Oni, et al., \Impact of impermeable shale streaks on production," in Nigeria Annual International Conference and Exhibition, Society of Petroleum Engineers, 2011. [35] U. Ahmed, S. Crary, G. Coates, et al., \Permeability estimation: the various sources and their interrelationships," Journal of Petroleum Technology, vol. 43, no. 05, pp. 578{587, 1991. [36] J.-P. Delhomme, \The quest for permeability evaluation in wireline logging," Aquifer Systems Management: Darcy's Legacy in a World of Impending Water Shortage: Selected Papers on Hydrogeology 10, p. 55, 2007. 99 [37] C. E. Rasmussen, \Gaussian processes for machine learning," 2006. [38] J.-S. Lim, \Reservoir permeability determination using articial neural network," J. Korean Soc. Geosyst. Eng, vol. 40, pp. 232{238, 2003. [39] F. Aminzadeh and P. De Groot, Neural networks and other soft computing tech- niques with applications in the oil industry. Eage Publications, 2006. [40] M. Nikravesh, L. A. Zadeh, and F. Aminzadeh, Soft computing and intelligent data analysis in oil exploration, vol. 51. Elsevier, 2003. [41] A. Ouenes, \Practical application of fuzzy logic and neural networks to fractured reservoir characterization," Computers & Geosciences, vol. 26, no. 8, pp. 953{962, 2000. [42] R. Gholami, A. Shahraki, and M. Jamali Paghaleh, \Prediction of hydrocarbon reservoirs permeability using support vector machine," Mathematical Problems in Engineering, vol. 2012, 2012. [43] V. De Oliveira, B. Kedem, and D. A. Short, \Bayesian prediction of transformed gaussian random elds," Journal of the American Statistical Association, vol. 92, no. 440, pp. 1422{1433, 1997. [44] B. Mat ern, Spatial variation, vol. 36. Springer Science & Business Media, 2013. [45] S. Lafon and A. B. Lee, \Diusion maps and coarse-graining: A unied framework for dimensionality reduction, graph partitioning, and data set parameterization," Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 28, no. 9, pp. 1393{1403, 2006. [46] R. R. Coifman and S. Lafon, \Diusion maps," Applied and computational harmonic analysis, vol. 21, no. 1, pp. 5{30, 2006. [47] A. Singer, R. Erban, I. G. Kevrekidis, and R. R. Coifman, \Detecting intrinsic slow variables in stochastic dynamical systems by anisotropic diusion maps," Proceed- ings of the National Academy of Sciences, vol. 106, no. 38, pp. 16090{16095, 2009. [48] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker, \Geometric diusions as a tool for harmonic analysis and structure denition of data: Diusion maps," Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 21, pp. 7426{7431, 2005. [49] W. J. Welch, R. J. Buck, J. Sacks, H. P. Wynn, T. J. Mitchell, and M. D. Morris, \Screening, predicting, and computer experiments," Technometrics, vol. 34, no. 1, pp. 15{25, 1992. [50] T. J. Santner, B. J. Williams, and W. I. Notz, The design and analysis of computer experiments. Springer Science & Business Media, 2013. 100 [51] G. E. Box and K. Wilson, \On the experimental attainment of optimum conditions," Journal of the Royal Statistical Society. Series B (Methodological), vol. 13, no. 1, pp. 1{45, 1951. [52] S. C. Ferreira, R. Bruns, H. Ferreira, G. Matos, J. David, G. Brandao, E. P. da Silva, L. Portugal, P. Dos Reis, A. Souza, et al., \Box-behnken design: An alternative for the optimization of analytical methods," Analytica chimica acta, vol. 597, no. 2, pp. 179{186, 2007. [53] W. G. Cochran and G. M. Cox, \Experimental designs .," 1957. [54] M. K. Deb, I. M. Babu ska, and J. T. Oden, \Solution of stochastic partial dierential equations using galerkin nite element techniques," Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 48, pp. 6359{6372, 2001. [55] R. G. Ghanem and P. D. Spanos, Stochastic nite elements: a spectral approach. Courier Corporation, 2003. [56] R. Tipireddy and R. Ghanem, \Basis adaptation in homogeneous chaos spaces," Journal of Computational Physics, vol. 259, pp. 304{317, 2014. [57] R. Courant and D. Hilbert, Methods of mathematical physics, vol. 1. CUP Archive, 1966. [58] U. EIA, \Annual energy outlook 2015: With projections to 2040. also see associated statistical tables," 2015. [59] C. D. Rudd, A. C. Long, K. N. Kendall, and C. Mangin, Liquid moulding technolo- gies: Resin transfer moulding, structural reaction injection moulding and related processing techniques. Elsevier, 1997. [60] K. Potter, Resin transfer moulding. Springer Science & Business Media, 2012. [61] A. Isayev, \Injection and compression molding fundamentals," Department of Poly- mer Engineering, 2012. [62] R. L. Iman and M. J. Shortencarier, \Fortran 77 program and user's guide for the generation of latin hypercube and random samples for use with computer models," tech. rep., Sandia National Labs., Albuquerque, NM (USA), 1984. [63] S. Sourour and M. Kamal, \Dierential scanning calorimetry of epoxy cure: isother- mal cure kinetics," Thermochimica Acta, vol. 14, no. 1, pp. 41{59, 1976. [64] E. GROUP et al., \Pam-rtm user's guide & tutorials version 2015," ESI Group, Paris, 2011. 101 Appendix A Attribute Value Attribute Value Spatial Coordinates (X,Y) 984406.156520 1256785.517478 Sand Type O API Well Number 608114027400 Water Depth (ft) 4289 Well Name 1 Proved Recoverable Oil (bbl) 10548216 Well Name Sux ST00BP00 Proved Recoverable Gas (Mcf) 18419100 Operator Number 2010 Proved Recoverable BOE (bbl) 13825635 Bottom Field Name Code GC608 Cumulative Oil Produced (bbl) 8542494 Spud Date 20000320 Cumulative Gas Produced (Mcf) 13455042 Bottom Lease Number G18402 Cumulative BOE Produced (bbl) 10936629 Continued on next page ... 102 Table 4 { continued from previous page Attribute Value Attribute Value RKB Elevation (ft) 82 Proved Remaining Recoverable Oil (bbl) 2005722 Borehole Total Measured Depth (ft) 14500 Proved Remaining Recoverable Gas (Mcf) 4964058 True Vertical Depth (ft) 14498 Proved Remaining Recoverable BOE (bbl) 2889006 Surface North/South Distance (ft) 7201 Unproved Remaining Recoverable Oil (bbl) 0 Surface North/South Code N Unproved Remaining Recoverable Gas (Mcf) 0 Surface East/West Distance (ft) 5313 Unproved Remaining Recoverable BOE (bbl) 0 Surface East/West Code W Discovered Oil (bbl) 10548216 Surface Area GC Discovered Gas (Mcf) 18419100 Surface Block 608 Discovered BOE (bbl) 13825635 Continued on next page ... 103 Table 4 { continued from previous page Attribute Value Attribute Value Bottom North/South Distance (ft) 7140 Subsea Depth (ft) 11624 Bottom North/South Code N Total Average Net Thickness (ft) 130.52 Bottom East/West Distance (ft) 5313 Total Area (acres) 249 Bottom East/West Code W Total Volume (acre-ft) 32499 Bottom Area GC Oil Average Net Thickness (ft) 159.41 Bottom Block 608 Oil Total Area (acres) 189 Total Depth Date 20000421 Oil Total Volume (acre-ft) 30129 Status Date 20000502 Gas Average Net Thickness (ft) 39.5 Type Code E Gas Total Area (acres) 60 District Code 2 Gas Total Volume (acre-ft) 2370 Status Code ST Drive Type PAR Depth of Water 4320 Reservoir Type U Continued on next page ... 104 Table 4 { continued from previous page Attribute Value Attribute Value Surface Longitude -90.180154 Porosity 0.32 Surface Latitude 27.364676 Water Saturation 0.22 Bottom Longitude -90.18015 Average Permeability (mD) 627.087 Bottom Latitude 27.364844 Initial Pressure (psi) 7008 Surface Lease Number G18402 Initial Temperature (F) 112 Complex ID Number Sand Pressure Gradient (psi per ft) 0.603 Structure Number Sand Temperature Gradient (F per 100 ft) 0.623 6.08114E+11 Initial Solution Gas-Oil Ratio (scf/stb) 500 Sand Sequence Number 417187 Yield (stb/MMcf) 19.231 Sand Name 1382 GC608 Proportion Oil (decimal) 0.927 Assessed? Y Gas-Oil Ratio (Mcf/bbl) 1.746 Continued on next page ... 105 Table 4 { continued from previous page Attribute Value Attribute Value Sand Discovery Date 4/21/2000 Gas Specic Gravity (decimal at 60 degrees F and 15.025 psi) 0.67 Sand Discovery Year 2000 Oil API Gravity (API Units) 30 Sand Discovery Date High 5/25/2007 Initial Gas Formation Volume Factor (scf/cf) 361.2 Sand Discovery Year High 2007 Initial Oil Formation Volume Factor (bbl/stb) 1.566 Field Name GC608 Recoverable Oil per acre-foot (bbl/acre-ft) 347.657 Field Class PDP Recoverable Gas per acre-foot (Mcf/acre-ft) 1616.152 Field Status A Oil in Place (bbl) 37255236 Field Structure Code F Gas in Place (Mcf) 8326695 Field Primary Trap Code Q Oil Recovery Factor (decimal) 0.28 Continued on next page ... 106 Table 4 { continued from previous page Attribute Value Attribute Value Field Secondary Trap Code Oil Reservoirs' Recoverable Oil (bbl) 10474557 Field Discovery Date 4/21/2000 Oil Reservoirs Recoverable Gas (Mcf) 14588820 Field Discovery Year 2000 Produced GOR for Oil Reservoirs (Mcf/stb) 1.409 EIA Identication Number 846608 Gas Recovery Factor (decimal) 0.46 Planning Area CGM Gas Reservoirs Recoverable Oil (bbl) 73659 Sand Name M SERIES Gas Reservoirs Recoverable Gas (Mcf) 3830280 Play Number 1382 Produced GOR for Gas Reservoirs (Mcf/stb) 52 Play Name MUU F2 Count of Non-associated Gas Reservoirs in the Sand 1 Continued on next page ... 107 Table 4 { continued from previous page Attribute Value Attribute Value Pool Name 1382 GC608 Count of Under-saturated Oil Reservoirs in the Sand 3 Chronozone MUU Count of Saturated Oil Reservoirs in the Sand 0 Play Type F2 Count of Total Reservoirs in the Sand 4 Proved or Unproved P Table 4: All attributes for a sand in GC608. 108

Abstract (if available)

Abstract This dissertation focuses on developing probabilistic models for predicting the behavior of physical systems in the presence of uncertainty. Novel contributions are made to distinct aspects of the problem, namely, 1) risk assessment and risk mitigation models for the ultra-deepsea drilling in the Gulf of Mexico, 2) regression models for intrinsic variables constrained onto a manifold embedded in an ambient space, and, 3) efficient isometry computation for the stochastic basis adaptation method to detect lower dimensional stochastic space adapted to quantities of interest.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

A polynomial chaos formalism for uncertainty budget assessment

PDF

Analytical and experimental studies in modeling and monitoring of uncertain nonlinear systems using data-driven reduced‐order models

PDF

Spectral optimization and uncertainty quantification in combustion modeling

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

Accurate and efficient uncertainty quantification of subsurface fluid flow via the probabilistic collocation method

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

$Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure$

PDF

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

PDF

Stochastic and multiscale models for urban and natural ecology

PDF

Model falsification: new insights, methods and applications

PDF

Uncertainty management for complex systems of systems: applications to the future smart grid

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures

PDF

Controlled substructure identification for shear structures

PDF

Comprehensive uncertainty quantification in composites manufacturing processes

Asset Metadata

Creator Thimmisetty, Charanraj (author)

Core Title Risk assessment, intrinsic interpolation and computationally efficient models for systems under uncertainty

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Publication Date 04/15/2016

Defense Date 03/25/2016

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag OAI-PMH Harvest,uncertainty quantification

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ghanem, Roger (committee chair), Aminzadeh, Fred (committee member), Masri, Sami (committee member)

Creator Email charanrajt@gmail.com,thimmise@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-228151

Unique identifier UC11277879

Identifier etd-Thimmisett-4251.pdf (filename),usctheses-c40-228151 (legacy record id)

Legacy Identifier etd-Thimmisett-4251.pdf

Dmrecord 228151

Document Type Dissertation

Format application/pdf (imt)

Rights Thimmisetty, Charanraj

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA