University of Southern California Dissertations and Theses

Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

(USC Thesis Other)

Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNCERTAINTY QUANTIFICA TION AND DA TA ASSIMILA TION VIA TRANSFORM PROCESS FOR STRONGLY NONLINEAR PROBLEMS by Qinzhuo Liao A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) December 2014 Copyright 2014 Qinzhuo Liao ii Acknowledgements I would like to express my sincerest gratitude to those who have contribution to my research work in this dissertation. I would like to thank my advisor, Dr. Dongxiao Zhang for his passion for research, his depth of knowledge, his innovative ideas and his drive for excellence in the research. I feel extremely fortunate and blessed to have the opportunity of having him as my advisor and mentor. I would like to express my appreciation to Dr. Iraj Ershaghi for his kindness and support during my graduate studies. I would like to thank other members of my dissertation committee, Dr. Roger Ghanem and Dr. Behnam Jafarpour for their time, expertise, and helpful comments. I would like to extend my thanks to Dr. Kristian Jessen, Dr. Fred Aminzadeh, Dr. Felipe de Barros, Idania Takimoto, Andy Chen and other faculties and staffs in USC, who helped me throughout the academic pursuit. I would also like to thank the former and current students and post-doctors, Dr. Heng Li, Dr. Haibin Chang, Dr. Lingzao Zeng, Dr. Liangsheng Shi, Dr. Weixuan Li, Heli Bao, Le Lu, and many others for their caring and help. Most importantly, I am deeply grateful to my parents and my whole family for their unconditional love, support and encouragement which gave me strength and motivation. iii Table of Contents Acknowledgements .................................................................................................................. ii List of Figures .......................................................................................................................... vi Abstract ................................................................................................................................... xv Chapter 1 Introduction ....................................................................................................... 1 1.1 Background ................................................................................................................ 1 1.1.1 Uncertainty Quantification .................................................................................. 1 1.1.2 Data Assimilation ................................................................................................ 3 1.2 Motivation .................................................................................................................. 4 1.2.1 Strong Nonlinearity in Uncertainty Quantification ............................................. 4 1.2.2 Strong Nonlinearity in Data Assimilation ........................................................... 9 1.3 Objective of This Study ............................................................................................ 14 1.4 Outline ...................................................................................................................... 16 Chapter 2 Governing Equations ...................................................................................... 17 2.1 Multiphase Flow Model ........................................................................................... 17 2.2 Solute Transport Model ............................................................................................ 19 Chapter 3 Review of Basic Concepts ............................................................................... 21 3.1 Basic Concepts in Uncertainty Quantification ......................................................... 21 3.1.1 Karhunen-Loeve Expansion (KLE) .................................................................. 21 3.1.2 Polynomial Chaos Expansion (PCE) ................................................................ 22 3.1.3 Lagrange Interpolation ...................................................................................... 23 3.1.4 Smolyak Sparse Grid Algorithm ....................................................................... 23 3.1.5 Probabilistic Collocation Method (PCM) ......................................................... 25 3.1.5.1 Pre-processing ........................................................................................... 26 3.1.5.2 Polynomial Construction ........................................................................... 26 3.1.5.3 Post-processing .......................................................................................... 27 3.1.5.4 Framework of PCM ................................................................................... 28 3.2 Basic Concepts in Data Assimilation ........................................................................ 29 3.2.1 Ensemble Kalman Filter (EnKF) ...................................................................... 29 3.2.1.1 Forecast Stage ............................................................................................ 29 3.2.1.2 Update Stage .............................................................................................. 29 Chapter 4 Location-based Transformed Probabilistic Collocation Method (xTPCM) .. ........................................................................................................................... 31 4.1 Theoretical Framework ............................................................................................ 32 iv 4.1.1 Pre-processing ................................................................................................... 32 4.1.2 Forward Transform from Response to Location ............................................... 32 4.1.3 Polynomial Construction ................................................................................... 33 4.1.4 Backward Transform from Location to Response............................................. 35 4.1.5 Post-processing ................................................................................................. 36 4.1.6 Framework of xTPCM ...................................................................................... 36 4.2 Case Studies .............................................................................................................. 37 4.2.1 Simple Illustrative Example .............................................................................. 37 4.2.2 One-dimensional Homogeneous Case .............................................................. 40 4.2.2.1 Multiphase Flow Test ................................................................................ 40 4.2.2.2 Solute Transport Test ................................................................................. 49 4.2.3 One-dimensional Heterogeneous Case ............................................................. 54 4.2.4 Two-dimensional Homogeneous Case .............................................................. 58 4.3 Discussion ................................................................................................................ 61 Chapter 5 Displacement-based Transformed Probabilistic Collocation Method (dTPCM) ........................................................................................................................... 65 5.1 Theoretical Framework ............................................................................................ 67 5.1.1 Pre-processing ................................................................................................... 67 5.1.2 Forward Transform from Response to Displacement ....................................... 68 5.1.3 Construct Polynomial ........................................................................................ 74 5.1.4 Backward Transform from Displacement to Response ..................................... 76 5.1.5 Post-processing ................................................................................................. 77 5.1.6 Framework of dTPCM ...................................................................................... 78 5.2 Case Studies .............................................................................................................. 79 5.2.1 One-dimensional Example with Single Random Variable ................................ 79 5.2.2 Two-dimensional Example with Single Random Variable ................................ 84 5.2.3 Two-dimensional Example with Multiple Random Variables .......................... 88 5.2.4 Two-dimensional Example with Random Field ................................................ 91 5.2.4.1 Effect of Interpolation Level ..................................................................... 96 5.2.4.2 Comparison of Computational Effort ........................................................ 99 5.2.4.3 Effect of Simulation Time ......................................................................... 99 5.2.5 Three-dimensional Example with Random Field............................................ 101 5.3 Discussion .............................................................................................................. 106 5.3.1 Sparse Grid in High Dimensional Problems ................................................... 106 5.3.2 Piecewise Polynomial Approximation ............................................................ 108 5.3.3 Limitation of the dTPCM and Possible Solution ............................................ 109 Chapter 6 Time-based Transformed Probabilistic Collocation Method (tTPCM) ... 111 6.1 Theoretical Framework .......................................................................................... 113 6.1.1 Pre-processing ................................................................................................. 113 6.1.2 Forward Transform From Response to Time .................................................. 113 6.1.3 Polynomial Construction ................................................................................. 115 6.1.4 Backward Transform from Time to Response ................................................ 117 6.1.5 Post-processing ............................................................................................... 117 6.1.6 Framework of tTPCM ..................................................................................... 118 6.2 Case Studies ............................................................................................................ 119 6.2.1 One-dimensional Illustrative Example ............................................................ 119 v 6.2.2 Two-dimensional Solute Transport Test .......................................................... 123 6.2.3 Three-dimensional Multiphase Flow Test ....................................................... 129 6.3 Discussion .............................................................................................................. 134 Chapter 7 Transformed Ensemble Kalman Filter (TEnKF) ....................................... 138 7.1 Theoretical Framework .......................................................................................... 139 7.1.1 Forecast Stage ................................................................................................. 139 7.1.2 Forward Transform Stage ................................................................................ 139 7.1.3 Update Stage ................................................................................................... 145 7.1.4 Backward Transform Stage ............................................................................. 146 7.1.5 Framework of TEnKF ..................................................................................... 147 7.2 Case Studies ............................................................................................................ 148 7.2.1 One-dimensional Homogeneous Case ............................................................ 148 7.2.2 Two-dimensional Heterogeneous Case ........................................................... 159 7.2.2.1 Implementation Details of TEnKF .......................................................... 160 7.2.2.2 Estimation of Saturation and Permeability with Material Balance Check .... 164 7.2.2.3 Match and Prediction of the Observation ................................................ 168 7.2.3 Three-dimensional PUNQ-S3 Case ................................................................ 169 7.2.3.1 Estimation of Saturation, Porosity and Permeability ............................... 171 7.2.3.2 Match and Prediction of the Observation ................................................ 177 7.2.3.3 Computational Cost Comparison. ........................................................... 179 7.3 Discussion .............................................................................................................. 181 Chapter 8 Conclusion and Future Work ....................................................................... 187 8.1 Conclusion .............................................................................................................. 187 8.2 Future Work ............................................................................................................ 190 Bibliography ......................................................................................................................... 192 List of Publications .............................................................................................................. 200 vi List of Figures Figure 1.1: Oscillation problem for different functions: (a) f 1 (x) = 1 / (1 + 25x 2 ); (b) f 2 (x) = exp (-15x 2 ); (c) f 3 (x) = heaviside (x), (i.e., f 3 (x) = 0, x < 0; f 3 (x) = 1, x > 0; f 3 (x) = 1/2, x = 0); (d) f 4 (x) = (1 + 8erf(x) ) / 2, where “erf” denotes error function). ..... 5 Figure 1.2: Challenge of the polynomial approximation for nonlinear problems: (a) concentration profile as a function of location for different parameter values (deterministic runs); (b) concentration as a function of parameter at a fixed location; (c) saturation profile as a function of location for different parameter values (deterministic runs); (d) saturation as a function of parameter at a fixed location.8 Figure 1.3: Illustration of non-physical problem due to nonlinearity of state vector and observation: (a) saturation before Kalman filter update; (b) saturation after standard Kalman filter update; (c) saturation after Kalman filter update with normal score transform; (d) saturation after Kalman filter update with error function transform. 13 Figure 2.1: Illustration of advection dominated effects: (a) in multiphase flow, the small capillarity leads to a steep saturation profile (discontinuous shock in the extreme case); (b) in solute transport, the small dispersivity leads to a steep concentration profile (Dirac delta function in the extreme case). .............................................. 19 Figure 4.1: Illustration of the difference between the traditional PCM and the xTPCM.37 Figure 4.2: Illustration of nonlinear/discontinuous effect: (a) saturation as a function of location with different random parameters; (b) saturation as a function of parameter at a fixed location (solid line for MC reference, dash-dotted line for 2 nd order polynomial, dashed line for 4 th order polynomial, similarly in Figure 4.2e); (c) histogram of saturation at a fixed location; (d) location as a function of saturation with different random parameters; (e) location as a function of parameter for a particular saturation; (f) histogram of location for a particular saturation. The log permeability ln(k) = 0.3 θ, where θ ~ N(0,1). ...................................................... 39 Figure 4.3: Saturation profiles: directly generated from deterministic model runs by (a) MC method, (b) 2 nd order PCM, and (c) 4 th order PCM; and randomly generated from polynomial approximation by (d) 2 nd order xTPCM, (e) 2 nd order PCM, and (f) 4 th order PCM. ...................................................................................................... 42 Figure 4.4: Saturation as function of random parameter at locations: (a) x = 0.2, (b) x = 0.3, and (c) x = 0.4; and corresponding PDF of saturation at these locations: (d) x = 0.2, (e) x = 0.3, and (f) x = 0.4. The log porosity ln( ϕ) = 0.2 + 0.316 θ, where θ ~ N(0,1). .................................................................................................................. 43 vii Figure 4.5: Location as function of random parameter at particular saturation values: (a) S w = 0.2; (b) S w = 0.5; and (c) S w = 0.8; and corresponding location PDF at these saturations: (d) S w = 0.2; (e) S w = 0.5; and (f) S w = 0.8. The log porosity ln( ϕ) = 0.2 + 0.316 θ, where θ ~ N(0,1). ................................................................................... 45 Figure 4.6: Statistical moments of saturation by the MC method and collocation methods: (a) mean; (b) variance; (c) third order moment; and (d) fourth order moment. ... 48 Figure 4.7: Statistical moments of saturation with different input variabilities, 2 ln φ σ = 0.01: (a) mean and (d) variance; 2 ln φ σ = 0.1: (b) mean and (e) variance; 2 ln φ σ = 1: (c) mean and (f) variance. .................................................................................... 49 Figure 4.8: Concentration profiles: directly generated from deterministic model runs by (a) MC method, (b) 2 nd order PCM, and (c) 4 th order PCM; and randomly generated from polynomial approximation by: (d) 2 nd order xTPCM, (e) 2 nd order PCM, and (f) 4 th order PCM. ...................................................................................................... 51 Figure 4.9: Transformed variables as function of random parameter: (a) scaling factor c max , (b) first location x 1 , and (c) second location x 2 ; and corresponding PDFs: (d) scaling factor c max , (e) first location x 1 , and (f) second location x 2 . The log conductivity ln(K) = 0.316 θ, where θ ~ N(0,1). .................................................. 52 Figure 4.10: Concentration as function of random parameter at locations: (a) x = 7.5, (b) x = 10, and (c) x = 12.5; and corresponding PDFs of concentration at these locations: (d) x = 7.5, (e) x = 10, and (f) x = 12.5. The log conductivity ln(K) = 0.316 θ, where θ ~ N(0,1). ............................................................................................................ 53 Figure 4.11: Moments of concentration with different dispersivities, α L = 0.01 m: (a) mean; (b) variance; and α L = 1 m: (c) mean; (d) variance. .................................. 54 Figure 4.12: Statistical moments of saturation with different number of random dimensions N and number of collocation points P = 2N + 1. Mean saturation: (a) N = 1, P = 3, (b) N = 3, P = 7, (c) N = 10, P = 21; and saturation variance: (d) N = 1, P = 3, (e) N = 3, P = 7, (f) N = 10, P = 21. ................................................................. 57 Figure 4.13: Errors of saturation moment in various number of retained random dimensions N: (a) mean error; and (b) variance error. ......................................... 57 Figure 4.14: Statistical moments including the mean saturation from: (a) MC method; (c) 2 nd order PCM; (e) 4 th order PCM; and (g) 2 nd order xTPCM; and the saturation variance from: (b) MC method; (d) 2 nd order PCM; (f) 4 th order PCM; and (h) 2 nd order xTPCM. ...................................................................................................... 60 Figure 4.15: Illustration of the difficulty for the xTPCM when the location is a non-monotonic function of the response.............................................................. 62 viii Figure 5.1: Illustration of digital image correlation for time varying image sequence: (a) reference image; (b) deformed image with displacement vector; (c) normalized cross-correlation coefficient. The displacement vector is determined through searching for the maximum of the correlation coefficient. .................................. 69 Figure 5.2: Illustration of digital image correlation for stochastic analysis: (a) reference image; (b) deformed image with displacement vector; (c) normalized cross-correlation coefficient. The displacement vector is determined through searching for the maximum of the correlation coefficient. .................................. 72 Figure 5.3: Illustration of scattered data interpolation by triangulation. ..................... 77 Figure 5.4: Illustration of the difference between the traditional PCM and the dTPCM.79 Figure 5.5: 1D solute transport and multiphase flow examples: (a) concentration profiles with different collocation points; (b) different orders of polynomial approximation to the concentration response at a fixed location; (c) PDF estimation of the concentration response at a fixed location; (d) saturation profiles with different collocation points; (e) different orders of polynomial approximation to the saturation response at a fixed location; (f) PDF estimation of the saturation response at a fixed location. ................................................................................................................ 80 Figure 5.6: Illustration of forward transform: (a) concentration profile as undeformed image; (b) concentration profile as deformed image; (c) determining the displacement by correlation coefficient. .............................................................. 81 Figure 5.7: Illustration of polynomial constructions of the displacement and the target location value. ...................................................................................................... 82 Figure 5.8: Illustration of backward transform: (a) and (d) displacement from polynomial interpolation; (b) and (e) target location value from polynomial interpolation; (c) and (f) response in meshed grid obtained from that in scattered grid. (a), (b) and (c) are for one realization; (d), (e) and (f) are for multiple realizations. ......................... 83 Figure 5.9: Final results by the dTPCM: (a) 20 realizations randomly generated by the MC method and the dTPCM using the same parameters; (b) true concentration response approximated by the dTPCM; (c) exact PDF approximated by the dTPCM. .............................................................................................................................. 84 Figure 5.10: Forward transform from response to displacement: (a) undeformed response image; (b) deformed response image; and (c) normalized cross-correlation coefficient. The displacement is determined when the correlation coefficient reaches its maximum......................................................................................................... 85 ix Figure 5.11: Construct polynomials of the displacement, the target location value and the response: (a) x-direction displacement as a function of parameter in the dTPCM (the solid line indicates the MC reference, the dashed line indicates the second order polynomial approximation; the triangles indicate the collocation points); (b) y-direction displacement as a function of parameter in the dTPCM; (c) target location value as a function of parameter in the dTPCM; and (d) fixed location response as a function of parameter in the PCM for comparison. ....................... 86 Figure 5.12: Backward transform from displacement to response: (a) and (d) displacement field 0 (,; ) dx t θ and target location value 0 (,; ) x d st θ generated by polynomial approximation; (b) and (e) scattered data 00 (,;) (,;) xx d sdt s t θ θ += from definition; (c) and (f) gridded data (, ; ) x st θ interpolated from 0 (,;) x sdt θ + . .............................................................................................................................. 88 Figure 5.13: Mean concentration <c> and concentration variance σ c 2 by the MC method, as well as the PCM and dTPCM with sparse grid level 1, 2 and 3. ..................... 90 Figure 5.14: Mean concentration <c> and concentration variance σ c 2 compared between the MC method and the PCM/dTPCM with sparse grid level 1, 2 and 3. ........... 91 Figure 5.15: Schematic diagram of model setup in the multiphase flow test. Lines A 1 -A 2 , B 1 -B 2 , and C 1 -C 2 indicate where the mean and variance are compared in detail. Points D 1 , D 2 and D 3 indicate where the PDFs are compared in detail. .............. 93 Figure 5.16: Illustration of motion analysis at 500 days: (a) log permeability with parameter 3 ξ ; (b) log permeability with parameter 5 ξ ; (c) water saturation with parameter 3 ξ as a deformed image; (d) water saturation with parameter 5 ξ as a deformed image; (e) water saturation with parameter 1 ξ as the undeformed image, and the displacement between 1 () ξ w S and 3 () ξ w S ; and (f) water saturation with parameter 1 ξ as the undeformed image, and the displacement between 1 () ξ w S and 5 () ξ w S . ................................................................................................................. 95 Figure 5.17: Moments of water saturation by different methods: (a) and (b) for the MC method; (c) and (d) for the PCM; and (e) and (f) for the dTPCM. (a), (c), and (e) are mean saturations; and (b), (d), and (f) are saturation variances. .......................... 97 Figure 5.18: Moments of water saturation along lines A 1 -A 2 (x + y = 2000 ft), B 1 -B 2 (x + y = 1600 ft), and C 1 -C 2 (x + y = 1200 ft). (a), (b), and (c) represent mean; and (d), (e), and (f) represent variance. .................................................................................... 98 Figure 5.19: PDFs of water saturation at different locations: (a) D 1 (x = 600 ft, y = 1400 ft); (b) D 2 (x = 640 ft, y = 1360 ft); and (c) D 3 (x = 800 ft, y = 1200 ft). ............. 98 Figure 5.20: Moments of water saturation along line A 1 -A 2 (x + y = 2000 ft) at four simulation times: (a), (b), (c) and (d) represent mean; and (e), (f), (g) and (h) represent variance. ............................................................................................. 101 x Figure 5.21: Schematic diagram of model setup in the solute transport test. Lines A 1 -A 2 , B 1 -B 2 , and C 1 -C 2 indicate where the simulation results are compared in detail.103 Figure 5.22: Moments of concentration using different methods: (a) and (b) for the MC method; (c) and (d) for the PCM; and (e) and (f) for the dTPCM. (a), (c), and (e) are mean concentrations; and (b), (d), and (f) are concentration variances. ............ 104 Figure 5.23: Moments of concentration along the lines A 1 -A 2 (x, y = 5 m, z = 2.5 m), B 1 -B 2 (y, x = 12 m, z = 2.5 m), and C 1 -C 2 (y, x = 14 m, z = 2.5 m). (a), (b), and (c) represent for mean concentrations; and (d), (e), and (f) for represent concentration variances. ........................................................................................................... 105 Figure 5.24: Mean and variance with different dimensions N: (a) and (b) N = 1; (b) and (e) N = 3; (c) and (f) N = 10. ................................................................................... 108 Figure 6.1: Illustration of the difference between the traditional PCM and the tTPCM.119 Figure 6.2: Response concentration profiles: (a) generated from model runs for the MC method, (b) generated from model runs for the PCM, (c) observed from another point of view for the tTPCM, (d) generated from polynomial approximation in the PCM, and (e) generated from polynomial approximation in the tTPCM. ......... 121 Figure 6.3: Relations between model response, arrival time and random parameter: (a) response as a function of parameter at the fixed location x = 0.3 and time t = 0.3; (b) PDF of the response; (c) arrival time as a function of parameter at the fixed location x = 0.3 and response c = 0.5; and (d) PDF of the arrival time. .......................... 122 Figure 6.4: Statistical moments of concentration by the MC method, the PCM and the tTPCM: (a) mean; (b) variance; (c) third-order moment; and (d) fourth-order moment. ............................................................................................................. 123 Figure 6.5: Schematic diagram of model setup in the 2D solute transport test. Lines A 1 -A 2 , B 1 -B 2 and C 1 -C 2 indicate where the simulation results are compared in detail. 124 Figure 6.6: Normalization process: (a) absolute (original) concentration profiles generated from model runs; and (b) relative concentration profiles after normalization. .................................................................................................... 125 Figure 6.7: Moments of concentration at ln 0.3 K σ = , including the mean concentration from: (a) the MC method, (c) the PCM, and (e) the tTPCM; and the concentration variance from: (b) the MC method, (d) the PCM, and (f) the tTPCM. .............. 127 Figure 6.8: Moments of concentration at ln 0.3 K σ = along the line A 1 -A 2 (y = 5 m), B 1 -B 2 (x = 18 m) and C 1 -C 2 (x = 22 m). (a), (b) and (c) represent mean concentrations; and (d), (e) and (f) represent concentration variances. ............. 128 Figure 6.9: Moments of concentration at ln 0.5 K σ = along the line A 1 -A 2 (y = 5 m), B 1 -B 2 (x = 18 m) and C 1 -C 2 (x = 22 m). (a), (b) and (c) represent mean concentrations; and (d), (e) and (f) represent concentration variances. ............. 129 xi Figure 6.10: Schematic diagram of model setup in the 3D multiphase flow test with one injection well and three production wells. ......................................................... 130 Figure 6.11: Moments of water saturation at the top, middle and bottom layers by the MC method, the PCM, and the tTPCM. (a), (c), and (e) describe mean saturations; and (b), (d), and (f) describe saturation variances. ................................................... 132 Figure 6.12: Statistical properties of the water saturation: (a) mean saturation along the line (x = y, z = 50 ft); (b) saturation variance along the line (x = y, z = 50 ft); (c) saturation PDF at the point (x = 310 ft, y = 310 ft, z = 50 ft). ............................ 133 Figure 6.13: Water cut at production well P1: (a) 100 realizations generated from model runs in the MC method; (b) 100 realizations generated from polynomial approximation in the PCM; (c) 100 realizations generated from polynomial approximation in the tTPCM; and (d) water cut PDF at time t = 400 day. ........ 133 Figure 7.1: Illustration of digital image correlation: (a) undeformed image; (b) deformed image with displacement vector; and (c) normalized cross-correlation coefficient. The displacement vector is determined through searching for the maximum of the correlation coefficient. ....................................................................................... 140 Figure 7.2: A simple example of forward transform: (a) undeformed image; (b)-(g) deformed images with different displacement vectors; and (h) normalized cross-correlation coefficient. The dashed square is the template window; the solid square is the target window; and the dash-dotted square is the search area. ..... 143 Figure 7.3: Illustration of scattered data interpolation: (a) the locations (A’, B’, and C’) from the previous updated stage are usually scattered, but we need the value of the query point (D’) in the meshed gridblocks; (b) find the triangle (A’, B’, and C’) that encloses the query point (D’) and calculate the value of the query point (D) through computing the weighted sum of values of the three vertices of the enclosing triangle (A, B, and C). ..................................................................................................... 147 Figure 7.4: Illustration of the differences between the standard EnKF and the transformed EnKF (TEnKF). ................................................................................................. 148 Figure 7.5: Illustration of the non-physical problem: (a) pressure before the Kalman update; (b) pressure after the standard Kalman update; (c) saturation before the Kalman update; (d) saturation after the standard Kalman update; (e) saturation after the Kalman update with normal score transform (NST); (f) saturation after the Kalman update with error function transform (EFT); (g) saturation after the Kalman update with restart; and (h) saturation after the Kalman update with iteration once. The thin curves are the ensemble, and the thick curve is the reference. ............ 152 Figure 7.6: Illustration of the nonlinearity/non-Gaussianity that causes the non-physical problem in saturation estimation: (a) the relation between porosity and pressure; (b) the relation between porosity and saturation; (c) porosity before the Kalman update; (d) pressure before the Kalman update; (e) saturation before the Kalman update; (f) PDF of porosity before the Kalman update; (g) PDF of pressure before the Kalman xii update; (h) PDF of saturation before the Kalman update; (i) porosity after the Kalman update; (j) pressure after the Kalman update (the blue cross indicates the pressure directly updated by the Kalman correction, and the red circle indicates the pressure calculated through rerunning the model using the updated porosity from the initial time; similar in (k)); (k) saturation after the Kalman update. (l) PDF of porosity after the Kalman update; (m) PDF of pressure after the Kalman update (the blue line indicates the pressure directly updated by Kalman correction, and the red line indicates the pressure calculated through rerunning the model using the updated porosity from initial time; similar in (n)); (n) PDF of saturation after the Kalman update. The black thick line indicates the true reference. .................................. 154 Figure 7.7: Illustration of the TEnKF: (a) forecast saturation from direct model runs; (b) forward transform for one of the realizations; (c) detail of forward transform with enlarged scale; (d) forecast displacement for one of the realizations; (e) forecast displacement by forward transform; (f) displacement at x 0 = 0.2 before the Kalman update; (g) displacement at x 0 = 0.2 after the Kalman update; (h) updated displacement; (i) updated displacement for one of the realizations; (j) detail of forward transform with enlarged scale; (k) (j) updated saturation for one of the realizations; (f) updated saturation by backward transform. The blue curves show the ensemble, the black curve is the reference, and the red curve is the saturation generated from the ensemble mean of porosity as the undeformed image. ....... 157 Figure 7.8: Material balance check for different methods. The cross represents the average, the rectangle represents the p10 and p90 (the values at 10% and 90% of the distribution, respectively), and the error bar represents the minimum and maximum values. ................................................................................................................ 159 Figure 7.9: Schematic diagram of the model setup in the 2D test. Line A 1 -A 2 indicates the places where the results are compared in detail. ................................................ 160 Figure 7.10: Sensitivity analysis for selecting the appropriate template size: (a) saturation generated from the ensemble mean as the undeformed image; (b)-(d) displacement fields calculated by forward transform using different template sizes; (e) saturation generated in one realization as the original deformed image; (f)-(h) reproduced deformed image by backward transform using different template sizes. The goal is to check if the reproduced deformed image can match the original deformed image.163 Figure 7.11: Determine the appropriate template size and search range: (a) relative error vs. half template size in 20 realizations (thin blue lines) and the mean error (thick black line); (b) histogram of displacement in the x-direction; (c) histogram of displacement in the y-direction. ......................................................................... 163 Figure 7.12: Water saturation updated at day 1,000 and 2,000 using different methods: the first row is the reference and mean saturation in the whole space at day 1,000; the second row is the saturation along line A 1 -A 2 at day 1,000; the third row is the reference and mean saturation in the whole space at day 2,000; and the fourth row is the saturation along line A 1 -A 2 at day 2,000. In the first and third rows, the xiii saturation values are indicated by color, and the white area in (b) and (j) indicate that the values are out of the range from 0 to 1 (negative in this case). In the second and fourth rows, the thin curves show the ensemble prediction, and the thick black curve is the reference. .................................................................................................. 165 Figure 7.13: Material balance check for standard EnKF and TEnKF in 40 realizations from the first to the last assimilation step. ......................................................... 166 Figure 7.14: Log permeability in the initial guess and after the first and final update: the first column is the true log permeability; the second column is the mean log permeability by standard EnKF; the third column is the mean log permeability by REnKF; and the fourth column is the mean log permeability by TEnKF. ......... 168 Figure 7.15: The production data of P2 estimated from the initial guess and from the final updated permeability using the EnKF, the REnKF, and the TEnKF. The first row is the oil production rate; and the second row is the water cut. The thin curves show the ensemble prediction, and the thick black curve is the reference. ....................... 169 Figure 7.16: Water saturation of layer 5 in each assimilation step: the first row is the true saturation, calculated from the exact permeability and porosity; the second row is the mean saturation, estimated using the EnKF; the third row is the mean saturation, estimated using the REnKF; the fourth row is the mean saturation, estimated using the TEnKF with 2D motion analysis for each layer; and the fifth row is the mean saturation, estimated using the TEnKF with 3D motion analysis. ..................... 173 Figure 7.17: Gas saturation of layer 1 in each assimilation step: the first row is the true saturation, calculated from the exact permeability and porosity; the second row is the mean saturation, estimated using the EnKF; the third row is the mean saturation, estimated using the REnKF; the fourth row is the mean saturation, estimated using the TEnKF with 2D motion analysis for each layer; and the fifth row is the mean saturation, estimated using the TEnKF with 3D motion analysis. ..................... 174 Figure 7.18: Porosity of layer 1 in each assimilation step: (a) is the true porosity. The first row except for (a) is the mean porosity, estimated using the EnKF; the second row is the mean porosity, estimated using the REnKF; the third row is the mean porosity, estimated using the TEnKF with 2D motion analysis for each layer; and the fourth row is the mean porosity, estimated using the TEnKF with 3D motion analysis.176 Figure 7.19: Horizontal log permeability of layer 3 in each assimilation step: (a) is the true log permeability. The first row except for (a) is the mean log permeability, estimated using the EnKF; the second row is the mean log permeability, estimated using the REnKF; the third row is the mean log permeability, estimated using the TEnKF with 2D motion analysis for each layer; and the fourth row is the mean log permeability, estimated using the TEnKF with 3D motion analysis. ................. 177 Figure 7.20: History matching and prediction of PRO-11: the first row is from the initial guess; the second row is from the EnKF; the third row is from the REnKF; the fourth row is from the TEnKF with 2D motion analysis for each layer; and the fifth xiv row is from the TEnKF with 3D motion analysis. The first column is the well bottom-hole pressure; the second column is the well oil production rate; the third column is the well gas-oil ratio; and the fourth column is the well water cut. .. 178 Figure 7.21: Computational cost comparison between different methods: (a) CPU time; and (b) CPU time ratio, defined as the ratio of the time in each method to the time in the EnKF. ........................................................................................................... 180 Figure 7.22: Limitation of the TEnKF and possible solution: (a) water saturation at layer 15 in one realization of the SPE9 case; (b) water saturation at layer 15 in another realization of the SPE9 case; (c) aperture problem in motion analysis where we cannot determine the true displacement; (d) possible solution to the aperture problem in the SPE9 case by focusing on the water flooding in the x-direction; (e) water saturation at layer 9 in one realization of the Brugge case; and (f) water saturation at layer 9 in another realization of the Brugge case. Since the geometric shape of saturation is complicated in (e) and (f), it is difficult to obtain the displacement field. ............................................................................................. 183 xv Abstract Based on the polynomial approximation, the traditional probabilistic collocation method (PCM) approximates the model output response, which is a function of the random input parameter, from the Eulerian point of view in specific location. In some cases especially when the advection dominates, the model response has a strongly nonlinear profile with discontinuous shock or large gradient. This nonlinearity in space domain is then translated to the nonlinearity in random parametric domain, which causes oscillation and inaccurate estimation using the traditional PCM. To address this issue, a new transformed PCM (TPCM) is developed in this study, where the model response is now represented by the an alternative variable (e.g., the location of particular response value, or displacement vector, or arrival time of particular response value), which is relatively linear to the random parameter with a smooth profile. This alternative variable is then approximated by the polynomial, from which a sufficient number of location samples are randomly generated and transformed back to obtain the response samples. Finally, the statistical moments and probability density functions of model response are estimated from these samples. The advantage of the TPCM is demonstrated through applications to multiphase flow and solute transport in porous media. It is shown that the TPCM achieves higher accuracy for the statistical properties than does the PCM, and produces more reasonable xvi realizations without oscillation, while the computational effort is greatly reduced compared to the direct sampling Monte Carlo method. The ensemble Kalman filter (EnKF) has been widely used for data assimilation. It is challenging, however, when the relation of state and observation is strongly nonlinear. For example, near the flooding front in an immiscible flow, directly updating the saturation using the EnKF may lead to non-physical results. One possible solution, which may be referred to as the restarted EnKF (REnKF), is to update the static state (e.g., permeability and porosity) and rerun the forward model from the initial time to obtain the updated dynamic state (e.g., pressure and saturation). However, it may become time consuming, especially when the number of assimilation steps is large. In this study, we develop a transformed EnKF (TEnKF), in which the state is represented by displacement as an alternative variable. The displacement is first transformed from the forecasted state, then updated, and finally transformed back to obtain the updated state. Since the relation between displacement and observation is relatively linear, this new method provides a physically meaningful updated state without re-solving the forward model. The TEnKF is tested in history matching of multiphase flow in a one-dimensional homogeneous medium, a two-dimensional heterogeneous reservoir, and a three-dimensional PUNQ-S3 model. The case studies show that the TEnKF produces physical results without the oscillation problem that occurs in the traditional EnKF, while the computational effort is reduced compared to the REnKF. 1 Chapter 1 Introduction 1.1 Background 1.1.1 Uncertainty Quantification There has been increasing interest in uncertainty quantification of geophysical analysis in recent years, which estimates the statistics of the system outputs, given the probabilistic description of the model characteristics and inputs [Zhang, 2002; Xiu, 2010; Le Maître and Knio, 2010]. Monte Carlo (MC) method is one of the most common approaches, in which the model inputs are randomly sampled and simulated to obtain the corresponding outputs that can be further analyzed statistically. However, it often requires a large number of realizations to reduce the sampling errors, thus it may become prohibitively expensive, especially for large-scale problems [Ballio and Guadagnini, 2004]. As an attractive alternative, the polynomial chaos expansion (PCE) techniques have found increased use in uncertainty quantification over the past decade. First introduced by Wiener [1938], the PCE constructs a model response surface by polynomial of uncertain parameters and offers an efficient way of including nonlinear effects in stochastic analysis [Ghanem and Spanos, 1991]. For parametric distributions other than 2 Gaussian distribution (e.g., Beta, Gamma, Uniform, etc.), Xiu and Karniadakis [2002] developed generalized polynomial chaos (gPC), where the stochastic solutions are expressed as different types of orthogonal polynomials to achieve better convergence. In general, the PCE techniques can be categorized as intrusive and non-intrusive approaches [Le Maître and Knio, 2010]. The former approaches have been termed intrusive, as they reformulate the governing equations that require new solvers. The best-known method in this group is the stochastic finite element method, relied on Galerkin projection reformulations [Ghanem and Spanos, 1991], and it has been applied to many areas including fluid flow in porous media [Ghanem, 1998; Deb et al, 2001]. However, in stochastic finite element method, the PCE coefficients are governed by a set of coupled equations, which is difficult to solve when the number of coefficients is large. Non-intrusive approaches employ deterministic sampling of the forward model run that can be treated as a black box. Tatang et al. [1997] developed the probabilistic collocation method (PCM), also referred to as the stochastic collocation method (SCM), which is successfully applied to uncertainty quantification [Webster et al., 1996; Pan et al., 1997; Tatang et al., 1997; Le Maître et al., 2002; Mathelin and Hussaini, 2003; Babuska et al., 2007; Xiu, 2007; Li and Zhang, 2007; Muller et al., 2011] and sensitivity analysis [Fajraoui et al., 2011; Oladyshkin et al., 2012; Formaggia et al., 2013; Ciriello et al., 2013; Dai et al., 2013]. In the PCM, the output random field is approximated either by Lagrange interpolation [Mathelin and Hussaini, 2003; Babuska et al., 2007], or by orthogonal basis function, whose coefficients are solved through pseudo-spectral projection [Le Maître et al., 2002; Xiu, 2007] or matrix inversion/least-squares regression 3 [Isukapalli et al., 1998; Li and Zhang, 2007]. Initially, the PCM is found to be quite promising for low-dimensional systems based on tensor product. Later, the choices of collocation points are extended to the Smolyak sparse grid [Xiu and Hesthaven, 2005] to gain efficiency for high-dimensional problems. In general, the PCM derives an uncoupled system; thus, it is highly-parallelizable and can be easily implemented with existing simulators [Shi et al., 2009]. 1.1.2 Data Assimilation There has been great progress in data assimilation within atmospheric, oceanographic and petroleum sciences during past few decades [Aanonsen et al., 2009]. In data assimilation, one aims to include the information from observation to estimate parameters and forecast state variables. The Kalman filter (KF) [Kalman, 1960] is an efficient recursive filter that estimates the state of a linear dynamical system from noisy measurements. It is extended to work with nonlinear models through the extended Kalman filter (EKF), which linearizes the model around the estimated mean of the state. However, the application of the EKF to highly nonlinear or large-scale problems is infeasible. The ensemble Kalman filter (EnKF), introduced by Evensen [1994], is widely used for data assimilation and history matching. It has also been successfully applied to petroleum engineering in recent years [Nævdal et al., 2002; Gu and Oliver, 2005; Jafarpour and McLaughlin, 2007; Chang et al., 2010]. The EnKF estimates the Kalman gain from 4 ensemble and thus is more robust for nonlinear models. The EnKF consists of two sequential stages in each assimilation step: the first is called the forecast stage, in which the state variables are advanced forward based on the model dynamics; the second is called update stage, in which the state variables are adjusted via assimilating the observations [Oliver et al., 2008]. Two extensive reviews of the history matching progress in reservoir engineering can be found in Aanonsen et al. [2009] and Oliver and Chen [2011]. 1.2 Motivation 1.2.1 Strong Nonlinearity in Uncertainty Quantification It is well recognized that the global PCE-based techniques have an inherent assumption of smoothness of the output response with respect to input parameter [Ghanem and Spanos, 1991; Le Maître and Knio, 2010]. In some cases where the relation is strongly nonlinear with large gradient or even discontinuity, the PCE approximation could be inaccurate and suffer from oscillation (i.e., Gibbs phenomenon [Gottlieb and Shu, 1997]). To better demonstrate the effect of nonlinearity, Figure 1.1 shows four mathematical models approximated by polynomials through nine interpolation nodes on the interval [ −1, 1] based on three different algorithms (i.e., Clenshaw-Curtis, Gaussian-Legendre, equispace). Figure 1.1a illustrate the well-known Runge phenomenon that a polynomial 5 interpolation with equispaced nodes fails to converge for the simple function f 1 (x) = 1 / (1 + 25x 2 ). Although the oscillation problem can be alleviated somehow by non-uniform nodes, the mismatch still exists. Figure 1.1b shows similar problems with unphysical negative values for a Gaussian-distribution-like function. Such phenomena are also observed in Figures 1.1c and 1.1d for discontinuous or steep functions. In sum, the strong nonlinearity/unsmoothness including shock front (strong discontinuity) and large gradient (weak discontinuity) leads to the Gibbs oscillation, which seriously affects the convergence rate of PCE-based methods. Figure 1.1: Oscillation problem for different functions: (a) f 1 (x) = 1 / (1 + 25x 2 ); (b) f 2 (x) = exp (-15x 2 ); (c) f 3 (x) = heaviside (x), (i.e., f 3 (x) = 0, x < 0; f 3 (x) = 1, x > 0; f 3 (x) = 1/2, x = 0); (d) f 4 (x) = (1 + 8erf(x) ) / 2, where “erf” denotes error function). 6 The problem of strongly nonlinearity such as oscillation and inaccuracy are also observed in physical models and jeopardize the collocation methods. For example, in a piston-like displacement of two-phase immiscible flow, the saturation at a given location and time can take on only two possible values (0 and 1) with a discontinuous shock at the flooding front. The probability density function (PDF) of the saturation near the front is bimodal that is strongly non-Gaussian, even when the input log porosity is Gaussian distributed [Gu and Oliver, 2006; Chen et al., 2009]. Under this condition, when the traditional PCM tries to approximate the saturations by polynomials of Gaussian random variables (as the input is Gaussian), it may produce inaccurate estimations of the PCE coefficients. Therefore, the realizations generated from the PCE may provide unphysical values outside the boundaries of 0 and 1, and the statistical moments and the PDFs may deviate from the true values. Another similar instance in solute transport is that, the point concentration has a bimodal PDF with one of the modes being zero and the other close to the concentration at the source, with small dispersion in early times [Dagan, 1982; Bellin et al., 1994]. This non-Gaussianity therefore brings about difficulty in predicting the concentration statistics by traditional collocation method. It was observed by Zhang et al. [2010] that the PCM may produce unphysical realizations with negative concentrations under certain circumstances. Actually such problems occur not only in non-intrusive collocation method, but also in intrusive Galerkin method [Tryoen et al., 2010] and moment equation method [Liu et al., 2007; Jarman and Tartakovsky, 2008]. Figure 1.2a shows five concentration profiles from deterministic model runs with regard to different random inputs θ, in the shape of steep “bell curves”. This nonlinearity 7 in the space domain therefore causes the nonlinearity in the random domain (solid line in Figure 1.2b), if we analyze the concentration c at the specific location x = 0.3. Since the random parameter θ follows a standard normal distribution, it is natural to use Gauss-Hermite abscissas as collocation points. As we can see in Figure 1.2b, the polynomials from curve fitting or interpolation lead to severe oscillation, which still exists even with an 8 th order polynomial construction. As for multiphase flow test, Figure 1.2c shows five saturation profiles in the shape of discontinuous shock waves. If we analyze the saturation S w at the specific location x = 0.3, a large interpolation error and the oscillation problem will also be observed as in Figure 1.2d. 8 Figure 1.2: Challenge of the polynomial approximation for nonlinear problems: (a) concentration profile as a function of location for different parameter values (deterministic runs); (b) concentration as a function of parameter at a fixed location; (c) saturation profile as a function of location for different parameter values (deterministic runs); (d) saturation as a function of parameter at a fixed location. A straightforward approach to alleviate this problem of strong nonlinearity is increasing the order of PCE, which requires more collocation points and model runs. Lin and Tartakovsky [2009] observed inaccurate concentration moments with small dispersivity, which they attributed to insufficient collocation points. Although they utilized the sparse grid algorithm and set a relatively large initial plume to reduce the 9 discontinuity, hundreds of collocation points were still required in a very high-order polynomial approach for only four random dimensions. In consequence, the traditional PCM may become too complicated and computationally unaffordable. Besides, there is no guarantee that high order polynomials will provide satisfactory results, as the oscillation problem does not disappear as the order increases. To address the limited utility of collocation methods when dealing with strong nonlinearity, Wan and Karniadakis [2005] formulated a multi-element gPC method to discretize the space of random inputs into multiple non-overlapping subspaces, which is combined with the collocation method later by Foo et al. [2008]. Le Maître et al. [2004] presented a wavelet-based method, which naturally leads to localized decompositions. Abgrall [2007] proposed a piecewise reconstruction approach via essentially non-oscillatory techniques. These methods suggest the possibility of a robust behavior, albeit at the expense of a slower rate of convergence, owing to more number of model runs in the subdomains. Recently, Poëtte et al. [2009] developed an entropy based method, which will bind the oscillation to a certain extent, but the overshoot still exists. 1.2.2 Strong Nonlinearity in Data Assimilation In the EnKF, if the relations between state variables and observations are relatively linear, the state variables can be updated successfully. If the relations are strongly nonlinear, however, the updated dynamic state (e.g., saturation and concentration) may be non-physical and inconsistent with the updated static state (e.g., permeability and 10 porosity). For example, in water flooding with weak capillarity, the saturation profile can be almost discontinuous near the flooding front. In some cases (e.g., the injection rate is fixed, the porosity is a random parameter, and the bottom-hole pressure is observed), the relation between the saturation and the observation is also discontinuous. In other words, the nonlinearity in the space domain is translated to the nonlinearity in the stochastic domain. In addition, the latter nonlinearity leads to a non-Gaussian (multimodal) distribution of saturation, even if the porosity is Gaussian-distributed. The strong nonlinearity and non-Gaussianity, hence, cause non-physical saturation values to be less than zero or greater than one, when the EnKF is used for updating saturation in the front area. Since the EnKF is a sequential method, the updated non-physical saturation is used as the initial condition for forecasting in the following time steps, which may induce unsatisfactory results for later data assimilation and prediction. This problem has drawn wide attention and been discussed in previous works. Some approaches are based on updating the static state only and resolving the forward model to obtain the updated dynamic state. The most thorough method among them is to rerun the forward model from the initial time, and may be referred to as the restarted EnKF (REnKF) [Gu and Oliver, 2007; Thulin et al., 2007; Li and Reynolds, 2009; Chen and Oliver, 2010; Wang et al., 2010; Song et al., 2014]. Although this method guarantees a physical dynamic state that is consistent with the static state, the computational effort may increase, especially when the number of assimilation steps increases. Wen and Chen [2006] proposed a confirming EnKF, in which the forward model is not rerun from the initial time, but rather from the previous time step. This method ensures consistency 11 between the static state and the dynamic state at the present time, but may also introduce inconsistency in previous times [Zafari and Reynolds, 2007; Song et al., 2014]. Gu and Oliver [2006] presented an iterative EnKF (IEnKF) [Liu and Oliver, 2005; Li and Reynolds, 2009; Lorentzen and Naevdal, 2011], in which an extra iteration step is added after rerunning the model from the previous time step. They reported that the state variables become more realistic, but the oscillation problem may still exist. Therefore, they suggest repeating the iteration step until the corrected state variables satisfy the physical bounds; however, the repetition process requires additional time. Li and Reynolds [2009] presented both the REnKF and the IEnKF algorithms, whose results were superior in comparison to standard EnKF. Krymskaya et al. [2009] proposed a method in which the mean of the estimated parameters at the end of the EnKF were used to initialize the ensemble for the next iteration. Song et al. [2014] evaluated the confirming EnKF and the REnKF, and proposed a modified REnKF to reduce computational cost by calculating only the mean simulation, not all the ensemble realizations, from initial time. Another class of approaches to alleviate the oscillation problem is based on variable transforms. One common transform is the normal score transform (NST) [Goovaerts, 1997; Bertino et al., 2003; Gu and Oliver, 2006; Zhou et al., 2011; Schöniger et al., 2012], which is constructed from the empirical cumulative density function of gridblock saturations. Another transform is the error function transform (EFT) [Jafarpour and McLaughlin, 2009], which is similar to the NSF, but from an analytical point of view. 12 Both methods ensure reasonable saturations inside the physical bounds, but the saturations may still oscillate spatially between high and low values. Figure 1.3a shows the saturation profiles before Kalman filter update, which are from the direct model simulation (we use Eclipse black-oil simulator here). Figure 1.3b shows the saturation profiles after Kalman filter update, with negative saturation values. Figure 1.3c shows the saturation profiles after Kalman filter update with normal score transform (NST). Although the NST successfully constraints the saturation to be non-negative. Some non-physical phenomena (e.g., low saturation value in upstream and high saturation value in downstream) are still observed. Similar problem (Figure 1.3d) exists in the saturation profiles after Kalman filter update with the error function transform (EFT), which is defined as S erf = erf -1 (2S − 1) and S = erf (S erf + 1) / 2. Actually, the EFT can be considered as a “global” NST, as the EFT transforms the state in various locations by a same analytical function, while the NST transforms the state in various locations by individual functions from the distribution of samples. 13 Figure 1.3: Illustration of non-physical problem due to nonlinearity of state vector and observation: (a) saturation before Kalman filter update; (b) saturation after standard Kalman filter update; (c) saturation after Kalman filter update with normal score transform; (d) saturation after Kalman filter update with error function transform. Another class of approaches is related to the idea of using alternative variables. Gu and Oliver [2006] presented an EnKF with changed variables, in which the location of the saturation front was used in the state variables instead of the actual saturation itself. This method improves the linearity of the relation between the state and observation, and works well for one-dimensional problems; however, it is unclear how to extend this method to two- or three-dimensional (2D or 3D) problems in which the front location is not well defined. Chen et al. [2009] introduced an EnKF with reparameterization, which used the arrival time of the saturation front as a state variable. This method can be applied 14 to 2D or 3D problems, but the simulator must be run beyond the current assimilation time to make a time window for interpolation. Butler and Juntunen [2012] proposed an idea related to principal analysis and spectral projection to construct a basis for state representation. They used a different basis to approximate solutions by constructing maps between the numerical basis and the representation basis, which was interpreted as a reparameterization of the state. Some other approaches with regard to taking the physical constraints into account [Phale and Oliver 2011] will not be discussed in detail here. 1.3 Objective of This Study In this research for uncertainty quantification, we propose a series of new transformed probabilistic collocation methods (TPCMs) to quantify the uncertainty accurately and efficiently under strongly nonlinear conditions. The original model response is transformed to a new variable, which is easier to analyze. This variable can be location, displacement or arrival time. Therefore, there are three approaches: location-based, displacement-based and time-based transformed probabilistic collocation method (xTPCM, dTPCM and tTPCM). The main difference between the traditional PCM and the new TPCM is that, the PCM approximates the model response directly, while the TPCM approximates an alternative variable that is connected to the model response through transform. Since the alternative variable is much more linear to the input parameter than the response, this transform process greatly improves the performance of collocation method. Case studies for multiphase flow and solute transport models are 15 conducted to demonstrate the advantage of TPCM over the PCM with higher accuracy and the same level of computational effort. While the results from the MC method are used as benchmarks to compare the effectiveness of the collocation methods. For data assimilation, we also propose a transformed ensemble Kalman filter (TEnKF) to assimilate the data accurately and efficiently under strongly nonlinear conditions. The main difference between the traditional EnKF and the TEnKF is that the former updates the state variables directly, while the latter updates the displacement as an alternative variable that is connected to the state through transform. Since the relation between displacement and observation is much more linear than the relation between state and observation, the transform process greatly improves the performance of the Kalman correction. The new TEnKF can be viewed as an extension of the EnKF with changed variables from saturation to location [Gu and Oliver 2006], as both methods are based on the linearity and Gaussianity of the transformed variables (i.e., displacement and location). The main difference is that the EnKF with changed variables from saturation to location, is limited to one-dimensional problems, since the locations are non-unique in 2D or 3D problems; whereas, the TEnKF can be applied to one- to three-dimensional problems, in which the displacement is either a scalar (in 1D) or a vector (in 2D or 3D). Through case studies of multiphase flow problems, we demonstrate that the TEnKF provides more physical results than does the traditional EnKF, and requires less computational effort than the REnKF, but has the same level of accuracy. 16 1.4 Outline This work is organized as below: in Chapter 2, we introduce the multiphase flow and solute transport models. in Chapter 3, we review some basic concepts including the traditional PCM and EnKF. In Chapter 4, we propose the location-based transformed PCM (xTPCM). In Chapter 5, we propose the displacement-based transformed PCM (dTPCM). In Chapter 6, we propose the time-based transformed PCM (tTPCM). In Chapter 7, we propose the transformed EnKF (TEnKF). Finally in Chapter 8, we provide conclusion and future work. 17 Chapter 2 Governing Equations In this study, the properties of the porous media, such as permeability or porosity, is treated either as homogeneous random constant or as heterogeneous random field. The saturation in multiphase flow or concentration in solute transport is considered as model response. Thus, the governing equation becomes stochastic partial differential equation whose solution is no longer deterministic but probabilistic characterized by statistical properties. 2.1 Multiphase Flow Model The oil/water two phase flow model for saturation tests is expressed by the following continuity equation as [Bear, 1974] () () (),, x x ri i ii i i kS p gz q i wo t κ ρφ μ ⎡⎤ ∂ ∇⋅ ∇ − ∇ + = = ⎢⎥ ∂ ⎣⎦ (2.1) where x and t denote the location and time, respectively, w and o denote the two phases (i.e., water and oil), φ is the porosity of the media, ,,, ii i i p S μ ρ are the viscosity, pressure, density and saturation of phase i, respectively, κ is the absolute permeability, k ri is the relative permeability of phase i, which is a function of S i , q i is the source or sink term, g is the gravitational acceleration, and z is the depth. Equation (2.1) is coupled with 18 1 () wo cw o w SS pSp p + = = − (2.2) where c p is the capillary pressure, which is a function of S w . The above equations can be simplified into a nonlinear advection-diffusion equation. When the capillary term is small compared to the advection term, the equation is hyperbolic and advection dominated [Zhang, 2002]. Figure 2.1a shows the effect of capillarity on the saturation profile in a simple 1D case. With large capillarity, the shape of saturation curve will be relatively smooth shown as the dashed line. While with small capillarity, the curve would indicate the movement of a persistent sharp front between the two phases as the solid line. It is pointed out that the capillarity may be neglected for studying oil/water displacement in oil reservoirs [Merle, 1981], which leads to a pure advection problem and usually called Buckley-Leverett displacement [Buckley and Leverett, 1942]. In the homogeneous case, analytical solutions to this problem are well-known by the method of characteristics with a discontinuous shock wave [Merle, 1981]. This discontinuity in space domain will be translated to the discontinuity in random domain, which makes it extremely difficult to quantify the saturation uncertainty by traditional collocation methods. Similar problems are observed in heterogeneous cases with spatially correlated random field, where the model is implemented with the commercial Eclipse black-oil simulator. 19 Figure 2.1: Illustration of advection dominated effects: (a) in multiphase flow, the small capillarity leads to a steep saturation profile (discontinuous shock in the extreme case); (b) in solute transport, the small dispersivity leads to a steep concentration profile (Dirac delta function in the extreme case). 2.2 Solute Transport Model The solute transport model for concentration tests under advection and dispersion is given as [Bear, 1974] () (, ) (, ) ( , ) (, ) ( , ) x xx vx x ij ct Dt c t t c t t ∂ =∇ ⋅ ∇ − ∇ ∂ (2.3) where x is the location, t is the time, c is the solute concentration, D ij is the hydrodynamic dispersion tensor, v is the pore water velocity, computed by (, ) ( , )/ vx u x tt φ = . And φ is the porosity, u is the Darcy velocity, computed by () () () ux x x Kh =−∇ , where () x K is the hydraulic conductivity and () x h is the hydraulic head. The hydrodynamic dispersion tensor is given by () (, ) (, ) vx vx ij ij L T T ij d ij vv DtD t ααα δ δ =− + + (2.4) 20 where L α is the longitudinal dispersivity, T α is the transverse dispersivity, (, ) vx t is the magnitude of the pore velocity, d D is the molecular diffusion coefficient and ij δ is the Kronecker operator ( 1, ij ij δ== ; and 0, ij ij δ = ≠ ). Figure 2.1b demonstrates the effect of dispersion by a synthetic test where a thin slug of pollutant is released initially. With large dispersivity, the concentration profile in the dashed line is relatively smooth. While with small dispersivity, the profile in the solid line has a large gradient. In the absence of dispersion, it would be like a Dirac delta function with strong discontinuity. This advection dominated phenomenon in space domain will lead to a strong nonlinearity in random domain, and thus strong non-Gaussianity for model response. As pointed out by Dagan [1982] and Bellin et al. [1994] that the point concentration has a bimodal distribution with small dispersion at early times. Zhang et al. [2010] reported that the traditional PCM may produce unphysical oscillatory realizations with negative concentrations under certain circumstances due to the nonlinearity and non-Gaussianity. Such problems can be addressed by the new TPCM, which will be discussed in detail in the following chapters. 21 Chapter 3 Review of Basic Concepts 3.1 Basic Concepts in Uncertainty Quantification 3.1.1 Karhunen-Loeve Expansion (KLE) Let the input m be a random function of space coordinate x and time coordinate t denoted by (, ; ) (, ) ( , ; ) xx x mt mt m t θ θ ′ =+ , where (, ) x mt is the mean and (, ; ) x mt θ ′ is the fluctuation. The covariance is described by 12 1 2 ( , ,) ( ,;) ( ,;) Cx x x x m tm t m t θ θ ′′ = , which can be decomposed as [Ghanem and Spanos, 1991] 12 1 2 1 (, ,) (,) ( ,) Cx x x x mnnn n tf tf t λ ∞ = = ∑ (3.1) where n λ and (, ) x n f t are the eigenvalues and eigenfunctions, respectively, and can be solved from the Fredholm equation of the second type. Thus the random field (, ) x m θ can be expressed by the Karhunen-Loeve expansion (KLE) in the truncated form as 1 (, ; ) (, ) ( , ) ( ) xx x N nn n n mt mt f t θ λξθ = =+ ∑ (3.2) where n ξ are the independent random variables. 22 3.1.2 Polynomial Chaos Expansion (PCE) Assume the model output s is treated as a random process (, ; ) x st θ . Since the process is usually non-Gaussian and the covariance structure of the output field is not known in advance, the output cannot be expanded using the KLE. Instead, the polynomial chaos expansion (PCE) as a more general representation for the random field may be used to express the output as [Ghanem and Spanos, 1991] 1 11 12 1 2 112 01 2 111 ˆ(, ; ) (, ) ( , ) ( ( )) (, ) ( ( ), ( )) ... xx x x i ii ii i i iii st t t t θα α ξ θ α ξ θ ξ θ ∞∞ === =+ Γ + Γ + ∑∑∑ (3.3) where 0 (, ) x t α and 12 ... (, ) x j ii i t α are deterministic coefficients, the basis 12 ( , ,..., ) j ji i i ξξξ Γ are a set of polynomial chaos with respect to the independent random variables 12 , ,..., j ii i ξξ ξ [Wiener, 1938]. For independent Gaussian random variables, 12 ( , ,..., ) j ji i i ξξξ Γ are the multidimensional Hermite polynomials of degree j [Ghanem and Spanos, 1991]. In case of other random distributions, generalized polynomial chaos (gPC) [Xiu and Karniadakis, 2002] can be used to represent the random field (e.g., Legendre-Chaos for uniform distribution, Jacobi-Chaos for Beta distribution). Equation (3.3) is usually truncated by a finite random dimension N and degree M as 1 ˆ(, ; ( )) (, ) ( ( )), ( )!/( ! !) x ξ x ξ Q ii i st a t Q N M NM θθ = =Ψ =+ ∑ (3.4) where (, ) x i at are the PCE coefficients, 11 ( , ,..., ) ξ N T ii i ξξ ξ = and (( )) ξ i θ Ψ are the orthogonal basis functions such that , ij ij δ Ψ Ψ= . In general, the convergence rate of the PCE depends on the smoothness of the function s in terms of random parameter ξ. 23 3.1.3 Lagrange Interpolation In the univariate case, let {} 1 ξ P i i = ∈ Ω be a set of prescribed nodes in the random space Ω, where P is the number of nodes. A Lagrange interpolation of the model response can be written as 1 ˆ(, ; ) (, ; ) ( ) x ξ x ξ ξ P ii i st st L = = ∑ (3.5) where {} 1 () ξ P i i L = are the corresponding Lagrange interpolation basis function, defined by 1 () ξ ξ ξ ξ ξ P j i j ij ji L = ≠ − = − ∏ (3.6) which satisfies () ,1 , ξ ij ij LijP δ =≤ ≤ . In the univariate case, the number of interpolation nodes is the same with the number of interpolation bases (also the number of orthogonal bases). Therefore, the interpolation approach is equivalent to a matrix inversion process, because of the uniqueness of the univariate interpolation structure [Xiu, 2010]. While in the multivariate case, the model response can be interpolated from tensor product construction [Babuska et al., 2007] or sparse grid construction [Xiu and Hesthaven, 2005]. 3.1.4 Smolyak Sparse Grid Algorithm The selection of collocation points is an important component of collocation method, especially in high dimensional problems. Although tensor product construction makes 24 mathematical analysis accessible [Babuska et al., 2007], the total number of nodes (M + 1) N grows exponentially for M-degree polynomial construction as the number of random dimensions N increases (i.e., the “curse of dimensionality”). To reduce the computational burden, the sparse grid strategy based on the Smolyak algorithm is employed by a linear combination of product formulas as [Xiu and Hesthaven, 2005] () () 1 1 1 (, ) 1 N qi i i qN i q N qN U U qi − −+≤ ≤ − ⎛⎞ Α= − ⋅ ⊗⊗ ⎜⎟ − ⎝⎠ ∑ (3.7) where U j is the one-dimensional interpolation formula for dimension j, N is the number of dimensions, the sparseness parameter q determines the order of the formula, and |i| = i 1 + i 2 + … + i N . Set k = q – N as the “level” of the Smolyak construction. To compute A(q, N), we only need to evaluate function on the sparse grid ( ) 1 11 1 (, ) N i i qN i q qN −+≤ ≤ = Θ× ×Θ ∪ H (3.8) where 1 j i Θ is the one-dimensional collocation point set for dimension index j, and (, ) qN H denotes the N-dimensional collocation point set. For functions in space { } :[ 1,1] | , , i MN Nj F f R f continuous i M j =− → ∂ ≤ ∀ , the interpolation error follows (2)( 1)1 , (, ) (log ) MMN NNl IAqN C P P − +−+ ∞ −≤⋅⋅ (3.9) where I N is the identity operator in a N-dimensional space, and P is the number of interpolation points. We can see that the convergence rate depends weakly on the dimension N but strongly on the smoothness index M. It is illustrated by Barthelmann et al. [2000] that the approximation result is poor for unsmooth function and completely useless for discontinuous function. In sum, the sparse grid strategy can somehow resolve 25 the problem termed as “curse of dimensionality”, but the problem of strong nonlinearity still exists and calls for proper solutions. 3.1.5 Probabilistic Collocation Method (PCM) The collocation method for solving differential equations can be derived from weighted residual method. Consider a stochastic partial differential equation (, ; ) (, ; ) xx st f t θ θ = L (3.10) where L is a differential operator, s is the model output, f is an known function and θ belongs to the event space Ω. Define a residual R as ˆ (, ; ) (, ; ) (, ; ) xx x Rtst ft θ θθ = − L (3.11) The notion in the weighted residual method is to solve equation (3.10) in a weak form by forcing the residual to zero in some average sense over the domain. That is ( ,; ) () () 0 x Rt d θω θ ρ θ θ Ω = ∫ (3.12) where ω is the weight function, and ρ is the probability density function. In the collocation method, the weight function is chosen as the Dirac delta function (( )) (( ) ( )) ξξ ξ jj ω θδ θ θ =− , which leads to an uncoupled system as (, ; ) 0 x ξ j Rt = (3.13) 26 3.1.5.1 Pre-processing In the KLE based PCM, the model input m is expanded by equation (3.2). To determine the PCE coefficients or to construct the Lagrange interpolation, we choose P sets of collocation points , 1,2,..., ξ j jP = for the random process of input, where 12 ( , ,..., ) ξ T N ξξ ξ = . By solving the model P times correspondingly, we get P sets of output (,; ), 1,2,..., x ξ j st j P = . 3.1.5.2 Polynomial Construction In the traditional PCM, the model response is approximated by the polynomial through two major approaches: the orthogonal basis approach [Le Maître et al., 2002; Xiu, 2007], and the Lagrange interpolation approach [Mathelin and Hussaini, 2003; Babuska et al., 2007]. In the former orthogonal basis approach, the model response s is approximated by orthogonal polynomials as equation (3.4). According to equation (3.13), we set ˆ(, ; ) (, ; ) x ξ x ξ jj st st = and obtain 1 ( , ) ( ) ( , ; ), 1, 2,..., x ξ x ξ Q iij j i at s t j P = Ψ= = ∑ (3.14) To determine the PCE coefficients ( , ), 1, 2,..., x i ati Q = , a pseudo-spectral projection process can be applied as [Xiu, 2007] 27 1 22 (, ; ) ( ) (, ; ) ( ) ( ) (, ; ), ( ) (, ) (), ( ) xξξ x ξξ ξξ xξξ x ξξ Q ji j j i j i i ii ii st w st d st at ρ = Ω Ψ Ψ Ψ == ≈ ΨΨ ΨΨ ∑ ∫ (3.15) where { } 1 , ξ Q jj j w = are a set of nodes and associated weights such that the summation approximates the integral. An alternative process is to treat equation (3.14) as a linear system in the matrix form [Li and Zhang, 2007] Za s = (3.16) where 12 1 2 ( , ,..., ) , ( , ,..., ) as TT QP aa a s s s == and ( ) Z ξ ji i j = Ψ is a PQ × Vandermonde-like matrix. For a well-posed problem, PQ ≥ is required. Thus the PCE coefficients could be attained by weighted least-squares regression as ( ) 1 a Z WZ WZ s TT − = (3.17) where W is a diagonal weight matrix with W jj = w j . It can be shown that the solutions by equation (3.15) and (3.17) are equivalent, with the same set of collocation points and orthogonal bases. While in the Lagrange interpolation approach, the model response is approximated by Lagrange polynomials as equation (3.5). 3.1.5.3 Post-processing In the traditional PCM, the moments of model response such as mean and variance can be directly calculated from the PCE coefficients obtained from equation (3.15) or (3.17) based on the orthogonality of basis function i Ψ as 28 222 1 2 , Q s ii i sa a σ = = =Ψ ∑ (3.18) An alternatively way is to first generate a sufficient number of response samples from equation (3.4) or (3.5), then estimating the moments from the ensemble as () () () 2 2 11 1 11 1 ,, ss s NN N nn is i i ii i ss s ss ssEss ss NN N σ == = ⎡⎤ == − − = − ⎣⎦ ∑∑ ∑ (3.19) where N s is the ensemble size. The advantage of this sampling approach is that we can produce the output realizations and obtain the PDFs in the process. 3.1.5.4 Framework of PCM The three major steps involved in the PCM: (1) pre-processing, i.e., representing the model input in terms of independent random parameters, selecting the collocation points and performing deterministic model runs at these points (section 3.1.5.1); (2) constructing polynomial of the model response, i.e., approximating the response by orthogonal basis function or Lagrange interpolating polynomial (section 3.1.5.2); (3) post-processing, i.e., estimating the statistical properties of the model response either from the PCE coefficients directly or from sampling (section 3.1.5.3). 29 3.2 Basic Concepts in Data Assimilation 3.2.1 Ensemble Kalman Filter (EnKF) The EnKF, which is essentially a Monte Carlo based Kalman filter, is usually implemented with the introduction of an augmented/jointed state vector [Oliver et al., 2008]. This augmented state s consists of the static state m and the dynamic state r in the form of [, ] sm r TTT = . In the EnKF, the augmented state is represented by an ensemble of realizations, involving a forecast stage and an update stage. 3.2.1.1 Forecast Stage The forecast stage is based on the solution of the dynamic equations for flow and transport in the reservoir as ,1, ( ), 1, 2,..., ss fu kj k j gj Ne − == (3.20) where g( x ) is the reservoir simulator; f represents “forecast”; u indicates “update”; k is the time step index; j is the ensemble member index; and Ne is the ensemble size. 3.2.1.2 Update Stage The update stage is the analysis or assimilation stage, in which the state variables are updated to honor the observations. The mean s f k and covariance k P can be estimated from the ensemble as 30 , 1 1 () ss s Ne f ff kk kj j E Ne = == ∑ (3.21) ()() ,, 1 1 1 Pssss Ne T ff ff kkjkkjk j Ne = =− − − ∑ (3.22) The Kalman gain is then calculated from the covariance as 1 () KPH HPH R TT kk k − =+ (3.23) where [, ] H0I = is the observation operator; and R is the covariance matrix of the observation errors. Each ensemble member is updated via ,, , , ( ), 1, 2,..., ss Kobs Hs uf f kj kj k k j k j j Ne =+ − = (3.24) where obs is the perturbed observation. 31 Chapter 4 Location-based Transformed Probabilistic Collocation Method (xTPCM) As we can see from Chapter 1 that the direct approximation of model response by the polynomial is not suggested for strongly nonlinear problems. A possible solution is to find an alternative variable to represent the model response, which is relatively linear to the random parameter. Intuitively inspired by the Lagrangian point of view, the location of particular response value is a good candidate for transform from the Eulerian coordinates, especially for the advection dominated problems. For 1D problem, it is actually a process of finding the iso-points. If the output value is a monotonic function of location in space domain, the approach is relatively simple (e.g., linear interpolation). Otherwise, there could be non-unique locations for a given response value, and these multiple locations can be analyzed individually. For 2D or 3D problems, the transform approach leads to searching for the iso-lines/contours or iso-surfaces, respectively. In some cases (e.g., homogenous media), it is possible to parameterize the response interfaces (iso-lines or iso-surfaces) by certain simple geometrical shapes. Unfortunately under general conditions, the interfaces are usually complicated and difficult to interpret. Possible solutions including motion analysis via digital image processing techniques will be presented in Chapter 5, and other solutions such as level set method are ongoing. 32 This chapter presents the work in published paper by Liao and Zhang [2013]. 4.1 Theoretical Framework The probabilistic collocation method with transform by location (xTPCM) can also be derived from the weighted residual method as equations (3.10)-(3.13). 4.1.1 Pre-processing In the KLE based xTPCM, the model input m is also expanded by equation (3.2). To determine the PCE coefficients or to construct the Lagrange interpolation, we choose P sets of collocation points , 1,2,..., ξ j jP = for the random process of input, where 12 ( , ,..., ) ξ T N ξξ ξ = . By solving the model P times correspondingly, we get P sets of output (,; ), 1,2,..., x ξ j st j P = . Actually, the collocation points in the xTPCM could be exactly the same with those in the PCM. Therefore, we can improve the existing PCM results by the xTPCM without additional model runs. For better comparison, we will use the same collocation points in the PCM and the xTPCM for the case studies. 4.1.2 Forward Transform from Response to Location One issue needs to be considered is that, if the minimum or maximum value of the output in the domain is changed, we may not be able to find the location of certain value 33 in some realizations. Therefore we should first normalize the outputs into the same scale, and then perform interpolation. Define a relative value after normalization as min max min min max (, ; ) (; ) (, ; ) (; ) ( ; ) (; ) min ( , ; ) (; ) max ( , ; ) x x x x x x r D D st s t st st s t st s t st s t θ θ θ θ θ θθ θθ ∈ ∈ − ≡ − = = (4.1) where the scaling factors s min and s max are the minimum and maximum value over the space domain D, respectively, which will be analyzed together with the location and used for generating samples in the following steps. In a word, the location x of any given response value s is now represented by a level set where the normalized model response s r (as a function of x) takes on the constant s at a certain time t as { } (, ; ) ( , ; ) xxx r st s t s θθ ≡ = (4.2) 4.1.3 Polynomial Construction In the xTPCM, the location in equation (4.2) as an alternative variable (instead of the model response in the PCM) is approximated by the polynomial. Similarly, there are two major approaches: the orthogonal basis approach and the Lagrange interpolation approach. In the former orthogonal basis approach, assume the location is represented by the PCE similar to equation (3.3) as a random process 1 11 12 1 2 112 01 2 111 ˆ(, ; ) (, ) ( , ) ( ( )) (, ) ( ( ), ( )) ... x ββ β i ii ii i i iii st st st st θξθ ξθξθ ∞∞ === =+ Γ + Γ + ∑∑∑ (4.3) and truncated by a finite number like equation (3.4) as 34 1 ˆ(, ; ( )) (, ) ( ( )) x ξ b ξ Q ii i st st θ θ = =Ψ ∑ (4.4) The PCE coefficients ( , ), 1,2,..., b i st i Q = are vectors since the location x is a vector in space domain (e.g., x = (x, y, z) in 3D space). According to equation (3.13), we set ˆ(, ; ) (, ; ) x ξ x ξ jj st st = and obtain 1 ( , ) ( ) ( , ; ), 1, 2,..., b ξ x ξ Q iij j i st st j P = Ψ= = ∑ (4.5) The coefficients ( , ) b i st can also be solved through pseudo-spectral projection as 1 22 (, ; ) ( ) (, ; ) () () (, ; ), ( ) (, ) (), ( ) x ξξ x ξξ ξξ x ξξ bx ξξ Q ij j i j i i ii ii st w st d st t ρ = Ω Ψ Ψ Ψ == ≈ ΨΨ ΨΨ ∑ ∫ (4.6) or weighted least-squares regression ( ) 1 Β ZWZ WZX TT − = (4.7) where 12 1 2 ( ), ( , ,..., ) , ( , ,..., ) Z ξΒ bb b X x x x TT ji i j Q P =Ψ = = . Similar to the location, the scaling factors s min and s max in equation (4.1) are approximated by the PCE as well min min, max max, 11 ˆˆ (; ( )) ( ) ( ( )), (; ( )) ( ) ( ( )) ξξ ξ ξ QQ ii ii ii st b t s t b t θ θθ θ == =Ψ = Ψ ∑∑ (4.8) By setting min min ˆ (; ) ( ; ) ξ ξ jj st st = and max max ˆ (; ) ( ; ) ξ ξ jj st st = , the PCE coefficients for scaling factors min, () i bt and max, () i bt can be solved through pseudo-spectral projection as 35 min max 11 min, max, 22 (; ) ( ) ( ; ) ( ) () , ( ) ξξ ξξ QQ ji j j ji j j jj ii ii st w s t w bt b t == ΨΨ ≈≈ ΨΨ ∑∑ (4.9) or weighted least-squares regression () ( ) 11 min min max max , b Z WZ WZ s b Z WZ WZ s TT TT −− == (4.10) While in the Lagrange interpolation approach, the location and scaling factors are approximated by Lagrange polynomials similar to equation (3.5) as 1 min min max max 11 ˆ(, ; ) (, ; ) () ˆˆ (; ) ( ; ) ( ), ( ; ) (; ) ( ) x ξ x ξξ ξξξ ξ ξ ξ P ii i PP ii ii ii st st L st st L s t s t L = == = == ∑ ∑∑ (4.11) In the orthogonal basis approach, after substituting equation (4.6) or (4.7) into equation (4.4) and generating a sufficient number of realizations, we may obtain the location samples ˆ(, ; ) x st θ . Similarly, we can get the scaling factor samples min ˆ (; ) st θ and max ˆ (; ) st θ from combining equation (4.8) and equation (4.9) or (4.10). While in the Lagrange interpolation approach, the samples of location and scaling factors can be generated through equation (4.11). 4.1.4 Backward Transform from Location to Response The location samples ˆ(, , ) x st θ may be transformed back to obtain the response samples over the domain ( , ; ) x st θ , by finding the output values in certain places (e.g., interpolation) as 36 { } ˆ (, , ) ( , , ) xx x st s st θθ ≡ = (4.12) Recall that the outputs are normalized, thus they have to be restored to the original state according to equation (4.1) as [ ] min max min ˆˆ ˆ ˆ ( ,;) (; ) ( ,; ) (;) (; ) xx s t st s t s t st θ θθ θ θ =+ ⋅ − (4.13) These responses then serve as the samples for statistical analysis of moments and PDFs. 4.1.5 Post-processing In the xTPCM, we may estimate the statistical moments similar to equation (3.19) in the PCM, but from the transformed samples in equation (4.13). It is further clarified that it may be impractical to calculate response moments without sampling process in the xTPCM, since we would obtain the location moments instead of the model response moments from equation (3.18). 4.1.6 Framework of xTPCM The five major steps involved in the xTPCM are: (1) pre-processing, i.e., representing the model input in terms of independent random parameters, selecting the collocation points and performing the deterministic model runs at these points (section 4.1.1); (2) transforming the model response to location of particular response value (along with the scaling factors, similarly hereinafter) (section 4.1.2); 37 (3) constructing polynomial of the location, i.e., approximating the location by orthogonal basis function or Lagrange interpolating polynomial, and generating location samples (section 4.1.3); (4) transforming the location samples back to the response samples (section 4.1.4); (5) post-processing, i.e., estimating the statistical properties of the model response from the samples (section 4.1.5). It is seen that the first and last steps are exactly alike in the PCM and xTPCM. While the differences are: the xTPCM has two more steps for transforms between the response and the location, and the target of polynomial construction is changed from the response to the location (Figure 4.1). Figure 4.1: Illustration of the difference between the traditional PCM and the xTPCM. 4.2 Case Studies 4.2.1 Simple Illustrative Example The effect of strong nonlinearity is first illustrated by a simple example. Consider a 1D piston-like displacement of oil by water through a homogeneous porous medium of 38 length L = 1 [L], with the viscosity μ w = μ o = 1 [ML -1 T -1 ], and the pressure difference between the entrance and exit is held constant as Δp = 1 [ML -1 T -2 ]. The model parameter log permeability follows a Gaussian distribution as ln(k) = 0.3 θ [L 2 ], where θ ~ N(0,1) a standard normal distribution. Thus we can obtain the shock front velocity v = k Δp / ( μL) [LT -1 ], and the shock front location x = vt [L]. The water saturation S w at time t = 0.4 [T] is observed as the model response s. Figure 4.2a shows five saturation profiles with discontinuous shock front in the space domain from different random parameters. The discontinuity in space domain therefore causes the discontinuity in random parametric space (Figure 4.2b), if we analyze the saturation S w at the specific location x = 0.5 [L]. This discontinuity thus results in a strong non-Gaussianity for the output, even with a Gaussian distributed random input. As we can see in Figure 4.2c that the histogram of the saturation indicates a bimodal distribution (ensemble size = 10000). 39 Figure 4.2: Illustration of nonlinear/discontinuous effect: (a) saturation as a function of location with different random parameters; (b) saturation as a function of parameter at a fixed location (solid line for MC reference, dash-dotted line for 2 nd order polynomial, dashed line for 4 th order polynomial, similarly in Figure 4.2e); (c) histogram of saturation at a fixed location; (d) location as a function of saturation with different random parameters; (e) location as a function of parameter for a particular saturation; (f) histogram of location for a particular saturation. The log permeability ln(k) = 0.3 θ, where θ ~ N(0,1). In this case, the direct approximation of the model response by polynomial in traditional PCM not only requires high order polynomial construction that increases the computational effort, but also results in oscillation that does not disappear as the order increases. In the PCM, based on the Gauss-Hermite quadrature, the collocation points are set to be the roots of Hermite polynomial. And for 2 nd and 4 th order polynomial, we need three and five collocation points, respectively. Figure 4.2b shows the true relation in solid line is approximated by the 2 nd order polynomial in dash-dotted line (through three collocation points (3,0, 3) − marked by triangles) and by the 4 th order polynomial in 40 dashed line (through five collocation points ( 5 10 , 5 10 ,0, 5 10 , 5 10 ) −+ −− − + marked by squares). It is clearly seen that the approximated results exceed the physical constraints of 0 and 1 by the definition of saturation. However, if we transform the saturation into the location of particular saturation value (Figure 4.2d), we could gain a relatively linear relation between the location and the random parameter (Figure 4.2e). It is seen that the 2 nd order polynomial approximation is already close enough to the reference, while the 4 th order polynomial almost overlaps the exact function. Figure 4.2f demonstrates the histogram of the location, which reveals a log normal distribution that is much easier to approximate than a bimodal distribution in Figure 4.2c. In summary, to overcome the difficulty in advection dominated problem, we use the location of particular response value as a new variable, since it increases the linearity of the relation between the input and the output. 4.2.2 One-dimensional Homogeneous Case 4.2.2.1 Multiphase Flow Test Consider an oil/water two-phase flow problem involving water flooding from one end of a 1D homogeneous aquifer with random porosity. The initial reservoir is saturated with oil. The absolute permeability is constant and the Corey-type relative permeabilities are quadratic functions of water saturation as 22 ,(1 ) rw w ro w kSk S ==− . The viscosity ratio is assumed to be 1. The injected volume is set to be 0.05 of the total volume. The gravity, 41 capillary pressure and liquid compressibility are neglected. In this standard Buckley-Leverett problem, let the log porosity be Gaussian distributed with a mean of 0.2 and variance of 0.1, as ln( ϕ) = 0.2 + 0.316 θ, where θ ~ N(0,1). The saturation as model response is solved analytically by the method of characteristics in the MC method, the PCM and the xTPCM. To obtain a converged MC reference, 10000 samples of the parameter θ are randomly generated to get the saturation realizations (Figure 4.3a). While in the collocation methods, the KLE is not needed since there is only one random parameter in this homogeneous case. The collocation points are selected as the roots of Hermite polynomial according to Gaussian quadrature. Figures 4.3b and 4.3c show the corresponding realizations, where the saturation profile consists of a discontinuous shock and a continuous rarefaction. 42 Figure 4.3: Saturation profiles: directly generated from deterministic model runs by (a) MC method, (b) 2 nd order PCM, and (c) 4 th order PCM; and randomly generated from polynomial approximation by (d) 2 nd order xTPCM, (e) 2 nd order PCM, and (f) 4 th order PCM. In the MC approach, statistical properties including moments and PDFs are calculated directly from the 10000 realizations. While in the PCM, we construct the polynomial from the realizations in Figures 4.3b or 4.3c to obtain the relation between the saturation and random parameter. These results are shown in Figure 4.4 at various places. It is seen in Figure 4.4a that the true relation (the MC reference) between the output saturation and the random parameter is strongly nonlinear at x = 0.2. While the 2 nd order PCM is actually trying to approximate this discontinuous function with a 2 nd order polynomial by fitting or interpolating the data at the three collocation points (3,0, 3) − . Therefore, it results in Gibbs oscillation due to the strong nonlinearity/discontinuity. And the 4 th order PCM is trying to approximate this function with a 4 th order polynomial through the data 43 at the five collocation points ( 5 10 , 5 10 ,0, 5 10 , 5 10 ) −+ −− − + . Unfortunately, the result does not show much improvement as the number of collocation points increases. From another point of view, Figure 4.4d shows that the true saturation PDF from the MC method is bimodal, but the 2 nd or 4 th order PCM fails to capture this non-Gaussian phenomenon by the polynomial of Gaussian variable. Similar phenomena can be observed at other places (x = 0.3 in Figures 4.4b and 4.4e and x = 0.4 in Figures 4.4c and 4.4f) with different relations and PDFs. Figure 4.4: Saturation as function of random parameter at locations: (a) x = 0.2, (b) x = 0.3, and (c) x = 0.4; and corresponding PDF of saturation at these locations: (d) x = 0.2, (e) x = 0.3, and (f) x = 0.4. The log porosity ln( ϕ) = 0.2 + 0.316 θ, where θ ~ N(0,1). Since the function in Figure 4.4a is discontinuous at only one place in random space, it is natural to think that changing the collocation points (e.g., adding more points near the 44 discontinuity, or moving the points next to the discontinuity) would help resolving the problem. However, it may not be practical in three aspects: (1) the position of discontinuity is not known beforehand and may be difficult to detect; (2) the position of discontinuity varies in different spatial locations as shown in Figures 4.4a, 4.4b and 4.4c, thus it is almost impossible to capture all discontinuities from one set of collocation points; (3) the oscillation occurs not only in discontinuous problems, but also in problems with steep profile or large gradient, but it is not clear how to change the collocation points in these cases. Similarly, the idea that decomposing the entire random space into a few smooth subdomains (e.g., on either side of the discontinuity) also suffers from such aspects. It is seen that the direct approximation of the saturation by polynomial is not encouraged from the above tests. Instead, we may approximate the location of particular saturation value, which can be obtained from the realizations in Figure 4.3 by linear interpolation. In this multiphase flow test, the boundaries of saturation are fixed as 0 and 1, thus the normalization process is not required. Figure 4.5 conveys the relations between the locations and the random parameters, and corresponding location PDFs. It is seen in Figure 4.5a that the location of particular response S w = 0.2 is much more linear to the random parameter than the response itself in Figure 4.4a. Therefore, the polynomial approximation from the 2 nd order xTPCM with three collocation points almost overlaps the MC reference. Consequently, the corresponding location PDF as shown in Figure 4.5d is closer to Gaussian and much easier to match than the response PDF in Figure 4.4d. And similar results can be observed in Figure 4.5 for S w = 0.5 and 0.8 as well. To insure 45 accuracy in forward transform, one hundred saturation levels uniformly distributed from 0 to 1 are used, while ten levels are enough for acceptable results based on our tests. Since the 2 nd order xTPCM needs only three realizations from the model runs, which are the same with the 2 nd order PCM (i.e., Figure 4.3b), the computational effort is greatly reduced compared to the MC method. The 4 th order xTPCM is not shown here as the 2 nd order xTPCM already performs well enough. Figure 4.5: Location as function of random parameter at particular saturation values: (a) S w = 0.2; (b) S w = 0.5; and (c) S w = 0.8; and corresponding location PDF at these saturations: (d) S w = 0.2; (e) S w = 0.5; and (f) S w = 0.8. The log porosity ln( ϕ) = 0.2 + 0.316 θ, where θ ~ N(0,1). In the PCM, a sufficient number (we use 10000 here) of saturation realizations can be randomly generated from the constructed polynomial approximation. This sampling process is actually very fast without deterministic model runs (with only addition, 46 subtraction and multiplication, but without division or matrix inversion), albeit a relatively large number. As for the xTPCM, the same number of location samples are generated and transformed back to obtain the saturation samples through linear interpolation. Figure 4.3e shows 50 realizations generated from the 2 nd order PCM by the same parameters as used in Figure 4.3a. It is seen that the results are unphysical in three aspects: (1) oscillatory saturation value outside the boundary of 0 and 1; (2) lower value at the upstream and higher value at the downstream; (3) zero output variability at about x > 0.55. We here explain the possible reasons one by one: for the first unphysical aspect, it is clear that the oscillation is caused by polynomial fitting/interpolation near discontinuity as shown in Figure 4.4; for the second aspect, it is because the saturation is analyzed individually at each location in the PCM, thus there is no guarantee that the saturation follows a decreasing trend from downstream to upstream; for the third aspect, the reason that the saturations in three realizations are all 0 at about x > 0.55 if we recall Figure 4.3b. While in the 4 th order PCM, although a higher order polynomial and more number of collocation points are used, unphysical problems are still observed in Figure 4.3f similar to Figure 4.3e. On the contrary, 50 saturation realizations from the 2 nd order xTPCM are much more reasonable as shown in Figure 4.3d. What’s more, generated by the same parameters, the 50 realizations from the xTPCM (Figure 4.3d) and those from the MC method (Figure 4.3a) are exactly alike, indicating a nice approximation of the response surface. To further illustrate the validity of the xTPCM, we also investigate the relation between the saturation and the parameter as well as the saturation PDFs as shown in Figure 4.4, which are obtained from the transformed saturation samples. Compared to 47 the result from the PCM, the result from the xTPCM almost overlaps with the MC reference, even for the complex discontinuous relation and bimodal PDF. Since the model uncertainties are usually quantified by the statistical moments, we compare up to fourth order moments of model response due to non-Gaussianity. As illustrated in Figure 4.6, the results by the PCM may deviate from those by the MC method, with weird discontinuous curves. What’s more, the moments stay at zero at places x > 0.55 for the 2 nd order PCM, because of zero output variabilities at these places as shown in Figure 4.3b. Such problems are also observed for the 4 th order PCM. On the other hand, the 2 nd order xTPCM clearly outperforms both the 2 nd and 4 th order PCM. Notice that the xTPCM with only 3 realizations generates almost the same results with the MC reference from 10000 realizations. For better comparison, we also investigated the effect of ensemble size in the MC method, and observed that the moments are not accurate enough until the ensemble size reaches about 1000 (not shown here due to limited space). 48 Figure 4.6: Statistical moments of saturation by the MC method and collocation methods: (a) mean; (b) variance; (c) third order moment; and (d) fourth order moment. Two additional cases with different log porosity variability are implemented, where the means and variances of the saturation are compared. When the variability is small ( 2 ln φ σ = 0.01), the shock front locations from different random parameters are close to each other as shown in Figures 4.7a and 4.7d. While the variability is large ( 2 ln φ σ = 1), a wide zone is affected by the shock front, as illustrated in Figures 4.7c and 4.7f. It can be concluded that the transform approach is more important with larger input variability. 49 Figure 4.7: Statistical moments of saturation with different input variabilities, 2 ln φ σ = 0.01: (a) mean and (d) variance; 2 ln φ σ = 0.1: (b) mean and (e) variance; 2 ln φ σ = 1: (c) mean and (f) variance. 4.2.2.2 Solute Transport Test In this section, we will apply the collocation methods to a solute transport model and analyze the concentration uncertainties. In the above examples, the responses are monotonic functions of locations, thus there is only one possible location for particular response value. To verify the xTPCM in examples with non-monotonic functions, we design a 1D solute transport test with instantaneous injection at the origin initially as 2 2 12 0 ,, (0) , ( ) , ( ,0) ( ), ( , ) 0, 0 LLL cc c Kh Dv D vv tx x x hhhL hcx c xc t t α φ δ ∂∂ ∂ ∂ =− = = ∂∂ ∂ ∂ == = ±∞=> (4.14) 50 where () x δ is the Dirac delta function. And the analytical solution is given as [Bear, 1974] 2 () 4 0 (, ) 2 L x vt D t L c cx t e Dt φπ − − = (4.15) The domain is 30 m long, where the two boundaries are imposed with constant heads (8 m and 5 m, respectively). The porosity is assumed constant as 0.1, the dispersivity α L is 0.1 m, the input log conductivity (in m/day) follows a normal distribution with a mean of 0 and variance of 0.1, as ln(K) = 0.316 θ, where θ ~ N(0,1). The concentration at time t = 10 days is investigated as model response. Similar to the saturation tests, 10000 realizations are carried out in the MC method as benchmarks for stochastic analysis. Figure 4.8 shows the concentration profile from direct model simulation. Compared to Figure 4.3, there is no discontinuous shock wave but continuous steep front. Notice that the maximum concentration in the domain varies among different realizations. To insure that there are two locations for any given concentration value, each concentration profile is normalized by its maximum value c max (c min is zero as a constant).Therefore the profile is rescaled in the range [0, 1] as a relative concentration c r in equation (4.1). Now for a certain concentration level, there are two locations, which can be approximated by the polynomial individually as x 1 and x 2 . Figure 4.9 shows the approximations of c max , x 1 and x 2 by the xTPCM in forms of max (; ) ct θ , 1 (, ; ) x ct θ and 2 (, ; ) x ct θ as equations (4.4) and (4.8). Here, the illustrated x 1 and x 2 are related to the relative concentration level c = 0.5, while other levels can be treated similarly. It is seen that these transformed variables are relatively linear to the random parameter and the corresponding PDFs are close to 51 Gaussian distribution. Therefore, the exact results from the MC approach are easily approximated by the xTPCM after transform, just like in Figure 4.5. Figure 4.10 shows the relations between the concentration and the random parameter at various locations and the corresponding PDFs, which is similar to Figure 4.4. The randomly generated realizations from the polynomial approximation by the PCM and the xTPCM are illustrated in Figures 4.8d, 4.8e and 4.8f, with similar observations from Figure 4.3. The unphysical problems are also observed for the PCM, while the xTPCM conveys reasonable results resembling the MC method. Figure 4.8: Concentration profiles: directly generated from deterministic model runs by (a) MC method, (b) 2 nd order PCM, and (c) 4 th order PCM; and randomly generated from polynomial approximation by: (d) 2 nd order xTPCM, (e) 2 nd order PCM, and (f) 4 th order PCM. 52 Figure 4.9: Transformed variables as function of random parameter: (a) scaling factor c max , (b) first location x 1 , and (c) second location x 2 ; and corresponding PDFs: (d) scaling factor c max , (e) first location x 1 , and (f) second location x 2 . The log conductivity ln(K) = 0.316 θ, where θ ~ N(0,1). 53 Figure 4.10: Concentration as function of random parameter at locations: (a) x = 7.5, (b) x = 10, and (c) x = 12.5; and corresponding PDFs of concentration at these locations: (d) x = 7.5, (e) x = 10, and (f) x = 12.5. The log conductivity ln(K) = 0.316 θ, where θ ~ N(0,1). Two additional cases with different dispersivities are implemented to test the effect of dispersion from the statistical moments. When the dispersivity is small ( α L = 0.01 m), the problem is advection dominated and the transform approach is a must, as shown in Figures 4.11a and 4.11d. While the dispersivity is large ( α L = 1 m), the concentration profile becomes relatively smooth, and the results by the PCM are somehow acceptable (still not as accurate as the xTPCM), as illustrated in Figures 4.11c and 4.11f. In general, the xTPCM clearly demonstrates its advantage over the PCM. 54 Figure 4.11: Moments of concentration with different dispersivities, α L = 0.01 m: (a) mean; (b) variance; and α L = 1 m: (c) mean; (d) variance. 4.2.3 One-dimensional Heterogeneous Case Let the model input be a spatially correlated random field in this heterogeneous case. Consider a water flooding problem with a 1D domain of 5000×100×100 ft (1 ft ≈ 0.305 m) divided into 100×1×1 grid blocks. The initial reservoir pressure is 4000 psi (1 psi ≈ 6895 Pa). There is one water injection well at one end with constant injection rate of 100 barrels per day, and one production well at the other end with constant pressure of 3000 psi. The absolute permeability is assumed constant as 20 mD (1 mD ≈ 9.87E −16 m 2 ) since saturation is not sensitive to it in this boundary condition. The irreducible water saturation (also the initial condition) S wi = 0.2, and the residual oil saturation S or = 0.3. 55 The relative permeability and capillary pressure follows the Brooks-Corey type model [Brooks and Corey, 1964] 22 1 31 2 ,(1 )(1 ), , 1 wwi rw e ro e e c b e e or wi SS kS k S S p pS S SS λλ λ ++ − − ==− − = = −− (4.16) where p b is bubbling pressure and set as 0 here for simplicity, λ is pore-size distribution index and set to be 2. Assume the log porosity to be a second-order stationary Gaussian random field with a mean of 0.1, and covariance 2 ln ln 1 2 exp( / ) Cxx φφ σ η =−− , where 2 ln 0.5 φ σ = and /0.2 L η = , where L is the domain length. The deterministic forward model run is implemented by the Eclipse simulator numerically. In the MC method, 1000 direct simulations are carried out. While in the collocation method, the log porosity random field is expressed by the KLE as equation (3.2) with truncated dimension N. The collocation points are selected based on level one Smolyak sparse grid. And the one-dimensional nodal set 1 i Θ is defined as the roots of Hermite polynomial with order 2i − 1 (e.g., 1 1 {0} Θ= , 2 1 {3,0, 3} Θ= − ) according to Gaussian quadrature. The Gauss-Hermite abscissas is said to be weakly nested (with 0 repeated), thus 12 11 Θ ⊂Θ and subsequently (, ) ( 1, ) NN N N ⊂+ HH , although in general 1 11 kk + Θ⊄Θ . Therefore, the number of collocation points is P = 2N + 1 as 1 23 1 23 21 [0,0,...,0] [0,0,..., 3] , [0,..., 3,0] ,..., [ 3,0,...,0] [0,0,..., 3] , [0,..., 3,0] ,..., [ 3,0,...,0] ξ ξξ ξ ξξ ξ T TT T N TT T NN N + ++ + = == = =− = − =− (4.17) 56 The polynomial approximating the model response is constructed with Lagrange interpolation approach (the orthogonal polynomial can be considered as well). Figure 4.12 shows that as N increases, the xTPCM will be more accurate. If N = 10 dimensions of KLE with about 94% energy 11 / N Ni i ii E λ λ ∞ == = ∑ ∑ are retained, the xTPCM with P = 21 times of model runs provide almost exactly the same results with the MC method from 1000 times of simulation. However, the traditional PCM even diverges with higher retained random dimension N and larger number of collocation points P. It is possibly because the polynomial fitting/interpolation process is more challenging for high dimensional problem with strong nonlinearity/discontinuity. To better compare the performs of these two collocation methods with respect to the random dimension, we define a root mean square error (RMSE) corresponding to an L2 error norm as 2 1 1 () b N estimation reference ii i b RMSE y y N = =− ∑ (4.18) where y denotes the target variable (i.e., mean saturation or saturation variance in this case), i is the grid block index, N b is the total number of blocks. In common senses, the result of collocation methods is expected to converge to the reference as the retained random dimension increases [Chang and Zhang, 2009]. It is seen in Figure 4.13 that the error from the xTPCM keeps decreasing, but the error from the PCM does not decline but in the opposite inclines. 57 Figure 4.12: Statistical moments of saturation with different number of random dimensions N and number of collocation points P = 2N + 1. Mean saturation: (a) N = 1, P = 3, (b) N = 3, P = 7, (c) N = 10, P = 21; and saturation variance: (d) N = 1, P = 3, (e) N = 3, P = 7, (f) N = 10, P = 21. Figure 4.13: Errors of saturation moment in various number of retained random dimensions N: (a) mean error; and (b) variance error. 58 4.2.4 Two-dimensional Homogeneous Case Consider a homogeneous domain of 30×3 m. The left and right sides are assigned as constant heads of 8 m and 5 m, respectively, while the top and bottom boundaries are impervious. If the instantaneous point source is released at the origin initially, the governing equation will be slightly different from equation (4.14) as 22 22 0 (, ,0) ( ) ( ) LT cc c c DD v tx y x cx y c x y δδ ∂∂ ∂ ∂ =+ − ∂∂ ∂ ∂ = (4.19) And the analytical solution is given as [Bear, 1974] 22 () 44 0 (, ) 4 LT xvt y DtDt LT c cx t e DD t πφ ⎡ ⎤ − −+ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ = (4.20) The conductivity is assumed constant as 1 m/day, and the random log porosity is Gaussian distributed with a mean of 0.1 and variance of 0.1. The longitudinal dispersivity α L is 0.1 m, and the transverse dispersivity α T is 0.01 m. Concentrations at time t = 10 days are compared between the MC approach and collocation methods. In the MC method, 1000 deterministic model runs are implemented. While in the PCM and xTPCM, since there is only one random parameter in this homogeneous case, the collocation points are chosen as the roots of Hermite polynomial similar to section 4.2.2. It is observed from equation (4.20) that for any concentration level, the locations satisfy the formula of an ellipse as 22 () 44 LT xvt y const Dt D t − += (4.21) 59 Thus, the process of transform from the data to the ellipse parameter is like this: (1) normalize the concentration by the maximum concentration c max as equation (4.1); (2) for given relative concentration value, find the locations by linear interpolation in Eulerian coordinates (here we use the function “contourc” in MATLAB); (3) solve for the ellipse parameter including the center positions x c , y c , semimajor axis r a and semiminor axis r b from the above locations by curve fitting. These parameters are then represented by the polynomial as max (; ) ct θ , (,;) c x ct θ , (,;) c y ct θ , (,;), a rct θ (,;) b rct θ whose coefficients are solved by the collocation method. After randomly generating a sufficient number of the ellipse parameters, these multiple levels of ellipses are then transformed back to obtain the concentrations in Eulerian coordinates for each realization by scattered data interpolation (here we use the function “scatteredInterpolant” in MATLAB). Finally, the statistical properties are estimated from these concentration samples. Figure 4.14 shows the mean and variance from these methods. It is seen that the PCM under or over estimates the moments, while the xTPCM shows a significant improvement and matches the solution from the MC simulation quiet well. The PDF results also confirm this conclusion, which are not shown here due to limited space. As for 3D homogeneous cases, similar approach can be carried out through replacing the ellipse by ellipsoid. 60 Figure 4.14: Statistical moments including the mean saturation from: (a) MC method; (c) 2 nd order PCM; (e) 4 th order PCM; and (g) 2 nd order xTPCM; and the saturation variance from: (b) MC method; (d) 2 nd order PCM; (f) 4 th order PCM; and (h) 2 nd order xTPCM. 61 4.3 Discussion In this research, the xTPCM is applied to 1D cases and a 2D homogeneous case. As for general 2D or 3D heterogeneous cases, how to parameterize the non-unique locations of a certain output level is an interesting topic. One solution is to utilize digital image processing techniques through motion analysis. The idea is to treat the data of each realization in the whole domain at a certain time as an image, which is analyzed as a function of random parameters in a stochastic domain, similar to the sequential image interpolation in a temporal domain. Other solutions, including the level set method with interface tracking (segmentation, description, and evolution), may also be tested. One possible limitation of the xTPCM is the “out of bounds” problem, where the location of a particular response value is outside of the spatial domain and, thus, difficult to obtain. One solution is expanding the domain to a larger area. Another practical solution is to use the “arrival time” of a particular response value instead of the “location”. That is, we transform the response to a time domain instead of a space domain, and run the model for a longer time if necessary. Another possible problem of the xTPCM is that the backward transform process may become difficult if the location is a non-monotonic function of the response, which is similar to the contradiction in the “method of characteristics”. Such an example is illustrated by the solid line in Figure 6.1, where the location x = x(s) is a non-monotonic function of the model response s. In this case, if we want to find the relation s = s(x), it becomes multi-valued resulting in a non-physical solution. However, this problem is 62 two-fold. First of all, it is unlikely to happen. Since the approximated realization x = x(s, θ*) is generated through polynomial construction from the other realizations x = x(s, θ j ), j = 1,2,...,P by deterministic model runs (which are physical), the x = x(s, θ*) is unlikely to be non-physical. At least, we didn’t observe such phenomenon in our tests. And even it happens, it could probably be eliminated by a higher order polynomial approximation with higher accuracy. Second, if the non-monotonic phenomenon is observed, which is similar to the contradiction in “method of characteristics”, it could be removed by the formation of a shock wave discontinuity (dashed line) in 1D cases. While in 2D or 3D cases, it still needs further investigation. These extensions and solutions are being investigated and will be discussed in the future. Figure 4.15: Illustration of the difficulty for the xTPCM when the location is a non-monotonic function of the response. In some previous tests, the PCM is successfully applied to geophysical models, including subsurface flow and transport, without reporting non-physical oscillation. One possible reason for this is that the model output is not sensitive to the random input in 63 certain boundary conditions. In general, if the hydraulic head/pressure is fixed on the boundary, the saturation or concentration is more sensitive to permeability than porosity; on the other hand, if the injection/production rate is fixed, the saturation or concentration is more sensitive to porosity than permeability. Therefore, under the condition of a constant injection rate and constant porosity, even the input permeability varies, the model output (e.g., saturation [Li et al., 2009] or concentration [Fajraoui et al., 2011]) does not change much; thus, no obvious oscillation was observed. Other possible reasons include relative small input variability [Li and Zhang, 2009] and weak discontinuity with large dispersivity [Shi et al., 2010]. It draws our attention that Gu and Oliver [2006] presented a similar idea of replacing the state variable by shock front location, which was applied to a 1D example by ensemble Kalman filter for data assimilation. We, however, introduced a transform process into the collocation method for uncertainty quantification. Moreover, we used the locations of multiple levels, instead of the shock front location alone, to gain robustness under all conditions from strong to weak discontinuity. In addition, we proposed a normalization step before interpolation for varying output boundaries, and also discussed possible extensions to 2D or 3D problem. It is also worth noting that our method is different from the combined methods of the polynomial chaos approach and Lagrangian particle simulation. Le Maître and Knio [2007] presented a scheme based on the incorporation of polynomial chaos into a particle approximation in Lagrangian coordinates, by using a unique set of particles to transport the stochastic modes of the solution. Muller et al. [2011] presented a so-called 64 “polynomial chaos/Monte Carlo method” for dispersion-free problems, by first analyzing the hydraulic head by the standard collocation method, and then applying the Monte Carlo method to the streamline technique for the transport problem. Those methods require either Lagrangian particle solvers or streamline simulators; our method, however, does not. The reason for this is that we used the location of particular model response values, instead of the position of particles. Thus, we utilized the advantage of representation in Lagrangian coordinates over that in Eulerian coordinates for advection-dominated problems through a transform approach. 65 Chapter 5 Displacement-based Transformed Probabilistic Collocation Method (dTPCM) As discussed in Chapter 1, the direct approximation of model response by the polynomial is not suggested due to strong nonlinearity. A possible solution is to find an alternative variable to represent the model response, which is relatively linear to the random parameter. The alternative variable could then be analyzed by traditional collocation method, in the form of a polynomial as a surrogate. To obtain the response statistics, a sufficient number of the samples of the alternative variable are generated through either projection or regression or interpolation, and then transformed back to obtain the response samples, from which the moments and distributions can be estimated readily. As we can see, the key is to find a suitable alternative variable which should be not only easily represented by the random parameter, but also relatively linear to the parameter, and it has to be conveniently transformed from and to the model response (i.e., extract the variable from the response, and reconstruct the response from the variable). Intuitively inspired by the Lagrangian point of view, the location of particular response value is a good candidate for transform from the Eulerian coordinates, especially for the advection dominated problems. However, the extension of this “location” concept to 2D or 3D general problems naturally leads to iso-line (contour) or iso-surface segmentation, description and evolution. Although the iso-line/iso-surface can be analyzed by some 66 methods (e.g., level set method), how to combine them with the collocation methods for uncertainty quantification is not straight forward. We find that another variable similar to location, the “displacement” vector, serves as a better option for 2D and 3D general problems, and thus will be adopted in this work. In this chapter, we propose a new displacement-based transformed probabilistic collocation method (dTPCM) to quantify the uncertainty accurately and efficiently in 2D and 3D general cases under strongly nonlinear conditions. In the dTPCM, the model response in space at a certain time is considered as an image, and the differences of the responses (with respect to different random parameters) are considered as shape changes between images. To obtain the displacement field from the response data, we introduce digital image processing techniques, specifically motion analysis. Therefore, the response field is now represented by a displacement field. Similar to the location in the xTPCM, the displacement in the dTPCM is relatively linear to the input parameter; thus, the transform process improves the performance of the collocation method. Furthermore, the displacement serves as a better alternative variable than the location, as the former can be easily extended to 2D and 3D general problems. The displacement is approximated by polynomial construction on collocation points, and the polynomial is then used to generate the displacement on any node in random parametric domain. Finally, the generated displacement is transformed back through scattered data interpolation to obtain the response. We conduct multiphase flow and solute transport tests to show the validity of the dTPCM through higher accuracy than the PCM, and less computational effort compared to the Monte Carlo (MC) method. 67 This chapter presents the work in paper which is under review by Liao and Zhang [2014a]. 5.1 Theoretical Framework The probabilistic collocation method with transform by displacement (dTPCM) can also be derived from the weighted residual method as equations (3.10)-(3.13). 5.1.1 Pre-processing In the KLE based dTPCM, the model input m is expanded by equation (3.2). In order to determine the model output statistics, we choose P sets of collocation points , 1,2,..., ξ j j P = for the random process of input, where 12 ( , ,..., ) ξ T N ξξ ξ = . By solving the model P times correspondingly, we get P sets of output ( , ; ), 1,2,..., x ξ j st j P = . Actually, the collocation points in the dTPCM could be exactly the same with those in the PCM. Therefore, we can improve the existing PCM results by the dTPCM without additional model runs. For better comparison, we will use the same collocation points in the PCM and the dTPCM for the case studies. 68 5.1.2 Forward Transform from Response to Displacement As discussed at the beginning of Chapter 5, the “displacement” is found to be a good candidate to represent the model response. But how to obtain the displacement from the response, that is to say, how to transform from the response to the displacement? This task can be achieved by digital image processing techniques, specifically motion analysis. The motion analysis, determining the movement of objects from a time varying image sequence, has been applied to a wide range of areas (e.g., image coding, particle image velocimetry, robot navigation, traffic monitoring, cloud and weather systems tracking, motion of biological cells) [Aggarwal and Nandhakumar, 1988; Vega-Riveros and Jabbour, 1989; Qin et al., 2007; Huang et al., 2009]. A common starting point for pixel motion analysis is to assume that pixel intensities are translated from the reference image at time t 0 (Figure 5.1a) to the deformed image at time t 1 (Figure 5.1b) as [Hild and Roux, 2006] 01 (, ) ( , ) xxd f tf t = + (5.1) where f (x,t) is the intensity as a function of location x and time t, and d is the displacement vector. It is noted that the “image” here is not limited to 2D frame (x = (x, y) T , d = (d x , d y ) T ), but also applied to 3D scene/vision (x = (x, y, z) T , d = (d x , d y , d z ) T ) as well. 69 Figure 5.1: Illustration of digital image correlation for time varying image sequence: (a) reference image; (b) deformed image with displacement vector; (c) normalized cross-correlation coefficient. The displacement vector is determined through searching for the maximum of the correlation coefficient. Among various approaches for motion analysis, the most widely used two are correlation-based approach and gradient-based approach [Barron et al., 1992; Jahne, 2005]. The former correlation-based approach, also referred to as digital image correlation (DIC), involves considering a small template window in the reference image (Figure 5.1a), and searching for the displacement which produces the “best match” among all possible target windows in the deformed image (Figure 5.1b). The target window is usually the same size as the template window, moving in the deformed image within a certain search range, based on the prior motion estimation. To describe the similarity of the template and target windows, a normalized cross-correlation coefficient is defined as [Lewis, 1995] 00 1 1 01 22 00 1 1 (,)(,) ( ,)( ,) (; , , ) (, ) ( , ) ( , ) ( , ) D DD ft f t f t f t d tt f tf t d f t f t d γ ∈ ∈∈ ⎡⎤⎡ ⎤ −+−+ ⎣⎦⎣ ⎦ = ⎡⎤ ⎡ ⎤ −+−+ ⎣⎦ ⎣ ⎦ ∫ ∫∫ x xx x x xd xd x dx xx x xd xd x (5.2) where D is the window size, and f denotes the mean value of f within this window. 70 Therefore, the displacement vector d at any given location x for time t 1 is determined by searching for the maximum of the correlation coefficient γ (Figure 5.1c). The displacement resolution of the DIC processing could be extended to the sub-pixel by interpolating the correlation map to increase accuracy [Schreier et al., 2000]. In this study, we use the MATLAB function “normxcorr2” in image processing toolbox to calculate the correlation coefficient γ in 2D space, and use the cubic spline interpolation for sub-pixel estimation with about 0.1 pixel precision. Another class of motion analysis is gradient-based approach. By assuming that the deformed image is well approximated by a first-order Taylor series, we may obtain a so-called gradient constraint equation, similar to the continuity equation that states the mass conservation in fluid mechanics [Jahne, 2005]. Usually, the gradient-based approach through directly solving a linear equation system is faster than the correlation-based approach by brute-force searching process. However, the gradient-based approach is more suitable for small displacement problems (best for <1 pixel), although some remedial techniques (e.g. coarse-to-fine strategy [Battiti et al., 1991]) can somehow relieve this constraint. In this study, we stick to the correlation-based approach for simplicity, while the idea of transform can be extended to other motion analysis approaches. After presenting the concept of motion analysis, we now demonstrate how to introduce it into the dTPCM. Motion analysis in digital image processing usually refers to determining the movement of an object in a time-varying image sequence (i.e., in temporal domain); whereas, motion analysis for stochastic analysis refers to determining 71 the movement of an object in parameter-varying images (i.e., in a random parametric domain). The latter is somewhat more complicated, since temporal domain is a one-dimensional domain, while parametric domain could be a multi-dimensional domain. In this work, the undeformed image, which can be recorded as 0 (, ; ) x ξ st , is set to be the model response (e.g., saturation or concentration) when the random parameter is located at the mean of the collocation points (e.g., 0 (0,0,...,0) ξ = for a multi-dimensional normal distribution with symmetric collocation point set). The deformed image, which can be recorded as ( , ; ), 1,2,..., x ξ j st j P = , is the model response when the random parameter is located at each collocation point. Note that we need an extra model run at 0 ξ , if it is not included in the collocation point set. Specifically, assume that pixel intensities are translated from the reference image (model response from random parameter 0 ξ ) (Figure 5.2a) to the deformed image (model response from random parameter , 1,2,..., j j P = ξ ) (Figure 5.2b) as: 0 (, , ) ( , , ) j ft f t = + x ξ xd ξ (5.3) To describe the similarity of the template and target windows, a normalized cross-correlation coefficient can be defined as: 00 0 22 00 (, , ) (, , ) ( , , ) ( , , ) (; , , , ) (, , ) (, , ) ( , , ) ( , , ) jj D j jj DD ft f t f t f t d t ftft d f t f t d γ ∈ ∈∈ ⎡⎤⎡ ⎤ −+−+ ⎣⎦⎣ ⎦ = ⎡⎤⎡ ⎤ −+−+ ⎣⎦⎣ ⎦ ∫ ∫∫ x xx x ξ x ξ xd ξ xd ξ x dx ξξ x ξ x ξ xxd ξ xd ξ x (5.4) 72 Therefore, the displacement vector d with respect to j ξ at any given location x and time t is determined by searching for the maximum of the correlation coefficient γ (Figure 5.2c). Figure 5.2: Illustration of digital image correlation for stochastic analysis: (a) reference image; (b) deformed image with displacement vector; (c) normalized cross-correlation coefficient. The displacement vector is determined through searching for the maximum of the correlation coefficient. From the undeformed image 0 (, ; ) x ξ st and deformed image (, ; ) x ξ j st , we can obtain the displacement vector field, which is recorded as ( , ; ) dx ξ j t . For example, if the template location in the undeformed image is x 0 , the “best match” target location in the deformed image will be x 0 + d. It is worth mentioning that the response value 0 (,;) xd ξ j st + at the target location in the deformed image may not be exactly the same as the response value 00 (,; ) x ξ st at the template location in the undeformed image. Thus, we need to obtain 0 (,;) xd ξ j st + , which can be recorded as 0 (,; ) x ξ dj st . The target location value 0 (,; ) x ξ dj st is acquired via gridded interpolation (we utilize linear interpolation in this work, i.e., bilinear/trilinear interpolation in the two- and 73 three-dimensional cases). The computational effort for motion analysis depends on the number of template locations, as we have to perform the DIC process once for each location. In this study, we choose the template locations to be the Eulerian grid points as in the model setup to keep the interpolation error in transform as the same level as the discretization error in simulation. For example, if a 2D space domain has 10 by 10 gridblocks, the template locations are also these gridblocks. The number of template locations can be further reduced if we only perform motion analysis on some selected sparse points. For example, we may set the template locations only on the gridblocks with even number index, which results in a 5 by 5 template locations in the previous example. The effect of the DIC process is mainly affected by two factors: the window size and the search range. We first set up a large enough search range and determine an appropriate window size by a sensitivity analysis. If the window size is too small, the template window cannot capture the feature around the template location, and the displacement will be affected by disturbance. On the contrary, if the window size is too large, the template window will contain locations far from the template location, and these locations are actually not moving similarly to the template location. This will make the displacement field too smooth, and also increase computational cost. The goal is to choose the best window size in this way: if the window size is appropriate, the displacement field should be accurate, and we should be able to reproduce the deformed image from the undeformed image and the displacement. In other words, the forward transform and backward transform should be a set of invertible transforms. Next, we determine the search range. In general, if the search range is smaller than the true 74 displacement, the estimated displacement will be bounded in the range and, thus, inaccurate. In contrast, if the search range is larger than the true displacement, the estimated displacement will be correct, but it will require additional time. Therefore, we should use the smallest search range that is just large enough. Essentially, the model response (,; ) x ξ j st is now represented by alternative variables, including the displacement 0 (,; ) dx ξ j t and the corresponding target location value 0 (,; ) x ξ dj st . 5.1.3 Construct Polynomial In the dTPCM, the displacement as an alternative variable (instead of the model response in the PCM) is approximated by the polynomial. Similarly, there are two major approaches: the orthogonal basis approach and the Lagrange interpolation approach. In the former orthogonal basis approach, assume the displacement is represented by the PCE similar to equation (3.3) as a random process 1 11 12 1 2 112 000 01 02 111 ˆ (,; ) (,) ( ,) (()) (,) ( (), ()) ... dx β x β x β x i ii ii i i iii tt t t θξθ ξθξθ ∞∞ === =+ Γ + Γ + ∑∑∑ (5.5) and truncated by a finite number like equation (3.4) as 00 1 ˆ (,;()) (,) (()) dx ξ bx ξ Q ii i tt θ θ = =Ψ ∑ (5.6) The PCE coefficients 0 ( , ), 1, 2,..., bx i ti Q = are vectors since the displacement d is a vector in space domain (e.g., d = (d x , d y , d z ) in 3D space). According to equation (3.13), we set 00 ˆ ( ,;) ( ,;) dx ξ dx ξ jj tt = and obtain 75 00 1 ( , ) ( ) ( , ; ), 1,2,..., bx ξ dx ξ Q iij j i ttj P = Ψ= = ∑ (5.7) The coefficients ( , ) b i st can also be solved through pseudo-spectral projection as 0 0 1 0 0 22 (,;) ( ) ( , ; ) () () (,;), () (,) (), ( ) dx ξξ dx ξξ ξξ dx ξξ bx ξξ Q ij j i j i i ii ii tw td t t ρ = Ω Ψ Ψ Ψ == ≈ ΨΨ ΨΨ ∑ ∫ (5.8) or weighted least-squares regression ( ) 1 Β ZWZ WZ D TT − = (5.9) where 12 12 ( ), ( , ,..., ) , ( , ,..., ) Z ξΒ bb b D dd d TT ji i j Q P =Ψ = = . Similar to the displacement, the target location value 0 (,; ) x ξ dj st can be approximated by the PCE as well 0,0 1 ˆ(,;()) (,) (()) x ξ x ξ Q dsii i st b t θ θ = =Ψ ∑ (5.10) By setting 00 ˆ(,; ) (,; ) x ξ x ξ dj d j st st = , the PCE coefficients for target location value can be solved through pseudo-spectral projection as 0 1 ,0 2 (,; ) ( ) ˆ (,) x ξξ x Q djijj j si i st w bt = Ψ ≈ Ψ ∑ (5.11) or weighted least-squares regression ( ) 1 b Z WZ WZ s TT s d − = (5.12) where ,1 ,2 , ,1 ,1 , ( ), ( , ,..., ) , ( , ,..., ) Z ξbs TT ji i j s s s s Q d d d d P bb b s s s =Ψ = = . While in the 76 Lagrange interpolation approach, the displacement and target location value are approximated by Lagrange polynomials similar to equation (3.5) as 00 0 0 11 ˆ ˆ (,;) (,; ) (), ( ,;) ( ,; ) () dx ξ dx ξξ x ξ x ξ ξ PP ii d d ii ii ttLst stL == == ∑∑ (5.13) In the orthogonal basis approach, after substituting equation (5.8) or (5.9) into equation (5.6) and generating a sufficient number of realizations, we may obtain the displacement samples 0 ˆ (,; ) dx t θ . Similarly, we can get the target location value samples 0 ˆ(,; ) x d st θ from combining equation (5.10) and equation (5.11) or (5.12). While in the Lagrange interpolation approach, the samples of displacement and target location value can be generated through equation (5.13). 5.1.4 Backward Transform from Displacement to Response After generating the displacement samples 0 ˆ (,; ) dx t θ and target location value samples 0 ˆ(,; ) x d st θ as discussed above, we readily obtained the response samples in a scatter plot 00 ˆ ˆ (,;) (,;) xd x d st s t θ θ += based on the definitions. Then we need to transform the scattered data 0 ˆ(,;) xd st θ + back to get the final response samples ˆ(, ; ) x st θ in original meshed grid. In this work, we use the “scatteredInterpolant” class in MATLAB for scattered data interpolation. The algorithm involves constructing an interpolating surface (Figure 5.3) by triangulating the points (A'-B'-C') and lifting the vertices by a magnitude into the 77 response space (A-B-C) orthogonal to the parameter space. And we choose the default linear interpolation method via calculating the weighted sum of values on the three vertices of the enclosing triangle (D' →D), while other interpolation methods such as nearest-neighbor could be applied as well. Though the illustration highlights 2-D interpolation, we can extend this technique to higher dimensions (actually we only need to handle up to 3-D problems in space). Figure 5.3: Illustration of scattered data interpolation by triangulation. 5.1.5 Post-processing In the dTPCM, we may estimate the statistical moments similar to equation (3.19) in the PCM, but from the back transformed samples obtained in section 5.1.4. It is further clarified that it may be impractical to calculate response moments without sampling process in the dTPCM, since we would obtain the displacement moments instead of the model response moments from equation (3.18). 78 5.1.6 Framework of dTPCM The five major steps involved in the dTPCM are: (1) pre-processing, i.e., representing the model input in terms of independent random parameters, selecting the collocation points and performing the deterministic model runs at these points (section 5.1.1); (2) transforming the model response to displacement (along with the target location value, similarly hereinafter) (section 5.1.2); (3) constructing polynomial of the displacement, i.e., approximating the displacement by orthogonal basis function or Lagrange interpolating polynomial, and generating displacement samples (section 5.1.3); (4) transforming the displacement samples back to the response samples (section 5.1.4); (5) post-processing, i.e., estimating the statistical properties of the model response from the samples (section 5.1.5). It is seen that the first and last steps are exactly alike in the PCM and dTPCM. While the differences are: the dTPCM has two more steps for transforms between the response and the displacement, and the target of polynomial construction is changed from the response to the displacement (Figure 5.4). 79 Figure 5.4: Illustration of the difference between the traditional PCM and the dTPCM. 5.2 Case Studies 5.2.1 One-dimensional Example with Single Random Variable Although the forward and backward transforms in the dTPCM are illustrated in 2D space as in Figures 5.2 and 5.3, we start with a simple 1D example to demonstrate the implementation details. Let an instantaneous point source released at the origin initially. Assume that the log conductivity θ = ln(k) is normally distributed, with a mean of 0 and a standard deviation of 0.3. Figure 5.5a shows the concentration profiles with different collocation points (e.g., θ = −1.732, 0 and 1.732 for 2 nd order approximation according to Gaussian quadrature rule for standard normal distribution). Figure 5.5b shows the true point concentration at x = 0.3 as a function of random parameter θ in black solid line. This exact response is approximated by the 2 nd , 4 th , 6 th and 8 th order polynomials. We can see that the result is not much improved as the polynomial order increases, due to the low regularity and unsmoothness. In consequence, the estimated PDFs are also deviated from the exact one as shown in Figure 5.5c. Similar phenomena are observed in multiphase 80 flow case, where the discontinuity in the space domain (Figure 5.5d) is translated to the discontinuity in the random parametric domain (Figure 5.5e), and the latter results in a bimodal distribution (Figure 5.5f). Figure 5.5: 1D solute transport and multiphase flow examples: (a) concentration profiles with different collocation points; (b) different orders of polynomial approximation to the concentration response at a fixed location; (c) PDF estimation of the concentration response at a fixed location; (d) saturation profiles with different collocation points; (e) different orders of polynomial approximation to the saturation response at a fixed location; (f) PDF estimation of the saturation response at a fixed location. To address the problem caused by the traditional PCM, we apply the proposed dTPCM in this example. First, we transform concentration to displacement. We consider the concentration profile with θ = 0 as the undeformed image (Figure 5.6a), and consider the 81 concentration profiles with other collocation points (e.g., θ = 1.732, Figure 5.6b) as the deformed images. For each template location (e.g., x 0 = 0.3), we perform motion analysis to obtain the displacement d(x 0 ; θ) by direct searching for the maximum of correlation coefficient γ(d) (Figure 5.6c). Meanwhile, the target location concentration c d (x 0 ; θ) with respect to the displacement d(x 0 ; θ) is recorded along with the displacement (Figure 5.6b). Figure 5.6: Illustration of forward transform: (a) concentration profile as undeformed image; (b) concentration profile as deformed image; (c) determining the displacement by correlation coefficient. Now the concentration is represented by the displacement (Figure 5.7a) and the target location value (Figure 5.7d). Next, we use polynomials to approximate the displacement (Figure 5.7b) and the target location value (Figure 5.7e). We can see that the 2 nd order polynomial approximation is already very accurate, while the 4 th order polynomial almost overlaps the true function. Therefore, the estimated PDFs are very close to the exact ones (Figures 5.7c and 5.7f). Compared to Figure 5.5, we can infer that the dTPCM improves the linearity of the input-output relation, hence the model response can now be approximated accurately using relatively low order polynomials to reduce computational effort. The idea of the dTPCM can be understood from another perspective: the model 82 response as one variable in the PCM (Figrue 5.4a) is now separated into displacement and target location value as two variables in the dTPCM (Figures 5.7a and 5.7d), taking advantage that these two variables are smoother than the original variable. Figure 5.7: Illustration of polynomial constructions of the displacement and the target location value. Since we have approximated the displacement and the target location value by polynomials, we can generate them with regard to any random parameter through polynomial interpolation (e.g., θ = 2, demonstrated by the blue solid lines in Figures 5.8a and 5.8b). Therefore, the concentration profile can be obtained from the target location value plus a shifting process related to the displacement, as depicted by the red dashed line in Figure 5.8b. Note that this concentration profile is defined in scattered grid (e.g., the template location x 0 = 0.3 plus the displacement d = 0.2464 results in a target location 83 x 1 = x 0 + d = 0.3 + 0.2464 = 0.5464 where the target location concentration c d = 1.392), we need to apply the scattered data interpolation (linear interpolation in this simple 1D case) to obtain the response in meshed grid, as revealed by the red solid line in Figure 5.8c. For any random parameter, we carry out the above backward transform process and generate multiple realizations of the model response (Figures 5.8d, 5.8e and 5.8f). Figure 5.8: Illustration of backward transform: (a) and (d) displacement from polynomial interpolation; (b) and (e) target location value from polynomial interpolation; (c) and (f) response in meshed grid obtained from that in scattered grid. (a), (b) and (c) are for one realization; (d), (e) and (f) are for multiple realizations. To test the performance of the dTPCM, we randomly generate 20 realizations from the MC method and the dTPCM using the same input parameters. Figure 5.9a shows that the realizations are identical between these two methods, indicating a good match. If we look into the point concentration as a function of random parameter, we can see that the dTPCM captures both the parameter-response relation (Figure 5.9b) and the response 84 PDF (Figure 5.9c) very well. Compared to Figure 5.5, we can also see that the PCM failed to approximate the MC results with even the 8 th order polynomial, on the contrary, the dTPCM estimate the results accurately with only the 2 nd order polynomial. Figure 5.9: Final results by the dTPCM: (a) 20 realizations randomly generated by the MC method and the dTPCM using the same parameters; (b) true concentration response approximated by the dTPCM; (c) exact PDF approximated by the dTPCM. 5.2.2 Two-dimensional Example with Single Random Variable The implementation of dTPCM is, perhaps, better explained by a solute transport example in a 2D homogeneous porous medium. If an instantaneous point source is released at the origin initially, the analytical solution is [Bear, 1974] 22 () 44 0 (, ) 4 LT xvt y DtDt LT c cx t e DD t πφ ⎡ ⎤ − −+ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ = (5.14) Assume that the log conductivity θ = ln(k) is normally distributed as the random variable, with a mean of 0 and a variance of 0.2. Figure 5.10a shows the response concentration as the undeformed image when the random parameter θ = 0. Figure 5.10b depicts one of the deformed images with regard to θ = 1.732. If we aim to determine the displacement at the 85 template location (x 0 = 0.3, y 0 = 0.05), we can calculate the normalized cross-correlation coefficient, as shown in Figure 5.10c. The correlation coefficient reaches its maximum when the displacement is (d x = 0.35, d y = 0.03); thus, the target location in the deformed image is (x 1 = x 0 + d x = 0.65, y 1 = y 0 + d y = 0.08). After we determine the displacement 00 0 0.3, 0.05, 1.732 (,; ) dx xy t θ θ == = , we perform data interpolation to obtain the response value at the target location 00 0 0.3, 0.05, 1.732 (,; ) x d xy st θ θ == = (i.e., the concentration at (x, y) = (0.65, 0.08) in Figure 10b). In commonly used Cartesian coordinates, the process is quite easy, and we use the linear interpolation here (i.e., bilinear/trilinear for 2D/3D case). Thus, the response (, ; ) x st θ is now represented by alternative variables, including the displacement 0 (,; ) dx t θ and the corresponding target location value 0 (,; ) x d st θ . Figure 5.10: Forward transform from response to displacement: (a) undeformed response image; (b) deformed response image; and (c) normalized cross-correlation coefficient. The displacement is determined when the correlation coefficient reaches its maximum. Next, we construct polynomial approximations for the displacement and the target location value. If we attend to the x-direction displacement 0 (,; ) x x dt θ as a function of parameter θ at the template location (x 0 = 0.3, y 0 = 0.05), we may find that d x is relatively 86 linear to θ (Figure 5.11a). Therefore, a second order polynomial is adequate by fitting/interpolating the three collocation points (i.e., θ = −1.732, 0, and 1.732), which are marked as circles. Similarly, the y-direction displacement 0 (,; ) x y dt θ and the target location value 0 (,; ) x d st θ can be well approximated by second order polynomials (Figures 5.11b and 5.11c). On the contrary, if we attend to the response ( , ; ) x st θ at the same location (x = 0.3, y = 0.05) in the traditional PCM, the relation between the parameter and the response is strongly nonlinear and thus cannot be properly captured by a second order polynomial (Figure 5.11d). These phenomena are similar to those observed in the previous 1D example (Figures 5.5 and 5.7). Figure 5.11: Construct polynomials of the displacement, the target location value and the response: (a) x-direction displacement as a function of parameter in the dTPCM (the solid line indicates the MC reference, the dashed line indicates the second order polynomial approximation; the triangles indicate the collocation points); (b) y-direction displacement as a function of parameter in the dTPCM; (c) target location value as a function of parameter in the dTPCM; and (d) fixed location response as a function of parameter in the PCM for comparison. 87 Then, we perform backward transform from the displacement to the response. The displacement field 0 (,; ) dx t θ and the corresponding target location value 0 (,; ) x d st θ , which are generated from polynomial interpolation with regard to parameter θ = 1, are shown in Figures 5.12a (top view) and 5.12d (side view). Therefore, we obtain the response value in scattered plot as 00 (,;) (,;) d st s t θ θ + = xd x (Figures 5.12b and 5.12e). Note that to perform data interpolation, we need to define the boundary values to avoid extrapolation. In this case, the boundary values are assigned to 0, assuming the domain is large enough and no solute reaches the boundary. Through scattered data interpolation, the final response ( , ; ) x st θ (in Cartesian coordinates/meshed grid) is presented in Figures 5.12c and 5.12f. Finally, we estimate the response statistics from the interpolated response samples. 88 Figure 5.12: Backward transform from displacement to response: (a) and (d) displacement field 0 (,; ) dx t θ and target location value 0 (,; ) x d st θ generated by polynomial approximation; (b) and (e) scattered data 00 ( ,;) ( ,;) xx d sdt s t θ θ + = from definition; (c) and (f) gridded data ( , ; ) x st θ interpolated from 0 (,;) x sdt θ + . 5.2.3 Two-dimensional Example with Multiple Random Variables In this example, we still analyze the solute transport model defined in equation (5.14), but consider four random parameters (i.e., conductivity k, porosity ϕ, longitudinal dispersivity α L and hydraulic head on boundary h). The domain size is 30m×10m, and each gridblock is 0.25m×0.25m. The top and bottom boundaries are impervious, while the left and right boundaries have constant head h 1 and h 2 , respectively. Set the transverse dispersivity α T = 0.5 α L , and h 2 = 0. The four random parameters are assumed to be 89 uniformly distributed independently (k ∈[2, 3], ϕ ∈[0.2, 0.3], α L ∈[0.1, 0.15], and h 1 ∈[4, 6]). Initially, the solute is released at cell (x = 5 m, y = 5 m) with unit concentration. The concentration is observed at time t = 5 day. We randomly generate 5,000 realizations in the MC method as benchmark, and use Smolyak sparse grid algorithm with Clenshaw-Curtis quadrature rule in the PCM and the dTPCM. Figure 5.13 shows the mean and variance of concentration. The first row shows the MC results. The second row shows the results by collocation methods with level 1 sparse grid (i.e., nine collocation points for four random variables). It can be seen that the PCM results have multiple peak values, which is inaccurate. On the other hand, the dTPCM results are much closer to the MC results. As the sparse grid level increases, the PCM results become better, but still different from the benchmark; whereas the dTPCM provides results almost identical to the benchmark. Figure 5.14 illustrates the mismatch of the data from the collocation methods and the MC simulation. As the level of sparse grid grows, the results are improved in both collocation methods. However, the concentration variance by level 3 PCM still has obvious deviation from the MC benchmark (Figure 5.14j). In contrast, the high R 2 values in the dTPCM indicate the good agreement with the benchmark results. 90 Figure 5.13: Mean concentration <c> and concentration variance σ c 2 by the MC method, as well as the PCM and dTPCM with sparse grid level 1, 2 and 3. 91 Figure 5.14: Mean concentration <c> and concentration variance σ c 2 compared between the MC method and the PCM/dTPCM with sparse grid level 1, 2 and 3. 5.2.4 Two-dimensional Example with Random Field Although the 1D and 2D examples with homogeneous media in previous sections 3.2 to 3.4 can be analyzed by the xTPCM, we still use these simple examples to show the implementation details and effects of the dTPCM clearly. In this section, we test the PCM and the dTPCM in an oil/water two-phase flow model with a heterogeneous medium. Consider a five-spot pattern consisting of 51×51×1 grid blocks, each measuring 92 40×40×40 ft 3 (1 ft = 0.305 m). There is one water injection well in the center with a constant injection rate of 1,000 barrels per day (1 barrel = 0.159 m 3 ) and four production wells in the corners with constant pressures of 2,000 psi (1 psi = 6895 Pa), as shown in Figure 5.15. The initial reservoir pressure is 4,000 psi, and the porosity is assumed to be constant at 0.2. The residual oil saturation is S or = 0.3, and the irreducible water saturation, as well as the initial condition, is S wi = 0.2. The relative permeability and the capillary pressure follow the Brooks-Corey type model as: 22 1 31 2 ,(1 )(1 ), , 1 wwi rw e ro e e c b e e or wi SS kS k S S p pS S SS λλ λ ++ − − ==− − = = −− (5.15) where p b = 100 psi is the bubbling pressure; and λ = 2.0 is the pore-size distribution index. The log permeability (in mD) is a second-order stationary Gaussian random field with a mean of 4.0, and a covariance of 2 ln ln 1 2 1 2 exp( ) kk x y Cxxyy σ ηη =−− −− , where ln 2.0 k σ = and the correlation lengths are 800 ft xy η η = = . The first N = 20 random dimensions in the KLE are retained. In the MC method, a relatively large ensemble size of 5,000 is used to reduce the sampling error. Since the traditional Gaussian quadrature nodes are unnested, we adopted the nested quadrature rules in Genz and Keister [1996]. Therefore, 41 and 801 collocation points are selected in the level 2 and 3 sparse grid algorithm, respectively. The level 4 sparse grid (over 10,000 collocation points are required for 20 random dimensions in this case) is not tested, since it requires more forward model runs than the MC method. This also infers the difficulty of using high order polynomial approximation in high dimensional problems. The saturation at 500 days is observed as model response. 93 Figure 5.15: Schematic diagram of model setup in the multiphase flow test. Lines A 1 -A 2 , B 1 -B 2 , and C 1 -C 2 indicate where the mean and variance are compared in detail. Points D 1 , D 2 and D 3 indicate where the PDFs are compared in detail. Figure 5.16 shows the motion analysis, in which we use the function “normxcorr2” in the MATLAB image processing toolbox to calculate the correlation coefficient γ, and the spline interpolation to estimate the sub-pixel displacement. When the input parameter is located at the mean of the collocation points (recorded as 0 ξ ), the log permeability field is homogeneous (ln(k) = 4.0 everywhere), and the corresponding output saturation as the undeformed image is symmetrical with respect to the lines x = 1,000 ft and y = 1,000 ft (Figure 5.16e). When the parameter is located at another collocation point (recorded as 3 ξ ), the log permeability value is higher on the right side (Figure 5.16a). Therefore, the corresponding output saturation as a deformed image shifts to the right (Figure 5.16c). From the above undeformed and deformed images, we can obtain the displacement 94 vector field (Figure 5.16e) using digital image correlation techniques. For another collocation point 5 ξ , the log permeability value is higher in the upper-left and lower-right quadrants (Figure 5.16b). Consequently, the plume moves faster in these two directions (Figure 5.16d), and the displacement vectors clearly show the characteristics (Figure 5.16f). 95 Figure 5.16: Illustration of motion analysis at 500 days: (a) log-permeability with parameter 3 ξ ; (b) log-permeability with parameter 5 ξ ; (c) water saturation with parameter 3 ξ as a deformed image; (d) water saturation with parameter 5 ξ as a deformed image; (e) water saturation with parameter 0 ξ as the undeformed image, and the displacement between 0 () w S ξ and 3 () w S ξ ; and (f) water saturation with parameter 0 ξ as the undeformed image, and the displacement between 0 () w S ξ and 5 () w S ξ . 96 5.2.4.1 Effect of Interpolation Level The first row in Figure 5.17 shows the mean saturation. We can see that the result from the level 2 dTPCM (Figure 5.17d) is very close to that from the MC method (Figure 5.17a) and is smoother than that from the level 2 PCM (Figure 5.17b). It can be seen that the PCM overestimates or underestimates the values. The second row in Figure 5.17 shows saturation variance with similar observations, where overestimation by the PCM is more obvious. Figure 5.18 depicts the saturation moments along the lines A 1 -A 2 , B 1 -B 2 , and C 1 -C 2 . Being consistent with the contours in Figure 5.17, the PCM results depart from the MC results; whereas, the dTPCM results match the MC results well. Although a higher level sparse grid (level 3 with 801 collocation points) is used in the PCM, the results are not improved, but become even worse in some places (larger deviation for mean and overestimation for variance). This is caused by interpolating strongly nonlinear/discontinuous functions in high dimensional problem. We will elaborate it in the discussion section. Figure 5.19 illustrates the saturation PDFs at three locations along the line A 1 -A 2 , which clearly explain the mismatch of the PCM and the MC method in the statistical moments. There are obvious non-physical values in the PDF in Figures 5.19a and 5.19b, since the traditional PCM does not consider physical bounds in polynomial approximation. We point out that the nontrivial value for level 3 PCM at S w = −1 in Figure 5.19a indicates that there are some estimated saturation less than −1. Similar phenomenon is observed at S w = 2 in Figure 5.19b. These are consistent with the results in Figures 5.18a and 5.18d. Note that in Figure 5.19c, the level 2 and 3 PCMs capture only one mode (S w = 0.7), but fail to capture the other mode (S w = 0.2). In general, the 97 PCM renders better results in location D 3 than in locations D 1 and D 2 . This proves that the PCM error occurs mainly near the flooding shock front, where the nonlinearity/discontinuity is strongest. On the contrary, the dTPCM captures the complex character of the bimodal distribution well in all positions, due to the linearity/smoothness of input parameter and displacement. Importantly, although the dTPCM runs the model only 41 times, it still produces results very close to those from the MC method with 5,000 model runs. Figure 5.17: Moments of water saturation by different methods: (a) and (b) for the MC method; (c) and (d) for the PCM; and (e) and (f) for the dTPCM. (a), (c), and (e) are mean saturations; and (b), (d), and (f) are saturation variances. 98 Figure 5.18: Moments of water saturation along lines A 1 -A 2 (x + y = 2000 ft), B 1 -B 2 (x + y = 1600 ft), and C 1 -C 2 (x + y = 1200 ft). (a), (b), and (c) represent mean; and (d), (e), and (f) represent variance. Figure 5.19: PDFs of water saturation at different locations: (a) D 1 (x = 600 ft, y = 1400 ft); (b) D 2 (x = 640 ft, y = 1360 ft); and (c) D 3 (x = 800 ft, y = 1200 ft). 99 5.2.4.2 Comparison of Computational Effort In the MC method, 19,140 seconds are required for 5,000 model runs. In the PCM method with level 2 sparse grid, 143 seconds are used for 41 model runs. In the dTPCM methods with level 2 sparse grid, additional 45 seconds are needed for forward transform and 36 seconds for backward transform; thus the total time is 224 seconds. We can see that the dTPCM costs about 50% more time than the PCM, but greatly improves the accuracy. In the PCM method with level 3 sparse grid, 3778 seconds are used for 801 model runs, but the results are still inaccurate. Furthermore, in the model runs, we use the Eclipse commercial software, which is based on FORTRAN and very fast. However, in the forward and backward transforms, we use MATLAB, which is an interpreted program and relatively slow. Therefore, the efficiency of the dTPCM can be further improved by rewriting the code in a compiled program (e.g., C or FORTRAN). 5.2.4.3 Effect of Simulation Time To compare the accuracy of the collocation methods in different simulation times, we run the model up to 1,000 days. In the dTPCM, since more water has been injected as time advances, the water zone is enlarged. In some realizations, the water zone has reached the domain boundary. In this case, the topology of the saturation image may be affected by the boundary. Although we can still apply motion analysis and obtain the displacement field, the displacement vector is bounded by the domain size. Therefore, the linearity/regularity of the parameter-displacement relation (as illustrated in Figures 5.7 and 5.11) becomes weak, and the accuracy may be reduced. 100 Figure 5.20 shows the saturation moments along line A 1 -A 2 at 250, 500, 750 and 1000 days. We can see that as the time advances, both the PCM and the dTPCM become less accurate, as the influence of the input random permeability field on the output saturation response is increased. However, the dTPCM still provides much better estimations than does the PCM. If the simulation time is further increased (e.g., up to 2000 days), the water breakthrough will occur in the production wells for most realizations. In that case, the topology of the saturation image changes dramatically, and the dTPCM will be no longer applicable. Another drawback of the dTPCM is that for each time step, the forward and backward transforms have to be implemented once to obtain the statistical properties. Therefore, if we aim for a model response in the whole space at a few time steps, the dTPCM will be suitable. However, if we want to obtain the response information in some locations at many time steps, the dTPCM may be inefficient. In such cases, arrival time instead of displacement could serve as another variable for transforms. 101 Figure 5.20: Moments of water saturation along line A 1 -A 2 (x + y = 2000 ft) at four simulation times: (a), (b), (c) and (d) represent mean; and (e), (f), (g) and (h) represent variance. 5.2.5 Three-dimensional Example with Random Field In this section, we will test the collocation methods in solute transport models and analyze the concentration uncertainties. Assume that the domain is divided into 121 × 41 × 21 grid blocks, each measuring 0.25 × 0.25 × 0.25 m (Figure 5.21). The flow and advection-dispersion equations are subject to the following boundary conditions: 102 0 at 0, 10, 0, 5, 10 at 0, and 5 at 30, 0 at all boundaries h yy z z n hx h x C n ∂ = == = = ∂ == = = ∂ = ∂ (5.16) An instantaneous source with unit concentration (1 g/L) is released in the cell (x = 5 m, y = 5 m, and z = 2.5 m) at the initial time. The porosity is assumed to be constant at 0.2. The log conductivity (in m/day) is represented by a second-order stationary Gaussian random field with a mean of 0, and a covariance of 2 ln ln 1 2 1 2 1 2 exp( / / / ) KKx y z Cxxyyzz σ ηη η =−− −− −− , where ln 0.5 K σ = and the correlation lengths are 5 xy z η ηη == = m. The longitudinal dispersivity L α is set at 0.1 m, and the transverse dispersivity T α is 0.1 of L α ; thus, the Peclet number Pe = L/ L α is O(10 2 ). In this case, the first N = 30 random dimensions in the KLE are retained. According to the level 2 Smolyak sparse grid algorithm, a total of 61 collocation points are selected for both the PCM and the dTPCM similar to those in section 3.5. The MC method with 5,000 realizations is carried out as a benchmark. The concentration at 10 days is observed and analyzed. 103 Figure 5.21: Schematic diagram of model setup in the solute transport test. Lines A 1 -A 2 , B 1 -B 2 , and C 1 -C 2 indicate where the simulation results are compared in detail. Figure 5.22 shows the moments of concentration compared among the MC method, the PCM, and the dTPCM. It can be seen that the mean concentration given by the MC method (Figure 5.22a) is very close to that given by the dTPCM (Figure 5.22e), but different from that given by the PCM (Figures 5.22c). Similar phenomena are observed for concentration variance (Figures 5.22b, 5.22d, and 5.22f). Figure 5.23 shows the concentration moments along the lines A 1 -A 2 , B 1 -B 2 , and C 1 -C 2 . It is clear that the dTPCM results match the MC results well; whereas, the PCM results deviate from both. More specifically, the PCM overestimates or underestimates the mean in various places, and overpredicts the variance for almost the entire domain. The negative values in the mean estimation by the PCM demonstrates the difficulty of approximating strongly nonlinear relations in high dimensional problems. On the other hand, the dTPCM using 104 61 collocation points produces results resembling those from the MC method using 5,000 samples, indicating a great improvement in efficiency. Figure 5.22: Moments of concentration using different methods: (a) and (b) for the MC method; (c) and (d) for the PCM; and (e) and (f) for the dTPCM. (a), (c), and (e) are mean concentrations; and (b), (d), and (f) are concentration variances. 105 Figure 5.23: Moments of concentration along the lines A 1 -A 2 (x, y = 5 m, z = 2.5 m), B 1 -B 2 (y, x = 12 m, z = 2.5 m), and C 1 -C 2 (y, x = 14 m, z = 2.5 m). (a), (b), and (c) represent for mean concentrations; and (d), (e), and (f) for represent concentration variances. In the MC method, the computational effort is 216914 seconds for 5,000 model runs. In the PCM, it is 1506 seconds for 61 model runs. In the dTPCM, additional 3342 and 2250 seconds are required for forward and backward transforms, respectively; thus the total is 7098 seconds. We can see that the computational cost in the dTPCM is approximately 5 times as much as that in the PCM. This ratio is larger than that in the previous 2D example (which is around 1.5), due to the complexity of transforms in 3D space. This indicates that the efficiency of dTPCM for extremely complex problems in 106 3D space requires further investigation. Nevertheless, in this example, the dTPCM is still very efficient compared to the MC method, as the former is 30 times faster than the latter. We also notice that the average time for each model run in the MC method, which is 216914/5000 = 43.3s, is longer than that in the collocation methods, which is 1506/61 = 24.7s. The reason is that the input parameters in the MC method are randomly generated; whereas those in the PCM and the dTPCM are selected on the collocation points. Therefore, the permeability fields in the MC method are more heterogeneous and it takes more time for simulator to converge. 5.3 Discussion 5.3.1 Sparse Grid in High Dimensional Problems When the PCM with sparse grid algorithm yields inaccurate results, it is natural to consider using a higher level approximation. In problems with only a few random dimensions, the above strategy is usually effective. However, in high dimensional problems, it may not work or even produce worse results. Here we use a simple mathematical formula to illustrate the above problem. Let θ i be uniform random variables, the response y along the 1D domain x is: 1 1, 1 () , , [ 1,1] 0, N ii i x pos y x pos xpos N θθ = ≤ ⎧ == ∈− ⎨ > ⎩ ∑ (5.17) This model can be considered as an imitation of 1D piston-like displacement in two-phase flow system. The pos denotes the shock position, and its mean and variance 107 are constant regardless of the number of random dimensions N. The response is approximated by the PCM with sparse grid algorithm and Clenshaw-Curtis quadrature nodes. Figures 5.24a and 5.24d show that when N = 1, the accuracy is increased as the sparse grid level increases. Figures 5.24b and 5.24e depict that when N = 3, no obvious improvement is observed when the sparse grid level increases. Figures 5.24c and 5.24f illustrate that when N = 10, the results with a higher level sparse grid are even worse than a lower level approximation (larger deviation in mean and more apparent overestimation in variance). These phenomena, which are similar to those in Figures 5.18 and 5.23 in previous examples, demonstrate the difficulty of sparse grid methodology in high dimensional problems with strong nonlinearity/low regularity. Formaggia et al. [2013] also reported that the sparse grid approximation may lose accuracy as the number of dimensions increases, and the extension to high-dimensional problems is still an open research topic. 108 Figure 5.24: Mean and variance with different dimensions N: (a) and (b) N = 1; (b) and (e) N = 3; (c) and (f) N = 10. 5.3.2 Piecewise Polynomial Approximation To cope with the Gibbs oscillation in polynomial approximation under discontinuities, the piecewise reconstruction approaches are proposed via (weighted) essentially non-oscillatory (ENO or WENO) techniques (Shu, 1998). These are high order accurate finite difference schemes designed for problems with piecewise smooth solutions containing discontinuities. They are later introduced into the non-intrusive PCM for uncertainty quantification (Abgrall, 2007; Barth, 2012; Witteveen and Iaccarino, 2012). However, the application of the piecewise approaches in stochastic analysis is more 109 challenging than that in physical problems due to the curse-of-dimensionality. In the conservation law based fluid flow problems, the model parameters are in space and time, thus there are at most four parametric dimensions (3D space and 1D time). On the contrary, in the stochastic analysis, the number of random dimensions is usually from O(1) to O(100). Since the piecewise reconstruction requires a relatively large number of nodes (especially near the shock) to capture the discontinuity even in 1D domain, the extension to high dimensional problems will result in a significant increase of the collocation points along with computational burden. We think that it is the reason that the uncertainty quantification examples using piecewise approximation are mostly for low random dimensional problems in current literature. We also point out that a similar idea of changing the collocation points (e.g., adding points near the discontinuity, or moving the points next to the discontinuity) would not help resolve the problem in general. The main reason is that the position of discontinuity varies in different spatial locations; thus, it is almost impossible to capture all discontinuities from one set of collocation points (Liao and Zhang, 2013). This leads to a dense mesh for collocation points which costs more simulation time. 5.3.3 Limitation of the dTPCM and Possible Solution One limitation of the dTPCM is what we refer to as the “out of bounds” problem. For example, the displacement generated from the polynomial approximation could be related to a target location that is out of the domain boundary. This is the case when boundaries 110 refer to the outer ones, and we can solve the problem by expanding the domain. However, there is another case when boundaries refer to the inner ones (e.g., injection or production well), and we cannot solve the problem by simply enlarging the area. One practical solution is to use the arrival time instead of the location or displacement as the alternative variable. Specifically, we transform the response to a time domain instead of a space domain. Although the time-based transform may also encounter the “out of bounds” problem in the time domain, this problem can be easily addressed by running the model longer. We will present the results in the next chapter. Another concern of the dTPCM is related to the isotropy for space, time and random variable. The traditional motion analysis is applied to time-varying images sequence, based on the assumption that the images in space are similar when the times are close. In this study, we adopted this idea and apply the motion analysis to parameter-varying images set, assuming that the images in space are similar when the random parameters are close. However, further analysis including the theoretical proof is required. Moreover, the above assumption may not hold in extremely heterogeneous media or with very large input variability. In such cases, the topology of the response image would change dramatically, and thus the images are no longer similar. 111 Chapter 6 Time-based Transformed Probabilistic Collocation Method (tTPCM) In Chapters 4 and 5, we present the xTPCM and dTPCM. However, these methods may encounter the “out of bound” problem. For example, the location or displacement generated from the polynomial approximation could be related to a target position that is out of the domain (outer boundary). In this case, we can solve the problem by enlarging the area. Another example is that the target position happens to be an injection or production well (inner boundary), and we cannot relieve the obstruction by simply expanding the domain due to the possible overlap in space. In addition, the above two methods are unable to tackle the variables that are only defined in certain specific places. For instance, if we aim for the production/injection rate in source or sink wells or the saturation/concentration at fixed positions, using the location- or displacement-based transform is not straightforward. Moreover, their computational costs grow with the number of time steps to be analyzed, since these methods investigate the output statistics at each time step individually. The stochastic model in the context of subsurface hydrology describes the spreading of the solutes by means of a Lagrangian formulation in some previous works [Dagan, 1982; Shapiro and Cvetkovic, 1988; Zhang and Neuman, 1995; Xin and Zhang, 1998; Zhang and Tchelepi, 1999]. A noteworthy simplification of the flow and transport problem may 112 be achieved by projecting the equation along a streamline and embedding all the heterogeneities of the transport properties within a single variable, the travel time. One implication of the time-based Lagrangian formulation is that the 3D Eulerian equations in spatial domain have been recast in a 1D form in temporal domain, because the original three-dimensional nature of the flow and transport is entirely encapsulated in the travel time. A similar idea is used here, but in a different scheme. The traditional stochastic Lagrangian approaches track the travel time of a large amount of solute particles, whereas we pursue the arrival time of certain concentration levels via interpolation. In this work, we use the arrival time instead of the location or displacement as the alternative variable to represent the response, and propose a time-based transformed PCM (tTPCM). Specifically, we transform from an output domain to a time domain instead of a space domain. Although the time-based transform may also encounter the “out of bounds” problem with regard to time, this problem can be easily addressed by running the model longer in the 1D time domain. Similar to the location/displacement in the xTPCM/dTPCM, the arrival time in the tTPCM is relatively linear to the input parameter; thus, the transform process improves the performance of the collocation method. Through solute transport and multiphase flow tests, we show that the tTPCM has higher accuracy than the PCM and requires lower computational burden compared to the direct-sampling MC method. This chapter presents partial works in paper which is in preparation by Liao and Zhang [2014b]. 113 6.1 Theoretical Framework Some basic concepts (e.g., Karhunen-Loeve expansion (KLE), polynomial chaos expansion (PCE), Lagrange interpolation and Smolyak sparse grid algorithm) have already been introduced in Chapter 3. In general, the collocation method can be derived from the weighted residual method as in equations (3.10)-(3.13). 6.1.1 Pre-processing In the tTPCM, the input random field can be represented by the KLE. To construct a polynomial approximation for the output response, we choose P sets of collocation points , 1,2,..., ξ j j P = for the random process of input, where 12 ( , ,..., ) ξ T N ξξ ξ = . By running the forward model P times correspondingly, we get P sets of output ( , ; ), 1,2,..., x ξ j st j P = . The collocation points in the tTPCM could be exactly the same as those in the PCM, which means that we can improve the existing PCM results by the tTPCM without additional model runs. For better comparison, we will use the same collocation points in the PCM and the tTPCM for the case studies. 6.1.2 Forward Transform From Response to Time From the direct model runs, we readily have the response values as functions of location, time and parameter (on the collocation points) as ( , ; ) x ξ j st . Usually, we show 114 these functions in the relation of response vs. location at a certain time, while we can also observe these functions in the relation of response vs. time at a fixed location. We then search for the arrival time of any given response value as ( , ; ) x ξ j ts using interpolation techniques. If there are more than one times are found (e.g., in the solute transport cases with an instantaneous release initially), we can analyze these times individually. One issue that needs to be considered is that, if the minimum or maximum values of the response in the time domain are changed, we may not be able to determine the arrival time of a certain response value in some realizations. Therefore, we should first normalize the responses into the same scale, and then perform interpolation. Define a relative value after normalization as: min max min min max (, ; ) ( ; ) (, ; ) (; ) ( ; ) (; ) min ( , ; ) (; ) max ( , ; ) x x xx x xx xx xx r D D st s st ss sst sst θ θ θ θ θ θθ θθ ∈ ∈ − ≡ − = = (6.1) where the scaling factors s min and s max are the minimum and maximum value over the time domain D, respectively. These factors will be analyzed together with the arrival time and used for generating samples in the following steps. Essentially, the arrival time t of any given response value s is now represented by a level set, where the normalized model response s r (as a function of t) takes on the constant s at a certain location x as: { } (, ; ) (, ; ) xx r ts ts t s θθ ≡ = (6.2) 115 6.1.3 Polynomial Construction For the tTPCM, the arrival time as an alternative variable (instead of the model response in the PCM) is approximated by the polynomial. Similarly, there are two major approaches to construct polynomials: the orthogonal basis approach and the Lagrange interpolation approach. In the former approach, assume that the arrival time is represented by the PCE similar to equation (3.3) as a random process: 1 11 12 1 2 112 01 2 111 ˆ (, ; ) (, ) ( , ) ( ( )) (, ) ( ( ), ( )) ... xx x x i ii ii i i iii ts s s s θβ β ξθ β ξθ ξ θ ∞∞ === = +Γ + Γ + ∑∑∑ (6.3) and truncated by a finite number similar to equation (3.4) as: 1 ˆ (, ; ( )) (, ) ( ( )) x ξ x ξ Q ii i ts b s θ θ = =Ψ ∑ (6.4) According to equation (3.13), we set ˆ (, ; ) (, ; ) x ξ x ξ jj ts t s = and obtain: 1 ( , ) ( ) ( , ; ), 1, 2,..., x ξ x ξ Q iij j i bs t s j P = Ψ= = ∑ (6.5) The coefficients ( , ) x i bs can be solved through a pseudo-spectral projection as: 1 2 (, ; ) ( ) (, ; ), ( ) (, ) (), ( ) x ξξ x ξξ x ξξ Q ij j j i i ii i ts w ts bs = Ψ Ψ =≈ ΨΨ Ψ ∑ (6.6) or a weighted least-squares regression as: ( ) 1 b Z WZ WZ t TT − = (6.7) where 12 1 2 ( ), ( , ,..., ) , ( , ,..., ) Z ξbt TT ji i j Q P bb b t t t =Ψ = = . 116 Similar to the location, the scaling factors s min and s max in equation (6.1) are also approximated by the PCE as: min min, max max, 11 ˆˆ ( ; ()) ( ) ( ()), ( ; ()) ( ) ( ()) x ξ x ξ x ξ x ξ QQ ii ii ii sb s b θ θθ θ == =Ψ = Ψ ∑∑ (6.8) By setting min min ˆ (; ) ( ; ) x ξ x ξ jj ss = and max max ˆ (; ) ( ; ) x ξ x ξ jj ss = , the PCE coefficients for scaling factors min, () x i b and max, () x i b can be solved through a pseudo-spectral projection as: min max 11 min, max, 22 (; ) ( ) ( ; ) ( ) () , ( ) xξξ xξξ xx QQ ji j j ji j j jj ii ii sw s w bb == ΨΨ ≈≈ ΨΨ ∑∑ (6.9) or a weighted least-squares regression: () ( ) 11 min min max max , b Z WZ WZ s b Z WZ WZ s TT TT −− == (6.10) In the Lagrange interpolation approach, however, the location and scaling factors are approximated by Lagrange polynomials in the univariate case as: 1 min min max max 11 ˆ (, ; ) (, ; ) ( ) ˆˆ (; ) ( ; ) ( ), ( ; ) (; ) ( ) x ξ x ξξ x ξ xξξ x ξ x ξ ξ P ii i PP ii ii ii ts t s L ss Ls s L = == = == ∑ ∑∑ (6.11) As for the multivariate case, the location and scaling factors can be interpolated from tensor product construction [Babuska et al., 2007] or sparse grid construction [Xiu and Hesthaven, 2005]. In the orthogonal basis approach, after substituting equation (6.6) or (6.7) into equation (6.4) and generating a sufficient number of realizations, we may obtain the time 117 samples ˆ (, ; ) x ts θ . Similarly, we can derive the scaling factor samples min ˆ (; ) st θ and max ˆ (; ) st θ from combining equation (6.8) and equation (6.9) or (6.10); whereas in the Lagrange interpolation approach, the samples of arrival time and scaling factors can be generated through tensor product or sparse grid based on equation (6.11). 6.1.4 Backward Transform from Time to Response The arrival time samples ˆ (, ; ) x ts θ may be transformed back to obtain the response samples (, ; ) x st θ , by searching for the response value in any given time (e.g., interpolation) as: { } ˆ (, , ) (, ; ) xx st st s t θθ ≡ = (6.12) Recall that the responses are normalized; thus, they have to be restored to the original state according to equation (6.1) as: [ ] min max min ˆˆ ˆ ˆ ( ,;) ( ; ) ( ,; ) ( ;) ( ; ) xx x x x st s s t s s θ θθ θ θ =+ ⋅ − (6.13) These responses then serve as the samples for statistical analysis. 6.1.5 Post-processing In the tTPCM, we estimate the statistical moments similar to equation (3.19), but from the transformed samples in section 6.1.4. Furthermore, calculating response moments 118 without the sampling process in the tTPCM may be impractical, since we would obtain the time moments instead of the response moments. 6.1.6 Framework of tTPCM The five major steps involved in the tTPCM are: (1) Pre-processing, i.e., representing the model input in terms of independent random parameters, selecting the collocation points and performing the deterministic model runs at these points (section 6.1.1); (2) Transforming the model response to the arrival time (along with the scaling factors) (section section 6.1.2); (3) Constructing polynomials of the arrival time, i.e., approximating the time by orthogonal basis functions or Lagrange interpolations, and generating time samples (section section 6.1.3); (4) Transforming the arrival time samples back to the response samples (section section 6.1.4); (5) Post-processing, i.e., estimating the statistical properties of the model response from the samples (section section 6.1.5). It can be seen that the first and last steps are exactly alike in the PCM and tTPCM. The differences, however, are: the tTPCM has two more steps for transforms between the response and the arrival time, and the target of the polynomial construction is changed from the response to the arrival time in the middle step (Figure 6.1). 119 Figure 6.1: Illustration of the difference between the traditional PCM and the tTPCM. 6.2 Case Studies We design a series of case studies with different model complexity (1D illustrative example with a single-input random variable, 2D solute transport and 3D multiphase flow examples with multiple-input random fields). For each deterministic model run, the MT3D or Eclipse code is used to numerically solve the differential equations in solute transport or multiphase flow example, respectively. We derive the statistical moments and PDFs of model responses using the PCM and the tTPCM. For the purpose of comparison, we conduct MC simulations with a large number of realizations. 6.2.1 One-dimensional Illustrative Example In this case, we consider a 1D homogeneous water-saturated porous medium, in which a steady state flow has been established. The domain length is L = 1 [L], and the difference of the hydraulic heads between the left and right boundaries is fixed at Δh = 120 0.1 [L]. The porosity is constant at ϕ = 0.1, and the conductivity is the only random parameter with a normal distribution as ln(K) = 0.3 θ, where θ ~ N(0,1). Thus, the pore velocity / vK h L φ =Δ is also random. If the contaminant is continuously released at a point source (x = 0), the analytical solution is well-known as: 0 (, ) exp 222 L LL c x vt xv x vt c x t erfc erfc D Dt Dt ⎡⎤ ⎛⎞ ⎛⎞ ⎛⎞ −+ =+⎢⎥ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎢⎥ ⎝⎠ ⎝⎠ ⎝⎠ ⎣⎦ (6.16) where () 0 0, t cc = is the source concentration and set at 1.0 [ML -3 ], LL Dv α = is the dispersion coefficient, and L α is the dispersivity and set at 0.005 [L]. The statistical properties of the solute concentration at time t = 0.3 [T] are analyzed. The MC method with 1,000 model runs (50 of them are shown in Figure 6.2a) serves as a benchmark. In the PCM, 3 model runs are conducted on the selected collocation points (Figure 6.2b). As for the tTPCM, the same collocation points are used, but observed in the time domain (Figure 6.2c) instead of the space domain. We perform the model runs up to time t = 1.0 [T] to create a long-enough time window, then we search for the arrival time of a given concentration value in the forward transform step. Since the concentration interval is fixed as [0, 1], the normalization process is not required in this simple case. 121 Figure 6.2: Response concentration profiles: (a) generated from model runs for the MC method, (b) generated from model runs for the PCM, (c) observed from another point of view for the tTPCM, (d) generated from polynomial approximation in the PCM, and (e) generated from polynomial approximation in the tTPCM. In the PCM, we directly analyze the model response by constructing the polynomial of random parameter. Figure 6.3a shows the concentration at x = 0.3 [L], t = 0.3 [T] as a function of parameter. The reference from the MC method (solid line) is strongly nonlinear, thus it cannot be approximated accurately by the second order polynomial (dashed line) through three collocation points (circles) in the PCM. Meanwhile, as shown in Figure 6.3b, the true non-Gaussian bimodal distribution from the MC method cannot be captured by the PCM as well. As for the tTPCM, we treat the arrival time as the polynomial of random parameter. Figure 6.3c reveals the time at x = 0.3 [L], c = 0.5 [ML -3 ] as a function of parameter. It can be seen that the reference from the MC method 122 (solid line) is very smooth and thus well approximated by the PCM (dash-dotted line) through three collocation points (triangles). Therefore, the arrival time PDF from the PCM resembles that from the MC method, as illustrated in Figure 6.3d. Figure 6.3: Relations between model response, arrival time and random parameter: (a) response as a function of parameter at the fixed location x = 0.3 and time t = 0.3; (b) PDF of the response; (c) arrival time as a function of parameter at the fixed location x = 0.3 and response c = 0.5; and (d) PDF of the arrival time. In the PCM, we can randomly generate concentration realizations from the polynomial approximation. Figure 6.2d depicts 50 of these realizations, where the non-physical problems with Gibbs oscillation are observed. As to the tTPCM, we generate arrival time samples from the polynomial construction and then transform them back to obtain the concentration realizations. Figure 6.2e shows 50 of these realizations, which are 123 physically meaningful and successfully reproduce the results in Figure 6.2a. Similar conclusions can be drawn from the concentration moments in Figure 6.4, which indicates that the tTPCM greatly improves the accuracy compared to the PCM. Figure 6.4: Statistical moments of concentration by the MC method, the PCM and the tTPCM: (a) mean; (b) variance; (c) third-order moment; and (d) fourth-order moment. 6.2.2 Two-dimensional Solute Transport Test In this section, we test the collocation methods in solute transport models and analyze the concentration uncertainties. In the 2D domain, the left and right sides are assigned as constant heads (20 m and 10 m, respectively), while the top and bottom boundaries are impervious (Figure 6.5). The domain is divided into 121 × 41 grid blocks with each block being 0.25 × 0.25 m. The porosity is assumed to be constant at 0.2. The log conductivity 124 (in m/day) is represented by a second-order stationary Gaussian random field with a mean of 0, and a covariance of 2 ln 1 2 ln 1 2 1 2 ( , ) exp( / / ) xx K Kx y Cxxyy σ ηη =−− −− , where ln 0.3 K σ = and 5 xy η η == m. The longitudinal dispersivity L α is set to be 0.1 m, and the transverse dispersivity T α is 0.1 of L α , thus the Peclet number Pe = L/ L α is O(10 2 ). There is an instantaneous source with unit concentration (1 g/L) released at the cells (x = 4.875-5.125 m and y = 3.875-6.125 m) at the initial time (Figure 6.5). Figure 6.5: Schematic diagram of model setup in the 2D solute transport test. Lines A 1 -A 2 , B 1 -B 2 and C 1 -C 2 indicate where the simulation results are compared in detail. In the MC method, 5,000 realizations are implemented by direct model runs as a reference. In the collocation methods, the first N = 20 random dimensions in the KLE are retained, and a total number of 2N + 1 = 41 collocation points are selected according to the level-one Smolyak sparse grid algorithm. In the forward transform step of the tTPCM, since the maximum concentration over the time domain varies with the random parameter (Figure 6.6a), we need to normalize the concentration to obtain relative values (Figure 125 6.6b). As there are two possible arrival times for a certain relative concentration level, we can analyze these two times individually. The concentration at time t = 10 day is observed, whereas in the tTPCM, the simulation time is t = 30 day to create a time window for interpolation. Figure 6.6: Normalization process: (a) absolute (original) concentration profiles generated from model runs; and (b) relative concentration profiles after normalization. The moments of concentration are compared in Figure 6.7. It can be seen that the mean concentration given by the tTPCM is smoother than that given by the PCM, and is closer to that of the MC method. The PCM overestimates the values near the position (x = 18 m, y = 5 m), while underestimates the values near the position (x = 22 m, y = 5 m). Similar phenomena are observed for concentration variance, where the PCM overestimates the values near the position (x = 22 m, y = 5 m). To further compare these results, Figure 6.8 shows the concentration moments along the profiles A 1 -A 2 (y = 5 m), B 1 -B 2 (x = 18 m) and C 1 -C 2 (x = 22 m). Being consistent with the contours in Figure 6.7, the PCM results deviate from the MC results, while the tTPCM matches the latter well. 126 Another case with a larger input variability ( ln 0.5 K σ = ) is also implemented, as shown in Figure 6.9. We can observe that the mismatch between the PCM and the MC method is more obvious, whereas the tTPCM solution is in good agreement with the reference. Note that the tTPCM runs the model only 41 times for a stochastic problem with 20 random dimensions, that is to say, it is an extremely low order polynomial approximation. However, by transform from the response to the arrival time, the tTPCM still produces results very close to those by the MC method with 5,000 model runs, indicating a significant improvement in the efficiency. We can use a higher-order polynomial approximation or a higher-level sparse gird construction if the input variability is even larger. It is worth mentioning that the dTPCM can be carried out as well, but the space domain needs to be expanded, since the solute may flow outside the right boundary in some realizations. 127 Figure 6.7: Moments of concentration at ln 0.3 K σ = , including the mean concentration from: (a) the MC method, (c) the PCM, and (e) the tTPCM; and the concentration variance from: (b) the MC method, (d) the PCM, and (f) the tTPCM. 128 Figure 6.8: Moments of concentration at ln 0.3 K σ = along the line A 1 -A 2 (y = 5 m), B 1 -B 2 (x = 18 m) and C 1 -C 2 (x = 22 m). (a), (b) and (c) represent mean concentrations; and (d), (e) and (f) represent concentration variances. 129 Figure 6.9: Moments of concentration at ln 0.5 K σ = along the line A 1 -A 2 (y = 5 m), B 1 -B 2 (x = 18 m) and C 1 -C 2 (x = 22 m). (a), (b) and (c) represent mean concentrations; and (d), (e) and (f) represent concentration variances. 6.2.3 Three-dimensional Multiphase Flow Test In this section, we test the PCM and the tTPCM in an oil/water two-phase flow model and focus on the saturation uncertainties. Consider a quarter of nine-spot pattern consisting of 25 × 25 × 5 grid blocks, each measuring 20 × 20 × 20 ft (1 ft = 0.305 m). There is one water injection well at the bottom layer in one corner with a constant injection rate of 1,000 barrels per day (1 barrel = 0.159 m 3 ) and three production wells at the top layer in other corners with constant pressures of 1,000 psi (1 psi = 6,895 Pa), as 130 shown in Figure 6.10. The initial reservoir pressure is 4,000 psi, and the porosity is constant at 0.2. The Corey-type relative permeability is assumed quadratic as 22 ,(1 ) rw w ro w kS k S ==− . The log permeability (in mD) is a second-order stationary Gaussian random field with a mean of 4.0, and a covariance of 2 ln 1 2 ln 1 2 1 2 1 2 ( , ) exp( ) xx kk x y z Cxxyyzz σ ηη η =−− −− −− , where ln 1.0 k σ = and / 0.4, , , ii Lixyz η == , where L i is the domain length in i-direction. The first N = 30 random dimensions in the KLE are retained. According to the level-one Smolyak sparse grid algorithm, a total of 2N + 1 = 61 collocation points are selected for both the PCM and the tTPCM. In the MC method, a relatively large ensemble size of 5,000 is used to reduce the sampling error. The saturation at time t = 400 day is observed, whereas in the tTPCM, the simulation time is t = 1,000 day to create a time window for interpolation. Figure 6.10: Schematic diagram of model setup in the 3D multiphase flow test with one injection well and three production wells. Figures 6.11 shows the saturation moments in the top, middle and bottom layers. We can see that the results from the tTPCM is very close to those from the MC method and 131 are smoother than those from the PCM. Figures 12a and 12b depict the saturation moments along the line (x = y, z = 50 ft). Being consistent with Figure 6.11, the PCM results depart from the MC results, while the tTPCM results match the MC results well. It is worth mentioning that the PCM error occurs mainly near the flooding front, where the nonlinearity/discontinuity is strongest. Figure 6.12c illustrates the saturation PDF at the point (x = 310 ft, y = 310 ft , z = 50 ft) with oscillatory values, which clearly explain the mismatch of the PCM and the MC method in the statistical moments. On the contrary, the tTPCM captures the complex character of bimodal distribution well. Importantly, although the tTPCM runs the model only 61 times, it still produces results very close to those from the MC method with 5,000 model runs. 132 Figure 6.11: Moments of water saturation at the top, middle and bottom layers by the MC method, the PCM, and the tTPCM. (a), (c), and (e) describe mean saturations; and (b), (d), and (f) describe saturation variances. 133 Figure 6.12: Statistical properties of the water saturation: (a) mean saturation along the line (x = y, z = 50 ft); (b) saturation variance along the line (x = y, z = 50 ft); (c) saturation PDF at the point (x = 310 ft, y = 310 ft, z = 50 ft). Figure 6.13a illustrates 100 realizations of water cut (the ratio of water rate and total rate) at production well P1 from direct model runs using the MC method. Figure 6.13b shows 100 realizations from polynomial approximation using the PCM. It can be seen that the water cut values may go beyond the physical range of [0, 1]. Figure 6.13c demonstrates the results using the tTPCM, which are very close to the reference. The water cut PDF at time t = 400 day is shown in Figure 6.13d, which proves the above statement from another aspect of view. Figure 6.13: Water cut at production well P1: (a) 100 realizations generated from model runs in the MC method; (b) 100 realizations generated from polynomial approximation in the PCM; (c) 100 realizations generated from polynomial approximation in the tTPCM; and (d) water cut PDF at time t = 400 day. 134 6.3 Discussion In this chapter, we propose the tTPCM, which is valid in 1D to 3D general cases. The idea of transform is initially inspired from the Lagrangian approach. In general, flow quantities in fluid mechanics depend upon both position and time. There are two approaches widely used to analyze flow and transport, either by observing the fluid properties at fixed positions, which yields an Eulerian representation, or by observing the trajectories of fluid particles, which yields a Lagrangian representation. For uncertainty quantification, one usually describes the model response at fixed positions in the traditional methods. While in our new methods, we focus on the location or time of given response values in the xTPCM or the tTPCM. It is noted that we do not really track the particles, but track the location or time by interpolation process; thus our methods are not limited to the deterministic solvers that are based on Lagrangian papproach. From this point of view, our methods are similar to the method of characteristics and the level set method. The Lagrangian scheme has some drawbacks, such as clustering of a large number of particles. The semi-Lagrangian scheme overcomes this difficulty by using a regular Eulerian mesh and tracking back the origin of each particle. We adopt the above idea and propose the dTPCM, which aims to find the displacement field between the deformed and undeformed images. Hence, we successfully extend the idea of tracking locations as in the xTPCM from 1D models to 2D and 3D models. To better understand the different levels of smoothness between parameter and response, we conduct a 1D solute transport test with instantaneous injection at the origin 135 initially, whose analytical solution is given as: 2 () (, ) exp 4 2 L L mxvt cx t Dt Dt φπ ⎡ ⎤ − =− ⎢ ⎥ ⎣ ⎦ (6.17) Substitute , / LL DvvKhL α φ ==Δ into equation (6.17), we obtain: 2 1 (, ) exp 42 4 4 LL L L mhtxxL cx t K LhtK ht K L φ αφ α α πα φ ⎡ ⎤ ⎛⎞ ⎛ ⎞ Δ =− +− ⎢ ⎥ ⎜⎟ ⎜ ⎟ Δ Δ ⎝⎠ ⎝ ⎠ ⎣ ⎦ (6.18) Consider the conductivity K as the input, thus equation (6.18) is in the form of: 1 1 2 12 34 () exp cK aK a K a a K − − ⎡ ⎤ =++ ⎣ ⎦ (6.19) If we approximate equation (6.19) using a Taylor series expansion, we can see that the smoothness between the input K and the output c mainly depends on the coefficients a 2 and a 4 , which clearly explains the reason that the input-output relation is strongly nonlinear especially when the dispersivity L α is small with a large Peclet number. On the other hand, in the xTPCM, we first normalize the concentration by a maximum value c max (c min = 0 is constant), and then set a relative concentration level c and search for the corresponding locations as: () 11 22 max 4ln () , ( , ) 4 L L ht c mht ct K xct K K LL ht L α φφ πα φ − Δ− ⎛⎞ Δ ==± ⎜⎟ Δ ⎝⎠ (6.20) Equation (6.20) describes functions of conductivity K in the forms of: 11 22 max 1 2 3 () , ( ) cK bK xK bK bK − ==± (6.21) It can be seen that the alternative variables c max and x in equation (6.21) are relatively 136 linear with regard to the input K. As for the dTPCM, the location and scaling factors are replaced by the displacement and target location value, which are also smooth with respect to the parameters. Similar smooth properties for the arrival time and scaling factors can be observed in the tTPCM, as well. The xTPCM and the dTPCM may encounter the “out of bound” problem. For example, in the 2D solute transport test in section 6.2.2, the solute may reach the boundary of the space domain (outer boundary) when the observation time is long enough. If we use the dTPCM in this case, we need to enlarge the domain area. Another example is the 3D multiphase flow test in section 6.2.3, where the breakthrough occurs in some of the production wells (inner boundary). In this case, we cannot accurately approximate the displacement whose corresponding target location is near the well using the xTPCM or the dTPCM, since the saturation shock front will stop at the position of well after breakthrough. However, the tTPCM does not suffer from the this impairment. In addition, the xTPCM and the dTPCM are only appropriate for the variables that are defined in the whole space domain (e.g., saturation and concentration), but not appropriate for the variables that are defined in some specific positions (e.g., water cut). On the other hand, the tTPCM is able to tackle both kinds of variables. Furthermore, since the xTPCM and the dTPCM (as well as the PCM) analyze the model response at a certain time, their computational efforts are proportional to the number of time steps. On the contrary, the tTPCM analyzes the model response along the time axis, thus it is convenient to obtain the results for many time steps. Last but not lease, it is more convenient to analysis the arrival time when the model response reaches a threshold using the tTPCM, which is 137 especially important in contaminant control. However, the simulator must be run beyond the current analyzing time in the tTPCM. Therefore, we recommend determining the conditions and requirements of the problem and then selecting the appropriate methods. It is worth noting that the idea of transform in this work is similar to preconditioning, which is a procedure that conditions a given problem into a form that is more suitable for calculation. As is well known, the preconditioning for linear system attempts to improve the spectral properties of the coefficient matrix. Whereas in this study, we transform the original model response to an alternative variable, in order to improve linearity. 138 Chapter 7 Transformed Ensemble Kalman Filter (TEnKF) In the traditional EnKF, we update the state vector in fixed locations from the Eulerian point of view (i.e., examining fluid motion that focuses on specific gridblocks through which the fluid quantity changes as time passes, which can be visualized by sitting on the bank of a river and watching the water pass a fixed location). However, in the TEnKF, we update the displacement from the Lagrangian point of view (i.e., examining fluid motion that focuses on an individual fluid parcel whose location changes as time passes, which can be visualized by sitting in a boat and drifting down a river). We think of the state in each gridblock as an individual fluid parcel and track its displacement, and treat the state in the whole space as an elastic membrane and record its deformation. Since the nonlinearity and non-Gaussianity is mainly observed in the advection-dominated dynamic state, including saturation and concentration, with non-physical oscillation and inconsistency problem, the transform process is only applied to these variables, while the rest of the state variables are treated as in the standard EnKF. Note that the transform process for the TEnKF in this chapter is similar to that for the dTPCM in Chapter 5. Anyhow, we still introduce the forward and backward transform processes again here for the completeness of the TEnKF, and for the convenience of readers who are only interested in data assimilation. 139 This chapter presents the work in paper which is under review by Liao and Zhang [2014c]. 7.1 Theoretical Framework 7.1.1 Forecast Stage The TEnKF has four stages for each time step. The first stage in the TEnKF is the same as in the EnKF ,1, ( ), 1, 2,..., ss fu kj k j gj Ne − == (7.1) 7.1.2 Forward Transform Stage The second stage is to transform the state to the displacement, which can be achieved using digital image processing techniques, specifically motion analysis. Motion analysis, i.e., determining the movement of objects from a time-varying image sequence, has been applied to a wide range of areas [Aggarwal and Nandhakumar, 1988; Vega-Riveros and Jabbour, 1989; Qin et al., 2007; Huang et al., 2009]. A common starting point for motion analysis is to assume that the pixel intensity f at location x in the undeformed image (Figure 7.1a) is the same as the intensity h at location x+d in the deformed image (Figure 7.1b) as [Hild and Roux, 2006] () ( ) xxd fh = + (7.2) 140 where d is the displacement vector. Note that the “image” here contains not only a 2D frame (x = (x, y) T , d = (d x , d y ) T ), but also a 3D scene/vision (x = (x, y, z) T , d = (d x , d y , d z ) T ), as well. Figure 7.1: Illustration of digital image correlation: (a) undeformed image; (b) deformed image with displacement vector; and (c) normalized cross-correlation coefficient. The displacement vector is determined through searching for the maximum of the correlation coefficient. Among various approaches for motion analysis, the two most widely used are the correlation-based approach and the gradient-based approach [Barron et al., 1992; Jahne, 2005]. The correlation-based approach, also referred to as the digital image correlation (DIC) processes, involves considering a small template window in the undeformed image (Figure 7.1a) and searching for the displacement, which produces the best match among all of the possible target windows in the deformed image (Figure 7.1b). The target window, moving in the deformed image within a certain search range based on the prior motion estimation, is usually the same size as the template window. To describe the similarity of the template and target windows, a normalized cross-correlation coefficient is defined as [Lewis, 1995] 141 22 () () ( ) ( ) (, ) () () ( ) ( ) D DD ff h h d f fd h h d γ ∈ ∈∈ ⎡⎤⎡ ⎤ −+−+ ⎣⎦⎣ ⎦ = ⎡⎤ ⎡ ⎤ −+−+ ⎣⎦ ⎣ ⎦ ∫ ∫∫ x xx x x xd xd x dx xx x xd xd x (7.3) where D is the window size; and f h denotes the mean value of f h within this window. The cross-correlation coefficient is within the range of [ −1, 1], according to the Cauchy-Schwarz inequality [Di Stefano and Mattoccia, 2003]. Then, the displacement vector d is determined by searching for the maximum of the correlation coefficient γ at any given location x (Figure 7.1c), and the resolution can be improved by using the sub-pixel interpolation [Schreier et al., 2000]. On the other hand, the gradient-based approach relies on a gradient constraint equation, which is similar to the continuity equation that states the mass conservation in fluid mechanics [Jahne, 2005]. Usually, the gradient-based approach is faster than the correlation-based approach; however, it is only suitable for small displacement problems, although some remedial techniques can somewhat relieve this limitation [Battiti et al., 1991]. In this study, we utilize the correlation-based approach, while other approaches of motion analysis can be applied, as well. To further illustrate the details of motion analysis, we show a simple example, as shown in Figure 7.2. Consider an undeformed image (Figure 7.2a) of size 11 by 11 pixels, the deformed image (Figure 7.2b) is generated through shifting the undeformed image to the right by 2 pixels and to the top by 1 pixel. To determine the displacement of the center gridblock (x = 6, y = 6), we set the half template size to be Nx = 1 and Ny = 1, which leads to a template window size of 2Nx + 1 = 3 by 2Ny + 1 = 3 pixels (shown as the 142 dashed square). The target window (shown as the solid square) has the same size as the template window. We then set the half search range to be Mx = 3 and My = 3, which leads to a search area of 2Mx + 1 = 7 by 2My + 1 = 7 pixels (shown as the dash-dotted square). Next, we move the target window in the search area (Figures 7.2b-g) and calculate the correlation coefficients, as in equation (7.3), with regard to different displacement vectors. Finally, we obtain a correlation coefficient matrix which is the same size as the search area (Figure 7.2h), and determine the true displacement (dx = 2, dy = 1) by directly searching for the maximum of the correlation coefficients. In this work, the template size is determined through a sensitivity analysis, in which we compare the results for different template sizes and select the template size with minimum error. The search range is set according to the observation of displacement in trial tests. We will demonstrate the details in the following 2D case. 143 Figure 7.2: A simple example of forward transform: (a) undeformed image; (b)-(g) deformed images with different displacement vectors; and (h) normalized cross-correlation coefficient. The dashed square is the template window; the solid square is the target window; and the dash-dotted square is the search area. After presenting the concept of motion analysis, we now demonstrate how to introduce it into the TEnKF. Motion analysis in original digital image processing usually refers to determining the movement of an object in a time-varying image sequence (i.e., in the temporal domain); whereas, the motion analysis for data assimilation refers to determining the movement of an object between ensemble members (i.e., in the stochastic domain). The latter one is somewhat more complicated, as the temporal domain is a one-dimensional space, while the stochastic domain could be a multi-dimensional space. In this work, the deformed image, which can be recorded as ( ), 1,2,..., x j sj Ne = , is the state vector in each ensemble member. The undeformed image, which can be recorded as 0 () x s , is set to be the state vector when the static state takes 144 the value of the ensemble average. Note that if the ensemble average is not included in the ensemble, an extra model run is necessary. In the case of multiple assimilation steps, the undeformed image is generated from the initial time using the ensemble average of the static state obtained in the previous assimilation step. For example, considering the assimilation step from t k −1 to t k , we calculate the ensemble average of the static state at t k −1 , and use it to run the forward model from t 0 to t k −1 to obtain the saturation at t k −1 as the undeformed image. The displacement field is obtained through the motion analysis of the undeformed image 0 () x f s and the deformed image () x f j s , and recorded as () dx f j . In other words, if the template location in the undeformed image is x with the corresponding value 0 () x f s , the best match target location in the deformed image will be xd f j + with the corresponding value () xd f f jj s + . This target location value () xd f f jj s + is acquired via gridded interpolation (we utilize linear interpolation in this work, i.e., bilinear/trilinear interpolation in the 2D or 3D cases) and will be used in the backward transform stage. The computational effort for motion analysis depends on the number of template locations to be analyzed, since we need to perform the DIC process at each location. In this study, the template locations are the same as the gridblocks in the simulator; for example, if a 2D case has 50 by 50 gridblocks, then the motion analysis is performed on all of these locations. However, if the motion analysis is only performed on some selected sparse points, the number of template locations can be reduced with increased efficiency. For instance, if only the even-numbered gridblocks are considered, there will be only 25 145 by 25 template locations to be calculated; however, this will reduce accuracy. Put simply, the state variable () x f j s is now represented by an alternative displacement variable () dx f j . 7.1.3 Update Stage The third stage is the analysis or data assimilation stage, in which the state variables are updated to honor the observations. Note that the transform process is only needed for the advection-dominated dynamic state (e.g., saturation or concentration), while the standard process is used for the static state (e.g., permeability or porosity) and other dynamic states (e.g., pressure). For instance, if an augmented state vector contains log permeability, porosity, pressure, and saturation as [ln , , , ] w sk pS TT T TT φ = in the EnKF, the transformed augmented state vector will contain log permeability, porosity, pressure, and displacement as [ln , , , ] sk pd TT T TT φ = in the TEnKF. The reason for this is that the strong nonlinearity/non-Gaussianity, which leads to non-physical updated results, is observed only in the advection-dominated dynamic state. On the other hand, the static state and other dynamic states can be updated using the standard EnKF. We will demonstrate the implementation details of the TEnKF in the following examples. Similar to the EnKF, the mean s f k and covariance k P can be estimated from the ensemble as , 1 1 () ss s Ne f ff kk kj j E Ne = == ∑ (7.4) 146 ()() ,, 1 1 1 Pssss Ne T ff ff kkjkkjk j Ne = =− − − ∑ (7.5) The Kalman gain is then calculated from the covariance as 1 () KPH HPH R TT kk k − =+ (7.6) where [, ] H0I = is the observation operator; and R is the covariance matrix of the observation errors. Each ensemble member is updated via ,, , , ( ), 1, 2,..., ss Kobs Hs uf f kj kj k k j k j j Ne =+ − = (7.7) where obs is the perturbed observation. The forms of the above equations for the TEnKF are the same as in the EnKF, but the advection-dominated dynamic state (e.g., saturation or concentration) is replaced by the displacement in the state vector. 7.1.4 Backward Transform Stage Because we focus on the location changes in the TEnKF, we update the displacement instead of state, and assume that the state value is invariable before and after the update as ,, ,, () () xd xd ff u u kj kj kj kj ss += + . Since we know the value of the state , u kj s at some scattered locations , xd u kj + , we need to transform it back to obtain the updated state in the original meshed gridblocks as , () x u kj s (Figure 7.3a). In this study, we call the “scatteredInterpolant” function in MATLAB for the scattered data interpolation (MATLAB scatteredinterpolant class, 2014). This algorithm involves constructing an interpolating surface via triangulating the points and lifting the vertices by a magnitude 147 into the response space, which is orthogonal to the parameter space (Figure 7.3b). We choose the default linear interpolation via calculating the weighted sum of the values on the three vertices of the enclosing triangle, while other options, such as the nearest-neighbor interpolation can be used, as well. Although the illustration highlights 2D interpolation, this technique can be applied to higher dimensions. For example, in 3D interpolation, tetrahedrons will be constructed instead of triangles. In sum, the displacement , () dx u kj is transformed back to obtain the state variable , () x u kj s . Figure 7.3: Illustration of scattered data interpolation: (a) the locations (A’, B’, and C’) from the previous updated stage are usually scattered, but we need the value of the query point (D’) in the meshed gridblocks; (b) find the triangle (A’, B’, and C’) that encloses the query point (D’) and calculate the value of the query point (D) through computing the weighted sum of values of the three vertices of the enclosing triangle (A, B, and C). 7.1.5 Framework of TEnKF In the EnKF, there are two stages (i.e., the forecast and the update stages); whereas, in the TEnKF, there are two more stages (i.e., the forward and backward transforms between 148 the state and displacement) (Figure 7.4). In some cases, especially for advection-dominated problems (e.g., water flooding), the relation between displacement and observation is usually more linear than the relation between state and observation. This linearity, hence, improves the performance of the Kalman update. Figure 7.4: Illustration of the differences between the standard EnKF and the transformed EnKF (TEnKF). 7.2 Case Studies 7.2.1 One-dimensional Homogeneous Case In the first case, we consider a one-dimensional piston-like displacement of oil by water through a homogeneous porous medium of length 1.0, with 100 gridblocks and a cross sectional area 1.0. We set a constant injection rate of 0.05, and an oil viscosity which is equal to the water viscosity as 1.0. Assume a Corey-type relative permeability function as k rw = S w 2 and k rw = (1 −S w ) 2 . The true porosity is 0.2, and 50 porosities are randomly generated following a Gaussian distribution, with a mean of 0.18 and a 149 standard deviation of 0.05 as the initial guess. Since we assumed a homogeneous medium, the porosity is constant in the entire space for each realization. As the boundary condition is related to a constant injection rate, the flooding front location is not sensitive to permeability, and absolute permeability is set to be 0.001. The injection well BHP is fixed at 100, and the production well BHP is observed at time 1.0. There is only one assimilation step in this simple test. In the Buckley-Leverett displacement, the forward model is solved analytically. Although the Kalman filter can be applied for this problem with an analytical solution, we utilize an ensemble treatment to demonstrate the behavior of the EnKF. Figure 7.5a shows the prior (thin blue curve) and the reference (thick black curve) pressure profiles before the Kalman update. Since the prior is very close to the reference before the update, the pressure is reasonable after the standard Kalman update, as shown in Figure 7.5b. Figure 7.5c shows the saturation profiles before the Kalman update. Figure 7.5d shows the profiles after the Kalman update with negative saturation values. This non-physical phenomenon is due to the strong nonlinearity between the saturation state and the pressure observation. Figure 7.5e illustrates the saturation profiles after the Kalman update with normal score transform (NST) [Goovaerts, 1997; Bertino et al., 2003; Gu and Oliver, 2006; Zhou et al., 2011; Schöniger et al., 2012]. The NST function can be derived through a graphical correspondence between the cumulative distributions of the original and standard normal variables. Let F(z) and G(y) be the stationary one-point cumulative density functions (CDFs) of the random function Z(u), and the standard normal random function Y(u) as 150 F(z) = Probability {Z(u) < z} and G(y) = Probability {Y(u) < y}. The NST allows one to go from the Z(u) with any cdf F(z) to the Y(u) with standard Gaussian cdf G(y) as Y(u) = G −1 [F(Z(u))] [Goovaerts, 1997]. Although the NST process successfully constrains the saturation to be non-negative, some non-physical phenomena (e.g., low saturation value in upstream and high saturation value in downstream) still exist [Gu and Oliver, 2006]. Similar problems are observed in the saturation profiles after the Kalman update with error function transform (EFT) [Jafarpour and McLaughlin, 2009], as shown in Figure 7.5f. The error function “erf” is a bijective function (one-to-one correspondence) mapping from ( −Inf, Inf) to ( −1, 1), which is related to the cumulative distribution of the standard normal distribution. In the Kalman update with EFT, we first transform from saturation to a pseudo-saturation as S erf = erf −1 (2S −1), then update the S erf using the standard Kalman update, and finally back-transform from pseudo-saturation to saturation as S = erf(S erf +1)/2. The EFT can be considered as a “global” and “analytical” NST, since the EFT transforms the state using a universal analytical function for various locations; whereas, the NST transforms the state using individual empirical function in each location from the ensemble. Figure 7.5g demonstrates the saturation profiles after the Kalman update with iteration [Gu and Oliver, 2006]. In the iterative method, the porosity (static state) is updated and simulated from the previous time step to obtain the corresponding saturation (dynamic state). More specifically, for the assimilation step from the previous time to the current time, we update the static parameter as in the standard Kalman correction, and go one step backward and start the simulation based on the updated static parameter and 151 previous time dynamic parameter to obtain the current time dynamic parameter. Then, the difference of the new computational data and the observation is used to update the state variables. If the state variables are still non-physical or inconsistent, the iteration step can be repeated until the corrected saturation profiles satisfy the physical bounds or the iteration exceeds a certain maximum number (we iterate once without repetition in this work). We can see that the results are improved compared to those from the standard Kalman update (Figure 7.5d), but the non-physical problem still occurs. Figure 7.5h shows the saturation profiles after the Kalman update with restart, in which the porosities are updated and then simulated from the initial time to obtain the corresponding saturations. These saturation values are physically meaningful, since they are generated from forward models. However, additional time is required to rerun the simulator, and the computational cost may increase as the number of time steps increases. 152 Figure 7.5: Illustration of the non-physical problem: (a) pressure before the Kalman update; (b) pressure after the standard Kalman update; (c) saturation before the Kalman update; (d) saturation after the standard Kalman update; (e) saturation after the Kalman update with normal score transform (NST); (f) saturation after the Kalman update with error function transform (EFT); (g) saturation after the Kalman update with restart; and (h) saturation after the Kalman update with iteration once. The thin curves are the ensemble, and the thick curve is the reference. To further analyze the oscillation phenomenon, we illustrate the relations between the static state, the dynamic state, and the observation. Figure 7.6a shows a nonlinear relation between porosity and pressure at x = 0.4. Figure 7.6b shows a discontinuous (also nonlinear) relation between porosity and saturation at x = 0.4. Figures 7.6c, 7.6d, and 7.6e illustrate the porosity, pressure, and saturation before the Kalman update, respectively. Figures 7.6f, 7.6g, and 7.6h illustrate the corresponding probability distribution functions (PDFs). Figures 7.6i, 7.6j, and 7.6k depict these state variables after the Kalman update. Figures 7.6l, 7.6m, and 7.6n depict the corresponding PDFs. We can see that the strong nonlinearity of static state and dynamic state (Figures 7.6a and 7.6b) is translated to the 153 strong nonlinearity of observation and dynamic state (Figures 7.6d and 7.6e). The former will bring about difficulty in uncertainty quantification or sensitivity analysis (forward modeling), while the latter will cause problems in data assimilation or history matching (inverse modeling). In this case, the pressure directly updated by the Kalman correction is illustrated as a blue cross in Figure 7.6j. We also show the pressure generated through rerunning the forward model using the updated porosity as a red circle in Figure 7.6j. We can see that the directly updated pressure is different from the pressure with the restart process; however, both of them are within a reasonable range and close to the exact value, which is shown as a thick black line in Figure 7.6j. This is because the standard deviation of pressure (1.04) is small compared to the mean pressure (73.46) before the update. Another reason is that the pressure distribution is approximately Gaussian (Figure 7.6g). On the other hand, the directly updated saturation (blue cross in Figure 7.6k), which is far from the saturation with restart process (red circle in Figure 7.6k), reveals obvious non-physical problem with negative values. This is because the standard deviation of saturation (0.3306) is large compared to the mean saturation (0.1938) before the update, and also due to additional non-Gaussianity, as well as nonlinearity (Figure 7.6h). 154 Figure 7.6: Illustration of the nonlinearity/non-Gaussianity that causes the non-physical problem in saturation estimation: (a) the relation between porosity and pressure; (b) the relation between porosity and saturation; (c) porosity before the Kalman update; (d) pressure before the Kalman update; (e) saturation before the Kalman update; (f) PDF of porosity before the Kalman update; (g) PDF of pressure before the Kalman update; (h) PDF of saturation before the Kalman update; (i) porosity after the Kalman update; 155 (Continued) (j) pressure after the Kalman update (the blue cross indicates the pressure directly updated by the Kalman correction, and the red circle indicates the pressure calculated through rerunning the model using the updated porosity from the initial time; similar in (k)); (k) saturation after the Kalman update. (l) PDF of porosity after the Kalman update; (m) PDF of pressure after the Kalman update (the blue line indicates the pressure directly updated by Kalman correction, and the red line indicates the pressure calculated through rerunning the model using the updated porosity from initial time; similar in (n)); (n) PDF of saturation after the Kalman update. The black thick line indicates the true reference. For a better comparison, we also test the proposed TEnKF in this example. The first forecast stage in the TEnKF (Figure 7.7a) is the same as in the EnKF (Figure 7.5c). The second forward transform stage is illustrated in Figure 7.7b, in which the prior saturation profile (blue curve) is considered as the deformed image, and the undeformed image is the saturation profile, in which the porosity takes the value of the ensemble average (red curve). Since the saturation value in the undeformed image is non-trivial in the area of x = 0 to 0.35, the displacement is required only in this range. In this one-dimensional problem, we use linear interpolation instead of the DIC to obtain the displacement, as the former approach is more convenient. The linear interpolation is based on creating a lookup table for saturation and location. For example, the displacement at x 0 = 0.2 can be obtained in this way: in the undeformed image, get the saturation value Sw = 0.794 at x 0 = 0.2; then, in the deformed image, check if 0.794 is between any two values in the adjacent gridblocks and find the gridblocks of x A = 0.35 and x B = 0.36 (Figure 7.7c); next, determine the location where the saturation is Sw C = 0.794 by linear interpolation as x C = (Sw C - Sw B ) / (Sw A - Sw B ) × (x A − x B ) and get x C = 0.358; finally, subtract the x 0 from x C and calculate the displacement as d| x0=0.2 = 0.358 − 0.2 = 0.158 (Figure 7.7d). The displacement at other locations and in other realizations is obtained similarly (Figure 156 7.7e). The third update stage is shown in Figures. 7.7f and 7.7g. It can be seen that the displacement (e.g., at x 0 = 0.2) is linear to the observation (Figure 7.7f) and, thus, is updated easily after the Kalman correction (Figure 7.7g). Similar phenomena are observed for various locations along the space domain (Figure 7.7h). To demonstrate the backward transform, we show one realization of the updated displacement in Figure 7.7i. The updated displacement at x 0 = 0.2 is d| x0=0.2 = 0.011; therefore, the target location after update is x 0 + d| x0=0.2 = 0.211 (marked as A in Figure 7.7j) with a saturation value of 0.794. Usually, we want the final results in meshed gridblocks. For example, to obtain the saturation at the location x = 0.21 (marked as B in Figure 7.7j), we need another updated displacement information at x 0 = 0.19 (marked as C in Figure 7.7j). The saturation can be calculated through interpolation as Sw B =(x B − x C ) / (x A − x C ) × (Sw A - Sw C ). The saturation at other locations is acquired similarly (Figure 7.7k). Figure 7.7l shows that the updated saturation (blue curves) are close to the reference (black curve), and almost the same as the results from the Kalman update with restart (Figure 7.5h). From the above test, we show that the TEnKF not only provides accurate and physically meaningful results, but also avoids rerunning the simulator. Note that the relation between the displacement and the observation before update is almost perfectly linear (Figure 7.7f), since the flow path is unique in this simple one-dimensional problem. As for general 2D or 3D problems, the relation may not be exactly linear, but still relatively linear. 157 Figure 7.7: Illustration of the TEnKF: (a) forecast saturation from direct model runs; (b) forward transform for one of the realizations; (c) detail of forward transform with enlarged scale; (d) forecast displacement for one of the realizations; (e) forecast displacement by forward transform; (f) displacement at x 0 = 0.2 before the Kalman update; (g) displacement at x 0 = 0.2 after the Kalman update; (h) updated displacement; (i) updated displacement for one of the realizations; (j) detail of forward transform with enlarged scale; (k) (j) updated saturation for one of the realizations; (f) updated saturation by backward transform. The blue curves show the ensemble, the black curve is the reference, and the red curve is the saturation generated from the ensemble mean of porosity as the undeformed image. Finally, we compare the material balance by checking the ratio of water-in-place and total injected water (Figure 7.8). We can see that the standard Kalman update cannot guarantee the material balance. After a simple truncation remedy by cutting off the unphysical saturation values, the result may be even worse. This is consistent with Figure 7.5d that if the negative saturation values are adjusted to be 0, then the total mass will be 158 greatly overestimated. The Kalman update with NST has the largest maximum value for material balance (Figure 7.8), because this method is based on “local” normal score transform and, thus, may have outliers (Figure 7.5e). The Kalman update with EFT has the best result for material balance in this test (Figure 7.8); the possible reason for this is that the error function transform considers the “global” saturation data in the entire space. However, the non-physical problem still exists between high/low values in the downstream/upstream (Figure 7.5f). The iterative approach and the proposed transformed approach have similar results for material balance, which are also close to that of the standard Kalman update without truncation (Figure 7.8). On the other hand, the transformed approach provides more reasonable saturation profiles (Figure 7.7l) than the other two (Figures. 7.5d and 7.5g), and it does not require rerunning the simulator. Unfortunately, none of these methods can guarantee a perfect material balance, which can only be achieved by a “thorough” restart process (i.e., rerun the forward model from the initial time using the updated static state). 159 Figure 7.8: Material balance check for different methods. The cross represents the average, the rectangle represents the p10 and p90 (the values at 10% and 90% of the distribution, respectively), and the error bar represents the minimum and maximum values. 7.2.2 Two-dimensional Heterogeneous Case In the second test, consider a five-spot pattern consisting of 51 × 51 × 1 gridblocks, with each block measuring 40 × 40 × 40 ft. There is one injection well in the center with a constant pressure of 5,000 psi and four production wells in the corners with constant pressures of 2,000 psi (Figure 7.9). The initial reservoir pressure is 4,000 psi, and the porosity is constant at 0.2. The relative permeability is k rw = S w 2 and k rw = (1 −S w ) 2 , and the capillary pressure is neglected. We assume the log permeability (in mD) to be a second-order stationary Gaussian random field with a mean of 4.0 and a covariance of 2 ln 1 1 2 2 ln 1 2 1 2 ( , ; , ) exp( ) kk x y Cxyx y x x y y ση η =−− −− , where 2 ln 1.0 k σ = is the variance and 800 ft x η = and 800 ft y η = are the correlation lengths in x and y directions, respectively. Based on the above geostatistics, the Karhunen-Loeve expansion [Ghanem 160 and Spanos, 1991] is used to randomly generate 40 realizations as the ensemble. The true reference is generated similarly, but with a lower mean value of 3.0. The water injection rate, oil production rate, and water cut are observed every 500 days for up to 2,000 days with 4 assimilation steps. Figure 7.9: Schematic diagram of the model setup in the 2D test. Line A 1 -A 2 indicates the places where the results are compared in detail. 7.2.2.1 Implementation Details of TEnKF In this section, we show how to perform the TEnKF. Figure 7.10a illustrates the saturation generated from the mean log permeability at day 500 (i.e., first assimilation step), which is considered as the undeformed image. Figure 7.10e reveals the saturation generated from one of the ensemble members at day 500, which is considered as the deformed image. In the forward transform stage, we need to obtain the displacement field from these two images. The effect of the DIC process is mainly affected by two factors: the template size and the search range. We first set up a large enough search range (we 161 choose 15 pixels for now, which means that the search area is 31 by 31) and determine the template size by a sensitivity analysis. If the template size is too small, the template window cannot capture the feature around the template location, and the displacement will be affected by disturbance. Figure 7.10b shows the displacement field when the half template size is set to 1 (i.e., the window is 3 by 3 pixels), from which we can see that the directions of the displacement vectors are random in some places. On the contrary, if the template size is too large, the template window will contain locations far from the template location, and these locations are actually not moving similarly to the template location. This will make the displacement field too smooth, and also increase computational cost. Figure 7.10d shows the displacement field when the half template size is set to 10 (i.e., the window is 21 by 21 pixels), from which we can see that the directions of the displacement vectors are similar to each other in the adjacent places. Figure 7.10c depicts the displacement when the half template size is set to 6 (i.e., the window is 13 by 13 pixels). The goal is to choose the best template size in this way: if the template size is appropriate, the displacement field should be accurate, and we should be able to reproduce the deformed image from the undeformed image and the displacement. In other words, the forward transform and backward transform should be a set of invertible transforms. Figures. 7.10f-7.10h demonstrate the reproduced images, which are all similar to the original image (Figure 7.10e). To quantitatively compare the results, we define a relative error as reproduced original original () () () x x xxx xx D D IdId d error Id d ∈ ∈ − = ∫ ∫ (7.8) 162 where Id is the deformed image; and D is the whole image domain. Figure 7.11a illustrates the error for 20 realizations randomly selected from the ensemble with regard to a half template size from 1 to 10. The thick black line clearly shows that when the half template size is too small (from 1 to 2) the mean error is large. Then, the error decreases and stays at a low level as the half template size grows (from 3 to 10). If the half template size is extremely large, the error will increase, which is not shown here. We observed that the relative error is not very sensitive to the half template size, as long as the latter is within a reasonable range (from 3 to 10). Considering both accuracy and efficiency, we use a half template size of 6 in this case. Next, we determine the search range. In general, if the search range is smaller than the true displacement, the estimated displacement will be bounded in the range and, thus, inaccurate. In contrast, if the search range is larger than the true displacement, the estimated displacement will be correct, but it will require additional time. Therefore, we should use the smallest search range that is just large enough. Figures. 7.11b and 7.11c show the estimated displacement in the x- and y-directions, respectively. We can see that most of the absolute values of the displacement are within 10 pixels and, thus, set the search range to 10 in both the x- and y-directions. After finishing this paper, we became aware of the recent work [Mariethoz et al., 2010], which bears similarities to the present approach in terms of window size and search range. Although their work focuses on direct sampling method to generate ensemble of stochastic realizations, some of their ideas (e.g., using a circle for template window instead of a rectangle) can be adopted in our study, which may be considered in future analysis. 163 Figure 7.10: Sensitivity analysis for selecting the appropriate template size: (a) saturation generated from the ensemble mean as the undeformed image; (b)-(d) displacement fields calculated by forward transform using different template sizes; (e) saturation generated in one realization as the original deformed image; (f)-(h) reproduced deformed image by backward transform using different template sizes. The goal is to check if the reproduced deformed image can match the original deformed image. Figure 7.11: Determine the appropriate template size and search range: (a) relative error vs. half template size in 20 realizations (thin blue lines) and the mean error (thick black line); (b) histogram of displacement in the x-direction; (c) histogram of displacement in the y-direction. 164 7.2.2.2 Estimation of Saturation and Permeability with Material Balance Check Figure 7.12 shows the updated saturation at day 1,000 and 2,000 (i.e., second and fourth assimilation step) using different methods. Figure 7.12a shows the reference saturation generated from the true model parameter. Figure 7.12b reveals non-physical phenomenon with values less than S wi (in white) using the standard EnKF. Figure 7.12e shows the detailed results with 40 realizations along the line A 1 -A 2 . Figure 7.12f depicts the results after a simple truncation remedy by forcing the values to be non-negative; however, the non-physical problem with high/low values in downstream/upstream still exists. Since the saturation value has to be within the physical range of [0, 1] as required in the simulator, the truncation remedy is performed for the standard EnKF in each assimilation step. In contrast, the results from the REnKF (Figures 7.12c and 7.12g) are physically meaningful, since they are generated by rerunning the forward model using the updated permeability from the initial time. The results from the TEnKF (Figures 7.10d and 7.10h) are also reasonable, but without additional model runs for each realization. Figures 7.12i-7.12p show saturation after the final (fourth) update at day 2,000 with similar conclusions. 165 Figure 7.12: Water saturation updated at day 1,000 and 2,000 using different methods: the first row is the reference and mean saturation in the whole space at day 1,000; the second row is the saturation along line A 1 -A 2 at day 1,000; the third row is the reference and mean saturation in the whole space at day 2,000; and the fourth row is the saturation along line A 1 -A 2 at day 2,000. In the first and third rows, the saturation values are indicated by color, and the white area in (b) and (j) indicate that the values are out of the range from 0 to 1 (negative in this case). In the second and fourth rows, the thin curves show the ensemble prediction, and the thick black curve is the reference. In order to check the material balance for standard EnKF and TEnKF in 40 realizations from the first to the last assimilation step, we compare the water-in-place to the total injected water (Figure 7.13). It can be seen that the standard EnKF overestimates the water-in-place (estimated from updated saturation with truncation remedy), because 166 all negative saturation values are set to be 0. It can also be seen that the ratios from the TEnKF are closer to 1 and more reasonable. The results from the REnKF are not shown, since the saturations are generated through rerunning the simulator from the initial time and, thus, the material balance is always satisfied. Figure 7.13: Material balance check for standard EnKF and TEnKF in 40 realizations from the first to the last assimilation step. Figure 7.14 illustrates the estimations of the permeability field. Figure 7.14a shows the true log permeability as a reference. Figures 7.14b-7.14d show the mean log permeability of the initial guess, which are the same in the EnKF, the REnKF, and the TEnKF. Figures 7.14e-7.14g show the mean log permeability after the first assimilation step, which are still identical in three methods, since the update procedures for log permeability are exactly the same. However, the update procedures for saturation are not 167 and, hence, affect the results in subsequent assimilation steps. Figures 7.14h-7.14j show the mean log permeability after the last assimilation step. We can see that the results from the REnKF and the TEnKF are close to the reference, while the result from the EnKF is significantly higher. We think that the reason is as follows: if the EnKF overestimates the water saturation in the current assimilation step, the water-in-place/oil-in-place is, thus, overestimated/underestimated; hence, the forecasted oil production rate in the next assimilation step will be underestimated (considering the relative permeability curve). To match the observed oil production rate that is higher than the forecasted value, the estimated permeability will be increased after the update. This is how the saturation (dynamic state) affects the permeability (static state) in this test. 168 Figure 7.14: Log permeability in the initial guess and after the first and final update: the first column is the true log permeability; the second column is the mean log permeability by standard EnKF; the third column is the mean log permeability by REnKF; and the fourth column is the mean log permeability by TEnKF. 7.2.2.3 Match and Prediction of the Observation The estimated production data of P2 from the initial guess and the updated permeability in the final assimilation step (at day 2,000) using the EnKF, the REnKF, and the TEnKF are plotted in Figure 7.15. The thin curves indicate the ensemble forecasts, and the thick black curve indicates the reference. We can see that the updated estimations in all three methods are improved compared to the initial guess (Figures 7.15a and 7.15e) 169 and match the reference during the assimilation period (day 0 to 2,000). Whereas, in the prediction period (day 2,000 to 8,000), the results from the REnKF (Figures 7.15c and 7.15g) and the TEnKF (Figures 7.15d and 7.15h) still match the reference well; however, those from the EnKF (Figures 7.15b and 7.15f) are deviated. Figure 7.15: The production data of P2 estimated from the initial guess and from the final updated permeability using the EnKF, the REnKF, and the TEnKF. The first row is the oil production rate; and the second row is the water cut. The thin curves show the ensemble prediction, and the thick black curve is the reference. 7.2.3 Three-dimensional PUNQ-S3 Case The PUNQ-S3 case study is a small-size reservoir engineering model set up by the PUNQ (Production forecasting with Uncertainty Quantification) project and commonly used for performing benchmarks [Floris et al., 2001; Gu and Oliver, 2005]. A full description of this case study can be found on the PUNQ-S3 webpage [PUNQ-S3 test 170 case, 2014]. The model contains 19×28×5 gridblocks, the dimensions of which in the x- and y-directions are 180×180 m 2 , and the thickness is defined in the online dataset. The field contains 6 production wells, but no injection wells, as it links to a strong aquifer. The true porosity and permeability data are given and used to produce the observation. We follow the procedure on the PUNQ-S3 webpage and generate 40 realizations of the static state including porosity, horizontal permeability, and vertical permeability using a sequential Gaussian simulation/co-simulation module in SGeMS (Stanford Geostatistical Modeling Software) [SGeMS webpage, 2014] based on Gaussian random fields. The original PUNQ-S3 project uses the production history of the first 8 years, including 1 year of well testing, 3 years of field shut-in, and 4 years of actual production. In our test, we skipped the first 4 years of testing or shut-in, and start actual production from the beginning, since our goal is to test the EnKF and the TEnKF, and compare them to the REnKF. The production rate is fixed at 150 m 3 /d, and all wells are shut-in at the first 2 weeks in each year, similar to the original model setup. The production data, including well bottom-hole pressure (WBHP), well oil production rate (WOPR), well gas-oil ratio (WGOR), and well water cut (WWCT) are assimilated once a year for up to 5 years; therefore, the total oil production is close to that in the original model setup, in which 8 years of production history are used, including 3 years of shut-in. In the EnKF, the water and gas saturation, as well as the pressure, are updated along with the porosity and permeability. In the REnKF, the pressure and saturation are produced using the updated porosity and permeability through rerunning the simulator from time 0. In the TEnKF, the pressure is updated along with the porosity and 171 permeability, while the water and gas saturations are updated separately. First, the saturation is transformed forward to the displacement; then, the displacement is updated; and, finally, the displacement is transformed back to obtain the saturation. The DIC process in the 3D case is similar to that in the 2D case, but involves displacement vectors in the x- y- and z-direction, and the template window and the search area will be cuboids instead of rectangles. We observed that in this PUNQ-S3 case, the geometric shape of the water or gas saturation is not similar in each layer. The reason is that in this model, the porosity and permeability are generated individually for each layer without any correlation. Therefore, aside from performing the full 3D displacement transform for the saturation in the whole space (referred to as TEnKF3d), we also test the 2D displacement transform for the saturation in each layer (referred to as TEnKF2d). Although the true movement of the saturation in 3D space cannot be captured by the horizontal direction only (since the buoyancy effect will drive the fluid in vertical direction), the saturation in each layer still has similarity in different realizations (e.g., if the fluid flows vertically from layer A to layer B, the size of fluid saturation in layer A/B will shrink/grow). 7.2.3.1 Estimation of Saturation, Porosity and Permeability Figure 7.16 shows water saturation in layer 5 from the initial time to the end of the assimilation steps. The true saturation generated using the exact porosity and permeability is illustrated in the first row. The mean saturation updated using the EnKF is illustrated in the second row. It can be seen that the saturation values are larger than 1 in 172 some places (as dark red). The updated saturations in the REnKF, the TEnKF2d and the TEnKF3d are illustrated in the third, fourth and fifth rows, respectively, with reasonable results. Figure 7.17 shows the truth and estimated mean gas saturation in layer 1. Similar to the water saturation, the gas saturation updated from the EnKF is larger than 0.8 (i.e., 1 −S wi , where S wi = 0.2) in some places (as dark red). Additionally, a saturation value significantly lower than the reference is observed near the gird block (x = 15, y = 16). 173 Figure 7.16: Water saturation of layer 5 in each assimilation step: the first row is the true saturation, calculated from the exact permeability and porosity; the second row is the mean saturation, estimated using the EnKF; the third row is the mean saturation, estimated using the REnKF; the fourth row is the mean saturation, estimated using the TEnKF with 2D motion analysis for each layer; and the fifth row is the mean saturation, estimated using the TEnKF with 3D motion analysis. 174 Figure 7.17: Gas saturation of layer 1 in each assimilation step: the first row is the true saturation, calculated from the exact permeability and porosity; the second row is the mean saturation, estimated using the EnKF; the third row is the mean saturation, estimated using the REnKF; the fourth row is the mean saturation, estimated using the TEnKF with 2D motion analysis for each layer; and the fifth row is the mean saturation, estimated using the TEnKF with 3D motion analysis. The true porosity in layer 1 is demonstrated in Figure 7.18a. The mean porosity updated using the EnKF in each assimilation step, is shown in Figures 7.18b-7.18g. We 175 can see that the EnKF overestimated the porosity in the bottom-left corner, and underestimated the porosity in the center area. In contrast, the REnKF, the TEnKF2d and the TEnKF3d captured the feature in general. Compared to the TEnKF3d, the TEnKF2d provided estimation closer to the true values. Similar conclusion can be drawn in Figure 7.19, which shows the horizontal log permeability (with base 10 to be consistent with the PUNQ webpage) in layer 3. It can be seen that the EnKF underestimated the log permeability in the bottom-left corner, and overestimated the log permeability near the position (x = 10, y = 20). Whereas the REnKF, the TEnKF2d and the TEnKF3d yielded correct characteristics (high or low value in different regions), albeit the TEnKF3d had a larger variance than the reference. 176 Figure 7.18: Porosity of layer 1 in each assimilation step: (a) is the true porosity. The first row except for (a) is the mean porosity, estimated using the EnKF; the second row is the mean porosity, estimated using the REnKF; the third row is the mean porosity, estimated using the TEnKF with 2D motion analysis for each layer; and the fourth row is the mean porosity, estimated using the TEnKF with 3D motion analysis. 177 Figure 7.19: Horizontal log permeability of layer 3 in each assimilation step: (a) is the true log permeability. The first row except for (a) is the mean log permeability, estimated using the EnKF; the second row is the mean log permeability, estimated using the REnKF; the third row is the mean log permeability, estimated using the TEnKF with 2D motion analysis for each layer; and the fourth row is the mean log permeability, estimated using the TEnKF with 3D motion analysis. 7.2.3.2 Match and Prediction of the Observation Figure 7.20 shows the history matching (from year 0 to 5) and prediction (from year 5 to 10) results for well PRO-11 generated from the initial guess (first row), and the final updated porosity and permeability. The thin curves indicate the ensemble forecasts, and the thick black curve indicates the reference. We can see that the EnKF (second row) 178 roughly matches the production history, but fails to predict the forecast period for WBHP, WOPR and WWCT, due to poor estimation in porosity and permeability. However, the REnKF, TEnKF2d and TEnKF3d all match the reference well. Figure 7.20: History matching and prediction of PRO-11: the first row is from the initial guess; the second row is from the EnKF; the third row is from the REnKF; the fourth row is from the TEnKF with 2D motion analysis for each layer; and the fifth row is from the TEnKF with 3D motion analysis. The first column is the well bottom-hole pressure; the second column is the well oil production rate; the third column is the well gas-oil ratio; and the fourth column is the well water cut. 179 7.2.3.3 Computational Cost Comparison. Figure 7.21a shows the CPU time as a function of assimilation step in all three methods. Figure 7.21b shows the time ratio, which utilizes the EnKF time as a baseline. We can see that the REnKF time is approximately twice as much as the EnKF time. The extra time in the REnKF is used for rerunning the simulator from the initial time to the current assimilation step. Although the length of the total simulation period increases as the number of assimilation step increases, the REnKF CPU time does not change much (although it increases slightly) in each step. We think that the reason for this is that, in this test, the partial differential equations in the simulator converge very fast (actually, the black-oil model is easy to solve); thus, the time step in the simulator can be very large. Therefore, the CPU time used for rerunning the simulator from the initial time to the current assimilation step is not much longer than the CPU time used for running the simulator from the previous assimilation step to the current assimilation step in the forecasting stage. However, it may not be true in other cases. For example, if the well control condition changes rapidly, or the fluid flow equation is coupled with heat transfer equation, the time step in the simulator has to be very small. In that case, the CPU time used for rerunning the simulator from time 0 will be much longer than that used in just running the simulator in one assimilation period, and the REnKF CPU time will be significantly increased. In general, if the number of assimilated time steps is N and the simulator running time in each assimilation step is t, then the total time of model runs in the EnKF is N×t, and the time in the REnKF is at least N×t + N×t = 2N×t and, at most, N×t + (N −1)×N/2×t = (N+1)×N/2×t, depending on the forward simulation performance 180 related to convergence speed. Figure 7.21b also shows that the time of TEnKF2d is approximately one and a half as much as the EnKF time, and the time of TEnKF3d is close to/slightly less than the REnKF time, since the displacement transform in 3D space costs more time than that in 2D space. We point out that in large scale practical cases, the efficiency of the TEnKF3d requires more tests and further investigation. The TEnKF time is expected to be further reduced if we utilize a more efficient compiler (e.g., C/C++ or Fortran) instead of the one that we currently use, i.e., MATLAB. In addition, since the TEnKF does not need to rerun the simulator, the computational cost will not be affected by the simulation period. Figure 7.21: Computational cost comparison between different methods: (a) CPU time; and (b) CPU time ratio, defined as the ratio of the time in each method to the time in the EnKF. 181 7.3 Discussion One of the problems in the TEnKF is that the forward transform from state to displacement may not be unique. We illustrate this using the ninth SPE comparative solution project (SPE9) in Figure 7.22. The SPE9 case was originally organized to investigate the complications for black-oil reservoir simulation in a heterogeneous permeability field (Killough 1995). The model has a dimension of 7200×7500×360 ft, which is represented by 24×25×15 gridblocks. The project involves water flooding from an injection well (located at gridblock x = 24, y = 25, z = 15) at the bottom of a reservoir, which is dipping in the x-direction at an angle of 10 degrees. We assume the log permeability (in mD) is a second-order stationary Gaussian random field, with a mean of 3.0 and a covariance of 2 ln 11 12 2 2 ln 1 2 1 2 1 2 ( , , ; , , ) exp( ) kk x yz Cxyzx y z x x y y z z σ ηη η =−− −− −− , where 2 ln 1.0 k σ = is the variance and η is the correlation length, which is set to be 0.4 of the domain size in each direction. In this study, 50 ensemble members are randomly generated and simulated up to 1 year. The ensemble mean is also simulated to generate the undeformed image in the TEnKF. Figure 7.22a shows the water saturation at layer 15 with regard to the ensemble mean. Figure 7.22b shows the saturation at the same layer for one of the ensemble members. The saturation profile has no variance in the y-direction, due to the dipping in the x-direction. If we apply the DIC process to determine the displacement between Figures 7.22a and 7.22b, we will meet the aperture problem. The aperture problem refers to the ambiguity of motion direction when looking at the visual 182 field through a small window or aperture, because a variety of motion directions can lead to identical observations in the aperture (Figure 7.22c). In some cases, enlarging the aperture can help to resolve this problem. Other solutions are mainly based on the idea of adding more information as additional constraints (e.g., Bruhn at el. 2005). In the SPE9 case, however, there is a simple and effective way to cope with the problem. Since the aperture problem is caused by the invariance of saturation in the y-direction, we can take advantage of the invariance and just consider the motion in the x-direction (Figure 7.22d). In other words, we not only reduce the motion analysis from the 3D to the 2D domain (i.e., determine the displacement in each layer, as in the PUNQ-S3 case), but further reduce the motion analysis from the 2D to the 1D (i.e., determine the displacement in each row). Another problem in the TEnKF is related to the existence of displacement. We illustrate this using the Brugge benchmark study, which tested the history matching and optimization methods in a closed-loop workflow [Peter et al., 2010]. The Brugge model, which was designed to mimic real field scenarios and consists of 139×48×9 gridblocks with 44,550 active cells, is by far the largest and most complex test case on closed-loop optimization [Chen and Oliver, 2010]. 104 realizations of geological models, including net-to-gross ratio, porosity and x-, y- and z-directional permeabilities, were provided by TNO as a representation of the prior knowledge. In this work, these realizations are simulated for 1 year to obtain the dynamic state. In the TEnKF, the ensemble mean is used to produce the undeformed image, and the 104 realizations are utilized to produce the deformed images. Figure 7.22e shows the water saturation at layer 9 from the 183 ensemble mean. Figure 7.22f shows the water saturation at layer 9 from one of the realizations. We can see that the geometric shapes of saturation are extremely complicated, with a large amount of irregular bends and breaks, and the images are very different, especially in the area marked by the red circle. These issues make it difficult to obtain the displacement field. Hence, the TEnKF may not work in such a case, at least not in a straightforward manner. Figure 7.22: Limitation of the TEnKF and possible solution: (a) water saturation at layer 15 in one realization of the SPE9 case; (b) water saturation at layer 15 in another realization of the SPE9 case; (c) aperture problem in motion analysis where we cannot determine the true displacement; (d) possible solution to the aperture problem in the SPE9 case by focusing on the water flooding in the x-direction; (e) water saturation at layer 9 in one realization of the Brugge case; and (f) water saturation at layer 9 in another realization of the Brugge case. Since the geometric shape of saturation is complicated in (e) and (f), it is difficult to obtain the displacement field. 184 In general, the forward transform from saturation to displacement may fail, if the topology of the saturation shape changes dramatically, due to following three main categories: (1) heterogeneity of porous media; (2) property difference of fluids; and (3) fluid separation/mixture. The first kind refers to significant changes of characteristics (e.g., permeability or porosity) of the porous media, which is illustrated in the Brugge case. The most common example in the second kind is viscous fingering, which means the formation of a unstable interface between two fluids that occurs when a less viscous fluid is injected displacing a more viscous one. The third kind is often related to boundary conditions or injection well locations. For example, if two water injection wells are very close (albeit uneconomical in practice), the initially two water zones will merge into one at a later time. From a stochastic point of view in random parametric domain, the number of water zones are either one or two in different realizations at a certain time. We think in above three categories of cases, the displacement fields are either difficult to obtain or not similar to each other in different realizations, thus the TEnKF is no longer valid. Note that the TEnKF may use more state variable than does the EnKF, since the displacement is a vector compared to the saturation that is a scalar. The increase of state variable numbers will have a great impact on the computation of the Kalman gain, especially for a large practical field model in 3D space. This will increases not only memory size but also computational burden during matrix multiplication. However, the displacement vector is only calculated on the locations where the saturation values are non-trivial, as stated in previous 1D example (other locations are not required since the saturation values are trivial, e.g., connate water saturation, and can be obtained easily 185 from interpolation in backward transform). The number of such locations are sometimes much smaller than the total number of gridblocks, especially in early times when the amount of injected water is small. We also point out that the sensitivity analysis for the parameter determination (i.e., template size and search range) will affect the computational burden as well. In this study, the sensitivity analysis are carried out in each assimilation step, while the updating process is waited before the optimal values are determined. To speed up the process, we may only perform the sensitivity analysis in the first assimilation step and use the same parameters in all steps, based on the observation that the accuracy of motion analysis is not very sensitive to the parameters as long as the latter ones are within reasonable ranges (as described in the 2D case). In the above test cases, the observations are collected at the wells as bottom-hole pressure and injection/production rates. If the saturation is also observed at the well, we can treat it as similar to other observations. It is possible that the observed saturation at the well is not honored in the TEnKF, since the TEnKF updates the displacement instead of the saturation. However, the REnKF has the same problem, because the REnKF updates the static state and obtains the dynamic state, including saturation, through a restart process. In general, this problem may be addressed using iteration to ensure that the observed data are honored, which requires further investigation. Actually, we do not recommend including the saturation in the observation because the point saturation (in a fixed location) is usually bimodal/non-Gaussian and nonlinear to the state variables (e.g., log permeability or porosity). 186 The method that we proposed does not necessarily require the DIC process in the forward transform stage. Any digital image processing techniques for motion analysis may be used in our approach, in which a transform of the state occurs before the update/analysis stage. In addition, the TEnKF may be applied to other physical models for sequential data assimilation. For example, in weather forecasting, the state in next time step is mainly affected by the state in current time step, whereas the input parameter is not as significant as in the multiphase flow model. In such cases, the accurate estimation of the state is more crucial to the successful history matching and prediction. Therefore, the TEnKF would be more meaningful in these problems. It is worth noting from a general point of view that, the idea of transform is not limited to state, but can be related to input parameter as well. Actually, some previous literatures have presented some work in this area, for example, the level-set method [Chang et al., 2010] and the discrete cosine transform approach [Jafarpour and McLaughlin]. 187 Chapter 8 Conclusion and Future Work 8.1 Conclusion For uncertainty quantification, this work presents three new methods, termed as location-based transformed probabilistic collocation method (xTPCM), displacement-based transformed probabilistic collocation method (dTPCM), and time-based transformed probabilistic collocation method (tTPCM) by adding a transform process into the probabilistic collocation method (PCM). In uncertainty quantification, the traditional PCM directly approximates the model responses by globally smooth polynomial basis functions. Therefore, it may produce unphysical oscillatory realizations as well as inaccurate estimation of the distributions and statistical moments in strongly nonlinear cases (e.g., discontinuity or large gradient), especially for advection dominated problems. While our new TPCM approaches address such a drawback via transform between the model response and the location. The traditional PCM and proposed TPCM have been tested in various examples of multiphase flow and transport in porous media, with different input variabilities and hydrodynamic dispersivities, and the Monte Carlo (MC) method serves as a benchmark. In summary, this work leads to the following conclusions: 188 (1) Through the geophysical examples, it is shown that with a small number of model runs, the TPCM can achieve a good agreement with the MC results from a large number of realizations. And the proposed TPCM clearly outperforms the PCM in accuracy for estimating the statistical moments and the probability density functions (PDFs) of the model outputs, as well as generating reasonable realizations without Gibbs oscillation. (2) Just like the traditional PCM, the proposed TPCM derives an uncoupled system, which is non-intrusive and highly parallelizable, hence it can be easily implemented with existing simulators. Furthermore, the TPCM could use the same collocation points as in the PCM without additional model runs, thus it is convenient to reanalyze and improve the existing PCM results by the TPCM with the transform process. (3) The traditional PCM provides acceptable results with only small input variability for advection dominated problems, but fails fast as the variability grows, while the TPCM is robust and effective from small to large variability. What’s more, the advantages of TPCM include the ability to handle discontinuous problem with low or even vanishing dispersivity (infinite Peclet number), and the flexibility to cope with complex and moving boundary problems or motions in large domains. (4) It is well known that the higher the retained random dimension, the better the results in collocation method. However, it is not the case when the problem is strongly nonlinear. Based on our numerical experiments, the results from the PCM could be even worse as the random dimension grows, due to the difficulty of fitting or interpolating the nonlinear response surface by polynomials in high dimensions. On the other hand, the TPCM error keeps decreasing, indicating a good convergence. 189 (5) It is natural to think that increasing the level of sparse grid approximation will improve the result in PCM. However, when the problem is high dimensional and strongly nonlinear, the result could be even worse as the level increases, according to our experience. (6) In this research, the TPCM only demonstrates its validity through the multiphase flow and solute transport examples, but it is expected to be also applicable to a wide range of physical models under nonlinear/discontinuous conditions (e.g., Burgers equation or Euler equation). Furthermore, the idea of transform could be extended to other stochastic methods (e.g., intrusive stochastic Galerkin method), and generalized polynomial chaos for various parametric distributions. For data assimilation, this work also presents a new method, which is referred to as the transformed ensemble Kalman filter (TEnKF), by adding a transform process in the traditional ensemble Kalman filter (EnKF). In inverse modeling, the EnKF directly updates the dynamic state along with the static state (as an augmented/joint state vector). Therefore, it may produce non-physical oscillatory realizations and inconsistent estimations (e.g., negative saturation values) in strongly nonlinear cases (e.g., discontinuity or large gradient), especially for advection-dominated problems. The restarted EnKF (REnKF), which updates the static state and reruns the forward model from the initial time to obtain the dynamic state, can generate physical values, but may require a large amount of time. Our new TEnKF addresses these problems via a transform between state and displacement. The case studies show that the TEnKF provides physical results similar to those from the REnKF, but reduces computational 190 effort. We also notice that the TEnKF may not be suitable for the cases in which the topology of the saturation shape changes dramatically, and its performance regarding feasibility and efficiency in large scale real applications requires more tests and further analysis. In this research, although the TEnKF only demonstrates its validity through the multiphase flow examples, it is expected to be applicable to a wide range of physical models under nonlinear/discontinuous conditions (e.g., solute transport model, Burgers equation, and Euler equation). Furthermore, the transform concept can be extended to other Kalman filter-based methods [Zhang et al., 2007; Zeng at al., 2011] and other data assimilation methods (e.g., Kalman smoother and particle filter). 8.2 Future Work The Bayesian inference provides a convenient framework to solve inverse problem. In this method, the parameters to be identified are treated as random variables. The prior knowledge, the system nonlinearity, and the measurement errors can be directly incorporated in the posterior PDF of the parameters. In most cases, the posterior PDF is known up to a multiplicative constant and difficult to sample directly. The Markov Chain Monte Carlo (MCMC) method is a powerful tool to generate samples from the posterior PDF. These samples provide the information of parameters conditioned to the measurement data. Since it requires a large number (usually in many thousands) of iterations of the solver, the MCMC method is usually regarded as a computationally 191 expensive method. To accelerate the MCMC sampling, instead of solving the original partial differential equations directly, a surrogate system can be constructed based on the interpolation to provide output realizations for given parameters. To cope with the system nonlinearity directly, the PCM based Bayesian method (more specifically, Markov Chain PCM (MCPCM)) has been proposed and shown to be more efficient than the standard MCMC. However, the accuracy and efficiency of the MCPCM depends on the accuracy of the surrogate system. In case of strong nonlinearity, the surrogate polynomial is poor, thus the posterior PDF from the MCPCM is inaccurate. This problem is similar to the nonlinear problem in uncertainty quantification. Therefore, we may introduce a transform process into the MCPCM to address this problem. 192 Bibliography Aanonsen, S. I., G. Nævdal, D. S. Oliver, A. C. Reynolds, and B. Vallès (2009), The ensemble Kalman filter in reservoir engineering ―a review. SPE J. 14, 393–412. Abgrall, R. (2007), A simple, flexible and generic deterministic approach to uncertainty quantifications in non linear problems: application to fluid flow problems, Rapport de recherche INRIA, 00325315. Aggarwal, J. K., and N. Nandhakumar (1998), On the computation of motion from sequences of images-A review, Proc. IEEE, 76, 917–935, doi:10.1109/5.5965. Babuska, I., F. Nobile and R. Tempone (2007), A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal., 45(3), 1005–1034. Barthelmann, V., E. Novak, and K. Ritter (2000), High dimensional polynomial interpolation on sparse grid, Adv. Comput. Math., 12, 273–288. Ballio, F., and A. Guadagnini (2004), Convergence assessment of numerical Monte Carlo simulations in groundwater hydrology. Water Resour. Res., 40, W04603, doi:10.1029/2003WR002876. Barron, J. L., D. J. Fleet, and S. S. Beauchemin (1994), Performance of optical flow techniques, Int. J. Comput. Vision, 12, 43–77. Battiti, R., E. Amaldi, and C. Koch (1991), Computing optical flow across multiple scales: an adaptive coarse-to-fine strategy, Int. J. Comput. Vision, 6, 133–145. Bear, J. (1972), Dynamics of Fluids in Porous Media, New York, Elsevier. Bellin, A., Y. Rubin, and A. Rinaldo (1994), Eulerian-Lagrangian approach for modeling of flow and transport in heterogeneous geological formations, Water Resour. Res., 30(11), 2913–2924, doi:10.1029/94WR01489. Bertino, L., G. Evensen, and H. Wackernagel (2003), Sequential data assimilation techniques in oceanography. International Statistical Review 71(2), 223–241. Bruhn, A., J. Weickert, and C. Schnörr (2005), Lucas/Kanade meets Horn/Schunck: combining local and global optic flow methods. International Journal of Computer Vision 61(3), 211–231. Brooks, R. H., and A. T. Corey (1964), Hydraulic properties of porous media, Hydrol. Pap. 3, Colorado State Univ., Fort Collins. Buckley, S. E., and M. C. Leverett (1942), Mechanism of fluid displacement in sands, Trans. AIME 146, 107–116. 193 Chang, H., and D. Zhang (2009), A comparative study of stochastic collocation methods for flow in spatially correlated random fields, Commun. Comput. Phys. 6, 509–535. Chang, H., D. Zhang, and Z. Lu (2010), History matching of facies distribution with the EnKF and level set parameterization, J. Comput. Phys., 229(20), 8011-8030. Chen, Y., D. S. Oliver, and D. Zhang (2009), Data assimilation for nonlinear problems by ensemble Kalman filter with reparameterization, J. Petrol. Sci. Eng. 66, 1–14. Chen, Y., and D. S. Oliver (2010), Ensemble-based closed-loop optimization applied to Brugge field. SPE Reserv. Evalu. Eng. 13(1), 56–71, doi:10.2118/118926-PA. Ciriello, V., V. Di Federico, M. Riva, F. Cadini, J. De Sanctis, E. Zio, and A. Guadagnini (2013), Polynomial chaos expansion for global sensitivity analysis applied to a model of radionuclide migration in a randomly heterogeneous aquifer. Stoch. Environ. Res. Risk. Assess. 27(4), 945–954. Dai, C., H. Li, and D. Zhang (2013), Efficient and accurate global sensitivity analysis for reservoir simulations using probabilistic collocation method, SPE J., doi:10.2118/167609-PA. Dagan, G. (1982), Stochastic modeling of groundwater flow by unconditional and conditional probabilities: 2. The solute transport, Water Resour. Res. 18(4), 835–848, doi:10.1029/WR018i004p00835. Deb, M. K., I. Babuska, and J. T. Oden (2001), Solution of stochastic partial differential equations using Galerkin finite element techniques, Comput. Methods Appl. Mech. Engrg. 190, 6359–6372. Di Stefano, L., and S. Mattoccia (2003), A sufficient condition based on the cauchy-schwarz inequality for efficient template matching. Paper presented at Image Processing International Conference IEEE. Evensen, G. (1994), Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics. J. of Geophysical Research 99(C5), 10143–10162. Fajraoui, N., F. Ramasomanana, A. Younes, T. A. Mara, P. Ackerer, and A. Guadagnini (2011), Use of global sensitivity analysis and polynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media, Water Resour. Res., 47, W02521, doi:10.1029/2010WR009639. Fiorotto, V., and E. Caroni (2002), Solute concentration statistics in heterogeneous aquifers for finite Peclet values, Transport Porous Med. 48, 331–351. Floris, F. J. T., M. D. Bush, M. Cuypers, F. Roggero, A-R Syversveen (2001), Methods for quantifying the uncertainty of production forecasts. Petroleum Geoscience 7, S87–S96. Foo, J., X. Wan, and G. Karniadakis (2008), The multi-element probabilistic collocation method: error analysis and simulation, J. Comput. Phys. 227, 9572–9595. 194 Formaggia, L., A. Guadagnini, I. Imperiali, V. Lever, G. Porta, M. Riva, A. Scotti, and L. Tamellini (2013), Global sensitivity analysis through polynomial chaos expansion of a basin-scale geochemical compaction model, Comput. Geosci. 17, 25–42, doi:10.1007/s10596-012-9311-5. Genz, A., and B. D. Keister (1996), Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight, J. Comput. Appl. Math. 71, 299–309. Ghanem, R. (1998), Scales of fluctuation and the propagation of uncertainty in random porous media, Water Resour. Res., 34(9), 2123–2136, doi:10.1029/98WR01573. Ghanem, R., and S. Spanos (1991), Stochastic Finite Element: A Spectral Approach, New York, Springer. Goovaerts, P. (1997), Geostatistics for Natural Resources Evaluation, New York, Oxford University Press. Gottlieb, D., and C. Shu (1997), On the Gibbs phenomenon and its resolution, SIAM Rev., 39(4), 644–668. Gu, Y ., and D. S. Oliver (2005), History matching of the PUNQ-S3 reservoir model using the ensemble kalman filter, SPE J. 10(2), 51–65, SPE-89942-PA. doi: 10.2118/89942-PA. Gu, Y., and D. S. Oliver (2006), The ensemble Kalman filter for continuous updating of reservoir simulation models, J. Energy Resour. Technol. 128(1), 79–87. Gu, Y., and D. S. Oliver (2007), An iterative ensemble Kalman filter for multiphase fluid flow data assimilation, SPE J. 12(4), 438–446, SPE-108438-PA, doi: 10.2118/108438-PA. Hild, F., and S. Roux (2006), Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties-a Review, Strain, 42, 69–80. doi:10.1111/j.1475-1305.2006.00258.x. Huang, J., X. Peng, L. Qin, T. Zhu, C. Xiong, Y. Zhang, and J. Fang (2009), Determination of cellular tractions on elastic substrate based on an integral Boussinesq solution, J. Boimech. Eng., 131, 061009, doi:10.1115/1.3118767. Isukapalli, S. S., A. Roy, and P. G. Georgopoulos (1998), Stochastic response surface methods (SRSMs) for uncertainty propagation: application to environmental and biological systems, Risk Anal. 18(3), 351–363. Jafarpour, B., and D. B. McLaughlin (2007), History matching with an ensemble Kalman filter and discrete cosine parameterization, Paper SPE 108761 presented at the SPE Annual Technical Conference and Exhibition, Anaheim, California, USA, 11–14 November, doi:10.2118/108761-MS. 195 Jafarpour, B., and D. B. McLaughlin (2009), Estimating channelized reservoir permeabilities with the ensemble Kalman filter: the importance of ensemble design, SPE J. 14(2), 374–388, SPE-108941-PA, doi:10.2118/108941-PA. Jahne, B. (2005), Digital image processing, Springer-Verlag, Berlin, Heidelberg. Jarman, K. D., and A. M. Tartakovsky (2008), Divergence of solutions to solute transport moment equations, Geophys. Res. Lett., 35, L15401, doi:10.1029/2008GL034495. Kalman, R. E. (1960), A new approach to linear filtering and prediction problems. Journal of Fluids Engineering, 82(1), 35–45. Killough, J. E. (1995) Ninth SPE comparative solution project: a reexamination of black-oil simulation, Paper SPE 29110 presented at the SPE Reservoir Simulation Symposium, San Antonio, Texas, USA, 12–15 February, doi:10.2118/29110-MS. Krymskaya, M. V., R. G. Hanea, M. Verlaan (2009), An iterative ensemble Kalman filter for reservoir engineering applications, Computers & Geosciences 13, 235–244. Le Maître, O., M. Reagan, H. Najm, R. Ghanem, and O. Knio (2002), A stochastic projection method for fluid flow: II. Random process, J. Comput. Phys. 181(1), 9–44. Le Maître, O., O. Knio, H. Najm, and R. Ghanem (2004), Uncertainty propagation using Wiener-Haar expansions, J. Comput. Phys. 197, 28–57. Le Maître, O., and O. Knio (2007), A stochastic particle-mesh scheme for uncertainty propagation in vortical flows, J. Comput. Phys. 226, 645–671. Le Maître, O., and O. Knio (2010), Spectral Methods for Uncertainty Quantification: with Applications to Computational Fluid Dynamics, Springer. Lewis, J. P. (1995), Fast Normalized Cross-Correlation, Vision interface, 10, 120–123. Li, G., and A. C. Reynolds (2009), Iterative ensemble Kalman filters for data assimilation, SPE J. 14(3), 496–505, SPE-109808-PA, doi:10.2118/109808-PA. Liu, N., and D. S. Oliver (2005), Ensemble Kalman filter for automatic history matching of geologic facies, Journal of Petroleum Science and Engineering 47(3), 147–161. Lorentzen, R. J., and G. Naevdal (2011), An iterative ensemble Kalman filter, IEEE Trans. Automat. Contr. 56,1990–1995. doi:10.1109/TAC.2011.2154430. Li, H., and D. Zhang (2007), Probabilistic collocation method for flow in porous media: comparisons with other stochastic methods, Water Resour. Res., 43, W09409, doi:10.1029/2006WR005673. Li, H., and D. Zhang (2009), Efficient and accurate quantification of uncertainty for multiphase flow with probabilistic collocation method, SPE J., 114802. Li, W., Z. Lu, and D. Zhang (2009), Stochastic analysis of unsaturated flow with probabilistic collocation method. Water Resour. Res., 45, W08425, doi:10.1029/2008WR007530. 196 Liao, Q., and D. Zhang (2013), Probabilistic collocation method for strongly nonlinear problems: 1. Transform by location, Water Resour. Res., 49(12), 7911–7928, doi:10.1002/2013WR014055. Liao, Q., and D. Zhang (2014a), Probabilistic collocation method for strongly nonlinear problems: 2. Transform by displacement, Water Resour. Res., under review. Liao, Q., and D. Zhang (2014b), Time-based transformed probabilistic collocation method: Framework and application to relative dispersion problems, in preparation. Liao, Q., and D. Zhang (2014c), Data Assimilation for Strongly Nonlinear Problems by Transformed Ensemble Kalman Filter, SPE J., under review. Lin, G., and A. M. Tartakovsky (2009), An efficient, high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media, Adv. Water Resour. 32(5), 712–722. Liu, G., Z. Lu, and D. Zhang (2007), Stochastic uncertainty analysis for solute transport in randomly heterogeneous media using a Karhunen-Loeve-based moment equation approach, Water Resour. Res., 43, W07427, doi:10.1029/2006WR005193. Mathelin, L., and M. Hussaini (2003), A stochastic collocation algorithm for uncertainty analysis, NASA Tech Rep, NASA/CR–2003–212153. MATLAB Scatteredinterpolant Class (http://www.mathworks.com/help/matlab/ref/scatteredinterpolant-class.html) Muller, F., P. Jenny, and D. W. Meyer (2011), Probabilistic collocation and lagrangian sampling for advective tracer transport in randomly heterogeneous porous media, Adv. Water Resour. 34, 1527–1538. Nævdal, G., T. Mannseth, and E. H. Vefring (2002), Near-well reservoir monitoring through ensemble Kalman filter, Paper SPE 75235 presented at the SPE/DOE Improved Oil Recovery Symposium, Tulsa, 13–17 April, doi:10.2118/75235-MS. Oliver, D. S., A. C. Reynolds, and N. Liu (2008), Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge: University Press. Oliver, D. S., and Y. Chen (2011), Recent progress on reservoir history matching: a review, Computational Geosciences, 15(1), 185–221. Oladyshkin, S., F. P. J. de Barros, and W. Nowak (2012), Global sensitivity analysis: a flexible and efficient framework with an example from stochastic hydrogeology. Adv. Water Resour., 37, 10–22. Pan, W., M. Tatang, G. McRae, and R. Prin (1997), Uncertainty analysis of direct radiative forcing by anthropogenic sulfate aerosols, J. Geophys. Res. 102, D18. Peters, L., R. J. Arts, G. K. Brouwer, C. R. Geel, S. Cullick, R. J. Lorentzen, Y. Chen, K. N. B. Dunlop, F. C. Vossepoel, R. Xu, P. Sarma, A. H. Alhutali, and A. C. Reynolds (2010), Results of the Brugge benchmark study for flooding optimization and history matching, SPE Reserv. Evalu. Eng. 13(3), 391–405. 197 Phale, H. A., and D. S. Oliver (2011), Data assimilation using the constrained ensemble Kalman filter, SPE J. 16(2), 331–342, SPE-125101-PA, doi:10.2118/125101-PA. Poëtte, G, B. Després and D. Lucor (2009), Uncertainty quantification for systems of conservation laws, J. Comput. Phys., 228, 2443–2467, doi:10.1016/j.jcp.2008.12.018. PUNQ-S3 test case, 2014, http://www3.imperial.ac.uk/earthscienceandengineering/research/perm/punq-s3model. Qin, L., J. Huang, C. Xiong, Y. Zhang, and J. Fang (2007), Dynamical stress characterization and energy evaluation of single cardiac myocyte actuating on flexible substrate. Biochem. Biophys. Res. Commun. 360, 352–356. Schöniger, A., W. Nowak, and H. J. Hendricks Franssen (2012), Parameter estimation by ensemble Kalman filters with transformed data: approach and application to hydraulic tomography, Water Resources Research, 48(4), W04502, doi:10.1029/2011WR010462. Schreier, H. W., J. R. Braasch, and M. A. Sutton (2000), Systematic errors in digital image correlation caused by intensity interpolation, Opt. Eng., 39, 2915–2921. SGeMS webpage 2014, http://sgems.sourceforge.net/. Shi, L., J. Yang, D. Zhang, and H. Li (2009), Probabilistic collocation method for unconfined flow in heterogeneous media, J. Hydrol., 365, 4–10. Shi, L., D. Zhang, L. Lin, and J. Yang (2010), A multiscale probabilistic collocation method for subsurface flow in heterogeneous media, Water Resour. Res., 46, W11562, doi:10.1029/2010WR009066. Shu, C. W. (1998), Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, Springer Berlin Heidelberg, 325–432. Song, X., L. Shi, , M. Ye, J. Yang, and I. M. Navon (2014), Numerical comparison of iterative ensemble Kalman filters for unsaturated flow inverse modeling. Vadose Zone Journal, 13(2). Tatang, M., A., W. Pan, R. G. Prinn, and G. J. McRae (1997), An efficient method for parametric uncertainty analysis of numerical geophysical models. J. Geophys. Res. 102(D18), 21925–21931. Thulin, K., G. Li, S. Aanonsen, and A. Reynolds (2007), Estimation of initial fluid contacts by assimilation of production data with EnKF, In: Proceedings of the SPE Annual Technical Conference and Exhibition, Anaheim, CA. 11–14 Nov. 2007, Soc. Pet. Eng., Richardson, TX. SPE 109975, doi:10.2118/109975-MS. Tryoen, J., O. Le Maître, M. Ndjinga, and A. Ern (2010), Intrusive Galerkin methods with upwinding for uncertain nonlinear hyperbolic systems, J. Comput. Phys., 229, 6485–6511, doi:10.1016/j.jcp.2010.05.007. Vega-Riveros, J. F., and K. Jabbour (1989), Review of motion analysis techniques, Proc. IEEE, 136, 397–404. 198 Wan, X., and G. E. Karniadakis (2005), An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. J. Computational Physics 209, 617–642. Wang, Y ., G. Li, and A. Reynolds (2010), Estimation of depths of fluid contacts by history matching using iterative ensemble-Kalman smoothers, SPE J. 15, 509–525. doi:10.2118/119056-PA. Webster, M., M. A. Tatang, and G. J. Mcrae (1996), Application of the probabilistic collocation method for an uncertainty analysis of a simple ocean model, MIT Jt. Program on the Sci. and Policy of Global Change Rep Ser 4, Mass. Inst. of Technol., Cambridge. Wen, X.-H., and W. H. Chen (2006), Real-time reservoir model updating using ensemble Kalman filter with confirming option. SPE J. 11(4), 431–442. SPE-92991-PA. doi: 10.2118/92991-PA. Wiener, N. (1938), The homogeneous chaos, Am. J. Math. 60, 897–936. Xiu, D. (2007), Efficient collocational approach for parametric uncertainty analysis, Commun. Comput. Phys 2(2), 293–309. Xiu, D. (2010), Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton Univ. Xiu, D., and G. E. Karniadakis (2002), Modeling uncertainty in steady state diffusion problems via generalized polynomial chaos, Comput. Methods Appl. Mech. Eng. 191(43), 4927–4948. Xiu, D., and J. S. Hesthaven (2005), High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput., 27(3), 1118-1139. Zafari, M., and A. Reynolds (2007), Assessing the uncertainty in reservoir description and performance predictions with the ensemble Kalman filter. SPE J. 12, 382–391. doi:10.2118/95750-PA. Zeng, L., L. Shi, D. Zhang, and L. Wu (2012), A sparse grid based Bayesian method for contaminant source identification, Adv. Water Resour., 37, 1–9, doi:10.1016/j.advwatres.2011.09.011. Zhang, D. (2002), Stochastic Methods for Flow in Porous Media: Coping with Uncertainties, Academic, San Diego, Calif. Zhang, D., L. Shi, H. Chang, and J. Yang (2010), A comparative study of numerical approaches to risk assessment, Stoch. Environ. Res. Risk. Assess. 24, 971–984. Zhang, G., D. Lu, M. Ye, M. Gunzburger, and C. Webster (2013), An adaptive sparse-grid high-order stochastic collocation method for Bayesian inference in groundwater reactive transport modeling, Water Resour. Res., 49, doi:10.1002/wrcr.20467. 199 Zhou, H., J. J. Gómez-Hernández, H. J. Hendricks Franssen, and L. Li (2011), An approach to handling non-gaussianity of parameters and state variables in ensemble Kalman filtering. Advances in Water Resources, 34(7), 844-864. 200 List of Publications Published Liao, Q., and D. Zhang (2013), Probabilistic collocation method for strongly nonlinear problems: 1. Transform by location, Water Resour. Res., 49(12), 7911–7928. Under review Liao, Q., and D. Zhang (2014), Probabilistic collocation method for strongly nonlinear problems: 2. Transform by displacement, Water Resour. Res., under review. Liao, Q., and D. Zhang (2014), Data Assimilation for Strongly Nonlinear Problems by Transformed Ensemble Kalman Filter, SPE J., under review. Liao, Q., and D. Zhang (2014), Constrained probabilistic collocation method for uncertainty quantification of geophysical models, Computat. Geosci., under review. In preparation Liao, Q., and D. Zhang (2014), Time-based transformed probabilistic collocation method: Framework and application to relative dispersion problems, in preparation.

Abstract (if available)

Abstract Based on the polynomial approximation, the traditional probabilistic collocation method (PCM) approximates the model output response, which is a function of the random input parameter, from the Eulerian point of view in specific location. In some cases especially when the advection dominates, the model response has a strongly nonlinear profile with discontinuous shock or large gradient. This nonlinearity in space domain is then translated to the nonlinearity in random parametric domain, which causes oscillation and inaccurate estimation using the traditional PCM. ❧ To address this issue, a new transformed PCM (TPCM) is developed in this study, where the model response is now represented by the an alternative variable (e.g., the location of particular response value, or displacement vector, or arrival time of particular response value), which is relatively linear to the random parameter with a smooth profile. This alternative variable is then approximated by the polynomial, from which a sufficient number of location samples are randomly generated and transformed back to obtain the response samples. Finally, the statistical moments and probability density functions of model response are estimated from these samples. ❧ The advantage of the TPCM is demonstrated through applications to multiphase flow and solute transport in porous media. It is shown that the TPCM achieves higher accuracy for the statistical properties than does the PCM, and produces more reasonable realizations without oscillation, while the computational effort is greatly reduced compared to the direct sampling Monte Carlo method. ❧ The ensemble Kalman filter (EnKF) has been widely used for data assimilation. It is challenging, however, when the relation of state and observation is strongly nonlinear. For example, near the flooding front in an immiscible flow, directly updating the saturation using the EnKF may lead to non-physical results. One possible solution, which may be referred to as the restarted EnKF (REnKF), is to update the static state (e.g., permeability and porosity) and rerun the forward model from the initial time to obtain the updated dynamic state (e.g., pressure and saturation). However, it may become time consuming, especially when the number of assimilation steps is large. ❧ In this study, we develop a transformed EnKF (TEnKF), in which the state is represented by displacement as an alternative variable. The displacement is first transformed from the forecasted state, then updated, and finally transformed back to obtain the updated state. Since the relation between displacement and observation is relatively linear, this new method provides a physically meaningful updated state without re-solving the forward model. The TEnKF is tested in history matching of multiphase flow in a one-dimensional homogeneous medium, a two-dimensional heterogeneous reservoir, and a three-dimensional PUNQ-S3 model. The case studies show that the TEnKF produces physical results without the oscillation problem that occurs in the traditional EnKF, while the computational effort is reduced compared to the REnKF.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

Deep learning architectures for characterization and forecasting of fluid flow in subsurface systems

PDF

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

PDF

Modeling and simulation of complex recovery processes

PDF

Robust real-time algorithms for processing data from oil and gas facilities

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation

Asset Metadata

Creator Liao, Qinzhuo (author)

Core Title Uncertainty quantification and data assimilation via transform process for strongly nonlinear problems

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Defense Date 07/28/2014

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

Tag collocation method,data assimilation,ensemble Kalman filter,OAI-PMH Harvest,strongly nonlinear problems,transform process,uncertainty quantification

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Zhang, Dongxiao (committee chair), Ershaghi, Iraj (committee member), Ghanem, Roger G. (committee member), Jafarpour, Behnam (committee member)

Creator Email liaoqz@gmail.com,qliao@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-469693

Unique identifier UC11287138

Identifier etd-LiaoQinzhu-2882.pdf (filename),usctheses-c3-469693 (legacy record id)

Legacy Identifier etd-LiaoQinzhu-2882.pdf

Dmrecord 469693

Document Type Dissertation

Format application/pdf (imt)

Rights Liao, Qinzhuo

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