University of Southern California Dissertations and Theses

Real-time reservoir characterization and optimization during immiscible displacement processes

(USC Thesis Other)

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content REAL-TIME RESERVOIR CHARACTERIZATION AND OPTIMIZATION DURING IMMISCIBLE DISPLACEMENT PROCESSES by Nelia Jafroodi 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 2009 Copyright 2009 Nelia Jafroodi Epigraph The splendor of life forever lies in wait about each one of us in all its fullness, but veiled from view, deep down, invisible, far off. It is there, though; not hostile, not reluctant, not deaf. If you summon it by the right word, by its right name, it will come. - Franz Kafka, Diaries ii Dedication To my mother and my father iii Acknowledgements I would like to thank most deeply, my parents, Fereshteh and Bahram, who have provided for and supported me generously far more than is typical of parents. I would also like to thank my brother Nima who has always encouraged me as long as I remember him. I wish them the best. I would like to thank the professors in my Ph.D guidance committee Drs. Iraj Ershaghi, Dongxiao Zhang, Hossein Hashemi, Kristian Jessen, and Katherine Shing. I would like to especially thank Professors Ershaghi, Zhang, and Hashemi for super- vising my Ph.D dissertation. I cannot easily express my gratitude to both Professor Ershaghi and Professor Zhang who advised me closely on my research and educa- tion. I would like to emphasize their role in my scientific development in sharing with me their knowledge and passion for research. I would like to thank Professor Ershaghi for providing for me a suitably free and supporting environment in which I could focus on my research. I am very grateful for his patience, guidance and encouragement. Numerous other professors in the Mork Family Petroleum Engi- neering department at USC have contributed to my education. That also goes for iv the professors and teachers at my former institute, Amir Kabir University of Tech- nology in Iran. I would like to thank them and also apologize that I cannot name them all in this short space. I would like to thank Chevron Cisoft fellowship that supported me for a portion of my studies at USC. Also, I would like to thank Chevron Corporation and Occi- dental Petroleum for providing me opportunities for summer internship that proved instrumental in my exposure to the industrial and real world aspects of petroleum engineering. In particular, I would like to thank Dr Kaveh Dehghani at Chevron who provided me with a chance to work on a project closely related to my research topic. I would like to thank my dearest and most inspiring friends: Farzaneh Vaziri who inspired in me a sense for culture and knowledge from early on and Farman Behboud who remains my most inspiring and at the same time, most strict teacher. I would also like to thank Parvin Moradi whose kind spirit is still a great encouragement to me to get out there and get things done. Same goes for all those little friends (all grown up now): Azadeh, Avideh, Amir-Ali, Shahab, and Arash. I wish them all the best of luck. I would also like to thank Faramarz Vaziri and Dr Behrooz Nahid who intro- duced me to the Petroleum Engineering program at USC which turned out to be a great place for me to learn and research. I would also like to thank Mali who generously took me in and supported me during my early months in Los Angeles. v Many thanks also to my fellow student friends at the Petroleum Department. They have all played a substantial role in my life at USC and helped me figure out so many things: Golnaz, Mohsen, Abdollah, Dalad, Reza, Farnaz, Parham, Hamid, Hamed, Tayeb, and Mohammad. Last, but far from least, I want to express my deep appreciation, gratitude, and love for my dear friend Kaveh who helped me tremendously with this dissertation despite my emotional outbursts, streaks of despair, and hours of solitude. vi Table of Contents Epigraph ii Dedication iii Acknowledgements iv List Of Tables x List Of Figures xi Abstract xxi Chapter 1: Introduction 1 1.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Research Objective . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . 6 Chapter 2: Reservoir Characterization Using Pressure Data 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Theory and Definition . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Estimation Approach . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.5.1 Front velocity . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5.2 Uniform Permeability profile . . . . . . . . . . . . . . . 19 2.5.3 Zones of different permeability profile . . . . . . . . . . . 20 2.5.4 An arbitrary Permeability Profile . . . . . . . . . . . . . 23 2.5.5 Non-linear Effects and Hysteresis . . . . . . . . . . . . . 28 2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 vii Chapter 3: Capacitance Resistance Modeling of Reservoirs un- der Waterflood 32 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.1 Capacitance model . . . . . . . . . . . . . . . . . . . . . 35 3.2.1.1 Time-variation of interwell connectivities and time lags . . . . . . . . . . . . . . . . . . . . . 40 3.2.1.2 Oil production model . . . . . . . . . . . . . . . 43 Chapter 4: Ensemble-based Characterization and Optimization 47 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1.1 Ensemble Kalman Filter . . . . . . . . . . . . . . . . . . 48 4.1.2 Production Optimization . . . . . . . . . . . . . . . . . . 53 4.2 Mathematical model and Methodology . . . . . . . . . . . . . . 55 4.2.1 Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2.2 Ensemble Kalman filter . . . . . . . . . . . . . . . . . . . 58 4.2.2.1 Ensemble representation of error statistics . . . 59 4.2.2.2 Ensemble representation of measurements . . . 60 4.2.2.3 Analysis equations . . . . . . . . . . . . . . . . 61 4.2.3 Ensemble Based Closed loop Production Optimization (EnOpt) . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Chapter 5: Numeric Studies with Synthetic Reservoirs 67 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.2 EnKF/CRM Setup . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.3 EnOpt/CRM additional Setup . . . . . . . . . . . . . . . . . . . 71 5.4 Introduction to case studies . . . . . . . . . . . . . . . . . . . . 72 5.5 Homogeneous Reservoir . . . . . . . . . . . . . . . . . . . . . . 77 5.5.1 3-Well Model (Case 1) . . . . . . . . . . . . . . . . . . . 77 5.5.2 Five Spot pattern (Case 2) . . . . . . . . . . . . . . . . . 95 5.6 Heterogeneous Reservoirs . . . . . . . . . . . . . . . . . . . . . . 108 5.6.1 Channelized Reservoir (Case 3) . . . . . . . . . . . . . . 108 5.6.2 Peripheral Pattern (Case 4) . . . . . . . . . . . . . . . . 123 5.7 Optimization Procedure . . . . . . . . . . . . . . . . . . . . . . 132 5.7.1 Streak Reservoir (Case 5) . . . . . . . . . . . . . . . . . 132 5.7.1.1 Scenario (a): EnKF with CRM framework . . . 133 5.7.1.2 Scenario (b): known geology . . . . . . . . . . . 147 5.7.1.3 Scenario (c): unknown geology . . . . . . . . . 149 5.7.2 Channelized Reservoir (Case 6) . . . . . . . . . . . . . . 158 5.7.2.1 Scenario (a): known geology . . . . . . . . . . . 159 5.7.2.2 Scenario (b): unknown geology . . . . . . . . . 165 viii 5.8 Application to Real field . . . . . . . . . . . . . . . . . . . . . . 173 5.8.1 Small pattern . . . . . . . . . . . . . . . . . . . . . . . . 173 5.8.2 Large Pattern . . . . . . . . . . . . . . . . . . . . . . . 180 Chapter 6: Summary and Conclusion 186 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 6.2 Future Research Proposals . . . . . . . . . . . . . . . . . . . . . 189 Glossary 191 References 194 ix List Of Tables 5.1 Synthetic Case 1 defining parameters. . . . . . . . . . . . . . . 77 5.2 Initial guess for synthetic Case 1 . . . . . . . . . . . . . . . . . 80 5.3 Final estimates (over 100 realizations) of all the unknowns in homogeneous synthetic Case 1 . . . . . . . . . . . . . . . . . . 82 5.4 Parameters a j and b j obtained from fitting the power law model 86 5.5 Final obtained values for all the unknowns in the EnKF-CRM methodology, synthetic Case 2 . . . . . . . . . . . . . . . . . . 102 5.6 Final obtained values for all the unknowns in the EnKF-CRM methodology, synthetic Case 3 . . . . . . . . . . . . . . . . . . 117 5.7 Various settings for synthetic Case 5. . . . . . . . . . . . . . . 134 5.8 Final estimates (over 100 realizations) of all the unknowns in synthetic Case 5; scenario (a) . . . . . . . . . . . . . . . . . . . 144 x List Of Figures 2.1 A one dimensional reservoir with varying lateral permeability profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Saturation profile (blue, solid) before breakthrough. Notice that it is far from a piston-like displacement (black dashed). . . . . 18 2.3 Front velocity vs injection rate . . . . . . . . . . . . . . . . . . 19 2.4 (a) Location of front x f as a function of time, (b),(c), and (d) denote ΔP , f w , and T (average transmissibility) as a function of x f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Error in approximating T eff and T ave for 4 cases: The simple case (blue) refers to a constant rate and a constant permeability profile. More detail is given in the text. . . . . . . . . . . . . . 21 2.6 Pressure drop as a function of the extrapolated front location for heterogeneous (non-uniform permeability) and homogeneous (uniform) cases. The non-uniform permeability is defined as 1, 0.1, and 0.25 resp. over the regions 0<x< 0.3, 0.3<x< 0.7, and 0.7<x< 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.7 Estimated profile from pressure derivative (blue, dashed) vs. ac- tual permeability profile (red, solid). Notice the shift and the asymptotic behavior for the piecewise constant profile. . . . . . 24 2.8 Estimated profile from pressure derivative (blue, dashed) vs ac- tual permeability profile (black, solid). Notice the shift for a smoothly varying permeability profile. . . . . . . . . . . . . . . 24 xi 2.9 The effect of viscosity on the shift between the estimated and ac- tual points of variation in the permeability profile. As expected from a friction factor, this effect is logarithmic in nature. d is the distance from the point of change in the slope of ΔP vs x f from the actual partition boundary as a function of the log( μo μw ) partition location is set at 0.2 and is given by K left = 1 on the left and K right = 0.1 on the right. . . . . . . . . . . . . . . . . 25 2.10 The effect of viscosity on the watercut f w as a function of front location x f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.11 Estimated Profile from Pressure Derivative for the case of vari- able injection rate (blue, dashed) and fixed injection rate (red, solid) vs Actual Permeability Profile (black, solid). . . . . . . . 26 2.12 Observables and variables of the smooth randomized permeabil- ity profile: Part (a) depicts the location of the flood front as a function of time which is expectedly linear (constant veloc- ity). Part (b) depicts the variations in the production pressure. Part (c) depicts the water cut as a function of the tracked front location. Part (d) depicts the underlying permeability profile. . 27 2.13 Estimated and actual inverse permeabilities of a randomized het- erogeneous reservoir. Notice the imprecise results still hold many clues to the actual permeability profile . . . . . . . . . . . . . . 28 2.14 Hysteresis curve depicting the nonlinearity of ΔP vs. q for the homogeneous case . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.1 Synthetic Case 1: Homogenous three well model . . . . . . . . 78 5.2 Injection rate used for synthetic Case 1 . . . . . . . . . . . . . 78 5.3 Observed total production rates for synthetic Case 1 . . . . . . 79 5.4 Observed oil production rates for synthetic Case 1 . . . . . . . 79 5.5 Matching total production rates for synthetic Case 1 . . . . . . 82 5.6 Matching oil production rates for synthetic Case 1 . . . . . . . 83 xii 5.7 Water saturation contour map at 4 different time: (a)t = 0; (b) t = 500 days; (c) t = 1500 days; (d) t = 3000 days, synthetic Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.8 {λ kj } variations over 100 realizations, synthetic Case 1 . . . . . 85 5.9 τ 1 and τ 2 variations over 100 realizations, synthetic Case 1 . . . 85 5.10 a j variations in a logarithmic scale over 100 realizations, syn- thetic Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.11 b j variations over 100 realizations, synthetic Case 1 . . . . . . . 87 5.12 cumulative water injected vs water oil ratio in a logarithmic scale, synthetic Case 1 . . . . . . . . . . . . . . . . . . . . . . . 88 5.13 RMSE variation of total production rates for various ensemble populations, synthetic Case 1 . . . . . . . . . . . . . . . . . . . 88 5.14 Spread variation of total production rates for various ensemble populations, synthetic Case 1 . . . . . . . . . . . . . . . . . . . 91 5.15 RMSE variation of total production rate for various σ k values, synthetic Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.16 Spread variation of total production rate for various σ k values in a logarithmic scale , synthetic Case 1 . . . . . . . . . . . . . 93 5.17 RMSE variation of total production rate for various σ obs values, synthetic Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.18 Spread variation of total production rate for various σ k values, synthetic Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.19 Synthetic Case 2: Homogeneous five spot pattern . . . . . . . . 96 5.20 Injection rates used for synthetic Case 2 . . . . . . . . . . . . . 96 5.21 Observed total production rates for synthetic Case 2 . . . . . . 97 5.22 Observed oil production rates for synthetic Case 2 . . . . . . . 97 xiii 5.23 Matching total production rates for synthetic Case 2 . . . . . . 99 5.24 Matching oil production rates for synthetic Case 2 . . . . . . . 99 5.25 Water saturation contour map at 4 different times: (a) t = 0; (b)t = 500 days; (c)t = 1500 days; (d)t = 3000 days, synthetic Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.26 {λ kj } distribution at 4 different time, synthetic Case 2. The horizontal axis is indexed by (j− 1)×m +k corresponding to {λ kj }. The maximum and minimum value ofλ for each case are depicted with error bars and the mean is depicted with circle marker. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.27 {λ kj } distribution, synthetic Case 2. Each arrow represents the mean estimate of{λ kj } at the final time of the simulation and the length of the arrows are linearly related to{λ kj }. . . . . . . 104 5.28 λ 23 variations over 100 realizations, synthetic Case 2 . . . . . . 105 5.29 λ 53 variations over 100 realizations, synthetic Case 2 . . . . . . 105 5.30 {τ j } distribution at 4 different times, synthetic Case 2. The horizontal axis represents the index j for{τ j }. The maximum and the minimum values of{τ j } are represented with error bars while the mean values are depicted with circle markers. . . . . 106 5.31 τ 1 variations over 100 realizations, synthetic Case 2 . . . . . . . 106 5.32 {λ kj } variations over 100 realizations, synthetic Case 2. . . . . 107 5.33 {τ j } variations over 100 realizations, synthetic Case 2 . . . . . 107 5.34 Synthetic Case 3: Heterogeneous channelized reservoir . . . . . 109 5.35 Well bottom hole pressure data, synthetic Case 3 . . . . . . . . 110 5.36 Observed total production rates for synthetic Case 3 . . . . . . 110 5.37 Observed oil production rates for synthetic Case 3 . . . . . . . 111 xiv 5.38 Matching total production rates for synthetic Case 3. The red curve shows the actual observed data and the black curves show the ensemble realizations. . . . . . . . . . . . . . . . . . . . . 112 5.39 Matching oil production rates for synthetic Case 3 . . . . . . . 113 5.40 Water saturation contour map at 4 different times: (a) t = 0; (b)t = 500 days; (c)t = 1500 days; (d)t = 3000 days, synthetic case 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.41 {λ kj } distribution at 4 different time, synthetic Case 3. The horizontal axis is indexed by (j− 1)×m +k corresponding to {λ kj }. The maximum and minimum value ofλ for each case are depicted with error bars and the mean is depicted with circle marker. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.42 {λ kj } distribution, synthetic Case 3. Each arrow represents the mean estimate of{λ kj } at the final time of the simulation and the length of the arrows are linearly related to{λ kj }. . . . . . . 115 5.43 λ 44 variations over 100 realizations, synthetic Case 3 . . . . . . 117 5.44 λ 34 variations over 100 realizations, synthetic Case 3 . . . . . . 118 5.45 {τ j } distribution at 4 different times, synthetic Case 3. The horizontal axis represents the index j for{τ j }. The maximum and the minimum values of{τ j } are represented with error bars while the mean values are depicted with circle markers. . . . . 118 5.46 τ 3 variations over 100 realizations, synthetic Case 3 . . . . . . . 119 5.47 {J j } distribution at 4 different times, synthetic Case 3. The horizontal axis represents the index j for{J j }. The maximum and the minimum values of{J j } are represented with error bars while the mean values are depicted with circle markers. . . . . 120 5.48 J 1 variations over 100 realizations, synthetic Case 3 . . . . . . . 121 5.49 {λ kj } variations over 100 realizations, synthetic Case 3 . . . . . 121 5.50 {τ j } variations over 100 realizations, synthetic Case 3 . . . . . 122 xv 5.51 {J j } variations over 100 realizations, synthetic Case 3 . . . . . 122 5.52 Synthetic Case 4: Heterogeneous peripheral pattern . . . . . . 124 5.53 Injection rates used for synthetic Case 4 . . . . . . . . . . . . . 125 5.54 Observed total production rates for synthetic Case 4 . . . . . . 125 5.55 Observed oil production rates for synthetic Case 4 . . . . . . . 126 5.56 Matching total production rates for synthetic Case 4 . . . . . . 127 5.57 Matching oil production rates for synthetic Case 4 . . . . . . . 128 5.58 RMSE variation of total production rate for various ensemble populations, synthetic Case 4 . . . . . . . . . . . . . . . . . . . 129 5.59 RMSE variation of total production rate for various σ k values, synthetic Case 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.60 RMSE variation of total production rate for various σ obs values, synthetic Case 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.61 Synthetic Case 5: Heterogenous streak model . . . . . . . . . . 133 5.62 Injection rates used for synthetic Case 5; scenario (a) . . . . . . 135 5.63 Bottom hole pressure used for synthetic Case 5; scenario (a) . . 135 5.64 Observed total production rate for synthetic Case 5; scenario (a) 136 5.65 Observed Oil production rate for synthetic Case 5; scenario (a) 136 5.66 Matching total production rates for synthetic Case 5; scenario (a) 137 5.67 Matching oil production rates for synthetic Case 5; scenario (a) 138 5.68 {λ kj } distribution at 4 different times: (a) t = 0 ; (b) t = 200 days; (c) t = 600 days; (d) t = 1000 days, synthetic Case 5; scenario (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 xvi 5.69 variation ofλ 11 estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. . . . . . . . . . . . . . . . . . . . . . . . 141 5.70 Water saturation contour maps at 4 different times; a) t = 0; b) t = 200 days; c) t = 600 days; d) t = 1000 days, synthetic Case 5; scenario (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.71 variation ofτ 3 estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. . . . . . . . . . . . . . . . . . . . . . . . 143 5.72 variation of{λ kj } estimates over time for synthetic Case 5; sce- nario (a). Black curves denote the ensemble members and red denotes the ensemble average. . . . . . . . . . . . . . . . . . . . 143 5.73 variation of{τ j } estimates over time for synthetic Case 5; sce- nario (a). Black curves denote the ensemble members and red denotes the ensemble average. . . . . . . . . . . . . . . . . . . . 144 5.74 variation of log(a j ) estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. . . . . . . . . . . . . . . . . 145 5.75 variation of{b j } estimates over time for synthetic Case 5; sce- nario (a). Black curves denote the ensemble members and red denotes the ensemble average. . . . . . . . . . . . . . . . . . . . 146 5.76 Injection rate vs. control step vs. iteration step in synthetic Case 5; scenario (b). The thick purple curve depicts the final optimized rate obtained for each injector. . . . . . . . . . . . . 148 5.77 Summary of optimized injection rate in synthetic Case 5; sce- nario (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 5.78 Results of oil and water production compared with the base case for synthetic Case 5; scenario (b). . . . . . . . . . . . . . . . . 150 xvii 5.79 {λ kj } distribution at 4 different data times for synthetic Case 5; scenario (c). The position (j− 1)×m +k on the horizontal axes denotes{λ kj }. The blue markers represent the mean val- ues of{λ kj } from the ensemble, while the standard deviation is represented with error bars. . . . . . . . . . . . . . . . . . . . 152 5.80 {τ j } distribution at 4 different times for synthetic Case 5; sce- nario (c). The position j on the horizontal axes corresponds to {τ j }. The blue markers represent the mean values of{τ j } from the ensemble, while the standard deviation are represented with error bars. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 5.81 Injection rate vs control step vs. iteration index for synthetic Case 5; scenario (c). The thick green curve denotes the final optimized injection rate. . . . . . . . . . . . . . . . . . . . . . . 154 5.82 Summary of optimized injection rate in synthetic Case 5; sce- nario (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 5.83 Comparison of the field production responses for different sce- narios: closed loop optimization (scenario (c), red), optimiza- tion with known geology (scenario (b), black), and the base case (green). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 5.84 Production match from the initial ensemble as a function of time; scenario (c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 5.85 Production match from the final ensemble as a function of time; scenario (c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 5.86 Injection rate vs. control step vs. iteration step in synthetic Case 6; scenario (a). The thick green curve depicts the final optimized rate obtained for each injector. . . . . . . . . . . . . 160 5.87 Summary of optimized injection rates in synthetic Case 6; sce- nario (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 5.88 Results of oil and water production compared with the base case for synthetic Case 6; scenario (a). The black curve shows the results of known geology scenario and the green curve shows the base case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 xviii 5.89 Sensitivity Analysis: Optimized profit vs control variable noise variance in Case 6; scenario (a) . . . . . . . . . . . . . . . . . 164 5.90 Sensitivity Analysis: Optimized profit vs control variable noise temporal correlation in Case 6; scenario (a) . . . . . . . . . . . 164 5.91 {λ kj } distribution at 4 different time, synthetic Case 6; scenario (b). The horizontal axis is indexed by (j− 1)×m +k corre- sponding to{λ kj }. The maximum and minimum value of λ for each case are depicted with error bars and the mean is depicted with circle marker. . . . . . . . . . . . . . . . . . . . . . . . . 166 5.92 {τ j } distribution at 4 different times, synthetic Case 6; scenario (b). The horizontal axis represents the index j for{τ j }. The maximum and the minimum values of{τ j } are represented with error bars while the mean values are depicted with circle markers. 167 5.93 Injection rate vs. control step vs. iteration step in synthetic Case 6; scenario (b). The red green curve depicts the final optimized rate obtained for each injector. . . . . . . . . . . . . . . . . . . 168 5.94 Summary of optimized injection rates in synthetic Case 6; sce- nario (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 5.95 Results of oil and water production compared with the base case for synthetic Case 6; scenario (b). The (black,red) curves show the results of (known geology,unknown) scenarios and the green curve shows the base case. . . . . . . . . . . . . . . . . . . . . 170 5.96 Production match from the initial ensemble as a function of time in Case 6; scenario (b) . . . . . . . . . . . . . . . . . . . . . . . 171 5.97 Production match from the final ensemble as a function of time in Case 6; scenario (b) . . . . . . . . . . . . . . . . . . . . . . . 172 5.98 Small Pattern, Real field . . . . . . . . . . . . . . . . . . . . . 174 5.99 Injection rates, Small Pattern . . . . . . . . . . . . . . . . . . . 174 5.100 Observed total production rates, Small Pattern . . . . . . . . . 175 5.101 Cumulative water injected-Cumulative production, Small Pattern 176 xix 5.102 Obtained match for the total production rates, Small Pattern . 178 5.103 RMSE variation of total production rate, Small Pattern . . . . 178 5.104{λ kj } distribution; Small Pattern. Arrow lengths represent the mean estimate of{λ kj } at the final time with a direction along (k,j). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 5.105 Maximum{λ kj } distribution; Small Pattern. Each arrow rep- resents the length and direction of the maximum{λ kj } at each injector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 5.106 Large pattern, Real field . . . . . . . . . . . . . . . . . . . . . . 180 5.107 Injection rates, Large Pattern . . . . . . . . . . . . . . . . . . . 181 5.108 Observed total production rates, Large Pattern . . . . . . . . . 181 5.109 Cumulative water injected-Cumulative production, Large Pattern 182 5.110 Obtained match for total production rates, Large Pattern . . . 183 5.111 RMSE variation of total production rate, Large Pattern . . . . 184 5.112{λ kj } distribution; Large Pattern. Arrow lengths represent the mean estimate of{λ kj } at the final time with a direction along (k,j). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 5.113 Maximum{λ kj } distribution; Large Pattern. Each arrow repre- sents the length and direction of the maximum{λ kj } at each injector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 xx Abstract Optimizing oil recovery subject to operational constraints requires accurate reser- voir models with high predictive capabilities. Specifically, during the immiscible injection processes such as water flooding or gas injection, changes occurring in the flow directions and proportions can affect the long term recovery factors dramati- cally. Monitoring such changes, as reflected in the rate and pressure measurements, in real time, can provide us with information that can prevent long term losses and help in optimal oil recovery. Such a monitoring platform is tied to a predic- tive model and should incorporates historic and the incoming data as it becomes available. In this dissertation, we used a continuous pressure interference testing to characterize a 1-D reservoir and subsequently obtained the permeability profile and were able to link the fluid front movement to the pressure difference varia- tions across the reservoir. To compensate for the typical scarcity of pressure data as opposed to the abundance of rate data we introduced a novel methodology for efficient characterization, prediction, and optimization of large water-flooding op- erations. This is based on ensemble based closed-loop production optimization xxi (EnOpt) and capacitance resistive model (CRM) as its linear underlying dynami- cal system where the injection rates are the driving force or input signal and the production rates are the dynamical variables or output signals. The production rate data are assimilated in real-time by an ensemble Kalman filter for characteri- zation of the reservoir. Simultaneously, the most up-to-date characterization of the reservoir produces optimal values to set the injection well rates to maximize the net present value of the reservoir. Basing the EnOpt method on CRM as opposed to reservoir simulation is computationally efficient even if limited geological data is available from numerous operating wells. furthermore, nonlinear effects associated with saturation movements (ignored in previous publish works) can be incorporated through evolving reservoir parameters. Synthetic and real field examples are used to demonstrate how EnOpt/CRM can match and predict oil and water production rates to probe successes and limitations of our methodology in terms of the reliabil- ity of characterization, improvement in optimization, and sensitivity to the choice of starting parameters. xxii Chapter 1: Introduction 1.1 Problem Statement Constructing a reservoir geologic model with accurate predictive capability is at the heart of reservoir management studies. Such a geological model can be used for sim- ulation and forecasting the reservoir behavior. An accurate geological description is typically based on estimates of the geological properties of the reservoir. Such estimates are typically approximate due to the lack of information underground. To develop a reliable accuracy to form a practical model of the reservoir, the static, dy- namic, and observed/measured data need to be integrated into the model. Static data for example may include information from conventional wire-line logs, core data, and seismic data. Dynamic data may include information from well test re- sults, tracer surveys, pressure data from sensors, temperature, and flow rate data. The term data assimilation refers to a sequential incorporation of the production data with an initially uncertain parameter model where the model uncertainty is 1 reduced with assimilation of every new piece of observed data. Automated assim- ilation of the observation data into a mathematical model of the reservoir is an increasingly important idea in reservoir management studies. Optimizing the economic efficiency of the reservoir is another important task in reservoir management. In practice, control of the reservoir oil output is constrained by logistic factors and costs. Also, estimating the net present value associated with productions in the future, requires a prediction of the reservoir outputs and ap- praisal of economic factors (e.g. discounting future cash flows by interest rates). Production optimization deals specifically with adjusting the available control pa- rameters of the reservoir within the operational constraints, to maximize the net present value of the reservoir. Inevitably, an optimization based on the future of the reservoir requires a reservoir model with good forecasting accuracy. The focus of this study, will (mostly) be the study of a prescribed technique that simultaneously decreases the uncertainty of a geologic model by assimilation of the incoming data and increases the economic efficiency of the reservoir and thus forms a practical solution to two distinct problems in reservoir management at once. Such a real-time reservoir management forms a closed loop production optimization procedure in which the continuous arrival of observed data from the reservoir is used to correct the geological description or the mathematical model of the reservoir and subsequently to adjust the control parameters to optimize the future production. Because at each point, the model describing the reservoir has 2 minimal uncertainty with respect to the available information, the optimization strategy based on the reservoir model is always up-to-date. Reservoir description studies typically require a detailed modeling scenario where large amounts of data are incorporated into the reservoir models in a form suitable for fluid simulators. Because the reservoir modeling in a fluid dynamical setting requires a large set of parameters, the processing and updating of the in- coming information is heavily time and memory consuming and the final results may still contain uncertain geological information that might even be irrelevant to the practical operations of the reservoir. In this study, we base our data assim- ilation on injection and production data as the most abundant source of data in waterflooding scenarios. We proceed by observing that in waterflooding projects, the reservoir heterogeneities and discontinuities play an important role in any opti- mization scenario for oil output maximization throughout the course of the reservoir life time. Defining a reliable notion of interwell connectivity among the injector- producer pairs can help locate the flow paths and barriers and in general provides an efficient approach to reservoir characterization. This paves the way towards reservoir optimization as the focus is shifted from the underlying geophysics to a practical input-output relationship. In fact as we shall describe in the capaci- tance resistance model, it is possible to approximately reduce the time-evolution of the reservoir to a linear response model in which various injection wells act as input signals and the production from the wells act as the output. This essentially 3 linear model is supplied with parameters describing interwell connectivities, time lags, and pressure-rate coefficients. Such linear models are highly desirable as they can be substituted for elaborate numerical simulation procedure based on a fluid- dynamically “exact” description of the reservoir. A reasoning based on the interwell connectivity is in fact a familiar concept. For example, consider in typical water- flooding scenarios where the injection wells are drilled near the production wells in order to maximize the oil recovery. Even in such cases, where in every subsection of the reservoir, the fluid flow is somewhat localized, the reservoir geologic distribu- tion may still significantly affect the results in a way that the interwell connectivity turn out to be uncorrelated to the well distances. In practice (and in certain phases of the reservoir life-time), well connectivities can even depend on the underlying reservoir geological constituency which are time-dependent. The aforementioned linear modeling of the reservoir, while typical in signal processing setting, may not be the most accurate when the reservoir saturation varies considerably over the time scale of observation. While a fundamental solution is to introduce non-linear equations into the model a practical solution is to allow for the parameters of the linear model to vary with time until they are stabilized. The evaluation of this premise is also a major theme in this thesis. 4 1.2 Research Objective Timely decision making plays an important role in reservoir management studies where evaluation of the reservoir performance in a short period of time is essential. Simple predictive models are typically based on a material balance equation with a simple geometry which is then turned into a linear dynamical system and used to evaluate the reservoir performance. Even with minimal dynamical information from well control variables such as injection or production rates, a fast and accu- rate description and time-propagation of the reservoir is possible. Such minimalistic models can be used prior to comprehensive simulation runs to form a preliminary idea about reservoir directionality, future reservoir performances, or even to kick start an optimization procedure in which the net present value of the reserves is being maximized. In this study, we introduce a new methodology that exploits one such linear parameterization, namely the capacitance resistive model (CRM) to- gether with an ensemble based closed loop production optimization (EnOpt, based on the ensemble Kalman filter, EnKF) in order to simultaneously characterize and optimize the reservoir of interest. We model the reservoir dynamics using CRM for- mulation with the CRM unknown parameters (such as interwell connectivities) ob- tained in turn using EnKF procedure. This methodology can also take into account variations of the CRM parameters in the dynamic updating algorithm (EnOpt) and effectively extends CRM to reservoir conditions where saturation changes are too 5 rapid or in waterflooding projects where the water breakthrough time has not been reached yet. The set of CRM parameters and a phenomenological set of parameters that describe the oil production, form the underlying geological model of the reservoir in our approach. At each time step, an ensemble of such geologic models are updated based on the available data and are immediately used for optimization of the injection rates that maximize the net present value of the reservoir. As described in the problem statement, our solution is a dynamical reservoir characterization and optimization scenario. For the synthetic numerical cases that we consider, we typically end up with a set of realizations (solutions) that lie statistically closer to the “true” reservoir geologic model with a far lower uncertainty. These solutions are further evaluated for predicting the reservoir output without any updating. We also interpret the final estimates of the CRM parameters in terms of coarse-grained reservoir characteristics such as directionality, channels and even try to map the flow barrier based on their time-variation history. In the optimization scenarios, we compare the final optimized well control variables against unoptimized scenarios with a focus on maximizing oil recovery. 1.3 Dissertation Outline • Chapter 1 is the current introduction which describes the motivation for ad- dressing the problem under investigation. 6 • Chapter 2 describes the pressure interference testing methodology where the pressure data is used to characterize the reservoir (as opposed to using rate data in CRM framework). This chapter contains a brief literature survey of interference testing, followed by the methodology and the formulation we have developed. The results of applying this method to a set of synthetic numerical cases are also provided in this chapter. • Chapter 3 describes CRM as means of quantifying the communication be- tween injector producer well pairs. A literature survey is followed by CRM methodology and formulations. While, CRM parameters were assumed to be constant in the previous CRM literature, we explain why CRM parameters may in fact vary. • Chapter 4 describes the ensemble based closed loop optimization technique (EnOpt) for modeling, updating and optimizing reservoir models. This chap- ter begins with a literature review of the Kalman filter(KF), ensemble Kalman filter (EnKF) in general, and ensemble based closed loop optimization tech- niques (EnOpt). The KF, EnKF and the EnOpt methodology and formu- lation are all described in this chapter in some detail in the context of our problem. • Chapter 5 describes how the CRM framework can be used with the EnKF and EnOpt methodology. We also describe how various initially uncertain estimates for the state and control vectors are constructed. We later apply 7 our technique to a set of synthetic cases and real field scenarios. This is where we provide the numerical results of our methodology. • Chapter 6 contains the summary and conclusions of this work. This chapter ends with some recommendations for future work. 8 Chapter 2: Reservoir Characterization Using Pressure Data 2.1 Introduction Single well and multiwell pressure transient tests have traditionally been used to characterize reservoir rock heterogeneities assuming single phase flow. Extensions to multi-phase flow were first addressed by the pioneering work of Perrine [Per56] and Martin [Mar59]. For single well tests, subsequent works by Ramsey [Ram75] and Kazemi et al. [KMJ72] focused on detecting movement of the injection front from single-well tests. Analysis of multi-well tests for characterizing two-phase flow has not received much attention mostly due to the fact that the process of data gathering has been limited to special circumstances of interference analysis when one or more wells can be shut down for observation purposes. In such cases, estimated values are the average values of transmissibility and storativity between active and responding wells at a given point in the life of the field. To extend the capability of multiwell tests in the tracking of the dynamic flood front, it is necessary 9 that we conduct a continuous monitoring of rate and pressure data. Thanks to the recent technological advances in smart wells and the ability for placing various sensitive sensors in the injection and production wells, there is now a potential paradigm for continuous front tracking via interwell tests. From the interpretation point of view, nonetheless, describing multiphase flow systems is complex because of the need to estimate in situ relative permeabil- ity characteristics that define the transport of the multiple phases (See [AL88], [Mar81], [Mar59], and [Rag89]). Relative permeabilities in immiscible systems are non-linear functions of saturations (See [CR56] and [Hup70]). Furthermore, while in the single-phase flow only a single equation is required to obtain the flow dy- namics, in the multiphase flow, multiple correlations are required to describe all phases (See [AS79]and [Mar81]). These equations, when combined with the rela- tive permeability functions, can in some cases become nonlinear and analytically intractable. A particular case of well test analysis on multiple wells is interference testing which produces an estimated permeability profile across the reservoir (See [ER73], [EO84], and [Ram75]). Once the effective directional reservoir permeability is mapped, such test can potentially help in reservoir fluid flow monitoring. In this study, we investigate continuous interference testing in a multi-phase flow (water-oil) for reproducing the permeability profile across a one dimensional reservoir. Essentially phenomenological and intuitive, our approach relates the immiscible flood front location or transmissibility variations across the reservoir to 10 the continued recording of pressures in production wells. In contrast to conventional interference testing, which requires shutting in some active wells and hence affecting the field operations, our approach shows how continuous monitoring of rate and pressure across the system can be used to locate the shock front. Also we show that based on the location of the shock front, a permeability profile can be obtained by monitoring the real-time pressure profile 2.2 Theory and Definition Consider a reservoir with a two-phase flow system such as oil and water as de- picted in Fig. 2.1. This reservoir is characterized by its permeability profile, initial saturation profile and boundary conditions. We assume that several injection and production wells are placed in this reservoir. The wells serve as flow sources and sinks with controllable rates. The injection and production fluid pressures can be monitored in these wells and will be the central observed quantities in our model. We describe the set of control, observed, and unknown variables including injection ratesq i (t) wherei corresponds to thei-th well, the pressure drops ΔP i and the wa- ter cutsf w,i ’s are the observed variables at each producing well. We assume no-flow boundaries for the reservoir. The fluid flow is then influenced by the permeability characteristics of the reservoir. In the general case, the initial saturations of oil and water can also influence the transmissibilities. The dynamics of saturation and 11 pressure are governed by the coupled partial differential equations of conservation of mass and the continuity equation: ∂(φρs α ) ∂t +∇.(ρ α v α ) =q α (2.1) where α =w,o. ρ α is density of phase α and v α is the velocity of phase α and q α is the corresponding source term. s α is the saturation of phaseα andφ is porosity. The velocities are related to the pressure gradients via Darcy’s law: v α = K×k rα μ α (∇P α −ρ α G) (2.2) where P α is the pressure of phase α and K is the matrix permeability and k rα is the relative permeability of individual phases and G is the acceleration due to gravity. In the set of differential equations (2.1,2.2) coefficients such as relative permeabilities depend explicitly on saturation or pressure. In certain cases, Eqs. (2.1,2.2) can be solved or characterized analytically but generally they are solved by numerical simulations. 2.3 Model To demonstrate our method clearly for a fundamental case where the given and unknown variables can be described in a simple manner, we use a one dimensional 12 Figure 2.1: A one dimensional reservoir with varying lateral permeability profile model for our fluid flow monitoring. Consider a two-phase flow condition in a one- dimensional reservoir of length L. The system is at 100% oil saturation at t = 0. The flow is initiated by injecting water at a given rate q at an injection well at x = 0. The rate q is fixed unless explicitly specified. The water flow pushes oil through the reservoir leading to a production well at x =L. We further assume a permeability profileK(x). Permeability is multiplied by the relative permeabilities of water and oil to produce transmissibility for each phase of the flow. We assume oil and water to have viscosity coefficients of μ o and μ w respectively. Relative permeability is here assumed to be related to the saturation of each phase based on the Corey type model. The pressure drop at the production well is defined by subtracting the initial pressure from the production pressure. We neglect the effects of gravity (horizontal flow) and the capillary pressure. In one dimension, the Buckley-Leverett solution [MM36] points to a flood front moving with a velocity given by the following equation: 13 x f = qt Aφ ∂f w ∂s w | sw =s f This model can be used to obtain a time-varying average transmissibility for the system formed by a Harmonic averaging along the various parts of the system. The change in the effective transmissibility due to the rearrangements of the regions of water and oil results in the change of pressure at the production well. Effec- tively, variations in pressure are related to the location of the fluid front and the permeability properties of the reservoir. We used a finite element numerical simulation to model this reservoir. Our code may be used with higher dimensional reservoirs given the earlier assumptions. All parameters and quantities reported in this paper are scaled by appropriate factors (some of which pertain to the numerical discretization) to be dimensionless and/or scale free. 2.4 Estimation Approach Our cases are based on estimating the permeability profile of a reservoir with injec- tion and production rate q and constant monitoring of ΔP . The variations in ΔP inversely relate to the variations in the effective transmissibility T eff defined as: T eff = q ΔP 14 We also define an average transmissibility T ave that we intend to use as a good approximation for T eff : T ave = L R L x=0 [T (x)] −1 dx where T (x) is the transmissibility at point x given by T (x) =K(x) k rw (x) μ w + k ro (x) μ o ! (2.3) and the relative permeabilities k rα are known functions of saturation at point x. The Corey type formulas relate the relative permeabilities to the saturations of the phases at each point: k rw (x) = s w (x)−s wi 1−s wi −s oi ! 2 (2.4) k ro (x) = 1− s w (x)−s wi 1−s wi −s oi ! 2 (2.5) where s w (x) is the water saturation at point x and s wi and s oi are initial water and residual oil saturations respectively. We base our estimates on the following assumption which we verify along with our results: T ave tT eff (2.6) 15 The saturation profile at any given time (or at any location of the front) can be used to calculate the relative permeability using a Corey type formula. Using the definition of average transmissibility we can then write: 1 T ave (x) = Z x=L x=0 1 K(x) k rw (x) μ w + k ro (x) μ o ! −1 dx Let us now assume that the saturation profile at any given time is close to a piston- like displacement. In a piston-like displacement, the reservoir has a piece-wise constant saturation profile which is constant on either side of the front at which it has a sharp drop [Fig2.2]. In a piston-like displacement, if the viscosities of oil and water are equal, the average transmissibility remains constant and as such will yield no information about permeability profile of the reservoir. While the saturation profile of the reservoir is different from piston-like in our simulations, we approximate the saturation profile with a linear drop from 1−s oi at x = 0 to s f at x = x f , followed by a constant saturation s wi from x = x f to x = L. The saturations f is the saturation at the front. Using this approximate saturation profile we can estimate the average transmissibility. Based on the relationship between the relative permeabilities and saturations given in Eqs. (2.4,2.5), let us define a phenomenological average transmissibility T phe corresponding to this saturation profile. In analogy to the Buckley-Leverett solution, one can show that the derivative of the inverse ofT phe with respect to the location of the front is given by: 16 d(1/T phe (x)) dx x=x f = (s wi +sro−1) 2 μwμo μw (s f +sor−1) 2 +(s f −s wi ) 2 μo − 1 K(x)| x=x f The expression appearing in the numerator can be simplified in the case where μ o = μ w = μ by using the Welge approach to obtain s f , where the numerator simplifies to 1 and we have: d(1/T phe (x)) dx x=x f = μ K(x)| x=x f (2.7) Eq.(2.7) demonstrates that the derivative of the pressure viewed as a function of the front locationx f , is a good estimate for 1/K. Therefore one may use the slope of the pressure curve to determine the permeability profile. The viscous nature of the flow does not allow it to immediately pick up the variations in the permeability profile and a ”lag” will be present. 2.5 Results In this section we report our results to form a set of tools that can be used to monitor the flood flow and relate the observed variations in the well pressures to the actual characteristics of a one dimensional reservoir. We study the above procedure in zones of uniform permeability, zones of different permeability profiles, varying injection rate with uniform permeability and varying injection rate with zones of different permeability. 17 Figure 2.2: Saturation profile (blue, solid) before breakthrough. Notice that it is far from a piston-like displacement (black dashed). 2.5.1 Front velocity We base our analysis on the tracking of the front locationx f as a linear function of time, i. e. we assume that the front has a constant velocity fixed by the injection rate q(t). This is consistent with the Buckley-Leverett solution as long as the saturation at the front is constant. We confirm this assumption in various settings with various permeability profiles and initial saturations. The location of the front in the simulation is obtained by numerically locating the position of the shock front in the saturation profile at the time of break-through and as a function of time. The position of the front is observed to be moving with a constant speed in various settings of the simulation such as permeability profiles, and constant injection rates. Furthermore, the velocity of the front is a linear function of q(t), even when q(t) 18 Figure 2.3: Front velocity vs injection rate is allowed to vary with time during the simulation. Fig. 2.3 depicts the linear relationship between v f , the front velocity with q for a case of homogeneous one dimensional reservoir with various values of fixed injection rateq. The fact that the speed of the front location is constant enables us to locate the position of the front based on the observed initial injection and final breakthrough times in a simple linear way. Once this is known, the front location can be used as a basic parameter against which the observed quantities are tabulated. 2.5.2 Uniform Permeability profile In this section we monitor the pressure drop at the production well in relations to the flood front movement when there is a uniform permeability profile. The results are displayed in Fig. 2.4. The movement of the front location is a linear function 19 of time and there is a sharp rise in the water cut as it reaches the production. This enables us to tabulate the observed values of pressure as a function of x f . Varying the initial oil and water saturation has no effect on the approximating effective transmissibility with average transmissibility: Furthermore, even a varying injection rate or variations in the permeability profile preserve the approximate equality of the average and effective transmissibilities. These results are displayed in Fig. 2.5. In this figure we plot the relative error in the above approximation Eq. (2.6) as a function of the time along the simulation. The zoned case refers to a constant rate and a zoned profiled where permeability is given by k left , k mid , k right in the zones given by x < 0.3, 0.3 < x < 0.7, and x > 0.7 respectively. The varying injection case (blue) refers to a case where the injection rate varies linearly from 0 to a maximum halfway through the simulation and then drops back to zero with the same rate. We observe that even at the sharp spikes due to the sudden changes in the permeability profile or injection rate, the relative error in approximatingT eff by T ave is at most 1.3%. 2.5.3 Zones of different permeability profile We pay special attention to the case where we define a permeability profile for the reservoir by including zones including zones of constant permeability. An example of such case is given in Fig. 2.6, where variations in pressure are observed as a function of the location of the front. These variations clearly coincide with the 20 Figure 2.4: (a) Location of front x f as a function of time, (b),(c), and (d) denote ΔP, f w , and T (average transmissibility) as a function of x f Figure 2.5: Error in approximating T eff and T ave for 4 cases: The simple case (blue) refers to a constant rate and a constant permeability profile. More detail is given in the text. 21 variations of the permeability profile. Furthermore, as the production phase kicks in, the dynamics becomes markedly different and we observe further discrepancies between the estimate permeability and the actual permeability. This is a less severe drawback for our estimates if the variations in the permeability profile are not too sharp. We have reported these estimates in Figures 2.7 and 2.8. As an impediment to flow, viscosity affects the lag in a logarithmic manner as observed in Fig. 2.9. Similarly the water cut is significantly affected by the viscosity factors. In our analysis, we have assumed that the water cut picks up only when the front is close to the producing well. Fig. 2.10 shows that this assumption breaks down when the viscosity ratio of oil and water increases. The estimation of the permeability profile in cases of high oil viscosity is difficult as the variations in pressure are lost due to dissipative effects. We also considered a case where the injection rate q is allowed to vary as a function of time during the movement of the front from injection to production. Fig. 2.11 demonstrates the results of estimating the permeability profile in contrast to the case of fixed injection rate. In this case, we assume that the variable injection rate increases linearly from 0 to a maximum injection rate and subsequently decreases at the same rate to 0 by the time the fluid front reaches the production well. We observe that the derivative of pressure, scaled by the injection rate at the corresponding time (transmissibility inverse) indeed approximate the actual permeability profile. This further confirms 22 Figure 2.6: Pressure drop as a function of the extrapolated front lo- cation for heterogeneous (non-uniform permeability) and homogeneous (uniform) cases. The non-uniform permeability is defined as 1, 0.1, and 0.25 resp. over the regions 0<x< 0.3, 0.3<x< 0.7, and 0.7<x< 1. the applicability of our scheme to real-world scenarios where a perfect control over the injection rate is not possible. 2.5.4 An arbitrary Permeability Profile In this section we consider a permeability profile, which is a smooth randomized function K(x) defined as K(x) =e u(x) 23 Figure 2.7: Estimated profile from pressure derivative (blue, dashed) vs. actual permeability profile (red, solid). Notice the shift and the asymptotic behavior for the piecewise constant profile. Figure 2.8: Estimated profile from pressure derivative (blue, dashed) vs actual permeability profile (black, solid). Notice the shift for a smoothly varying permeability profile. 24 Figure 2.9: The effect of viscosity on the shift between the estimated and actual points of variation in the permeability profile. As expected from a friction factor, this effect is logarithmic in nature. d is the distance from the point of change in the slope of ΔP vs x f from the actual partition boundary as a function of the log( μo μw ) partition location is set at 0.2 and is given by K left = 1 on the left and K right = 0.1 on the right. Figure 2.10: The effect of viscosity on the watercut f w as a function of front location x f 25 Figure 2.11: Estimated Profile from Pressure Derivative for the case of variable injection rate (blue, dashed) and fixed injection rate (red, solid) vs Actual Permeability Profile (black, solid). where u(x) is a random number with a Gaussian distribution with mean of zero and standard deviation equal to 1. The exponential function is used to make the distribution positive. We consider a system with the initial oil and water saturations ofs wi =s oi = 0. For simplicity, we consider the case where the viscosities of oils and water are the same and equal to 1. The main observables and variables of this case are included in Fig. 2.12. In Fig. 2.13, we use the slope of the pressure curve as a function of the front to estimate the permeability profile. Note that while sharp fluctuations of permeability make this estimate quite imprecise. The main features of the permeability profile are, however, captured in the estimated permeability. 26 Figure 2.12: Observables and variables of the smooth randomized per- meability profile: Part (a) depicts the location of the flood front as a function of time which is expectedly linear (constant velocity). Part (b) depicts the variations in the production pressure. Part (c) depicts the water cut as a function of the tracked front location. Part (d) depicts the underlying permeability profile. 27 Figure 2.13: Estimated and actual inverse permeabilities of a randomized heterogeneous reservoir. Notice the imprecise results still hold many clues to the actual permeability profile 2.5.5 Non-linear Effects and Hysteresis Most of our simulations results suggest that in the one dimensional case, a simple effective transmissibility can be defined which is equal to the average transmissi- bility defined by harmonic averaging. Nonetheless, the actual relationship between the pressure change and the injection rate (in analogy with voltage and current in electric circuits) can be non-linear as transmissibility depends implicitly on the injection rate via saturation dependent relative permeabilities. To observe such nonlinear effects, we obtain a hysteresis curve for our system by varying the injec- tion rate from 0 to some maximum value and then bringing it back to 0 and record the variations in pressure. For a linear system, a graph of ΔP vs. q could be linear but a non-linear system could display a hysteresis curve. Fig. 2.14 shows that 28 Figure 2.14: Hysteresis curve depicting the nonlinearity of ΔP vs. q for the homogeneous case indeed even the case of a homogeneous permeability profile exhibits nonlinearity and hysteresis. Such nonlinearities might be utilized as a memory effect to further assist us in profiling the reservoir but in general they make the task of production control further difficult even in simple one dimensional. 2.6 Conclusions One dimensional case studies were conducted to monitor the saturation front move- ment with various lateral permeability profiles. We monitor the production well pressure changes and relate that to the movement of the immiscible flood front movement. The derivative of pressure at the production well as a function of the front location (tested in various permeability profiles) is used to obtain the actual 29 permeability profile. We introduce an effective transmissibility in the system to ob- tain a dynamical handle on the flood front movement as the movement of the front causes the effective transmissibility to vary. In order to confirm this scheme, we run a numerical simulation on reservoirs with various permeability profiles, including a random permeability profile. Derivatives of the pressure variations are compared with the permeability profile and indicate that the pressure derivative (inverse ef- fective transmissibility derivatives) can be used to estimate the permeability profile. We support this idea with a phenomenological model of the saturation profile along the reservoir. The major disagreement between our estimate and the actual per- meability profile is characterized with a shift in the location of the permeability variations. This is mainly due to the dissipation and viscous nature of flow result- ing in loss of information. This is confirmed by our sensitivity studies as the shift is logarithmically related to the viscosity factors. Our results are best covered in the case when there is a constant flow along the reservoir. We further confirm that the front velocity is linearly related to the injection rate even in cases of varying injection rate. We further confirm our estimation method when the injection rate is variable, to allow for uncertainties in the injection rate in our estimates. We also consider the linearity of the relationship between the pressure drop and the injection ratio. We also considered the variations of the water cut at the production well with different oil viscosities as a function of the front location. As expected the water cut only picks up around the time of the water breakthrough. The rise in 30 the water cut appears earlier when the viscosity factors are higher and clearly hints at the inadequacy of our phenomenological estimate of the permeability profile at high viscosity ratios. Nonetheless, our simulation results strongly suggest that one dimensional two-phase reservoir with a heterogeneous structure can be linked to the Buckley-Leverett model of one dimensional reservoir with a uniform permeability. This result, along with the connection between the observed pressure data and the permeability profile can enable us to construct an approximate description of the reservoir structure by simply observing the pressure variations, the water cut, and the injection rate. 31 Chapter 3: Capacitance Resistance Modeling of Reservoirs under Waterflood In this chapter, we provide some background on capacitance resistive model (CRM) including a concise literature review. We then provide a brief mathematical descrip- tion and the theoretical methodology for the CRM method. This methodology forms the foundation of the implementations of Chap. 5. 3.1 Introduction Injection and production rates are the most abundant data available for character- izing reservoir heterogeneity and inter-well connectivity. Injection and production rates data obtained from the reservoir can be used to characterize the parameters describing the dynamics of the reservoir and the interwell connectivity. In princi- ple, such information can then be used to optimize oil recovery by varying injection patterns, in-fill drilling, and re-completion of wells. In particular, the effect of injectors on producers in the reservoir has been stud- ied extensively. The stringent computational resource requirements of large scale 32 reservoir simulation highlight the need for simplifying/approximating the relation- ship between production and injection rates. We can identify the work by Heffer et al. [HFMK97], where the Spearman rank correlations were used to obtain an injector-producer relationship. This relationship was subsequently linked to ge- omechanics. Spearman rank correlation coefficient is a non parametric measure of correlation and assesses how well an arbitrary monotonic function could describe the relationship between two variables, without making other assumptions about the particular nature of the relationship between the variables. The work by Re- funjol [Ref96] also used the Spearman rank analysis to determine the reservoir flow trend. Time lags were then used to find extreme coefficient values in order to relate the injectors to their surrounding producers. Sant’Anna Pizarro [SP98] tested the validity of the Spearman rank analysis using numerical simulation and character- ized the advantages and limitations associated with it. Other approaches include the work by Panda and Chopra [PC98] where artificial neural networks were used to identify the injector-producer relationship and the work by Soeriawinata and Kelkar [SK99] that employed a statistical approach to relate each injectors to its surrounding producers and the concept of superposition was introduced to allow for constructive and destructive interferences between injector effects. 33 Albertoni and Lake [AL03] used a linear multivariate regression technique (MLR) to estimate the interwell connectivities (weights) between pairs of injection- production wells. The linear model coefficients obtained from MLR approach out- line a response relationship between each set of injection-production wells and cor- respondingly a signal processing approach was used to introduce a time lag in well pairs. The interwell connectivity coefficients were subsequently interpreted in the work by Gentil [Gen05], in terms of reservoir transmissibility. A comprehensive capacitance resistive model (CRM) was developed by Yousef et al. [YGJL06] to include generic mathematical features of resistance and capacitance. Explicitly, two parameters were quantified for each injector-producer pair: a weight that quantifies the connectivity between wells and a time constant that signifies a degree of stor- age for each well. Sayarpour et al. [SZKL07] used the analytical solutions of the CRM differential equations in a characterization/optimization procedure based on the available production data. More recently, Liang et al. [LWE + 07] also presented a methodology for optimizing the oil production by adjusting the water injection rates using a capacitance model (CRM). Somewhat more mathematically advanced methods include the work by Thiele and Batycky [TB03] where a streamline based work flow was used to obtain the effective flow units in the reservoir. Lee et al. [LON + 09] used a novel data mining method to characterize the flow units between the injections and the production wells using manual perturbations in injection rates. A wavelet based approach was 34 used to design perturbation in the injection rates in order to observe the variations in the production rates. In this work, we use CRM as the underlying dynamics for the EnOpt methodol- ogy (See Chapter 4) to allow for rapid updating and optimization of waterflooding systems. More specifically we use the framework outlined by Yousef et al. [YGJL06] (the reservoir as a system that converts injection signals to production response) but we allow the system to have time-varying parameters to characterize non-linearity. We note that while in the standard CRM description, the parameters are constant, in principle they must be allowed to vary in order to capture the possibly evolving pathways in the reservoir in the course of its dynamics as the underlying physical quantities such as the transmissibility may vary with time. We will now provide a concise mathematical description of the CRM model of injection-production rela- tionship. A more elaborate description may be found in Yousef et al. [YGJL06]. 3.2 Mathematical Model 3.2.1 Capacitance model CRM is a total mass balance equation with compressibility. We assume there are m injectors acting simultaneously on n producers. The governing material balance 35 for the control volume around each producer is given by the following differential equation: q j (t) = m X k=1 λ kj i k (t)−c t V p d¯ p j(t) dt , (3.1) where c t is the total compressibility around the control volume, V p is the pore vol- ume around each producer, ¯ p j is the average pressure in the volume drained by producer j∈{1,...,n}, i k is the injection rate of each injector k∈{1,...,m}, and q j is the total production rate of producer j and{λ kj } defines the connec- tivity between (k,j) pairs of injector k and producer j, given by Sayarpour et al. [SZKL07]: λ kj = q kj (t) i k (t) , (3.2) where q kj is the production at producer j due to injector k. In the absence of injected fluid (i k (t) = 0), the system undergoes pure depletion (See [WL03]). Gentil [Gen05] interprets{λ kj } directly in terms of the transmissibility of the reservoir by λ kj = T kj P n j=1 T kj . (3.3) where T kj is the effective transmissibility between injector k and producer j. The effective transmissibility T kj is a function of the reservoir’s static characteristics but may also depend on the reservoir state through the dependence of say, the permeability at a given point on the saturation at that point. This will effectively 36 require us to allow for{λ kj } to vary over time as we shall describe later. We remark in passing that, in essence, Eq. (3.1) states that at any time, the total production of the producer j is determined by the change in the average pressure around the producer j and the flow share of the producer j from the total injection. This statement is valid under the assumptions of small and fixed compressibility and that no fluid flows in or out other than the injection and production wells. Eq. 3.1 can be interpreted as an interference testing where the rates are used rather than the pressures. In fact we may base Eq. 3.1 more directly on the rates if we use a productivity index model: q j =J j ( ¯ p j −p wf,j ) (3.4) where p wf,j and J j are the flowing bottom-hole pressure and productivity index of the producer j, respectively. Eliminating the average reservoir pressure ¯ p j in Eqs. (3.1) and (3.4) leads to the basic differential equation of the CRM model: q j (t) = m X k=1 λ kj i k (t)−τ j dq j (t) dt −τ j J j dp wf,j (t) dt (3.5) where we have defined the time constant{τ j } of the volume drained by each pro- ducer asτ j = ctV P J j . Notice that we have m×n parameters for{λ kj },n parameters 37 for{τ j }, and n parameters for{J j }. The solution to Eq. (3.5) can be obtained analytically: q j (t) = A z }| { q j (t 0 )e −(t−t 0 )/τ j + B z }| { e − t τ j Z t t 0 e ξ τ j 1 τ j m X k=1 λ kj i k (ξ)dξ − C z }| { e − t τ j Z t t 0 e ξ τ j J j dp wf,j dξ dξ (3.6) wheret 0 is the initial time from which the reservoir is monitored andξ is a dummy integration variable. Eq. 3.6 states that at any time the total production of each producer consists of three parts, A,B, and C as highlighted. Part A is the decaying remnant of the initial production rate. Part B is the accumulated response from the injection signal history. Part C denotes the effect of changing bottom hole pressure on total production rates. Eq. 3.6 states that the total production at time t is a function of primary production and the injection history as well as pressure history between t and t o . Let us assume that we have only one injector and one producer with a fixed bottom hole pressure. This would remove part C from Eq. 3.6. The contribution from the past injection history would then be controlled by the damping factor τ. Effectively for τ < 1 (time unit), the producers output will coincide the injection variations: injection signal will reach the producer instantly. Thus, if there were no dissipation in the reservoir, τ would be small and regardless of the actual distance 38 between injection and production the injection signal would arrive at the production response immediately. For τ > 1 , the injection step change makes a damped impression on the production response and the accumulated effect from previous time steps also contributes to the production response at time t. For τ 1 there is more attenuation and considerable delay in production response. Physically, the parameter τ describes the fluid pressure dissipation factor from injector to the producer. Unbalanced capacitance model (UCM) The formulation just given, technically describes the balanced capacitance model (BCM) in which the field-wide sum of injection rates is approximately equal to the sum of field-wide production rates. In case of any significant difference between the two sums, we need to describe the system by the unbalanced capacitance model (UCM). Yousef et al. [YGJL06] implement UCM by adding a constant rate q c,j to the solution of the differential equation of the CRM model, Eq. 3.6. In this work, we introduce this extra parameter as a (time-dependent) parameter along with the other parameters [{λ kj },{τ j },{J j }]. The formulation for the UCM is thus adjusted to: 39 q j (t) = q j (t 0 )e −(t−t 0 )/τ j +e − t τ j Z t to e ξ τ j 1 τ j m X k=1 λ kj i k (ξ)dξ − e − t τ j Z t to e ξ τ j J j dp wf,j dξ dξ +q c,j (3.7) UCM should be used when the waterfall system is naturally unbalanced or when a part of the field with flows across the boundaries is being analyzed. In fact, in the real field scenarios we end up using Eq. 3.7 [Chapter 5]. 3.2.1.1 Time-variation of interwell connectivities and time lags Previous studies on CRM are based on the premise that the interwell connectivities and time lags are constant with time. In particular, Yousef et al. [YGJL06] used numerical simulation of a homogeneous reservoir to verify that{λ kj } do not vary with varying injection rates. This was followed by Gentil [Gen05] by interpreting the interwell connectivities physically and relating them to the transmissibility regions between each pair of wells as given in Eq. (3.3). The effective transmissibilities T kj between injector k and producer j can be described in terms of the following effective geological parameters: rock permeability K, cross section A, viscosity μ, distance L, and relative permeability K r . One may thus write ([Gen05]): λ kj = ( KrK μL A) kj P n j=1 ( KrK μL A) kj (3.8) 40 CRM models typically assume that: (i) the reservoir is mature, i.e. the reservoir has an almost uniform and slowly varying saturation. This allows us to ignore the dependence of relative permeability on saturation (changes) to be dropped from Eq. (3.8). (ii) the fluid properties (μ o ,μ w ) do not change with pressure (time). If these assumptions are satisfied then the connectivities{λ kj } can be simplified to λ kj = ( 1 L ) kj P n j=1 ( 1 L ) kj (3.9) Eq. (3.9) characterizes a situation in which{λ kj } depend on injection-production pair distances and do not change with time. In contrast, distinct injection rate patterns applied at the early stages of the reservoir production or over very long time scales may change the saturation distribution of the reservoir drastically. By relaxing the assumption (i), we allow the fluid properties, in particular the effective transmissibilities, to vary along with saturation movements and thus with time. Adding new injectors or producers can also affect the saturation movements and subsequently the transmissibilities or even{λ kj } directly. A similar situation occurs for the time lags{τ j } : Yousef et. al [YGJL06] formulates{τ j } in terms of the total compressibility C t , pore volume injected V p , and the productivity index J j : τ j = C t ×V P J j (3.10) 41 The total compressibility C t can be expanded as C t =C o S o +C w S w +C r in which C o , C w , and C r stand for oil, water, and rock compressibilities. In the development of CRM models, the fluid compressibilitiesC w andC o are assumed to be small and constant. This effectively removes the dependence of C t on the fluid properties including the saturations. Furthermore, the productivity index (J j ) and the pore volume V p are assumed to be constant. In reality, C t may depend on sat- uration distribution and subsequently on time. For instance, C o orC w may not be small or a gas phase might contribute significantly to the total compressibility. Sim- ilarly,J j andV p are sensitive to the saturation distribution and adding/elimination of wells. The solution given by Eq. (3.6) for the CRM dynamical system assumes that the parameters{λ kj },{τ j }, and{J j } are constant. In this work, we model the possible variations of the parameters with time through a dynamical system with piecewise constant parameters. In such a case, Eq. (3.6) can be written as: q j (t n ) = q j (t n−1 )e −(tn−t n−1 )/τ j (tn) + e − tn τ j (tn) Z tn t n−1 e ξ τ j (ξ) 1 τ j (ξ) m X k=1 λ kj (ξ)i k (ξ)dξ − e − tn τ j (tn) Z tn t n−1 e ξ τ j (ξ) J j (ξ) dp wf,j dξ dξ (3.11) 42 whereτ j (t n ) andλ kj (t n ) are assumed to be constant in each time interval [t n−1 ,t n ]. This solution leads to a time-discretization in which in every discrete time segment the parameters are assumed to be constant. This extends the applicability of CRM to immature reservoirs and provides more quantitative information on the reservoir. As we shall see this is provided by integration of the observed data in the real time from the EnKF procedure (Chap. 4) that provides the best current match consistent with the geologic/fluid properties of the reservoir progressively in time. The unbalanced capacitance model (UCM) is similarly modified as follows: q j (t n ) = q j (t n−1 )e −(tn−t n−1 )/τ j (tn) + e − tn τ j (tn) Z tn t n−1 e ξ τ j (ξ) 1 τ j (ξ) m X k=1 λ kj (ξ)i k (ξ)dξ − e − tn τ j (tn) Z tn t n−1 e ξ τ j (ξ) J j (ξ) dp wf,j dξ dξ +q c,j (t n ) (3.12) 3.2.1.2 Oil production model The capacitance model describes the total fluid dynamics between injection and production wells through the reservoir and does not differentiate between oil and water as the constituents of the fluid. Evaluation of oil production as a function of time is an important factor in oil rate prediction. Oil cut (ratio of oil to total production rate) varies during the production life time according to the changes in the recovery process. In the initial stages of the recovery process, when the field 43 production rates mainly consist of oil production, the oil cut is close to 1.0 and nearly constant. When secondary recovery process begins, E. g. waterflooding, the oil saturation decreases and can reach the level of residual oil saturation (S orw ) and brings the oil cut monotonically to a minimum. Immiscible oil fractional flow models are either based on saturation front movements or empirical fractional flow models. Buckley-Leverett based fractional flow models can for example be used to obtain the dependence of the oil cut on saturation distribution variations. Knowledge of the saturation variations as a function of time, residual water and oil saturations, and relative permeability function parameters are necessary in using this method which make it inherently difficult to embed such a procedure with the CRM framework esp. for multiwell systems. In contrast to such analytic formulations, we use an empirical oil fractional flow model that requires less computational effort and can be incorporated into the CRM framework more handily. Various empirical models are used to obtain a relationship between oil produc- tion at the wells and various fluid properties [See [CH59], [Tim71], [EO78]]. These fractional flow models are either a function of cumulative oil production or average reservoir saturations which are both unknown during the forward estimation. In this study, we use a power law model (following Gentil [Gen05] and Liang et al. [LWE + 07]) to describe the relationship between the cumulative water injected toward the j-th production well (W j ) and the ratio of water to oil production q w,j /q o,j =F wo,j : 44 F wo,j (t) =a j W j (t n ) b j (3.13) The total water injected (W j ) is given by summing over all the water injected from all the injectors and thus defined as: W j (t n ) = Z tn t 0 m X k=1 λ kj (ξ)i k (ξ)dξ (3.14) Combining with Eq. (3.11) for the total fluid rate q j (t), Eq. (3.13) can be used to obtain the oil production rate (using q w,j +q o,j =q j ): q o,j (t n ) = q j (t n ) 1 +a j W j (t n ) b j (3.15) The parametersa j andb j (another set of 2n parameters) may be obtained by fitting to the power law model in Eq. (3.13). In fact the linear section of water oil ratio vs cumulative water injected in a logarithmic scale after 50% watercut, turns out to be the most suitable point for obtaining these parameters and is discussed in details in [EA84]. In our setting, the variable W j (t) implicitly depends on a priori unknown pa- rameters{λ k,j }. We will thus incorporate Eqs. (3.14) and (3.15) along with Eq. (3.11) to define the dynamical system describing the oil rates of the wells along with the total production rates. This becomes especially useful when the CRM model is 45 used to optimize the net present value of the reservoir that differentiates between oil and water production rates. 46 Chapter 4: Ensemble-based Characterization and Optimization In this chapter, we first provide a brief literature review and background on ensem- ble Kalman filter (EnKF) methodology with a focus on reservoir characterization studies. We then focus on ensemble based closed-loop optimization (EnOpt). We provide the mathematical terminology and a basic description of Kalman filter (KF), EnKF and EnOpt procedures in the context of reservoir characterization. 4.1 Introduction History matching in reservoir simulation concerns adjusting geological parameters such as permeabilities, porosities, and saturations in a way that the difference between the calculated and actual measurements is minimized. History matching is now automated to a large extent and can be performed during a multi-phase flow simulation by applying on-the-fly closed-loop numerical optimization. This automated process is our main focus in this chapter and throughout this thesis in general. 47 4.1.1 Ensemble Kalman Filter Kalman filter (KF) is a mathematical parameter estimation model that can be used for updating the state model of a linear dynamical system by assimilating the observation data (as an output of the dynamical system)[Gel02]. KF procedure is most efficient when a small number of variables (and parameters) is considered. The fact that KF assumes a linear transient matrix (from the linear dynamical system) to update the error covariance matrix prevents direct application of KF to nonlinear systems. Extended Kalman filter (EKF) is a generalized form of the Kalman filter in which the non-linearities of the dynamical system describing the evolution are quantified perturbatively by linearizing the differential equation by using the multi- variable high order derivatives. This can improve the applicability of the Kalman filter to nonlinear problems but may fail when the number of degrees of freedom becomes large since the number of parameters grows quickly with the degree of nonlinearity of the dynamical system. Some of the redundancy introduced into the statistical description of the parameters can be removed by finding an adaptive effective state space. An analysis of the EKF procedure with a large number of parameters appears in the work by Graham [Gra02]. The ensemble Kalman filter (EnKF) is another generalization of the original Kalman filter that can be used with nonlinear systems [VH01]. In essence, EnKF is a Monte Carlo method in which an ensemble of reservoir state models is generated 48 initially and evolves in time. At each estimation step, the parameters in each ensemble member model are updated by comparison with the observation data in such a way that the statistical properties of the ensemble match the uncertainty at that step and the conditions of the observables. The overall evolution of the ensemble effectively replaces the analytic procedure for updating correlation matrix in the original KF formulation. EnKF procedure typically starts with an initial ensemble sampled from some prior probability density function (PDF) of the model state. Each member of the ensemble contains a model state vector with the static (parameters) and dynamic variables along with the observation values. For example in petroleum reservoir problems, the constant values of permeability and porosity can be static variables whereas the values of saturation and pressure are dynamic variables. Observation values can refer to any value obtained from the field at a given (real) time and as such they act as controllers or probability conditioners. These controllers could be well rates, bottom hole pressures, water-oil ratios, or gas-oil ratios. The observation values are kept in the state vector in order to relate the model state to the observed values from the reservoir directly. The ensemble is propagated forwarded in time (forecast step) at every time step of a “simulation” of the dynamical system and is updated (analysis step) each time a new observation data becomes available. Static variables do not change with time in the forecasting step but get updated in the analysis step after each data assimilation. The updating component of EnKF involves only the first two statistical moments of the ensemble 49 and thus ignores any departures from the Gaussian distribution in the model state variables. EnKF was applied by Evenson [Eve03] to a fluid dynamical problem with a high degree of nonlinearity and a large number of parameters. Similarly, the flow equa- tions describing petroleum reservoirs, are typically non-linear partial differential equations and involve a large number of parameters (i.e. proportional to the geo- metric size of the reservoir or the simulation model approximating it). In the work by Naevdal and Mannseth [NMV02], EnKF was indeed used for updating the static parameters near well reservoir models with promising results. This was followed in Naevdal and Johnsen, [NJAV05] where a constructive model for permeability was presented and the EnKF procedure was shown to successfully reduce reservoir uncertainties resulting in a successful history matching. Gu and Oliver [GO05] used EnKF for parameter estimation with a fairly small ensemble size and observed issues related to overshooting problem of permeability and porosity. Gao et al. [GZR06] pointed out the similarities between the EnKF and randomized maximum likelihood in reservoir parameter estimation problems. Liu and Oliver [LO05] used EnKF for facies estimation in reservoir models. In particular, the reservoir used in their work, included large nonlinearities in the permeability and porosity distributions based on a multi-modal PDF for the petrophysics parameters. They represented the fa- cies distribution as a two-normal Gaussian field adapted to be solved using EnKF. They showed that, EnKF can fail when used with multi-modal distributions and 50 there might be no systematic way to apply EnKF for such problems. Wen and Chen [WC05] used EnKF for permeability estimation in a 2D reservoir problem where the effects of the ensemble size on the parameter estimation problem was an- alyzed. Zafari and Reynolds [ZR05] investigated the validity of the linear updating algorithm of EnKF with simple but highly nonlinear models and also reported that EnKF cannot be used with the resulting multi-modal distributions where simple average and variance do not effectively characterize the distribution. Their results nonetheless showed a clear improvement as the non-linearity of the model was re- duced. Other applications of EnKF for parameter estimation can be found in Refs. [LNL03], [Kiv03], [AH04], [AHEM05], and [MSGH05] to name a few. EnKF has recently gained great popularity due to its simplicity in implemen- tation and the flexibility in handling errors in initialization and evolution of the model state. Furthermore, the “Monte Carlo” nature of EnKF makes it directly suitable for running on parallelized computing resources. Utilizing EnKF in high dimensional problems is not straightforward as EnKF requires an ensemble of model states that approximates a posterior PDF and ob- taining an accurate PDF in high dimensional problems may produce undesired redundancy and extreme uncertainties. By nature, EnKF locates the “mean” so- lution as opposed to the “mode” and is limited in a search space spanned by the initial ensemble. Even more problematic are emerging non-Gaussian distributions in the course of the dynamics. More explicitly, the forecasting step in the EnKF 51 propagates the probability distribution fully where as the update part solely deals with the mean and variance of the model state. The Gaussian assumption makes the update part relatively easy to perform numerically but can carry large statisti- cal mismatches when the model variables are non-Gaussian. An important example of a non-Gaussian behavior is the variation in water saturation at the fluid front. Several methods have been introduced to deal with the non-Gaussianity problem. For example, a re-parametrization of saturation values near the front was used by Chen et al [COZ09], where the emergence time of the front at a specific satura- tion value is designated as a new variable and is updated along with the rest of the model variables. The saturation values near the front are then continuously dropped from the model state to ensure that the parameters in the model state remain approximately Gaussian. It was shown that using the time of front arrival instead of saturations near the front area improves EnKF considerably. Focusing on the front saturations, Oliver and Gu [OG06] proposed a new method of reducing the non-Gaussian distribution of the saturation values at the front area by utiliz- ing the normal score transformation of the saturation values at the margins. Even after the back-transformation, the saturation values showed large variations. This required using the location of the front instead of the front saturations in the up- date step. This approach works well with EnKF in one dimension but is yet to be fully implemented for two and three dimensional problems where the location of the front is difficult to define/find. 52 4.1.2 Production Optimization In production optimization of reservoirs, we focus on maximizing ultimate recovery and financial gain from the reservoir by assigning different control patterns for the operating wells. This needs a reliable simulation model with an accurate forecasting ability. Geological models come with high uncertainties that reduce the accuracy of the optimization procedure unless they are combined with a parameter estimation technique that also reduces the uncertainties on the fly. Such “closed-loop produc- tion optimization” combines a mathematical procedure for production optimization with a data assimilation technique in order to both characterize and optimize the output from the reservoir of interest. Typically the work flow of a closed loop pro- duction optimization is as follows: first an initial geological model along with an initial production strategy is built based on initial knowledge of the reservoir. Sub- sequently, at each data assimilation step (with arrival of actual production data) the reservoir geologic parameters are updated using the EnKF procedure for ex- ample. Production strategy is then optimized using the newly updated reservoir parameters. This set of update/optimize scenario is the basis of the ensemble based closed loop production optimization. Gradient based approaches are commonly used to solve optimization problem where the direction of search for a local minimum is obtained by calculating the gradient of the objective function (NPV in this case) with respect to the control variables (injection rates). Owing to the nonlinearity of multiphase flow problems, 53 the gradient based method has to be done in an iterative form and requires extensive computational resources. As an approximation, an ensemble based method may be used to compute the gradient through a Monte Carlo evaluation of the gradient over an ensemble of states corresponding to the uncertainties in the description of the model. Any optimization technique can in principle be combined with an incremen- tal history matching procedure to build a closed loop production optimization as addressed in several works such as those of Brouwer et al. [BNJ + 04], Sarma et al. [SDA05], and Wang et al. [WLR07]. Brouwer et al. [BNJ + 04] used an ad- joint method for optimization problem and used an EnKF for model updating. Sarma et al. [SDA05] used an adjoint method for both model updating and pro- duction optimization. Wang et al. [WLR07] used EnKF with three optimization methods on a small reservoir and concluded that the steepest ascent method with gradients obtained from numerical perturbation leads to better results compared to those obtained using an ensemble. Chen et al. [COZ08] proposed the concept of an ensemble based closed loop optimization (EnOpt) in which an ensemble of model system is used for both characterization and production optimization and can be used with essentially any reservoir simulation/modeling mechanism. They concluded that EnOpt has two distinct features. First, the search direction used in optimization scenario can be obtained using an ensemble. Second, EnOpt focuses on maximizing the expectation of the NPV instead of maximizing the NPV over 54 a single reservoir model. This leads to the robustness of EnOpt with respect to uncertainties of the reservoir geologic model. 4.2 Mathematical model and Methodology 4.2.1 Kalman filter KF is the starting point for understanding EnKF. A good description of KF is given in Ref. [BH92]. KF has been used for parameter estimation problems that involve a time evolution or more generally a dynamical system in the mathematical sense [Gel02]. The dynamical system is described by a state vector that contains static (parameters) and dynamical variables that encompass all degrees of freedom in the problem. The components of the variables are forecasted at each time step by employing a prescribed dynamics of the system and are subsequently updated when the observation data becomes available. The state vector in general contains variables that are not known precisely and need to be estimated. We also typically include the set of observables in the state vector for convenience. The state vector y at each time consists of N y model variables m (e. g., porosity, permeability, saturations, pressures, etc. at the grid coordinates) andN d observation variablesd, (e.g., well production rates, well bottom hole pressures, water-oil ratios, etc). Thus at each time step the state vector y can be represented as: 55 y =      m d      for simplicity we assume that the observation variables d can be obtained from the model state by a linear projection: d obs =G(m true ) + =Hy true + whereG is the so called sensitivity matrix and is the vector of measurement errors with a known covariance given by: E[ T ] =C D . The matrixH∈R N d ×Ny is the measurement operator that projects the state vector y to the observation valuesd. In many cases measurements are not linearly related to the components of the state vector, yet by including the observables along with the variables in the state vector we can still conveniently relate the observables to the state vector by the matrixH. Since observation values are part of state vector, the matrix H will only consist of 0s and 1s. More specifically: H = [O|I], 56 whereO denotes anN d ×(N y −N d ) matrix with all zero entries andI is theN d ×N d identity matrix. Let y p denote the prior estimated vector for the variables. We define an objective function S(y) that gives us the best estimate once minimized. S(y) is written as: S(y) = 1 2 (Hy−d obs ) T C −1 D (Hy−d obs ) + 1 2 (y−y p ) T C −1 y (y−y p ) in which C y is the covariance matrix for the variables: C y =      C M C M G T GC M GC M G T      where C M is the covariance matrix of the model variable m. The best estimate for y, obtained from minimizing S(y) is given by: hyi = y p +K(d obs −Hy p ) K = C y H T (HC y H T +C D ) −1 where K is the so called Kalman gain matrix. The covariance matrix C M 0 after each assimilation step is given by 57 C M 0 = (C −1 M +G T C −1 D G) −1 =C M −C M G T (C D +GC M G T ) −1 GC M The covariance matrixC M 0 may become too large for nonlinear problems and hence lead to increased uncertainty. Furthermore, the sensitivity matrix G may not be easily available as a linear function or maybe hard to formulate analytically. In practice, these problems prevent Kalman filter to be used in nonlinear regimes. 4.2.2 Ensemble Kalman filter Ensemble Kalman filter (EnKF) is essentially a Monte Carlo generalization of KF. An ensemble of realizations is used to describe the probability distribution of the model state. Similar to KF, the model state contains both the static and dynamic variables as well as the observation values. In Ensemble Kalman filter, the model state is updated based on a weighted covariance function between the measurements and the evaluated values of the observables from the model variables at each time of measurement. Algorithmically, EnKF consists of two parts: a forecast step and an analysis step. In the forecast step, the model variables in for all ensemble members are propagated forward in time based on the model dynamics through simulation. The update step modifies the state variables in all ensemble members based on comparison with the observation values obtained at each assimilation time with a weight given by the Kalman gain matrix K. 58 4.2.2.1 Ensemble representation of error statistics Consider the matrix Y that consists of columns of variables of each member j of the ensemble y j (t i )∈R n at the time t i . As described in section 4.2.1, the state vectory j (t i ) is a combination of the dynamic variablesη j (t i )∈R nη and static ones α j (t i )∈R nα where n =n η +n α total number of variables. Thus we can represent the matrix Y as follows: Y i =Y (t i ) =      η 1 (t i ) ... η Ne (t i ) α 1 (t i ) ... α Ne (t i )      where N e is the total number of realizations or the ensemble population. The ensemble mean is stored in each column of the matrixhY (t i )i: hY (t i )i =Y (t i )1 Ne where 1 Ne ∈R Ne×Ne is a matrix with all elements given by 1 Ne and matrix product is used. The perturbation matrix Y 0 is then defined as Y 0 (t i ) =Y (t i )−hY (t i )i =Y (t i )(I− 1 Ne ). Thus the ensemble covariance C Y (t i )∈R n×n can be defined as follows C Y (t i ) = Y 0 (t i )Y 0 (t i ) T N e − 1 . 59 The ensemble covariance denotes the smoothness and variation of the ensemble members. 4.2.2.2 Ensemble representation of measurements We define the vector measurement d∈R m where m is the number of total mea- surements. We perturb the N e measurement vectors d j as follows d j =d + j and define the matrix of measurements D: D = (d 1 ,d 2 ,...,d Ne )∈R m×Ne . The ensemble of measurement perturbations with mean zero is stored in the matrix E: E = ( 1 , 2 ,..., Ne )∈R m×Ne The covariance of the measurements errors is given by C D = EE T N e − 1 60 4.2.2.3 Analysis equations We can write the analysis (update) step for the ensembley a j (j denotes the ensemble member index between 1 and N e ) as follows y a j =y j +C Y H T (HC Y H T +C D ) −1 D 0 j (4.1) where the ensemble innovation vector D 0 j is given by D 0 j =d j −Hy j whereH project the model state on to the measurements. If there is no such linear relationship between the measurement values and the ensemble variables we can modify the state vector so that it also contains a set of calculated (by any means) observation values as well as we did in Section 4.2.1. We have to take special care when we prepare the initial representation of the ensemble. Since the estimates we get at the end (after each update of EnKF) will be contained in the space spanned by the initial ensemble members. This amounts to setting up the initial ensemble in such a way that the range of the initial ensemble parameters contain the true values. This can be solved if we increase the ensemble size and pick up our initial ensemble carefully. 61 4.2.3 Ensemble Based Closed loop Production Optimization (EnOpt) In this section, we briefly describe the formulation of the EnOpt algorithm used in our study. A more elaborate description can be found in Chen et al. [COZ08]. Let y denote the state vector of the system as before. For example, for the injection and production wells in our study, the variable y at each time contains (i) static parameters{λ kj },{τ j },{J j },{a j } and{b j }, (ii) the dynamical variables: accumulated injected water for each well{W j }, and (iii) the set of observables: total production for each well{q j } and oil production rate{q o,j }. An ensemble of state vector variables is stored in Y = (y 1 ,y 2 ,y 3 ,...,y Ne ). The EnKF procedure provides the updating procedure for the ensemble of states, given the observation data d obs : y u j =y f j +C Y H T (HC Y H T +C D ) −1 (d obs,j −Hy f j ) (4.2) where j is the ensemble member index. The variables d obs,j are the perturbed observations for each ensemble member. The ensembleY u = (y u 1 ,··· ,y u Ne ) indicates the updated ensemble while Y f = (y f 1 ,··· ,y f Ne ) indicates the forecasted ensemble variables from the reservoir simulation. Matrix C D indicates the covariance of the measurement noise and C Y is the covariance matrix for the state vector variables 62 andhY f i is the mean state vector as defined earlier. The matrix H projects the observed variables from the state vector. For optimization, we need to store an ensemble of control variables: let x be a vector describing the control variables that contains all well constrains (at all control steps) x = [x 1 ,x 2 ,x 3 ,...,x Nx ] where N x is the total number of control variables (e.g. number of wells× number of control steps). The net present value (NPV) is the objective function for the optimization in this thesis and is given as a function of the control vector x and the reservoir state y: g(x,y) = Nt X i=1 v o Q o,t | x,y −v w Q w,t | x,y (1 +r τ ) t i /τ (4.3) wherei is the time step index,N t is the number of time steps, andr τ is the discount rate in terms of time spanτ. The constantsv o andv w are the unit prices of oil sale and water disposal cost respectively andQ o,t | x,y andQ w,t | x,y are total oil and water produced from all the wells at t for the control x and state variable y. The control variables are optimized within the uncertainties of the model by averaging: g Y (x) = 1 N e Ne X j=1 g(x,y j ) (4.4) 63 The subscriptY indicates that NPV is based on the updated model obtained from EnKF. The set of model state vectors Y does not change during the optimization part where NPV is optimized based on the control variablesx. At each optimization step an updated control variable x l is produced via x l+1 = 1 α C x C x,g(x) +x l (4.5) where l denotes the optimization iteration step, α l is a tuning parameter which determines the step size,C x,g(x) is the cross covariance between the control variables and g(x): C x,g(x) t 1 N e − 1 Ne X j=1 (x l,j −hx l i) (g(x l,j ,y j )−<g(x l ,y)>) (4.6) wherex l.j is an ensemble of perturbed values of the control variablex of sizeN e with an averagehx l i and prior covariance C x . We use a temporally correlated Gaussian random field with zero mean as the perturbation for x l . In the notation we just introduced, the EnOpt procedure can be represented in a simplified form as the following algorithm. The stopping criteria typically includes a maximum optimizations stepl max or unsuccessful search for the tuning parameter α: 1. Initialize the state vector ensemble Y and the initial guess for the control x. Set t = 0. 64 2. Propagate the state vector one step in time using a dynamical equation (Eq. (3.11)) and the current value of x and set t =t + 1. 3. Use EnKF (Eq. (4.2)) to update Y . 4. Start optimization at l = 1. If this is the first time step (t = 1) then x 1,j (control variables at iteration 1 over ensemble members) would be generated in two steps. First, a mean control is sampled from the uniform distribution with the upper and lower limits equal to the maximum and minimum possible well constraint for each realization of each well. Second, a temporally correlated Gaussian random field with zero mean is generated for each realization of each well and added to the mean control. We then define x 1 = 1 Ne P Ne j=1 x 1,j . If this not the first time step (t6= 1) thenx 1 =x and the temporally correlated Gaussian random field with zero mean is generated to each realization of each well and added to the control variables x 1 for form x 1,j . 5. If this is not the first iteration (l6= 1) then the temporally correlated Gaussian random field with zero mean is generated for each realization of each well and will be added to control variables x l . 6. Calculate g(x l,j ,y j ) [Eq. 4.3] by running the CRM formulation from the current time to the end time of the simulation. y j are member of the Y ensemble. 65 7. Calculate the cross covariance C x,g(x) using Eq. 4.6 we then compute the updated control variables x l+1 using Eq. 4.5. 8. Obtain the objective function g(x l+1 ) using Eq. 4.4 which only needs N e simulation runs. 9. If g(x l+1 ) > g(x l ), let x l = x l+1 and l = l + 1. If not keep x l , increase α, calculate x l+1 using Eq. 4.5 and repeat from step 8. 10. If a stopping criteria is reached , setx =x l and exit the loop of optimization. We then repeat from step 2 up to the end of data assimilation time. otherwise repeat from step 5. The arbitrary parameters in the optimization procedure such as the stopping cri- teria or the initial value of α are adjusted initially at the actual implementation of the simulation/optimization procedure. The stopping criteria used in our work is (i) the relative increase of the net present value had been less than 0.01%, (ii) the change in the control variables at two previous times had been less than 1%, and (iii) we were not forced to increase the tuning parameter α more than twice. 66 Chapter 5: Numeric Studies with Synthetic Reservoirs 5.1 Introduction In this chapter, we report the numeric implementations of our methodology. We first discuss the setting and notation for our implementation more in terms of the actual reservoir observables and the CRM parameters. We then provide the details of our case studies that cover a few some interesting examples which serve to evaluate the effectiveness of our procedure. In these case studies, with the exception of the cases 7 and 8, the role of the true reservoir is played by a numerical reservoir simulator (Computer Modeling Group). 5.2 EnKF/CRM Setup Each state vector y consists of the following variables: y = [{λ ij },{τ j },{J j },{a j },{b j },{W j },{q o,j },{q t,j }] † (5.1) 67 where the index i refers to the injection wells and j refers to the production well indices. We assume that there are N e members in the ensemble. The CRM and phenomenological parameters are all flattened to form a one-dimensional vector. The dynamical and observable variables{W j },{q o,j },{q t,j } denote the cumulative injected water, oil production and total production for the production wells. The total number of variables in state vectory is equal tom×n+7×n in whichm is the number of injectors and n is the number of producers. For the case of unbalanced capacitance model an extra parameterq c,j (inequality term) has to be added to Eq. 5.1. We first describe how the guess distribution for the state vector at the initial time (t =t o ) is prepared. Each element iny is taken from an independent Gaussian distribution with a prescribed average and standard deviation. Further knowledge of interdependence or correlation among the parameters may require generation of a specific distribution with a more complex correlation matrix as opposed to the simple standard deviation. Mathematically this is equivalent to assuming an already diagonal correlation matrix for the initial distribution of the variables. The initial estimated averages for the static variables{λ kj } are each taken to be proportional to the inverse distance between the wellsk andj as expected from 68 the interpretation of{λ ij (t 0 )} in terms of the effective transmissibilities. Similarly the initial guess for the static variables{τ j } is chosen as τ j (t 0 )∝ m X k=1 1 λ kj (t 0 ) the coefficient of proportionality is deduced from the time scales of the reservoir. This inverse relationship betweenτ j (t n ) andλ kj (t n ) is discussed in detail in Yousef et al. [YGJL06] and is based on fluid dynamical considerations. Nonetheless, in our numerical study, we opted not to use the distance based initial estimate for λ kj (t 0 ) and τ j (t 0 ) in order not to oversimplify the CRM characterization problem for our demonstration purpose. The average guess for{λ kj } is thus taken to be the naive value of 1 n , in which n is the number of production wells. The{τ j } parameters assumed to be between 1− 10 days for different synthetic examples. The initial average values fora j andb j andJ j are based on the prior knowledge of the reservoir or an initial curve fitting. These are described for each synthetic numerical study later in this Chapter. The initial guess for the observable variables q t and q o are set at the initial total and oil production since these quantities are known initially. The initial guess for the dynamic variableW j is equal to the initial production rate. The standard deviations of the variables are set at a percentage of the mean. The standard deviation of the initial static variables in our numerical study for example is taken to be 50% of the average values for each element. This might be adjusted for the individual parameters to ensure that the range of initial variables is 69 large enough to contain “true parameters” . We then initialize the state matrixY as an ensemble ofN e instances ofy k variables sampled from the statistical distribution of the initial guess. Y = [y 1 ,y 2 ,···y Ne ] (5.2) The rest of the setup follows the standard EnKF procedure. At each time, we propagate each ensemble member (y j ) from the current time to the next available data assimilation step t i using Eq. (3.11), Eqs. (3.14), and (3.15) and compare the results of the forecast with perturbed real observable data. We subsequently use Eq. (4.2) for updating each member of ensemble before proceeding to the next step. Notice that the matrix H is a linear projection (onto q j and q o,j subspace of the dynamical variables) and produces the observations from each model state vector. The updated values of the state vector y j are then used to propagate the realizations forward in time up until the next data assimilation step. At each time, the average over the state vectors of the ensemble describes our best estimate of the system including the CRM and phenomenological parameters. The variance over the ensemble measures the current uncertainty/error in describing the system. 70 5.3 EnOpt/CRM additional Setup For optimization, we have parameterize the injection rates as a control vectorx. In the EnOpt/CRM setup at each step, once the CRM state vectors in the ensemble y k (Eq. 4.2 ) are determined, we optimize x, through the algorithm in Chapter 4, Sec. 4.2.3. The control vector is based on a discretization of the “life time” of the reservoir into n opt steps. We assume that the injection rate is constant during each such step and only apply the optimization procedure at the beginning of each optimization step. In the actual implementation of the algorithm, only the forthcoming “future” values of the injection rate are fixed so that the past history is shared and unaltered. This amounts to a control vector that hasm×n opt elements for them injection wells. The optimization algorithm in turn requires an ensemble of control vectors X: X = [x 1 ,x 2 ,x 3 ,··· ,x Ne ] where N e is the number of realizations (taken to be the same as the ensemble size for Y ). We assume that the only operational constraints are on the injection rates and thus on the vector x itself. We also assume that there is a maximum and minimum operating injection rate associated with each well. At the beginning of the simulation, the control vector needs to be initialized: we pick the mean control from the uniform distribution bounded by the minimum and 71 maximum possible well constraint for each realizations of each well. At later stages of the optimization algorithm a temporally correlated Gaussian random field with zero mean is added to the mean control vector. By using proper regularization, we explicitly enforce the well constraints on the perturbed control vector. Calculation of the objective function for each control/state vector pair requires running the CRM model from the current time to the supposed final time of the simulation. This in general, is the most time-consuming numerical step of the procedure. Once the control vector is determined at the end of the optimization algorithm, we use it as the “true” driving term for the CRM ensemble and the actual reservoir and continue with characterization and propagation of the system until the next control step comes up. The results of this combined procedure for various synthetic and real cases appears next. 5.4 Introduction to case studies We shall cover a series of synthetic and real cases for which we implement EnKF and/or EnOpt in the CRM framework. We use these cases to analyze the effective- ness and the sensitivity of our numerical procedure for a variety of situations. As we shall see, while more complex reservoirs with high degree of heterogeneities and rapid uneven saturation movements might not immediately satisfy the geological requirements for the CRM model in theory, it is still possible to at least partially 72 characterize and even optimize the oil production in such cases using EnKF and EnOpt procedures. A brief summary of the cases studied follows: Case 1 This case demonstrates a three well model in a homogeneous set up. In this case the transmissibility distribution does not vary with time due to the symmetry of the pattern. We model this reservoir within the CRM framework and obtain the unknown CRM parameters using the EnKF methodology. Detailed explanation of the implementation and the initialization setup are given. Almost perfect matches for both the total and oil production and the final values for all the unknown parameters are obtained in comparison to a fluid dynamical simulation. Since no variations in the transmissibility distribution occur, the CRM parameters are saturation independent and remain constant. The simplicity of this model on the other hand, allows us to do a set of sensitivity analysis with respect to the initial model estimate, ensemble size, and observation noise. The root mean square error (RMSE) and spread of the results are compared for different settings and an optimal setting is identified. Case 2 This case describes a homogeneous reservoir in which the transmissibility distribu- tion varies drastically because of the larger number of wells in contrast to Case 1. 73 Due to the variation in the transmissibility distribution across the reservoir, CRM parameters are shown to vary with time and are shown to be in accordance with saturation movements. The variation in CRM parameters can indeed be captured rapidly using EnKF. We then performed a forecast scenario in which we used the final obtained EnKF/CRM parameters to propagate the system. Due to the time- dependence of the parameters using fixed values expectedly does not give us perfect matches in the forecast scenario in the long run. A series of update-forecast cases should be implemented if near-perfect matches are desired. Case 3 This case describes a channelized reservoir in which besides varying the injection rates, we manually vary the production bottom hole pressures. This case evaluates the ability of EnKF/CRM framework to match for both the total and oil produc- tion rates when the productivity index is included to captures the effects of the variation in the bottom hole pressure. The variation of the unknown CRM param- eters with time were even more pronounced in this case due to the changes in both the injection rates and producer bottom hole pressures. Such changes were shown to be in accordance to the saturation movements and can be captured by EnKF. Approximate location of the high permeability channel was obtained through the final obtained values of{λ kj } and{τ j }. The final values for all the unknowns are then used for forecast scenarios. 74 Case 4 This case describes a peripheral water flood pattern with 12 injectors and 9 pro- ducers. A high permeability channel is also included in this case. This case tests the ability of the EnKF/CRM framework to match for the total and oil production rates where we have different and numerous well patterns. This complex scenario is used to test the variations of all the unknown CRM framework with time and with respect to the saturation movements. The final values of the unknowns especially {λ kj } and{τ j } can be used to approximately locate the channel in this scenario. Case 5 This case is a heterogeneous reservoir with two streaks of higher permeability re- gions. In this case, we introduce EnOpt/CRM framework in which we apply both the characterization along with optimization scheme . We proceed with two differ- ent scenarios. In the first scenario, we assume that the reservoir geologic features are already known and proceed only with optimizing the injection rates in order to increase the net present value of the reservoir without using a CRM model and the obtained injection rates were applied to the reservoir simulator directly. In the second scenario, we assume that the reservoir geologic features are unknown and we proceed with obtaining both the reservoir parameters and optimizing the injection 75 rates in order to increase the NPV. We will show that in both scenarios the re- sults of the cumulative production increases and the water production dramatically decreases compared to a base neutral set of controls. Case 6 This case is a channelized reservoir in which we implement Enopt/CRM techniques. This example contains a sensitivity analysis with respect to parameters associated with closed loop production optimization. The set of parameters used in the sensi- tivity analysis for this case are the spread of control noise and its temporal corre- lation. We also repeat the procedures of case 5 for this case as well. Cases 7 and 8 These cases are taken from two different sections of a field in California. We test the ability of EnKF/CRM method for matching the total production rates when we select a small and large pattern from the actual real field data. The estimated {λ kj } parameters are in general accordance with the fracture orientation in this field. 76 permeability 40md porosity 0.18 dimensions 40× 40× 6ft 3 grid 31× 31× 5 oil-water mobility ratio 1 oil compressibility 5× 10 −6 1 psi water compressibility 10 −6 1 psi rock compressibility 10 −6 1 psi bottom hole pressure 250psi Table 5.1: Synthetic Case 1 defining parameters. 5.5 Homogeneous Reservoir 5.5.1 3-Well Model (Case 1) The synthetic field in this case consists of 1 injector and 2 producers over a 2D square area depicted in Fig. 5.1. The system is homogeneous and uniform with the details defining this field given in Table 5.1. The injection well is controlled by a water injection rate and the production wells are controlled by fixed bottom hole pressures of 250 psi. The actual reservoir response is simulated using the Computer Modeling Group (CMG) reservoir simulator. Fig. 5.2 depicts the applied injection rate that provides the time-dependent driving term in the CRM dynamical system. Figs. 5.3 and 5.4 depict the total production rates and oil rates for all production wells respectively and form the observable data from the reservoir simulator. Notice the similarities between the two production rates due to the homogeneous characteristics of the system causing a symmetric production pattern. 77 Figure 5.1: Synthetic Case 1: Homogenous three well model Figure 5.2: Injection rate used for synthetic Case 1 78 Figure 5.3: Observed total production rates for synthetic Case 1 Figure 5.4: Observed oil production rates for synthetic Case 1 79 inter-well connectivities λ kj 0.5 time-lag τ j 2days a j 10 −8 b j 1.3 relative standard deviation 50% initial total production q j (t 0 ) 3.26× 10 3 initial oil production q o,j (t 0 ) 2.32× 10 3 Table 5.2: Initial guess for synthetic Case 1 We start with an initial “guess” distribution of the static variables as given in Table 5.2. The relative noise amplitude used to perturb the actual observation data is 10%. With all the elements in place, we apply EnKF/CRM progressively in time. To emphasize the large spread of our initial estimate, we simply propagate the CRM formulation up to 500 days without utilizing EnKF. After this period, we switch on EnKF until day 2500. We revert back to pure propagation until the final day of 3000 using the resulting estimate for the state vector to validate the forecasting abilities of our estimate for the CRM parameters. Fig. 5.5 depicts the results of EnKF matching of the total production rate. The red curve shows the actual observed data and the black curves depict the results of 100 ensemble realizations. Clearly, the ensemble members converge to the actual observation data once EnKF is switched on after day 500. Also note that after day 2500 when EnKF is switched off, the state vectors of the ensemble still closely follow the observed data until day 3000. Fig. 5.6 shows similar trends for the oil production rates. 80 We next focus on the connection between the time variation of{λ kj } and the saturation movements. Fig. 5.7 depicts the saturation contour maps for the initial time (day 0), time at which we start the EnKF (day 500), mid-time (day 1500), and the final time of the simulation (day 3000). The homogeneous nature of the system forces the front movement to be symmetric in shape and size in the X-Y plane. In order to relate the interwell connectivities to the saturation movements, we similarly depict the{λ kj } distribution obtained from EnKF-CRM in Table (5.3) which displays (mean values of ) the final estimate for the unknown parameters in the EnKF/CRM formulation for the homogeneous system. The estimated value of {λ kj } is very close to 0.5 and clearly indicate that the production wells receive the same portions from the injection well. This is expected from the symmetric position of the production wells with respect to the injection well (and the reservoir). We do expect variations in transmissibility due to relative permeability variations based on the saturation but this does not yet affect the values of{λ kj } due to symmetric positions of the production wells. We report the variation of λ 11 and λ 12 in the reservoir in Fig. 5.8 over 100 realizations where the average estimate for λ 11 and λ 12 is highlighted. Most prominent feature in this period is the reduction in the uncertainty (spread) of λ 11 and λ 12 among the ensemble members immediately at the time of applying EnKF. In particular, notice that the values of λ 11 and λ 12 almost instantly converge to a final constant value of 0.5 after day 500. 81 P 1 (j = 1) P 2 (j = 2) λ 1j 0.5 0.5 τ j (days) 3.44 3.36 a j 0.19× 10 −6 0.16× 10 −6 b j 1.16 1.17 Table 5.3: Final estimates (over 100 realizations) of all the unknowns in homogeneous synthetic Case 1 Figure 5.5: Matching total production rates for synthetic Case 1 82 Figure 5.6: Matching oil production rates for synthetic Case 1 In contrast with{λ kj }, the values for{τ j } depend more critically on the injection rates and thus are expected to vary with time initially. Fig. 5.9 shows the variation of τ 1 and τ 2 over 100 realizations. Notice that only after day 1700, the values of {τ j } do converge to a constant values. Incidentally, the time at which the{τ j } values stabilize in this simple case is close to the water breakthrough time. The final mean values of the parameter τ j are reported in Table (5.3). Along with the unknown parameters of the CRM formulation, we have also estimated the phenomenological parametersa j andb j that characterize the oil pro- duction rate. The estimates for these phenomenological parameters approach a constant value in the later stages of the reservoir depletion as expected from the phenomenological model. In fact, Eq. (3.13) is best applicable only after water breakthrough. Figs. 5.10 and 5.11 depict the variation of the estimates for{a j } 83 Figure 5.7: Water saturation contour map at 4 different time: (a) t = 0; (b) t = 500 days; (c) t = 1500 days; (d) t = 3000 days, synthetic Case 1 84 Figure 5.8: {λ kj } variations over 100 realizations, synthetic Case 1 Figure 5.9: τ 1 and τ 2 variations over 100 realizations, synthetic Case 1 85 P 1 (j = 1) P 2 (j = 2) a 1j 0.11× 10 −6 0.11× 10 −6 b j 1.2 1.2 Table 5.4: Parameters a j and b j obtained from fitting the power law model and{b j }. The results clearly depict that only after day 1500 the estimates for the parametersa j andb j do converge to the constant values. To validate that our phe- nomenological parameters do correspond to the reservoir characteristics, we plot the cumulative water injected vs the water oil ratio in a logarithm scale for each production well. The slope of the curve beyond the 50% water cut point is used to obtain the parameters a j and b j (Fig. 5.12). Table (5.4) contains the results of this fitting procedure. This model provides a constant set of a j and b j that can only match the oil production at the late stages of reservoir depletion. The fitting procedure gives values for a j and b j that are close to the values obtained from the EnKF/CRM procedure. In our opinion, EnKF/CRM algorithm can estimate the phenomenological parameters (as time-varying parameters) even at the early stages of the reservoir and lead to an accurate match of the oil production rate from early on [See Fig. 5.6]. Sensitivity Analysis In this section, we provide a sensitivity analysis of the efficiency of our EnKF/CRM scheme as a function of the simulation constraints and the initial state distribution. Results of a sensitivity analysis in principle allow us to optimize the arbitrary 86 Figure 5.10: a j variations in a logarithmic scale over 100 realizations, synthetic Case 1 Figure 5.11: b j variations over 100 realizations, synthetic Case 1 87 Figure 5.12: cumulative water injected vs water oil ratio in a logarithmic scale, synthetic Case 1 Figure 5.13: RMSE variation of total production rates for various en- semble populations, synthetic Case 1 88 parameters in the simulation so that the highest accuracy for a given model may be obtained in shortest time by efficient use of computer memory and also provides an insight into how the accuracy of EnKF depends on various statistical and physical assumptions inherent in our analysis. We particularly focus on ensemble population size, initial model estimates, and observable noise. Two quantities form the basis of our sensitivity analysis: (i) The (relative) root mean square error (RMSE) directly describes the error in the observable quantities from the ensemble estimate: RMSE = v u u u t 1 n n X j=1 " hq j (t)i−q true j (t) q true j (t) # 2 (5.3) where q true j (t) is the actual observable data at time t from producer j which may be obtained from the true field or a simulated scenario andhq j (t)i represents the average over the realizations at time t for producer j and is defined by: hq j (t)i = 1 N e Ne X m=1 q m j (t) (5.4) q m j (t) represents observable data at time t from producer j in the m-th member of the ensemble (ii) The ensemble spread σ describes the uncertainty range in the observable quantities from our ensemble values : σ 2 = 1 N e − 1 n X j=1 " Ne X m=1 h q m j (t)−hq j (t)i i 2 # (5.5) 89 Fig. 5.13 depicts the variation of RMSE of the total production rates as a func- tion of simulation time for different ensemble population sizesN e , varying between 20 and 200. The rest of the parameters in this simulation are the same as in Section 5.5.1. Notice that RMSE is plotted starting with the day at which EnKF proce- dure begins (day 500). Clearly, RMSE decreases with increasing the ensemble size. Notice also that when number of ensemble members are too small (N e = 20) the results are not satisfactory. This is naturally expected in any Monte Carlo based numerical method as increasing the ensemble size allows for smoother and finer variations among the ensemble members. This is a useful fact to consider when confronted with a large initial uncertainty as the effective size of the possible state space is enlarged along to allow for a broader range of parameters. Nonetheless, there is a limit to the accuracy obtainable by increasing the ensemble size as the state space is not practically infinite. This can be observed in Fig. 5.13 by not- ing the decreasing difference between ensemble error results as the ensemble size is increased, in particular notice that RMSE values for N e = 100 and N e = 200 are practically of the same order. Using a larger ensemble size also requires higher computational power and computer memory and leads to longer running times. We have thus fixed the ensemble size atN e = 100 which provides the right balance of accuracy and speed. The ensemble spread as a function of time with different ensembles sizes is depicted in Fig. 5.14 where the effective state space shrinks as a function of time for various values of the ensemble size. 90 Figure 5.14: Spread variation of total production rates for various en- semble populations, synthetic Case 1 The efficiency of any parameter estimation scheme can depend critically on the initial estimate for the parameters. In our simulation, the initial estimate is de- scribed by a normal distribution of{λ kj },{τ j },{a j }, and{b j } each characterized by an average and a (relative) standard deviation. We first vary the relative stan- dard deviation of the static variables (σ k ) that describes the initial range of the ensemble estimate. In Fig. 5.15, we graph the error as a function of time for differ- ent values of the initial standard deviation of the ensemble. Note that for certain values of σ k the space span by the initial estimate ensemble is simply too narrow [e.g. σ k = 0.01] to include the true field values . As we increase σ k the error shrinks, which guarantees that the parameter space is large and thus more likely includes the true parameters. For this reason, care must be taken for the choice of 91 Figure 5.15: RMSE variation of total production rate for various σ k values, synthetic Case 1 the initial estimate for the parameter space. Just as a narrow guess might miss the mark with EnKF, a guess that is too broad will not yield meaningful results either. It turns out that a relative standard deviation of σ k = 0.5 provides a reasonable initial spread. Similarly, the spread of the ensemble is depicted in Fig. 5.16 with various values for σ k . Increasing σ k values will increase the spread as expected. Finally, we consider the sensitivity of EnKF/CRM procedure to the observation noise (σ obs ). Figs. 5.17 and 5.18 depict the results of RMSE and spread variations of the total production rates with time for σ obs ∈{0.05, 0.1, 0.5}. In contrast to the ensemble standard deviation σ k , increasing σ obs results in the error increas- ing. Effectively, a noisy observation data is associated with the uncertainty of the measurements. A small perturbation on observed values is required for ENKF to 92 Figure 5.16: Spread variation of total production rate for various σ k values in a logarithmic scale , synthetic Case 1 provide a constant small variance on the updated ensemble parameters when the procedure converges to the actual parameters. A larger perturbation generated by higher σ obs values on the other hand may cause the solution to diverge from the true values. To summarize, we applied EnKF/CRM methodology to a simple 3-well model in a homogeneous system. We showed that the interwell connectivities{λ kj } are not saturation dependent and remain fixed due to the symmetry. On the other hand, the time lag parameter{τ j } does vary with the injection rate variations. Near perfect matches for the total and oil production model obtained and the estimated parameters were successfully evaluated for forecasting the system. We also analyzed the sensitivity of the results with respect to the ensemble size, the 93 Figure 5.17: RMSE variation of total production rate for various σ obs values, synthetic Case 1 Figure 5.18: Spread variation of total production rate for various σ k values, synthetic Case 1 94 initial model estimate, and the observation noise. We conclude that increasing ensemble sizes will decrease the error of our matching algorithm up to a point. Similarly, increasing the initial model uncertainty decreases error to a point. On the other hand, increasing the observation noise increases the error. 5.5.2 Five Spot pattern (Case 2) Synthetic case 2 consists of 5 injectors and 4 producers over a homogeneous 2D square area as depicted in Fig .5.19. The grid block sizes, dimensions, and the fluid properties are the same as Case 1 in Sec. 5.5.1. We particularly focus on the variation of{λ kj } and{τ j } in this homogeneous system where the transmissibility distribution can vary with time. According to Eq. 2.3 the variation in transmissibil- ity is mostly due to the saturation changes across the reservoir. Consecutively, the estimate for the parameters{λ kj } and{τ j } are directly correlated to the saturation movement. The actual injection and the production responses are depicted in Figs. 5.20, 5.21, and 5.22. Notice the large fluctuations of the injection wells in this setup. The similarities between the production rates for all the producers are due to the homogeneous nature of the system. The injection rates are the time-dependent driving terms in the EnKF/CRM framework and are the input to our dynamic algorithm. 95 Figure 5.19: Synthetic Case 2: Homogeneous five spot pattern Figure 5.20: Injection rates used for synthetic Case 2 96 Figure 5.21: Observed total production rates for synthetic Case 2 Figure 5.22: Observed oil production rates for synthetic Case 2 97 We have chosen the naïve initial guess ofλ kj (t 0 ) = 1 4 in our numerical procedure. The initial guess for τ j (t o ) is 10 days. The initial average values of a j and b j are based on the prior knowledge of the reservoir at 10 −14 and 1 [in units of days and bbl/day] respectively. The relative standard deviation of the initial static variables is equal to 50%. The initial guess for the dynamical variable q j (total production) and q o,j (oil production) are equal to the actual values at t = 0. The noise used to perturb the actual observation data is 10%. The number of realizations used in this study is 100. Similar to the setup for case 1 [Sec. 5.5.1], we apply EnKF/CRM progressively in time. To emphasize the large spread of our initial estimate, we simply propagate the CRM formulation up to 500 days without utilizing EnKF. After this period, we switch on EnKF up to day 2500. We revert back to pure propagation until day 3000 using the resulting estimate for the state vector to validate the forecasting abilities of our estimate for the CRM parameters. Fig. 5.23 depicts the results of EnKF matching of the total production rate. The red curve shows the actual observed data and the black curves are the results of 100 ensemble realizations. The ensemble members converge to the actual observation data once EnKF is switched on after day 500. Also note that once EnKF is switched off after day 2500, the state vectors of the ensemble still very closely follow the observed data until day 3000. Fig. 5.24 shows similar trends for the oil production rates. 98 Figure 5.23: Matching total production rates for synthetic Case 2 Figure 5.24: Matching oil production rates for synthetic Case 2 99 We next focus on the variations in{λ kj } and the saturation distribution. Fig. 5.25 shows the saturation contour maps for the initial time (day 0), the time at which we start the EnKF (day 500), mid-time (day 1500), and the final time of the simulation (day 3000). Injection rates are different (randomly chosen) for each in- jector and cause the saturation fronts to propagate in different directions across the reservoir. This irregular saturation movement causes a variation in the transmissi- bility and thus the estimate for the parameter{λ kj } and{τ j }. In order to relate the interwell connectivities to the saturation movements, we similarly depict the {λ kj } distribution obtained from EnKF-CRM methodology in Fig. 5.26. Clearly the initial realizations (t = 0) do not carry much information about the system and display a broad variation. The variation among the ensemble realizations is reduced once we implement EnKF (t = 500). As we go forward in time (t = 3000), the un- certainty associated to{λ kj } values is systematically reduced. Notice also the last plot showing the final estimated{λ kj } values. For example, we can identify index 9 on the horizontal axis with the largest{λ kj } value corresponding to (k,j) = (4, 2). The index 15 is the next largest value corresponding to (k,j) = (3, 3). These{λ kj } values are consistent with saturation movements around these wells. Fig. 5.27 shows the mean values of the final estimate of{λ kj } values in the EnKF/CRM for- mulation at the end of the simulation. Note that nearby well pairs expectedly have larger{λ kj } compared to distant wells. These results confirm our initial conjecture that in case of large variations in reservoir saturation, a time-varying (possibly on 100 a phase by phase basis) set of parameters for CRM model can better model the reservoir at the cost of adding to the complexity of the model. To explore this premise even further, we follow the variation of λ 23 and λ 53 in the reservoir. Figs. 5.28 and 5.29 show the time-variation of λ 23 and λ 53 over 100 realizations where the average estimate for λ 23 andλ 53 is also highlighted. As noted, the data for the first 500 days depict the results of evolving the initial estimate ensemble while the EnKF procedure is applied from day 500 up to day 2500. Most prominent feature in this period is the reduction in the uncertainty (spread) of λ 23 and λ 53 among the ensemble members. Finally we switch off EnKF at day 2500 and use the up- to-date estimate for CRM parameters at day 2500 to forecast the behavior of the reservoir up to day 3000. The ensemble parameters are constant during this final period. The evolving nature of the CRM parameters (through their dependence on saturations) suggests that a final estimate for the parameters might not be reliable beyond a certain time for forecasting unless the reservoir saturation features have been stabilized. In particular, notice that the values of λ 23 and λ 53 converge to a final constant value on and after day 2300. In fact, following Fig. 5.25 for satura- tions, we observe that λ 23 and λ 53 converge only after the saturation distribution has been more or less stabilized. We similarly depict the{τ j } parameters obtained from EnKF-CRM methodology in Fig. 5.30. Clearly, the uncertainty among{τ j } parameters is reduced systematically. It is evident from Fig. 5.30 that τ 1 and τ 2 have the smallest values corresponding to the production wells 1 and 2 possibly due 101 P (j = 1) P (j = 2) P (j = 3) P (j = 4) λ 1j (k = 1) 0.32 0.25 0.23 0.18 λ 2j (k = 2) 0.34 0.17 0.31 0.17 λ 3j (k = 3) 0.20 0.26 0.21 0.31 λ 4j (k = 4) 0.18 0.44 0.12 0.24 λ 5j (k = 5) 0.16 0.18 0.37 0.27 τ j (days) 0.7 0.1 7.6 7.4 a j 0.05× 10 −10 0.06× 10 −10 0.0009× 10 −10 0.43× 10 −10 b j 1.96 1.91 2.21 1.7 Table 5.5: Final obtained values for all the unknowns in the EnKF-CRM methodology, synthetic Case 2 to the fact that these production wells are located near injectors 1 and 2 with higher water injection rates. A similar pattern is also observed in Fig. 5.31 depicting the variation of τ 4 over 100 realizations. Table 5.5 depicts the final values of all the unknown in the EnKF-CRM methodology. Figs. 5.32 and 5.33 depict the variation of the estimates for{λ kj } and{τ j } over time respectively for all injectors k and producers j. These graphs, together with the saturation distribution evolution in Fig. 5.25 affirm our intuition that the variations in injection and production patterns can influence the front movement and through the implied changes in the saturation lead to a varying set of CRM parameters. In general, we expect the variations in the estimates of{λ kj } and{τ j } to provide an insight into fluid front movement in reservoirs. In this section we applied the EnKF/CRM methodology to a 5× 4 well model in a homogeneous system. This is an example of a system in which transmissibility distribution does change with time due to the large number of wells and uneven and moving saturation distribution. We showed that all the estimated CRM parameters 102 Figure 5.25: Water saturation contour map at 4 different times: (a)t = 0; (b) t = 500 days; (c) t = 1500 days; (d) t = 3000 days, synthetic Case 2 103 Figure 5.26:{λ kj } distribution at 4 different time, synthetic Case 2. The horizontal axis is indexed by (j− 1)×m +k corresponding to{λ kj }. The maximum and minimum value of λ for each case are depicted with error bars and the mean is depicted with circle marker. Figure 5.27:{λ kj } distribution, synthetic Case 2. Each arrow represents the mean estimate of{λ kj } at the final time of the simulation and the length of the arrows are linearly related to{λ kj }. 104 Figure 5.28: λ 23 variations over 100 realizations, synthetic Case 2 Figure 5.29: λ 53 variations over 100 realizations, synthetic Case 2 105 Figure 5.30:{τ j } distribution at 4 different times, synthetic Case 2. The horizontal axis represents the index j for{τ j }. The maximum and the minimum values of{τ j } are represented with error bars while the mean values are depicted with circle markers. Figure 5.31: τ 1 variations over 100 realizations, synthetic Case 2 106 Figure 5.32: {λ kj } variations over 100 realizations, synthetic Case 2. Figure 5.33: {τ j } variations over 100 realizations, synthetic Case 2 107 were effectively saturation dependent and thus vary with time. Such changes can be fully captured using EnKF as our dynamic estimation procedure. We also evaluated the ability of our estimated EnKF/CRM parameters for forecasting the system. Due to the time-dependence of the parameters we may use the forecast for a short period of time unless the values are stabilized. If a longer forecast scenario is desired, we need to combine the update and forecast procedures. The range of our forecast period typically depends on the variation time scale of our the CRM parameters. We also showed that the estimates for the CRM parameters, even while being time- dependent still captured the underlying geological characteristics of the reservoir to some extent. 5.6 Heterogeneous Reservoirs 5.6.1 Channelized Reservoir (Case 3) Synthetic case 3 consists of 5 injectors and 4 producers over a 2D square area depicted in Fig. 5.34. The grid block sizes and the fluid properties are the same as in Sec. 5.5.1. There is a high permeability channel connecting injection 1 to production 4 withk = 500 md. Permeability in every other grid block is uniformly set to k = 10 md. The injection wells are controlled by water injection rates given in Sec. 5.5.2. Production wells are controlled by bottom hole pressures that vary between 200 and 300 psi (depicted in Fig. 5.35). We have generated a variable 108 Figure 5.34: Synthetic Case 3: Heterogeneous channelized reservoir bottom hole pressure data to fully validate the ability of EnKF to match for all the unknown parameters including{J j } in CRM framework. With fixed bottom hole pressures, the parameters{J j } become redundant and we already exploited this fact in Sec. 5.5.2 where{J i } were removed from CRM control variables altogether unlike the current case. The observed total production and oil rates are shown in figs. 5.36 and 5.37 respectively. Notice that production 4 has a higher total and oil production rates. This is due to the high permeability channel passing through producer 4. The sharp jump in the response of producer 4 is due to the variation of bottom hole pressure. We start with an initial guess distribution of the static variables: λ kj (t 0 ) = 1 4 and τ j (t o ) = 10 days. The initial average values for a j , b j and J j are based on a prior 109 Figure 5.35: Well bottom hole pressure data, synthetic Case 3 Figure 5.36: Observed total production rates for synthetic Case 3 110 Figure 5.37: Observed oil production rates for synthetic Case 3 knowledge of the reservoir which are 10 −14 , 1, 1 bbl day·psi [a j is given in the appropriate units of day, bbl/day and psi] respectively. The relative standard deviation of the initial static variables is set to 50%. The initial guess for the dynamic variable q j (total production) and q o,j (oil production) is set to the actual (real) initial values without uncertainty. The relative noise amplitude used to perturb the observation data for EnKF procedure is 10%. As before, the CRM formulation is run up to 500 days without utilizing EnKF. After this period, we switch on EnKF up to day 2500. We revert back to pure propagation until day 3000 using the resulting estimate for the state vector to validate the forecasting abilities of our estimate for the CRM parameters. Fig. 5.38 depicts the results of EnKF matching for the total production rate. The ensemble members converge almost instantly to the actual observation data 111 Figure 5.38: Matching total production rates for synthetic Case 3. The red curve shows the actual observed data and the black curves show the ensemble realizations. once EnKF is switched on after day 500. After day 2500 when EnKF is switched off, the state vectors of the ensemble still very closely follow the observed data until day 3000 that confirms the predictive power of our estimates. Fig. 5.39 shows similar trends for the oil production rates. We consider the time variation of{λ kj } and the saturation movements for this case as well. Fig. 5.40 shows the saturation contour maps at various times. The saturation movements results in a change in transmissibility and subsequently in the estimates for{λ kj } and{τ j }. To pinpoint this relationship, we similarly dis- play the{λ kj } distribution obtained from EnKF-CRM methodology in Fig. 5.41. 112 Figure 5.39: Matching oil production rates for synthetic Case 3 Clearly the initial realizations (t = 0) do not carry any characteristics of the sys- tem. The variation among ensemble realizations is reduced when EnKF is started (t = 500). Towards the end of the simulation, the uncertainty of{λ kj } over the ensemble members is systematically reduced. In Fig. 5.41, we can identify indices 16 and 18 on the horizontal axis with the largest λ values corresponding to the (k,j) = (1, 4) and (k,j) = (3, 4) respectively. These{λ kj } values are consistent with saturation movements around these wells: Inj1 and Prod4 are exactly located inside the channel and Inj3 is closer to Prod4 compared to Inj1. Fig. 5.42 shows the average values of the final estimate of{λ kj } values in the EnKF/CRM formula- tion at the final time of the simulation. Notice that all injectors show higher{λ kj } values toward the production 4 situated inside the channel. Even injection 2 which 113 Figure 5.40: Water saturation contour map at 4 different times: (a)t = 0; (b) t = 500 days; (c) t = 1500 days; (d) t = 3000 days, synthetic case 3 is immediately outside the higher permeability channel shows a large connectivity to production 4. This is also evident from the saturation movements 5.40. These estimates confirm the overall suitability of the CRM parameter modeling and EnKF estimation procedure for this model. Moving beyond the static parameters, we consider the variation of λ 14 and λ 34 in particular. Figs. 5.43 and 5.44 depict the time-variation of λ 44 and λ 34 over 100 realizations where the average estimates for λ 44 and λ 34 are highlighted. The 114 Figure 5.41:{λ kj } distribution at 4 different time, synthetic Case 3. The horizontal axis is indexed by (j− 1)×m +k corresponding to{λ kj }. The maximum and minimum value of λ for each case are depicted with error bars and the mean is depicted with circle marker. Figure 5.42:{λ kj } distribution, synthetic Case 3. Each arrow represents the mean estimate of{λ kj } at the final time of the simulation and the length of the arrows are linearly related to{λ kj }. 115 most prominent feature in this period is the reduction in the uncertainty (spread) of λ 44 and λ 34 among the ensemble members. Notice that the ensemble parame- ters are constant during the 2500-3000 days where the ensemble members evolve freely without parameter adjustment. The evolving nature of the CRM parameters (through their dependence on saturations) suggests that a final estimate for the parameters might not be reliable beyond a certain time for forecasting unless the reservoir saturation features have been stabilized. In particular, notice that the values of λ 44 and λ 34 converge to a final constant value on and after day 2300. In comparison with Fig. 5.40 it is understood thatλ 44 andλ 34 converge only after the saturation distribution has been stabilized. Any abrupt change in the injection and production patterns imply a sudden variation of{λ kj } and can directly influence the front movement (and the saturation distribution evolution) which is an impor- tant feature of a dynamical characterization of the reservoir. We similarly depicted the{τ j } parameters obtained from EnKF-CRM methodology in Fig. 5.45. Clearly, the uncertainty among{τ j } parameters is reduced systematically. Note thatτ 4 has the smallest value as Prod4 is inside the channel and has a direct connection to all injectors and the smallest time lag. A similar pattern is also displayed in Fig. 5.46 for τ 3 . The final obtained values for{λ kj } and{τ j } have been depicted in Table 5.6. We finally consider{J j } parameters. Fig. 5.47 shows the average values of obtained{J j } parameters from EnKF-CRM methodology at initial time (t=0), 116 P (j = 1) P (j = 2) P (j = 3) P (j = 4) λ 1j (k = 1) 0.03 0.06 0.004 0.9 λ 2j (k = 2) 0.18 0.02 0.16 0.63 λ 3j (k = 3) 0.001 0.066 0.008 0.92 λ 4j (k = 4) 0.1 0.02 0.1 0.76 λ 5j (k = 5) 0.01 0.011 0.19 0.77 τ j (days) 28.3 4.3 10.43 3.5 J j 0.44 0.01 0.29 12.3 a j 2.5× 10 −13 2.06× 10 −9 1.3× 10 −11 7.5× 10 −10 b j 2.05 0.25 1.7 1.22 Table 5.6: Final obtained values for all the unknowns in the EnKF-CRM methodology, synthetic Case 3 Figure 5.43: λ 44 variations over 100 realizations, synthetic Case 3 117 Figure 5.44: λ 34 variations over 100 realizations, synthetic Case 3 Figure 5.45:{τ j } distribution at 4 different times, synthetic Case 3. The horizontal axis represents the index j for{τ j }. The maximum and the minimum values of{τ j } are represented with error bars while the mean values are depicted with circle markers. 118 Figure 5.46: τ 3 variations over 100 realizations, synthetic Case 3 EnKF start day (500), a middle time (1500), and the final time of the simulation (3000). Fig. 5.47 shows a systematic reduction in uncertainty of{J j } parameters as the EnKF procedure is implemented. We notice that Prod4 has the largest productivity index. This is due to the high permeability channel passing through this producer. Productivity index as defined in Eq. 3.4 is the ratio of the variations in the total production rate to the pressure variations near the producers. Assuming a constant value for this parameter in a characterization step may lead to large error. As with the rest of the parameters, we relax the CRM formulation to account for possible changes in this parameter over time. To explore this premise, we depict the variation of J 1 over 100 realizations in Fig. 5.48. It is clearly evident that the estimate for this parameters changes with time until it stabilizes towards the end 119 Figure 5.47:{J j } distribution at 4 different times, synthetic Case 3. The horizontal axis represents the index j for{J j }. The maximum and the minimum values of{J j } are represented with error bars while the mean values are depicted with circle markers. of simulation around day 2000. Compared with variations of{λ kj } and{τ j }, this parameter shows less and/or smoother variations with time. We reiterate our expectation that the variations in the estimates of{λ kj } and {τ j } and{J j } may provide an insight into fluid front movement in reservoirs. Figs. 5.49 and 5.50 depict the variation of the estimate for{λ kj } and{τ j } and{J j } over time respectively for all injectors k and producers j. Table 5.6 depicts the final values of the initially unknown parameters in the EnKF-CRM methodology implemented in this case. In this section, we applied the EnKF/CRM methodology to a channelized sys- tem. In this system we allowed the production well bottom hole pressures to also 120 Figure 5.48: J 1 variations over 100 realizations, synthetic Case 3 Figure 5.49: {λ kj } variations over 100 realizations, synthetic Case 3 121 Figure 5.50: {τ j } variations over 100 realizations, synthetic Case 3 Figure 5.51: {J j } variations over 100 realizations, synthetic Case 3 122 vary and affect the dynamics. This introduces the additional parameters into the CRM formulation, namely the productivity indices. The variation of the CRM pa- rameters with time are more pronounced in this case due to the fluctuations of both the injection rates and producers bottom hole pressures as well as the existence of a preferred direction of transmission in the reservoir. Expectedly, the variations in the CRM parameters are correlated with the saturation movements and can be cap- tured by the EnKF methodology. Approximate location of the high permeability channel can be verified by the final estimates of the CRM parameters. 5.6.2 Peripheral Pattern (Case 4) Synthetic case 4 consists of 12 injectors and 9 producers in a peripheral water flood pattern over a 2D square area depicted in Fig. 5.52. The grid block sizes and the fluid properties are the same as in Sec. 5.5.1. This is a channelized reservoir with a permeability of 100 md in the channel and 20 md in other grid blocks. The injection wells are controlled by water injection rates and are depicted in fig 5.53. Notice the high fluctuations of the injection wells: the injection rates are produced as a random walk with a mean value of 500 bbl/day and the standard deviation of 250 bbl/day. Production wells are controlled by fixed bottom hole pressures of 250 psi. The observed total production and oil rates are shown in Figs. 5.54 and 5.55 respectively. Notice that Prod 4 located inside the channel has a higher total and oil production rates compared with the rest of the producers. In this 123 Figure 5.52: Synthetic Case 4: Heterogeneous peripheral pattern section, we focus on the forecasting properties of our method, given that a possible change in transmissibility across the reservoir can change the parameters in the CRM methodology (e.g. λ kj , τ j ,...) as a function of time; such variations will reduce the prediction power of the CRM methodology. In order to have better and more reliable prediction properties, one needs a model that describes the variation in the parameters or alternate between prediction and characterization stages. This predict-update scenario leads to an expectedly better results at the cost of more frequent measurement data. We start with an initial guess distribution of the static variables: λ kj (t 0 ) = 1 9 and τ j (t o ) = 1 day. The initial average values for a j , b j , and J j are based on a prior knowledge of the reservoir given by (a j ,b j ,J j ) = (10 −5 , 1, 1( bbl day×psi )). The 124 Figure 5.53: Injection rates used for synthetic Case 4 Figure 5.54: Observed total production rates for synthetic Case 4 125 Figure 5.55: Observed oil production rates for synthetic Case 4 relative standard deviation of the initial static variables is set to 50%. The initial guess for the dynamic variable q j (total production) and q o,j (oil production) is set to the actual (real) initial values without uncertainty. The relative noise amplitude used to perturb the observation data for EnKF procedure is 10%. In comparison with the previous cases, we increase the duration of the simulation to 4000 days to increase the prediction period. The CRM formulation is then run for the ensemble members up to day 500 without EnKF. After this period, we switch on EnKF up to day 2500. We revert back to pure propagation until day 4000 using the resulting estimate for the state vector to validate the forecasting abilities of our estimate for the CRM parameters. Fig. 5.56 depicts the results of EnKF matching for the total production rate. The ensemble members converge almost instantly to the actual observation data 126 Figure 5.56: Matching total production rates for synthetic Case 4 once EnKF is switched on after day 500. Note that after day 2500 when EnKF is switched off the state vectors of the ensemble still very closely follow the observed data until day 4000 that confirms the predictive power of our estimates. Fig. 5.57 shows similar trends for the oil production rates. Sensitivity Analysis In this section, we provide another sensitivity analysis for the predictive efficiency as a function of the simulation constraints and the initial state distribution. We focus on initial model estimates, observable noise, and ensemble population size. For each case a plot of RMSE defined in Eq. 5.3 is used to test the predictive power of our proposed method. For all settings considered here, the error after day 2500 127 Figure 5.57: Matching oil production rates for synthetic Case 4 (forecasting period) will gradually increase. This suggests that an update-forecast scheme should be used if reliable results for the prediction scenario is needed. Fig. 5.58 depicts the variation of RMSE of total production rate as a function of simulation time for different ensemble populations N e , varying between 20 and 200. The other parameters of our model in this simulation are the same as before in case 4. Notice that RMSE is plotted starting at day 500 when EnKF procedure kicks in. Clearly, RMSE decreases with increasing the ensemble size. Similar to the previous sensitivity analysis, when the ensemble population is too small (N e = 20), the results are not satisfactory as expected in Monte Carlo methods. Increasing the ensemble size allows for smoother and finer variations among the ensemble mem- bers. Notice however that practically the improvement by increasing the ensemble 128 Figure 5.58: RMSE variation of total production rate for various ensem- ble populations, synthetic Case 4 population is limited. This is observed in Fig. 5.58 where the RMSE values for N e = 100 and N e = 200 are practically of the same order. We fixed the ensemble size at N e = 100 as an optimal value to avoid using excessive running times and computer memory. In our simulation, the initial estimate is described by a normal distribution of {λ kj },{τ j },{a j }, and{b j } distribution characterized by an average and a propor- tional standard deviation σ k . We first vary the relative standard deviation of the static variables (σ k ) that describes the initial range of the ensemble estimate. In Fig. 5.59, we graph the error as a function of time for different values of the initial standard deviation of the ensemble. Notice that for certain values of σ k the space span by the initial estimate ensemble is too narrow [e. g. σ k = 0.01] to include the true field values. As we increase σ k the error shrinks for the updating period 129 Figure 5.59: RMSE variation of total production rate for various σ k values, synthetic Case 4 t∈ (500, 2500) which guarantees that the estimate space of parameters is larger and may more likely include the true solution. For this reason, care must be taken for the choice of the initial estimate for the parameter space. Just as a narrow guess might miss the mark with EnKF, a guess that is too broad will not yield meaningful results either. In contrast, in the forecasting period t∈ (2500, 4000) the error will gradually increase for all the standard deviations. Finally, we consider the sensitivity of EnKF/CRM procedure to the obser- vation noise (σ obs ). Fig 5.60 depicts the results of RMSE with time for σ obs ∈ {0.05, 0.1, 0.5}. In contrast to the previous cases, increasingσ obs results in the error increasing. Effectively, a noisy observation data is associated with the uncertainty of the measurements. A small perturbation on the observed values is required for 130 Figure 5.60: RMSE variation of total production rate for various σ obs values, synthetic Case 4 EnKF to provide a constant small variance on the updated ensemble parameters when the procedure converges to the actual parameters. A larger perturbation generated by higher σ obs values on the other hand may cause the solution to di- verge from the true values. Notice the error in the forecasting period increases dramatically with increasing σ obs . In this section we applied the EnKF/CRM to a peripheral water flood with a high permeability channel. This case evaluates the EnKF/CRM in a situation with a large number of wells in a complex pattern. Close matches for the total and oil production model were obtained. We also analyzed the forecast period for sensitivity to the numerical parameters of the procedure: a plot of the forecast error was produced for different ensemble populations, initial uncertainties, and 131 observation noise amplitudes. The sensitivity results suggest that there are practical optimal values for the numerical parameters of the procedure. 5.7 Optimization Procedure 5.7.1 Streak Reservoir (Case 5) In this section, we evaluate EnOpt/CRM procedure for the synthetic reservoir shown in Fig. 5.61. We proceed with three different scenarios. In scenario (a), we use EnKF/CRM to validate the ability of the CRM modeling technique to obtain the reservoir directionality. In scenario (b), we use EnOpt/CRM with the assump- tion that all the reservoir geologic parameters are known and thus focus solely on optimizing injection rates in order to increase the NPV. This case tests the abil- ity of the EnOpt approach to increase the sweep efficiency when CRM parameters are known. In scenario (c), we apply EnOpt/CRM for both characterization and optimization. The synthetic field 5 consists of 5 injectors and 4 producers over a 2D square area depicted in Fig. 5.61. There are two streaks of higher permeability regions. The permeability is 5 md everywhere except for the streaks where it is 1000 md for the injector 1-producer 1 area and 500 md for the injector 3-producer 4 area. The porosity is uniform throughout the whole region and is set to 0.18. The synthetic field dimensions and the grid sizes are 40× 40× 6 ft and 31×31×5, respectively. 132 Figure 5.61: Synthetic Case 5: Heterogenous streak model I01, . . . , I05 denote the injection wells, while P01,. . . ,P04 denote the production wells. The oil water mobility ratio is equal to 1 and the oil, water, and rock compressibilities are 5×10 −6 , 1×10 −6 , and 1×10 −6 1 psi , respectively. The injection wells are controlled by water injection rates. The reservoir description is given in Table 5.7. The actual reservoir response is simulated by using Computer Modeling Group (CMG) reservoir simulator. 5.7.1.1 Scenario (a): EnKF with CRM framework Fig. 5.62 and 5.63 depict the actual injection rates of all the injection wells and the bottom hole pressures for all the production wells respectively and provide the time- dependent driving terms in the CRM dynamical systems. Notice that injectors 1 133 Parameter Value oil price v o 70$/bbl water disposal cost v w 5$/bbl. total time 1000 days control steps 20 discount rate r 0.1 per year number of realizations N e 100 Table 5.7: Various settings for synthetic Case 5. and 2 have higher water injection rates and that the bottom hole pressures among the production wells varies between 200 psi and 300 psi. Figs. 5.64 and 5.65 depict the total rates and oil rates for all the production wells respectively and form the observable data obtained from the reservoir simulator. Clearly, producer 1 and 4 that are near the two streaks of higher permeability regions have higher values of total and oil production rates. We have generated a variable bottom hole pressure to fully validate the ability of the EnKF to match for all the unknown parameters including{J j } in the CRM framework. We reiterate that with fixed bottom hole pressures, the parameters{J j } are redundant and we exploit this fact in the optimization section to remove bottom hole pressure from the CRM control variables altogether for simplicity. We propagate the CRM based ensemble up to 500 days without EnKF. After this period, we switch on EnKF up to day 2500. We revert back to pure propagation until day 3000 with the resulting estimate for the state vector so that the forecasting abilities of our estimate can be validated. This procedure is similar to cases 1 through 4. 134 Figure 5.62: Injection rates used for synthetic Case 5; scenario (a) Figure 5.63: Bottom hole pressure used for synthetic Case 5; scenario (a) 135 Figure 5.64: Observed total production rate for synthetic Case 5; sce- nario (a) Figure 5.65: Observed Oil production rate for synthetic Case 5; scenario (a) 136 Figure 5.66: Matching total production rates for synthetic Case 5; sce- nario (a) Fig. 5.66 depicts the results of EnKF matching for the total production rate. The red curve denotes the actual observed data and the black curves refer to the ensemble trajectories. Clearly, the ensemble members converge to the actual ob- servation data once EnKF is switched on after day 500. We also note that after day 2500 when EnKF is switched off the state vectors of the ensemble still closely follow the observed data until day 3000. Fig. 5.67 shows similar trends for the oil production rates. We consider the correlation between the saturation movements and the time variation of{λ kj } as usual. Fig. 5.70 shows saturation contour maps for the initial time (day 0), EnKF starting day (500), mid-time (1500), and the final time of the simulation (3000). Initially, the saturations for all the grid blocks are uniform 137 Figure 5.67: Matching oil production rates for synthetic Case 5; scenario (a) throughout the system at 20%. At day 500, five distinct regions can be identified in the model. Two elongated streaks of higher permeability occur between I01 and P01 and also between I03 and P04. Due to these higher streaks, I01 and I03 support only P01 and P04. Notice that water front reaches P01 earlier than P04. Also notice the extended high-water region around I02 supporting both P01 and P03. This region is larger than the high-water region near I04 and I05 as injector I02 starts with a higher rate. Similarly another extended region appears around I05 that mostly supports P03 and P04 and another surrounding I04 that supports P02 and P04. Notice that the area surrounding I01 remains largely un-swept due to the higher permeability streak. 138 Figure 5.68:{λ kj } distribution at 4 different times: (a) t = 0 ; (b) t = 200 days; (c) t = 600 days; (d) t = 1000 days, synthetic Case 5; scenario (a) In order to relate the interwell connectivities to the saturation movements we similarly depicted the{λ kj } distribution obtained from EnKF-CRM methodology in Fig. 5.68. Notice that the higher permeability regions connecting I01 with P01 and I03 with P04 do absorb large amounts of water that translates to higher values (longer arrows) in λ 11 and λ 34 distributions. The higher values for λ 54 and λ 21 in the regions connecting the corresponding wells are also consistent with variations in the water saturation in these regions. This confirms our initial conjecture that in the limit of large variations in reservoir saturation, a time-varying (possibly on a phase by phase basis) set of parameters for CRM model can better model the reservoir at the cost of adding to the complexity of the model. In particular, we follow the variation of λ 11 in the reservoir. 139 Fig. 5.69 shows the time-variation ofλ 11 over 100 realizations where the average estimate for λ 11 is also highlighted. The EnKF procedure is applied in the 500- 2500 days interval and the most prominent feature in this period is the reduction in the uncertainty (spread) of λ 11 among the ensemble members. The ensemble parameters are constant before and after this period. The evolving nature of the CRM parameters (through their dependence on saturations) suggests that a final estimate for the parameters might not be reliable beyond a certain time for fore- casting unless the reservoir saturation features have been stabilized. In particular, notice that the values of λ 11 converge to a final constant value on and after day 2100. In comparison with Fig. 5.70 for saturations, λ 11 converges only after the saturation distribution has been stabilized. Notice also the sharp jump in λ 11 at around day 1600 related to the sharp decrease in the injector 1 and 3 patterns at the time. We can also identify a similar sharp jump in the producer 1 and 4 as well. Such abrupt changes in the injection and production patterns imply a sudden variation of{λ kj }. The relationship between variations inλ and the injection rates is mediated by the saturation distribution evolution and is ultimately related to the front movement as an important feature of a dynamical characterization of the reservoir. A similar structure is also observed in Fig. 5.71 for τ 3 . Compared with Fig. 5.69 for λ 11 , the variations in this graph are smoother. Figs. 5.72 and 5.73 depict the variation of the estimate for{λ kj } and{τ j } over time respectively for all injectors k and producers j. 140 Figure 5.69: variation of λ 11 estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. Along with the CRM parameters, we have also estimated the phenomenological parametersa j andb j that characterize the oil production. The graphs in Figs. 5.74 and 5.75 depict the variation of the estimates for{a j } and{b j } over time. Table 5.8 displays the final results for all the variable in the CRM framework. Notice that the values of τ 1 and τ 4 are the smallest. This is due to the higher permeability regions around producers 1 and 4 that shortens the fluid travel time along the streaks. The fact that the estimate for τ 1 is less than the estimate for τ 4 is caused by the relative difference in the permeability values around production wells 1 and 4. The results for{J j } similarly depict that producers 1 and 4 have the largest well productivity indexes as expected. 141 Figure 5.70: Water saturation contour maps at 4 different times; a) t = 0; b) t = 200 days; c) t = 600 days; d) t = 1000 days, synthetic Case 5; scenario (a) 142 Figure 5.71: variation of τ 3 estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. Figure 5.72: variation of {λ kj } estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. 143 Figure 5.73: variation of{τ j } estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. P 1 (j = 1) P 2 (j = 2) P 3 (j = 3) P 4 (j = 4) λ 1j (k = 1) 0.9 0.008 0.003 0.08 λ 2j (k = 2) 0.72 0.01 0.16 0.09 λ 3j (k = 3) 0.28 0.001 0.0001 0.7 λ 4j (k = 4) 0.0004 0.18 0.02 0.79 λ 5j (k = 5) 0.04 0.004 0.2 0.75 τ j (day) 7.3 19.5 14.6 12.2 J j (RB/(day·psi)) 2.2 0.01 0.19 2.4 a j (1× 10 −7 ) 0.06 0.00001 0.0005 0.13 b j 1.4 2.2 1.8 1.3 Table 5.8: Final estimates (over 100 realizations) of all the unknowns in synthetic Case 5; scenario (a) 144 Figure 5.74: variation of log(a j ) estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. 145 Figure 5.75: variation of {b j } estimates over time for synthetic Case 5; scenario (a). Black curves denote the ensemble members and red denotes the ensemble average. 146 5.7.1.2 Scenario (b): known geology The accuracy of the matches in scenario (a) (and the other cases considered in this chapter) suggests that EnKF-CRM combination may successfully estimate the unknown parameters of CRM for the reservoir studied in many cases. In scenario (b), in contrast, we assume that the CRM parameters are already characterized precisely and proceed solely with the problem of optimization. ACMG simulator is used to propagate the data forward in time as a substitute for the real dynamics at every step using the control parameters at the time. The production wells are controlled by a constant bottom hole pressure of 250 psi. This renders the parameters{J i } redundant and removes the bottom hole pressure from the control variables. As in every realistic situation, the main challenge of the optimization procedure is to incorporate constraints without which the optimized solutions would be trivial. In this case, we constrain the total water injection rate ( P m k=1 i k (t)) at 5000 rb/day at all times. The fact that the CRM parameters are all known eliminates step 3 from the implementation of EnOpt/CRM in Sec. 4.2.3 so that the state vectors{y j } are all identical and known at each time. Each evaluation of the objective function in Eq. (4.4) requires a single simulation model run to compute g(x l+1 ,y ref ). This is performed through the CRM while the actual injection rates are applied to the reservoir (being simulated). Fig. 5.76 depicts the injection rates vs. the control steps for various optimization iteration steps. Notice in particular the trends associated with each injection rate. 147 Figure 5.76: Injection rate vs. control step vs. iteration step in synthetic Case 5; scenario (b). The thick purple curve depicts the final optimized rate obtained for each injector. Clearly, the injection wells 1 and 3 show a decreasing trend whereas the injection wells 2, 4 and 5 show an increasing trend. The results for the optimized injection rates are shown in Fig. 5.77. They suggest settingi 1 (t) andi 3 (t) lower than the rest of the injection rates. This is motivated by the higher permeability regions between Inj 1-Prod 1 and Inj 3-Prod 4. In order to get a higher NPV, the injection rates near the higher permeability regions need to be reduced to allow for a more even sweep of the reservoir. This can in fact be validated by the saturation movements (5.70) and the un-swept area around injector 1. Fig. 5.78 shows the results for the optimized total cumulative oil and water productions (top and bottom left, resp.) and oil and water production rates (top 148 Figure 5.77: Summary of optimized injection rate in synthetic Case 5; scenario (b). and bottom right) for the entire field. The results are contrasted with a base case in which all injection rates are fixed at 1000 rb/day. The results clearly demonstrate a higher cumulative oil production and a lower water production after optimization. 5.7.1.3 Scenario (c): unknown geology In this scenario, the geologic model is a priori given with a large uncertainty. We seek not only to characterize the reservoir but also to optimize its NPV over a period of 1000 days. After day 1000, we propagate the system (with the finalized injection rates) to evaluate the forecasting ability of the estimated parameters. Fig. 5.79 depicts the results for the distribution of{λ kj } at 4 different times through the propagation of the system. Notice the initial range of the distribu- tion: for example for λ 11 the initial range is in (0, 0.95). Clearly, the initial{λ kj } 149 Figure 5.78: Results of oil and water production compared with the base case for synthetic Case 5; scenario (b). 150 distribution does not carry any of the characteristics of the true field but already at 200 days (top right plot), most of the characteristics of the geological model is captured in{λ kj } values: λ 11 andλ 34 (points 1 and 18 on the horizontal axis) have the highest values as expected from the high permeability streaks that connect the corresponding well pairs. This distinction is highlighted even more so, towards the end of the simulation (bottom left, bottom right). In fact by the end of the run, we may even deduce λ 11 >λ 34 as suggested by the relative values of the permeability at the streaks connecting the corresponding pairs. The distribution ofτ is similarly plotted in Fig. 5.80. We have chosen the range of the initial distribution of{τ j } based on a prior estimate of the model. As we progress towards later times, we notice smaller estimates forτ 1 andτ 4 corresponding to higher permeability streaks. Similarly, τ 2 and τ 3 have the highest values of time lag. Both Figs. 5.79 and 5.80 clearly depict the decreasing trend in the uncertainty of the main CRM parameters τ and λ in this scenario. Fig. 5.81 depicts the variation of injection rates over time inside the optimiza- tion loop of the procedure right after the second control step has been reached in propagation of the reservoir. Notice the decreasing/increasing trends in various in- jection wells. These trends are in general agreement with the geological description of the reservoir, as required, and are decreased/increased for an even sweep of the reservoir and higher NPVs. Fig. 5.82 depicts the final optimized results in this scenario over time. 151 Figure 5.79: {λ kj } distribution at 4 different data times for synthetic Case 5; scenario (c). The position (j− 1)×m +k on the horizontal axes denotes{λ kj }. The blue markers represent the mean values of{λ kj } from the ensemble, while the standard deviation is represented with error bars. 152 Figure 5.80: {τ j } distribution at 4 different times for synthetic Case 5; scenario (c). The position j on the horizontal axes corresponds to{τ j }. The blue markers represent the mean values of{τ j } from the ensemble, while the standard deviation are represented with error bars. 153 Figure 5.81: Injection rate vs control step vs. iteration index for syn- thetic Case 5; scenario (c). The thick green curve denotes the final optimized injection rate. Figure 5.82: Summary of optimized injection rate in synthetic Case 5; scenario (c). 154 Fig. 5.83 depicts the data for total oil and water productions during the opti- mization procedure for scenarios (b) [known geology] and (c) [unknown geology], along with the base case defined in Sec. 5.7.1.2 with the constant injection rates. We observe that the cumulative oil production in scenario (c) is significantly higher than the base case and the production results for scenarios (b) and (c) are close, which suggests that the characterization goal is also achieved in scenario (c). We further examine the quality of the characterized CRM parameters{λ kj },{τ j },{a j }, and{b j } for forecasting the reservoir in the forecast phase (after day 1000) using the optimized controls with the initial distribution of these parameters and the final distribution of these values obtained in our procedure. Figs. 5.84 and 5.85 depict the forecasted observables for the initial and final estimates of the CRM parameters, respectively. The red curves denote the actual production data from reservoir simulator while the black curves denote the state vectors propagated using CRM. Clearly, the responses with the initial distribution show a high uncertainty. The final estimates on the other hand have a much smaller spread and provide a better match for the actual reservoir production data. This is yet another indicator that the characterization phase was successful in scenario (c) and could be used to forecast the reservoir for a significant time period. In this section, we implemented EnOpt/CRM methodology for a system with two streaks of higher permeability regions. In implementing this method, we fol- lowed two different scenarios. In the first scenario, we assumed that the reservoir 155 Figure 5.83: Comparison of the field production responses for different scenarios: closed loop optimization (scenario (c), red), optimization with known geology (scenario (b), black), and the base case (green). 156 Figure 5.84: Production match from the initial ensemble as a function of time; scenario (c) Figure 5.85: Production match from the final ensemble as a function of time; scenario (c) 157 parameters are known and focused on the optimization of the injection patterns for increasing the net present value of the reservoir. The results show that cumu- lative oil production is higher and the cumulative water production is lower than the case where the injection rates are equally distributed among wells. In the sec- ond scenario, we assumed that the reservoir geologic parameters are unknown and proceeded with both the characterization of the reservoir and optimization of the injection rates to increase the net present value. The final obtained CRM param- eters indicate the location of the two streak in this case. The results from the optimization procedure show an increased net present values. The results obtained in this section verify the ability of the EnOpt/CRM framework in order to both characterize and optimize the reservoir. 5.7.2 Channelized Reservoir (Case 6) In this section, we repeat EnOpt/CRM method for a synthetic reservoir shown in 5.6.1 and focus on the optimization of injection rates. We proceed with sce- narios involving partial and accurate knowledge of the reservoir. In scenario (a), EnOpt/CRM framework is used along with a fully known set of reservoir geologic parameters and the focus is on optimization of injection rates in order to increase the NPV. This scenario tests the ability of the EnOpt approach to increase the sweep efficiency when the CRM parameters are known. In scenario (b), we apply EnOpt/CRM for both characterization and optimization. The reservoir geologic 158 feature and fluid properties along with optimization settings are the same as in Sec. 5.6.1 and 5.7 unless otherwise specified. This case contains a set of sensitivity analysis with respect to the parameters involved in optimization set up. 5.7.2.1 Scenario (a): known geology Similar to the previous case (case 5, scenario b) we assume that the CRM param- eters are already characterized precisely and proceed solely with the problem of optimization. The fact that the reservoir dynamics is known, eliminates step 3 from the implementation of EnOpt in Sec. 4.2.3 so that the state vectors{y j } are all identical and known at each time. Each evaluation of the objective function Eq. (4.4) requires a single simulation model to compute g(x l+1 ,y ref ). Fig. 5.86 depicts the injection rates vs the control step for various optimization iteration steps. Notice in particular the trends associated with each injection rate. At the end of iterations, we are left with a set of controls for the reservoir for optimal NPV. For example, consider I01 located inside a high permeability channel. The optimized injection rate for this well (thick green curve in top-left subplot in Fig. 5.86) suggests reducing the injection rate initially, presumably until all production rates are saturated and then increasing it afterward for a maximum output. I03 located near the channel also shows a similar pattern. Summary of optimized injection rates appears in Fig. 5.87. Notice the high injection pattern associated with both I02 and I05 initially. These injectors are located far from the high 159 Figure 5.86: Injection rate vs. control step vs. iteration step in synthetic Case 6; scenario (a). The thick green curve depicts the final optimized rate obtained for each injector. permeability channel and can be used to sweep large parts of the reservoir before the production wells are saturated. The decreasing trends associated with I02 and I05 are in accordance with increasing trend for injectors I01 and I03. It should be noted that due to non-uniqueness in both characterization and optimization, it may be possible to obtain other control trends with possibly higher NPV in other runs of the procedure. Fig. 5.88 shows the optimized total cumulative oil and water productions (top and bottom left, resp.) and oil and water production rates (top and bottom right) for the entire field. The optimized results are clearly contrasted with a base case in which all injection rates are fixed at 1000 rb/day. The results demonstrate a 160 Figure 5.87: Summary of optimized injection rates in synthetic Case 6; scenario (a). higher cumulative oil production and a lower water production after optimization resulting in a higher NPV. Sensitivity analysis We analyze the effects of the parameters specific to the implementation of the EnOpt procedure namely the variance (s C ) and temporal correlation (γ) of the noise added to the control variables inside the optimization algorithm loop. See Sec. 4.2.3. In particular let i j (k) denote the injection control variable for the j-th injection well index and k-th time optimization step. Let i max and i min denote the 161 Figure 5.88: Results of oil and water production compared with the base case for synthetic Case 6; scenario (a). The black curve shows the results of known geology scenario and the green curve shows the base case. 162 maximum and minimum well constraints respectively. The distribution of the noisy control variables i j (k) is a Gaussian distribution with the following properties: hi j (k)i = r∈ [i min ,i max ] chosen uniformly randomly hi j 1 (k 1 )i j 2 (k 2 )i−hi j 1 (k 1 )ihi j 2 (k 2 )i = s c δ j 1 ,j 2 e −γ(k 2 −k 1 ) 2 (i max −i min ) where δ j 1 ,j 2 is the Kronecker’s delta (δ ij = 1 if and only if i =j and 0 otherwise). Notice that the added noise is independent for well indices and correlated for time indices. Figs. 5.89 and 5.90 depict the final NPV value from optimization with respect tos C andγ. Notice that albeit the improvement is relatively small (1%), choosing an optimal value fors C (around 5×10 4 in this case) andγ (around 10 −1 ) can lead to a better performance of EnOpt procedure. Clearly choosing a zero spread (s C = 0) will lead to no sampling in the control variable space and conversely a spread that is too large will simply break down the perturbative assumptions in the gradient search method. Thus an optimal value for s C must exist. A temporally smooth perturbation (γ = 0) leads to a wide sampling of the control variables, while a completely temporally uncorrelated (γ =∞) set of control variables also breaks down the perturbative assumptions of the gradient and thus an optimal value of γ must exist. 163 Figure 5.89: Sensitivity Analysis: Optimized profit vs control variable noise variance in Case 6; scenario (a) Figure 5.90: Sensitivity Analysis: Optimized profit vs control variable noise temporal correlation in Case 6; scenario (a) 164 5.7.2.2 Scenario (b): unknown geology In this scenario, the geologic model is a priori given with a large uncertainty through a wide initial guess of CRM parameters. A CMG simulator is used to provide the actual observed data incorporating the injection rates output by the optimization procedure. The initial estimates and various settings specific to the EnKF/CRM procedure are identical to case 3 in Sec. 5.6.1. Fig. 5.91 depicts the results for the distribution of{λ kj } at 4 different times through the propagation of the system. Notice the initial range of the distribution: for example for λ 11 the initial range is in (0, 0.6). Clearly, the initial{λ kj } distri- bution does not carry any of the characteristics of the true field but already at 200 days (top right plot), most of the characteristics of the geologic model is captured in{λ kj } values. The final distribution is similar to Sec. 5.6.1 e. g. λ 16 ···λ 20 that connect the injectors to Prod. 4 have the highest values. This is predicted as Prod. 4 is inside the high permeability channel. The distribution of τ is similarly plotted in Fig. 5.92 and is numerically similar to the results of Sec. 5.6.1 (case 3); e.g. τ 4 has the lowest time lag as Prod. 4 is inside the high permeability channel. Both Figs. 5.91 and 5.92 clearly depict the decreasing trend in the uncertainty of the main CRM parameters τ and λ in this scenario. Fig. 5.93 depicts the variation of injection rates over time inside the optimiza- tion loop of the procedure right after the second control step has been reached in 165 Figure 5.91:{λ kj } distribution at 4 different time, synthetic Case 6; sce- nario (b). The horizontal axis is indexed by (j−1)×m+k corresponding to {λ kj }. The maximum and minimum value of λ for each case are depicted with error bars and the mean is depicted with circle marker. 166 Figure 5.92: {τ j } distribution at 4 different times, synthetic Case 6; scenario (b). The horizontal axis represents the index j for{τ j }. The maximum and the minimum values of{τ j } are represented with error bars while the mean values are depicted with circle markers. 167 Figure 5.93: Injection rate vs. control step vs. iteration step in synthetic Case 6; scenario (b). The red green curve depicts the final optimized rate obtained for each injector. propagation of the reservoir. Fig. 5.94 depicts the final optimized results in this sce- nario over time. Notice the decreasing/increasing trends in various injection wells in general agreement with the geological description of the reservoir as required to be reduced/increased for an even sweep of the reservoir and higher NPVs. Notice the similar trends associated with final optimized rates obtained in Figs 5.94 for the high uncertainty and in Fig. 5.87 obtained for the known geology case (low uncertainty). Fig. 5.95 depicts the data for total oil and water productions during the op- timization procedures for known geology case (a) and unknown geology case (b) along with the base case defined in Sec. 5.7.1.2 with constant injection rates. We 168 Figure 5.94: Summary of optimized injection rates in synthetic Case 6; scenario (b). observe that the optimized cumulative oil production is significantly higher than the base case with scenarios (a) and (b) being close; this implies that the charac- terization goal is also achieved in scenario (b) as far as the optimization of NPV is concerned. We further examine the quality of the characterized CRM parameters{λ kj }, {τ j },{a j },{b j } for replicating the reservoir behavior in a forecast procedure using the optimized injection pattern (See Fig. 5.94) from day 1 with the initial distri- bution of these parameters and the final distribution of these values obtained in our procedure. Figs. 5.96 and 5.97 depict the forecasted observables for the initial and final estimates of the CRM parameters, respectively. The red curves denote the actual production data from reservoir simulator while the black curves denote 169 Figure 5.95: Results of oil and water production compared with the base case for synthetic Case 6; scenario (b). The (black,red) curves show the results of (known geology,unknown) scenarios and the green curve shows the base case. 170 Figure 5.96: Production match from the initial ensemble as a function of time in Case 6; scenario (b) the state vectors propagated using CRM. Clearly, the responses with the initial distribution show a high uncertainty. The final estimates on the other hand have a much smaller spread and provide a better match for the actual reservoir production data. This is yet another indicator that the characterization phase was successful in scenario (b). In this section, we have implemented the full EnOpt/CRM methodology for a channelized system. We followed two different scenarios. In the first scenario, we assume that the reservoir geological parameters are known and focused on the optimization of the injection patterns in order to increase the net present value of the reservoir. This scenario resulted in higher cumulative oil production rates and lower cumulative water production than a base neutral case where the injection 171 Figure 5.97: Production match from the final ensemble as a function of time in Case 6; scenario (b) rates are equally distributed among wells. In the second scenario, we assume that the reservoir geologic parameters are unknown and focus on both characterization of the reservoir and the optimization of the injection rates in order to increase the net present value. The final obtained CRM parameter clearly identify the direction of the channel in this case with respect to the wells. The results of the optimization in this scenario also increases the net present value of the reservoir compared to the base case. The accuracy of the results obtained in this section testify the ability of the EnOpt/CRM framework in order to both characterize and optimize the reservoir. 172 5.8 Application to Real field In this section, we apply the EnKF/CRM procedure to a real reservoir setting in California. We focus on two parts (small, large) in a section of the reservoir which is under waterflooding. We thus analyze EnKF/CRM method for matching total production rates and evaluate the forecasting power of our estimated characteriza- tion parameter of the reservoir. We also show that the obtained values for{λ kj } are in general agreement with the fracture orientation in this field. 5.8.1 Small pattern We first apply EnKF/CRM to a small pattern area consists of 6 injectors and 6 producers as shown in Fig. 5.98. The blue (green) markers depict the injection (production) wells. The wells have a line drive pattern in this small section. Note that each injection well may have 1 to 3 perforation interwalls. Each perforation interwall is accounted as a single injection well with an associated index in{λ kj } values. In this way, we distinguish between the perforation interwalls and their effects on the producers. Section 1 contains 6 injection wells with 2 perforations. This allows us to model the EnKF/CRM method with 12 separate injectors. Fig. 5.99 depicts the actual injection rates. In processing the actual rates, we found several dates in which the data was either unavailable or unacceptable. We filtered the missing/unacceptable data by a moving average method with a period of a month. The actual production responses have been depicted in Fig. 5.100. 173 Figure 5.98: Small Pattern, Real field Figure 5.99: Injection rates, Small Pattern 174 Figure 5.100: Observed total production rates, Small Pattern The CRM dynamical equations used in previous cases throughout this chapter were for a balanced capacitance model in which the total amount of water injected were equal to the total production rates. In the case of an unbalanced capacitance model we need to add another (initially unknown) parameter to the CRM equations. This lack of mass balance can be due to various reasons and in particular in this case is due to the fact that we consider only part of the pattern to work with and ignore the rest of the field. In general the injectors outside the selected area need to be considered as well if close matches are desired. Fig. 5.101 depicts the results of the cumulative water and cumulative total productions for section 1 along with the total unbalanced difference in the rates. Previous works on the CRM subject assumed that this difference is fixed with time [[YGJL06],[SZKL07]]. The EnKF/CRM procedure also allows us to consider a difference (q C ) that varies 175 Figure 5.101: Cumulative water injected-Cumulative production, Small Pattern (slowly) with time . We thus use Eq. 3.7 for propagating the ensemble members in time. We start with an initial guess distribution of the static variables: λ kj (t 0 ) = 1 6 , τ j (t o ) = 1 day and q c = 10 bbl/day. The relative standard deviation of the initial static variables is set to 50%. The initial guess for the dynamic variable q j (total production) is set to the actual (real) initial values without uncertainty. The relative noise used to perturb the observation data for EnKF procedure is 10%. With all the elements in place, we applied EnKF/CRM progressively in time in accordance with the above described set of driving terms and observed quantities. To emphasize the large spread of our initial estimate, we propagate the CRM models up to 200 days without utilizing EnKF. After this period, we switch on EnKF up to day 1000. We revert back to pure propagation until day 1200 with the resulting 176 estimate for the state vector so that the forecasting abilities of our estimate can be validated. Fig. 5.102 depict the results of EnKF matching for the total production rate. The red curve shows the actual observed data and the black curves depict the ensemble realizations. The ensemble members provide the actual observation data once EnKF is switched on after day 200. After day 1000 once EnKF is switched off, the state vectors of the ensemble still follow the observed data until day 1200. To make this numerically precise, we have considered the relative RMSE (5.3) that measures the deviation of the observed quantities of the ensemble from the actually observed data. Fig. 5.103 depicts RMSE of the total production rates as a function of time. Clearly the error increases initially until day 200 on which we start the EnKF/CRM procedure. Subsequently, the error is systematically reduced until day 1000 at 6%. After day 1000, the forecast error is slightly larger at 18%. In order to relate the interwell connectivities to the field fracture orientation, we have depicted the{λ kj } distribution from the EnKF-CRM methodology in Fig. 5.104 where the position of the wells in the reservoir is approximately identified. Each arrow’s length is proportional to the{λ kj } and points along the (k,j) line. To have a better understanding of the reservoir directionality, we have identified the vectors corresponding to the maximum{λ kj } values for each injector in Fig. 5.105. The direction of the arrows suggests a trend from the northeast to southwest in accordance with field fracture orientation estimate from geological measurements. 177 Figure 5.102: Obtained match for the total production rates, Small Pattern Figure 5.103: RMSE variation of total production rate, Small Pattern 178 Figure 5.104: {λ kj } distribution; Small Pattern. Arrow lengths represent the mean estimate of{λ kj } at the final time with a direction along (k,j). Figure 5.105: Maximum{λ kj } distribution; Small Pattern. Each arrow represents the length and direction of the maximum{λ kj } at each injector 179 Figure 5.106: Large pattern, Real field 5.8.2 Large Pattern In this section we consider a larger section consists of 19 injectors and 19 producers as shown in Fig. 5.106 . The blue (green) markers depict injection (production) wells. The wells have a line drive pattern and each injection well has 2 perforation interwalls. This forces us to use 38 separate injectors in the EnKF/CRM method- ology. Fig. 5.107 depicts the actual injection rates. We have removed and filtered the missing/unacceptable data using a moving average filter with a period of one month. The actual production response is depicted in Fig. 5.108. Similar to the smaller pattern, we have an unbalanced capacitance model in which the total amount of water injected were not equal to the total production rates. Fig. 5.109 shows cumulative water and cumulative total productions for the larger pattern along with the total difference in the rates. To account for 180 Figure 5.107: Injection rates, Large Pattern Figure 5.108: Observed total production rates, Large Pattern 181 Figure 5.109: Cumulative water injected-Cumulative production, Large Pattern the difference we have introduced a (possibly time-dependent) extra parameter q c according to Eq. 3.12. We start with an initial guess distribution of the static variables: λ kj (t 0 ) = 1 19 , τ j (t o ) = 1 day and q c = 100 bbl/day. The relative standard deviation of the initial static variables is set to 50%. The initial guess for the dynamic variable q j (total production) is set to the actual (real) initial values without uncertainty. The relative noise used to perturb the observation data for EnKF procedure is 10%. As in the case of the small pattern, The EnKF procedure is performed between the days 200 and 1000 and the final estimate is used to propagate the ensemble in time up to day 1200. Fig. 5.110 depicts the results of the matching of the total production rate. The red curves denote the actual observed data and the black curves show the ensemble realizations. The ensemble members converge to the 182 Figure 5.110: Obtained match for total production rates, Large Pattern actual observation data once EnKF is switched on after day 200. After day 1000, once EnKF is switched off, the state vectors still follow the observed data until day 1200. The relative RMSE of total production rate as a function of time is depicted in Fig. 5.111 and provides a numerical measure of the match quality. Clearly, the error increases until day 200 and once we start the EnKF procedure it is reduced until day 1000 at 5%. Subsequently, the forecast error is larger at 25%. We have depicted the geometric distribution of{λ kj } from our final estimate in Fig. 5.112. To have a clear view of the reservoir directionality, we have identified the maximum{λ kj } for each injector in Fig. 5.113. The direction of arrows suggest a trend from the northeast to southwest part in agreement with the field fracture orientation. 183 Figure 5.111: RMSE variation of total production rate, Large Pattern Figure 5.112: {λ kj } distribution; Large Pattern. Arrow lengths represent the mean estimate of{λ kj } at the final time with a direction along (k,j). 184 Figure 5.113: Maximum{λ kj } distribution; Large Pattern. Each arrow represents the length and direction of the maximum{λ kj } at each injector In this section we implemented the EnKF/CRM procedure for an actual reser- voir. Reliable matches and forecasts for the total production rates were obtained despite the fact that the reservoir under consideration is practically unbalanced. Furthermore, the final obtained well connectivities{λ kj } are in accordance with field’s fracture orientation. The results confirm the efficiency of the EnKF/CRM model in real world scenarios involving large number of injection and production wells. 185 Chapter 6: Summary and Conclusion In this chapter, we conclude with a summary of the results of this dissertation. A list of possible future directions and suggestions for research also appears. 6.1 Summary A new method is presented for assisting in characterization and modeling of the immiscible displacement problem in heterogeneous reservoirs that also allows for a simultaneous optimization of net present value of the reserves by controlling the injected fluid rates. This method advances the previous techniques by allowing the capacitance resistance and oil production models to allow for time-varying parame- ters to achieve a higher predictive power. Furthermore, this method requires far less computational resources than statistical methods based on the fluid dynamics of the reservoir although we may still infer some large scale geologic features and the directionality of the reservoir through the correlations between various connectivity parameters. 186 A pressure interference testing problem was the starting point of this disser- tation. We focused on the monitoring of a moving interface in immiscible floods through continuous pressure interference testing. A correspondence between the location of the front and variations in pressure was presented for monitoring satu- ration fronts in an extended case of a classical 1-D Buckley Leverett setting with various lateral permeability profiles. The flooding in our system was assumed to be at a constant injection rate but with the opportunity to continuously monitor in- jection pressure, fractional composition of produced fluid and production pressure. Lateral changes in effective transmissibilities between the injector and the produc- ers caused by movement of immiscible fluid-fluid interface and detected from the solution of interference pressure equation provide an added dimension for tracking the saturation fronts. The effect of the rock heterogeneities is filtered out by an analysis of the incremental changes observed with respect to a base measurements at the early stages of the flood. Our method effectively allows us to connect the flood front speed with pressure derivatives in a 1-D reservoir. To compensate for the typical scarcity of the pressure data we subsequently shifted our focus to using actual injection and production rates to characterize complex reservoir environ- ments. The abundance of production data and the wide range of injection patterns in the history of real reservoirs provides an opportunity to use numerical estimation methods to characterize the reservoir as a dynamical system. we used a conceptual simplification of the reservoir from a complex fluid flow problem into a linear model 187 with few variables, namely the capacitance resistive model (CRM). This approx- imation is not perfect: non-homogeneous saturation distributions may contribute to a non-linear behavior. Nonetheless, CRM provides an efficient and practically objective tool for describing the reservoir input-output relationship, and nonlinear effects can to some extent be absorbed into time-dependent parameters. Using CRM as the driving dynamical system of the reservoir, we implemented the ensem- ble Kalman filter (EnKF) procedure for estimating the CRM parameters (interwell connectivity, time lag, etc.) and also obtained a time-dependent estimate for the parameters to accommodate non-linearities in the input-output relationships. We showed through synthetic and real case studies that the uncertainty in estimation of interwell connectivities and time lags can be systematically reduced with the EnKF procedure. The real time nature of the EnKF procedure in particular allows for such time-dependent effective parameters to be characterized easily. Indeed in some of the synthetic cases that we studied, the CRM parameters turned out to have a time-dependent estimate as expected from the dependence of transmissibility on the time-variations of the saturation in the reservoir. The final stabilized estimates for the CRM parameters in all cases provided a reliable forecast of the reservoir’s future production. We also identify the limitations of our scheme in cases where the number of wells makes the numeric computations intractable or when due to the co-linearity of and therefore the redundancy of the data, the unknown parameters cannot be identified. 188 As a final utility, we employed the forecasting ability of the CRM description in conjunction with a phenomenological model of oil production, to optimize the net present value of the reservoir by controlling injection patterns using the so called EnOpt optimization algorithm which benefits from the compactness of the CRM formulation. The optimized controls yielded a higher cumulative oil production and a lower cumulative water production compared with a neutral base case. 6.2 Future Research Proposals • The pressure interference testing analysis in Chapter 2 is directly applicable to a 1-D reservoir only. Yet one may envisage a scenario in which an ar- ray of probing wells are employed to provide a relationship between pressure variations and location of the front in a more complex reservoir. In such a scenario, a direct analytic expression for the reservoir evolution does not ex- ist and one must rely on grid-based simulations for describing the reservoir dynamics. This can be done for example through the EnKF procedure for estimating the grid parameters. The actual parametrization of the fluid front in 2D is the main challenge in such a description. • It is possible to adjust the CRM linear dynamical systems by non-linear per- turbation defined by introducing new sets of parameters. We may be able to probe such extra parameters provided that we can physically motivate all 189 such non-linear terms. Such a non-linear description, while still approximate, may provide an even better approximation of the reservoir dynamics. • With the exception of certain symmetric cases, it is possible to visually re- late the stabilization of the effective time-dependent CRM parameters to the stabilization of the reservoir saturation or the so called water break-through time. An interesting problem is to numerically establish the connection at a deeper level such that the time-variations of the effective CRM parameters directly probe the saturation movement or stabilization. • An even more ambitious proposal that somewhat encompasses the previous points is to introduce a new set of variables to each injection/production well that characterizes “local saturation” . One may then directly relate such a lo- cal saturation with the mass balance associated with the oil and water phases separately. Such an enlarged state vector may also replace the phenomeno- logical model used oil/water production ratios in conjunction with the CRM model thorough this thesis. 190 Glossary α Phase: oil o, water w, or gas g. α l Tuning parameter. ΔP i Pressure drop in the i-th well. Vector of measurement errors. λ ij Interwell-connectivity between injection well i and production well j. μ α Viscosity of phase α. φ Porosity. ρ α Density of phase α. τ j Time lag of the j-th production well at time t. ξ Dummy integration variable. A Cross-sectional area. a,b Phenomenological parameters for matching the oil production rate c o Oil compressibility. c r Rock compressibility. c t Total compressibility around the control volume. c w Water compressibility. C D Covariance of measurement noise. C M Covariance matrix of the model variables. C x,g(x) Cross covariance between the control variables x and the net present value g(x). C Y Ensemble covariance matrix. 191 D Matrix of measurement. d Observation variables. f w,i Water cut ratio in the i-th well. G Acceleration due to gravity. H Projection matrix. I Identity matrix. i, j, k Well indices. i k (t) Total injection rate in the k-th injection well at time t. J j Productivity index of the j-th production well. K Matrix permeability. k rα Relative permeability of phase α. L One dimensional reservoir length. l Optimization iteration index. m Model variables. N e Total number of realizations. N t Total number of simulation. P α Pressure of phase α. p j (t) Average reservoir pressure around producerj at timet. ¯ p j denotes its spatial average. p wf,j Flowing bottom hole pressure of j-th production well. q j (t) Total rate in the j-th well at time t. q j (t 0 ) Total rate in the j-th well at time t 0 . Q o Cumulative oil production. Q w Cumulative water production. r τ Discount rate. s α Saturation of phase α. 192 T ave Average transmissibility. T eff Effective transmissibility. T phe Phenomenological transmissibility. T kj Effective transmissibility between injection well k and producing well j. v α Velocity of phase α. V p Pore volume around each producer. v o Oil sale price. v w Cost of water disposal. W j (t n ) Cumulative water injected toward j-th production well at time t n . x Vector of control variables. x f Front location. y State vector which consists of static,dynamic and observable variables. Y 0 Perturbation matrix state. y f Forecasted ensemble member. y p Prior estimated vector of variables. y u Updated ensemble member. 193 References [AH04] J. D Annan and J. C Hargreaves. Efficient parameter estimation for a highly chaotic system. ellus. Series A, Dynamic meteorology and oceanography, 56(5):520–526, 2004. [AHEM05] J. D. Annan, J. C. Hargreaves, N. R. Edwards, and R. Marsh. Param- eter estimation in an intermediate complexity earth system model using an ensemble Kalman filter. Ocean Modelling, 8(1-2):135–154, 2005. [AL88] C. Ayan and WJ Lee. The effects of multiphase flow on the interpreta- tion of pressure buildup tests. SPE formation evaluation, 3(2):459–466, 1988. [AL03] A. Albertoni and L. W. Lake. Inferring interwell connectivity only from well rate fluctuations in waterfloods. SPE Reservoir Evaluation and Engineering, 6(1):6–16, 2003. [AS79] K. Aziz and A. Settari. Petroleum reservoir simulation. Chapman & Hall, 1979. [BH92] R. G. Brown and P. Y. C. Hwang. Introduction to Random Signals and Applied Kalman Filtering. 1992, volume 261. Wiley, 1992. [BNJ + 04] D. R. Brouwer, G. Naevdal, J. D. Jansen, E. H. Vefring, and C. van Kruijsdijk. Improved reservoir management through optimal control and continuous model updating. SPE Journal, 90149:27–29, 2004. [CH59] B. C. Craft and M. F. Hawkins. Applied petroleum reservoir engineering. Englewood Cliffs New York, 1959. [COZ08] Y. Chen, D. S. Oliver, and D. Zhang. Efficient ensemble based closed- loop production optimization. In SPE/DOE Improved Oil Recovery Symposieum, 2008, number SPE 112873, 2008. [COZ09] Y. Chen, D. S. Oliver, and D. Zhang. Data assimilation for nonlinear problems by ensemble Kalman filter with reparameterization. Journal of Petroleum Science and Engineering, 66(1-2):1–14, 2009. 194 [CR56] A. T. Corey and C. H. Rathjens. Effect of stratification on relative permeability. Trans., AIME, 207:358–360, 1956. [EA84] I. Ershaghi and D. Abdassah. A prediction technique for immiscible processes using field performance data. Journal of petroleum technology, 36(4):664–670, 1984. [EO78] I. Ershaghi and O. Omoregie. A Method for Extrapolation of Cut vs Recovery Curves. In Journal of Petroleum Technology Forum, February 1978, 1978. [EO84] M. J. Economices and D. O. Ogbe. Single-well and multiwell pressure interference analysis. In SPE California Regional Meeting, 1984, 1984. [ER73] R. C. Earlougher and H. J. Ramey. Interference analysis in bounded systems. Journal of Canadian Petroleum Technology, Oct.–Dec.:33–45, 1973. [Eve03] G. Evensen. The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean Dynamics, 53(4):343–367, 2003. [Gel02] A. Gelb. Applied optimal estimation. MIT press, 2002. [Gen05] P. Gentil. The use of multilinear regression models in patterned water- floods: Physical meaning of the regression coefficients. Master’s thesis, The university of Texas at Austin, 2005. [GO05] Y. Gu and D. Oliver. History matching of the PUNQ-S3 reservoir model using the ensemble Kalman filter. SPE Journal, 10(2):217–224, 2005. [Gra02] W. D. Graham. Estimation and prediction of hydrogeochemical parame- ters using Extended Kalman filtering, pages 327–363. American Society of Civil Engineers Press, 2002. [GZR06] G. Gao, M. Zafari, and A. Reynolds. Quantifying uncertainty for the PUNQ-S3 problem in a Bayesian setting with RML and EnKF. SPE Journal, 11(4):506–515, 2006. [HFMK97] K. J. Heffer, R. J. Fox, C. A. McGill, and N. C. Koutsabeloulis. Novel techniques show links between reservoir flow directionality, Earth stress, fault structure and geomechanical changes in mature waterfloods. SPE Journal, 2(2):91–98, 1997. [Hup70] JD Huppler. Numerical investigation of the effects of core hetero- geneities on waterflood relative permeabilities. SPE Journal, 10(4):381– 392, 1970. 195 [Kiv03] G. A. Kivman. Sequential parameter estimation for stochastic systems. Nonlinear Processes in Geophysics, 10(3):253–260, 2003. [KMJ72] H. Kazemi, L. S Merrill, and J. R. Jargon. Problems in interpretation of pressure fall off tests in reservoirs with and without fluid banks. Journal of Petroleum Technology, Sept.:1147–1156, 1972. [LNL03] R. J. Lorentzen, G. Nævdal, and A. Lage. Tuning of parameters in a two phase flow model using an ensemble Kalman filter. International Journal of Multiphase Flow, 29(8):1283–1309, 2003. [LO05] N. Liu and D. S. Oliver. Ensemble Kalman filter for automatic history matching of geologic facies. Journal of Petroleum Science and Engineer- ing, 47(3-4):147–161, 2005. [LON + 09] K. H. Lee, A. Ortega, A. Nejad, N. Jafroodi, and I. Ershaghi. A Novel Method for Mapping Fractures and High Permeability Channels in Wa- terfloods Using Injection and Production Rates. In SPE Western Re- gional Meeting, 2009. [LWE + 07] X. Liang, D. B. Weber, T. F. Edgar, L. W. Lake, M. Sayarpour, and A. Al-Yousef. Optimization of oil production based on a capacitance model of production and injection rates. In SPE Hydrocarbon Economics and Evaluation Symposium, Dallas, 2007, number SPE 107713, 2007. [Mar59] J. C. Martin. Simplified equations of flow in gas drive reservoirs and the theoretical foundation of multiphase pressure buildup analysis. Trans. AIME, 216:309–311, 1959. [Mar81] C. Marle. Multiphase flow in porous media. Editions Technip, 1981. [MM36] M. Muskat and M.W. Meres. The flow of heterogeneous fluids through porous media. Journal of Applied Physics, 7:346, 1936. [MSGH05] H. Moradkhani, S. Sorooshian, H. V. Gupta, and P. R. Houser. Dual state-parameter estimation of hydrological models using ensem- ble Kalman filter. Advances in Water Resources, 28(2):135–147, 2005. [NJAV05] G. Nævdal, L. M. Johnsen, S. I. Aanonsen, and E. H. Vefring. Reser- voir monitoring and continuous model updating using ensemble Kalman filter. SPE Journal, 10(1):66–74, 2005. [NMV02] G. Nævdal, T. Mannseth, and E. H. Vefring. Near-well reservoir mon- itoring through ensemble Kalman filter. SPE Journal, 75235:13–17, 2002. 196 [OG06] D. S. Oliver and Y. Gu. Iterative Ensemble Kalman Filter for Quan- tifying Uncertainty in Nonlinear Problems. In 2006 Western Pacific Geophysics Meeting, 2006. [PC98] M. N. Panda and A. K. Chopra. An integrated approach to estimate well interactions. In SPE India Oil and Gas Conference and Exhibition, 1998. [Per56] R. L. Perrine. Analysis of pressure buildup curves. Drilling and Pro- duction Practices, n. s.:482, 1956. [Rag89] R. Raghavan. Well test analysis for multiphase flow. SPE formation evaluation, 4(4):585–594, 1989. [Ram75] H. J. Ramey. Interference analysis for anisotropic formations-a case history. Journal of Petroleum Technology, 27(10):1290–1298, 1975. [Ref96] B.T. Refunjol. Reservoir characterization of North Buck Draw field based on tracer response and production/injection analysis. Master’s thesis, University of Texas at Austin, 1996. [SDA05] P. Sarma, L. J. Durlofsky, and K. Aziz. Efficient closed-loop production optimization under uncertainty, SPE paper 94241, presented at the SPE Europec. In EAGE Annual Conference, Madrid, Spain, pages 13–16, 2005. [SK99] T. Soeriawinata and M. Kelkar. Reservoir management using produc- tion data. In SPE Mid-Continent Operations Symposium, 1999. [SP98] J. O. Sant’Anna Pizarro. Estimating injectivity and lateral autocorrela- tion in heterogeneous media. PhD thesis, University of Texas at Austin, 1998. [SZKL07] M. Sayarpour, E. Zuluaga, C. S. Kabir, and L. W. Lake. The use of capacitance resistive models for rapid estimation of waterflood per- formance and optimization. In SPE Annual Technical Conference and Exhibition, 2007, number SPE 110081, 2007. [TB03] M. Thiele and R. Batycky. Water injection optimization using a stream- line based workflow. In SPE Annual Technical Conference and Exhibi- tion, 2003. [Tim71] E. H. Timmerman. Predict performance of water floods graphically. Petroleum Engineering, 43:12, 1971. 197 [VH01] M. Verlaan and A. W. Heemink. Nonlinearity in data assimilation ap- plications: a practical method for analysis. Monthly Weather Review, 129(6):1578–1589, 2001. [WC05] X. H. Wen and W. H. Chen. Real-time reservoir model updating using ensemble Kalman filter. In SPE reservoir simulation symposium, 2005. [WL03] M. P. Walsh and L. W. Lake. A generalized approach to primary hydro- carbon recovery. Elsevier Science, 2003. [WLR07] C. Wang, G. Li, and A. C. Reynolds. Production Optimization in the Closed-loop Reservoir Management. University of Tulsa, 2007. [YGJL06] A. A. Yousef, P. Gentil, J. L. Jensen, and L. W. Lake. A Capacitance Model To Infer Interwell Connectivity From Production and Injection Rate Fluctuations. SPE Reservoir Evaluation & Engineering, 6:630– 646, 2006. [ZR05] M. Zafari and A. Reynolds. Assessing the uncertainty in reservoir de- scription and performance predictions with the ensemble Kalman filter. In SPE Annual Technical Conference and Exhibition, 2005. 198

Abstract (if available)

Abstract Optimizing oil recovery subject to operational constraints requires accurate reservoir models with high predictive capabilities. Specifically, during the immiscible injection processes such as water flooding or gas injection, changes occurring in the flow directions and proportions can affect the long term recovery factors dramatically. Monitoring such changes, as reflected in the rate and pressure measurements, in real time, can provide us with information that can prevent long term losses and help in optimal oil recovery. Such a monitoring platform is tied to a predictive model and should incorporates historic and the incoming data as it becomes available. In this dissertation, we used a continuous pressure interference testing to characterize a 1-D reservoir and subsequently obtained the permeability profile and were able to link the fluid front movement to the pressure difference variations across the reservoir. To compensate for the typical scarcity of pressure data as opposed to the abundance of rate data we introduced a novel methodology for efficient characterization, prediction, and optimization of large water-flooding operations. This is based on ensemble based closed-loop production optimization (EnOpt) and capacitance resistive model (CRM) as its linear underlying dynamicalsystem where the injection rates are the driving force or input signal and the production rates are the dynamical variables or output signals. The production rate data are assimilated in real-time by an ensemble Kalman filter for characterization of the reservoir.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

Flood front tracking and continuous recording of time lag in immiscible displacement

$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

Modeling and simulation of complex recovery processes

PDF

A method for characterizing reservoir heterogeneity using continuous measurement of tracer data

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

Parameter estimation to infer injector-producer relationships in oil fields: from hybrid constrained nonlinear optimization to inference in probabilistic graphical model

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

PDF

Investigating statistical modeling approaches for reservoir characterization in waterfloods from rates fluctuations

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

PDF

Steam drive: Its extension to thin oil sands and reservoirs containing residual saturation of high gravity crude

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

A study of diffusive mass transfer in tight dual-porosity systems (unconventional)

PDF

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

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

Asset Metadata

Creator Jafroodi, Nelia (author)

Core Title Real-time reservoir characterization and optimization during immiscible displacement processes

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering (Smart Oilfield Technologies)

Publication Date 11/11/2009

Defense Date 10/16/2009

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

Tag capacitance resistive modeling,ensemble Kalman filter,history matching,OAI-PMH Harvest,Petroleum Engineering,reservoir characterization

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ershaghi, Iraj (committee chair), Hashemi, Hossein (committee member), Zhang, Dongxiao (committee member)

Creator Email jafroodi@usc.edu,neliajafroodi@gmail.com

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

Unique identifier UC1131686

Identifier etd-jafroodi-3302 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-273804 (legacy record id),usctheses-m2731 (legacy record id)

Legacy Identifier etd-jafroodi-3302.pdf

Dmrecord 273804

Document Type Dissertation

Rights Jafroodi, Nelia

Type texts

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

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu