University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content HISTORY MATCHING FOR FRACTURED RESERVOIRS WITH HOUGH-TRANSFORM-BASED PARAMETERIZATION by Le Lu A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) December 2015 Copyright 2015 Le Lu Acknowledgements I would like to express the deepest appreciation to my advisor, Professor Dongxiao Zhang, for making this dissertation possible. I am thankful for his guidance and continuous support, and for his patience, motivation, and immense knowledge. My sincere thanks also goes to the professors on my dissertation committee: Dr. Iraj Ershaghi and Dr. Roger Ghanem, as well as Dr. Felipe de Barros, Dr. Behnam Jafarpour and Dr. Kristian Jessen, who served on my qualifying exam committee, for their insightful comments and encouragement. Special thanks to the fellow students and post-doctors in our research groups at USC, Heli Bao, Dr. Haibin Chang, Dr. Parham Ghods, Dr. Hamid Reza Jahangiri, Dr. Heng Li, Dr. Weixuan Li, Dr. Qinzhuo Liao, Dr. Liangsheng Shi, Dr. Lingzao Zeng, for the discussions that enlightened me about my research, and for the many unforgettable happy moments together. I would like to acknowledge the nancial supports provided by US National Science Foundation, National Natural Science Foundation of China, the Major Projects of the Ministry of Science and Technology of China, and Mork Family Department of Chemical Engineering and Materials Science. Finally, I would like to thank my wife Linzi Zhang and my parents, for their unwavering love and support. 1 Contents Acknowledgements 1 List of Tables 5 List of Figures 6 Abstract 16 1 Introduction 18 1.1 Numerical Simulation of Fluid Flow in Fractured Reservoirs . . . . . 18 1.2 History Matching and Reservoir Characterization for Fractured Reser- voirs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2 Ensemble Based Methods for History Matching 24 2.1 Ensemble Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2 Ensemble Smoother . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3 Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4 Covariance Localiztion in Ensemble Kalman Filters and Smoothers . 32 2.5 Assimilation of Time-Lapse Seismic Data . . . . . . . . . . . . . . . . 35 3 Hough-Transform-Based Parameterization 39 3.1 Hough Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2 Hough-Transform-Based Parameterization for Fractured Reservoirs . 43 3.3 Implementation of Hough-Transform-Based Parameterization in Ensemble- Based History Matching . . . . . . . . . . . . . . . . . . . . . . . . . 52 2 3.4 Integration of Prior Information . . . . . . . . . . . . . . . . . . . . . 56 3.5 Covariance Localization . . . . . . . . . . . . . . . . . . . . . . . . . 57 4 History Matching for a Synthetic Fractured Reservoir Model 61 4.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.2 History Matching Methods . . . . . . . . . . . . . . . . . . . . . . . . 64 4.3 Case 1: Eectiveness of The Hough-Transform-Based Parameterization in a Simple Reference Field . . . . . . . . . . . . . . . . . . . . . . . 65 4.4 Case 2: Improvements from covariance localization . . . . . . . . . . 72 4.5 Case 3: History Matching for Reference Field with Two Fractures . . 74 4.6 Case 4: Integration of Prior Information . . . . . . . . . . . . . . . . 77 4.7 Case 5: History Matching for a Reference Field with Fractures Away from Wells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5 Assimilation of Time-lapse Seismic Data with Hough-Transform- Based Parameterization 88 5.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.2 History Matching Methods . . . . . . . . . . . . . . . . . . . . . . . . 90 5.3 Case 1: Assimilation of Time-Lapse Seismic Data in a Simple Reference Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.4 Case 2: Assimilation of Time-Lapse Seismic Data for Reference Field with Two Fractures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6 Characterization of Structural Control in Geothermal Reservoirs 107 6.1 Summary of Soultz-sous-For^ ets EGS site . . . . . . . . . . . . . . . . 108 6.2 Case 1: History Matching for a Tracer Test at Soultz-sous-For^ ets . . . 108 6.3 Case 2: Characterization of Large-Scale Fault Network at Soultz-sous- For^ ets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.3.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . 115 3 6.3.2 Scenario 1: History Matching with Prior Information at Pro- ducing Wells . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.3.3 Scenario 2: History Matching with Prior Information for All Faults . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 7 Conclusions and Future Work 132 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 Nomenclature 136 Bibliography 140 4 List of Tables 3.1 Accumulator array for the binary image in Fig. 1. Only part of the array within 2 [0:1; 1:5] and 2 [0:0; 35:0] is provided. The largest entry in the array occurs at (0:9; 20:0), corresponding to the straight line in Fig. 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.1 Parameters used in the rock physics model for seismic impedance . . 91 6.1 Parameters used in the geothermal simulations . . . . . . . . . . . . . 109 5 List of Figures 2.1 Seismic impedance vs. water saturation for various porosity values, assuming gas saturation is zero and the other parameters are constant. 38 3.1 Angle-radius parameterization of a straight line. is the angle of its normal, and is its algebraic distance from the origin o. . . . . . . . 40 3.2 (a) Pixels in a binary image are transformed into (b) sinusoidal curves in the Hough space by Eq. 3.3. The curves belonging to co-linear pixels will cross at a point in the Hough space which represents the straight line through the pixels. . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3 The Hough-transform-based parameterization for a fracture. is the angle of its normal; is its algebraic distance from the origin o; L is its length; and D is the axial displacement of its midpoint from its perpendicular foot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4 The relationship between the approximation of facies map error g FE and Hough space resolution N H , for a 40 40 ow domain. N H is the resolution of a N H N H Hough space. . . . . . . . . . . . . . . . . . 48 3.5 The relationship between the approximation of facies map error g FE and Hough space resolutionN H , for (a)60 60 ow domain and (b)80 80 ow domain. N H is the resolution of a N H N H Hough space. . . . . 48 3.6 The relationship between the approximation of facies map error g FE and Hough space resolutionN H , for (a)20 40 ow domain and (b)40 20 ow domain. N H is the resolution of a N H N H Hough space. . . . . 49 6 3.7 The relationship between the approximation of facies map error g FE and Hough space resolutionN H , for (a)20 40 ow domain and (b)40 20 ow domain. N H is the resolution of a N H N H Hough space. . . . . 50 3.8 The linear relationship between the recommended Hough space reso- lution and ow domain resolution. N H is the resolution of aN H N H Hough space, recommended for a N 1 N 2 ow domain. . . . . . . . . 50 3.9 The Hough space is discretized into a 100 100 Cartesian grid. The black colored grid blocks are active during the history matching. . . . 52 3.10 The work ow of history matching using the iterative smoother LM- EnRML and Hough-transform-based parameterization. . . . . . . . . 55 3.11 Distance between a Hough space grid (;) and a data point (x o ;y o ) in ow domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.12 The ow domaine is discretized into 41 41 regular grid. Nine ob- servation wells marked as black dots on this map provide observation data to be assimilated to update the Hough space parameters. . . . . 59 3.13 The localization function maps for the Hough space parameters to the observations at wells (a)O1, (b)O2, (c)O3, (d)O4, (e)O5, (f)O6, (g)O7, (h)O8, and (i)O9, respectively. . . . . . . . . . . . . . . . . . . . . . . 60 4.1 The ow domain and well locations of the synthetic cases. Oil rate ob- servations are available at all eight producers denoted by lled circles, and bottom hole pressure is available at the injector in the center of the domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.2 The Hough space is discretized into a 100 100 Cartesian grid. The black colored grid blocks are active during the history matching. . . . 63 4.3 Case 1: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix; the producer P3 and the injector I are connected by a fracture. . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.4 Case 1: Reduction in the mean of the data mismatch with iteration. . 66 7 4.5 Probability maps of fractures for Case 1: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 17 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.6 Facies maps of three randomly selected realizations for Case 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 69 4.7 Mean Hough function maps for Case 1: (a) at the initial stage, (b) after the 5 t h iteration step, and (c) after the 17 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.8 The predictions from all the ensemble members before and after the history matching for Case 1: oil rate at (a)P 2 and (c)P 3; water cut at (b) P 2 and (d) P 3; (e) the bottom hole pressure at the injector. The black colored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. The vertical axes in the updated predictions in (a), (c), and (e) are adjusted to show graphs clearly. . . . . . . . . . . . . . . . . . . . . . 71 4.9 Case 2: Reduction in the mean of the data mismatch with iteration. . 73 4.10 Probability maps of fractures for Case 2: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 24 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.11 Mean Hough function maps for Case 2: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 24 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 8 4.12 Case 3: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix; the producer P 2 is connected to P 4 and P 5 by two fractures, respectively. . . . . . . . . . . . . . . . . . . . . 75 4.13 Case 3: Reduction in the mean of the data mismatch with iteration. . 75 4.14 Probability maps of fractures for Case 3: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 20 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.15 Mean Hough function maps for Case 3: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 20 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.16 Facies maps of three randomly selected realizations for Case 3: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 78 4.17 The predictions from all the ensemble members before and after the history matching for Case 3: oil rate at (a) P 2; (c) P 5, and (e) P 7, as well as water cut at (b) P 2, (d) P 5, and (f) P 7. The black col- ored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. The ver- tical axes in the updated predictions in (a), (c), and (e) are adjusted to show graphs clearly. . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.18 Case 4: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix. There is a fracture in NE-SW direction that passes throughP 3, and another fracture in NW-SE direction that does not pass through any wells but connects the former one. . . . . . 80 4.19 Case 4: Reduction in the mean of the data mismatch with iteration. . 81 9 4.20 Probability maps of fractures for Case 4: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 16 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.21 Mean Hough function maps for Case 4: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 16 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.22 Facies maps of three randomly selected realizations for Case 4: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 83 4.23 The predictions from all the ensemble members before and after the history matching for Case 4: oil rate at (a) P 2, (c) P 3, and (e) P 7, as well as water cut at (b) P 2, (d) P 3, and (f) P 7. The black col- ored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. The ver- tical axes in the updated predictions in (a), (c), and (e) are adjusted to show graphs clearly. . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.24 Case 5: The reference facies map. Black grid blocks are fracture facies, and the rest are rock matrix. . . . . . . . . . . . . . . . . . . . . . . . 86 4.25 Probability maps after updating for Case 5 from three separate EnKF assimilation runs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.1 The ow domain and well locations of the synthetic cases. Oil rate ob- servations are available at all eight producers denoted by lled circles, and bottom hole pressure is available at the injector in the center of the domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 10 5.2 Case 1: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix. . . . . . . . . . . . . . . . . . . . . . . . 92 5.3 Case 1: Reduction in the mean of the data mismatch with iteration. . 92 5.4 Probability maps of fractures for Case 1: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 13 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.5 Facies maps of three randomly selected realizations for Case 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 95 5.6 Mean Hough function maps for Case 1: (a) at the initial stage, (b) after the 5 t h iteration step, and (c) after the 13 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.7 The predictions from all the ensemble members before and after the history matching for Case 1: oil rate at (a)P 2 and (c)P 3; water cut at (b) P 2 and (d) P 3; (e) the bottom hole pressure at the injector. The black colored sections of the lines are the period (rst 20 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 30 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. . 97 5.8 Mean impedance dierence maps at the end of the 20 th day for Case 1: (a) the observations, (b) at the initial stage, and (c) after the 13 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . 98 5.9 Mean impedance maps at the end of the 20 th day for Case 1: (a) the observations, (b) at the initial stage, and (c) after the 13 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 11 5.10 Case 2: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix. . . . . . . . . . . . . . . . . . . . . . . . 99 5.11 Case 2: Reduction in the mean of the data mismatch with iteration. . 100 5.12 Probability maps of fractures for Case 2: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 15 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.13 Mean Hough function maps for Case 2: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 15 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.14 Facies maps of three randomly selected realizations for Case 2: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.15 The predictions from all the ensemble members before and after the history matching for Case 2: oil rate at (a) P 2 and (c) P 3, as well as water cut at (b) P 2 and (d) P 3; (e) the bottom hole pressure at the injector. The black colored sections of the lines are the period (rst 20 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 30 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.16 Mean impedance dierence maps at the end of the 20 th day for Case 2: (a) the observations, (b) at the initial stage, and (c) after the 15 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. . . . . . . . . . . . . . . . . . . . . . . . . 105 6.1 Deep borehole conguration at Soultz EGS site, based on the sketch from Gessner et al. (2009). The fault is shown as the blue ellipse, and the trajectories of the wells are shown as black curves. . . . . . . . . 109 12 6.2 The ow domain and well locations of Case 1. Color of the grid blocks indicates the facies distribution in the reservoir. . . . . . . . . . . . . 110 6.3 Case 1: Reduction in the mean of the data mismatch with iteration. . 111 6.4 Probability maps of fractures for Case 1: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 9 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.5 Facies maps of three randomly selected realizations for Case 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 113 6.6 The predictions of tracer concentration from all the ensemble members before and after the history matching for Case 1: tracer concentration at (a)(b) GPK2 and at (c)(d) GPK4. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.7 The predictions from all the ensemble members before and after the history matching for Case 1: water rate at (a) GPK2, (c) GPK4, and (e) GPK3, as well as bottomhole pressure at (b) GPK2, (d) GPK4, and (f) GPK3. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . 116 6.8 The ow domain and well locations of Case 2. Color of the grid blocks indicates the facies distribution in the reservoir. Names fo the wells are marked as red, and names of the faults are marked as blue. . . . . 117 6.9 Scenario 1: Reduction in the mean of the data mismatch with iteration.119 6.10 Probability maps of fractures for Scenario 1: (a) at the initial stage; (b) after the 1 st iteration step; and (c) after the 5 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . . . . . 119 13 6.11 Facies maps of three randomly selected realizations for Scenario 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 1 st iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 121 6.12 The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 1: tracer concen- tration at (a)(b) P 1 and at (c)(d) P 2. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.13 The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 1: tracer concen- tration at (a)(b) P 3 and at (c)(d) P 4. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.14 The predictions from all the ensemble members before and after the history matching for Scenario 1: water rate at (a) P 1, (c) P 4, and (e) I1, as well as bottomhole pressure at (b) P 1, (d) P 4, and (f) I1. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . 124 6.15 Scenario 2: Reduction in the mean of the data mismatch with iteration.125 6.16 Probability maps of fractures for Scenario 2: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 12 th iteration step when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 126 6.17 Facies maps of three randomly selected realizations for Scenario 2: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. . . . . . . . . . . . . . . . . . . . . . . 127 14 6.18 The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 2: tracer concen- tration at (a)(b) P 1 and at (c)(d) P 2. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.19 The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 2: tracer concen- tration at (a)(b) P 3 and at (c)(d) P 4. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.20 The predictions from all the ensemble members before and after the history matching for Scenario 2: water rate at (a) P 1, (c) P 4, and (e) I1, as well as bottomhole pressure at (b) P 1, (d) P 4, and (f) I1. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. . . . . . . . . . . . . . 131 15 Abstract Natural or hydraulically induced fractures provide the primary ow capacity in many reservoir, such as tight sand and shale gas reservoirs, where the permeability of the rock matrix is very low compared with the conventional reservoirs. Hence, the success of production in fractured reservoirs is signicantly dependent on the knowledge of the spatial distribution of fractures. However, it is a challenge to estimate fracture distri- bution by conventional history matching methods because of the complex geology of such reservoirs. Although there have been great progress in assisted history match- ing techniques during the last couple of decades, estimating fracture distributions in fractured reservoirs is still inecient due to the strong heterogeneity and spatial discontinuity of model parameters. It is well known that the ensemble-based history matching methods, such as the Ensemble Kalman Filter (EnKF), performs best for models with Gaussian distributed parameters, since two-point statistics is used to represent the relationship between the model parameters and the observations. The performance of the EnKF can be far from optimum for the non-Gaussian distributed parameters, such as facies indicator, eective permeability, and porosity. Although the geometric shapes of fractures may be properly generated at the initial step, the geometrical shape and spatial continuity of the features are dicult to preserve after updating, resulting in geologically unrealistic fracture distribution maps. In this research, we develop an assisted history matching method for fractured reservoirs with a Hough-transform-based parameterization. The facies maps of frac- tured reservoirs are parameterized into Hough function elds in a discrete Hough space, while each grid block in the Hough space represents a fracture dened by its two Cartesian coordinates: angle of its normal and its algebraic distance from 16 the origin in the ow domain. The length and axial position of the fractures are dened by two additional parameters nested on the same grid. The Hough function value of each grid block in the Hough space is used as the indicator of the existence of the fracture in the facies map. Facies maps can be uniquely generated from the Hough space parameters by an inverse transform, which guarantees that no matter how the Hough space parameters are modied, the fractures in facies maps are always geologically realistic. In assisted history matching with Hough-transform-based parameterization, the parameter elds in the Hough space, instead of the facies maps, are updated with production history and seismic data. Before each uid ow simulation during the history matching, an inverse transform is performed to generate facies maps for the reservoir simulator. Point-wise prior information such as known fractures discovered from well log data, as well as the statistics of regional fracture orientations, can be honored by the inverse transform. The eectiveness of this method is investigated with two dimensional synthetic water ooding examples. Seismic data is also assimilated in some of the examples to improve the characterization of fractures away from wells. The fracture distributions in the reference maps are identied by this method, and updated models are capable to provide improved predictions for prolonged periods of production. A few realistic examples based on the tracer tests at Soultz-sous-For^ ets geothermal reservoir are used to demonstrate the performance of proposed method in characterizing complex geological structure with fault network. 17 Chapter 1 Introduction Natural or hydraulically induced fractures and faults provide the primary ow ca- pacity in many reservoirs, especially in tight sand and shale gas reservoirs, where the permeability of the rock matrix is very low compared with the conventional reser- voirs. In such strongly heterogeneous and anisotropic reservoirs, the ow pattern and production behavior are signicantly in uenced by the distribution of fractures. Thus a good estimation of the fracture distribution is essential to achieve commercial production rate and reduce the risk of early water breakthrough. 1.1 Numerical Simulation of Fluid Flow in Frac- tured Reservoirs Numerical simulation of uid ow in fractured reservoirs is challenging, because of the high contrast between the ow behavior in the highly conductive fractures and that in the tight rock matrix. The strong heterogeneity of rock properties and geometrical complexity of fracture network make it dicult for the numerical models to accurately characterize the fractured reservoir. Several numerical simulation techniques have been developed for fractured reservoirs: discrete fracture network (DFN) models, dual continuum models, and facies models. Discrete fracture network models explicitly characterize the individual fractures 18 with their orientation, length, aperture, location, etc.. Data from dierent sources, such as well logs, seismic surveys, outcrops, production tests, are integrated to gener- ate the unstructured fracture networks, consisting of discrete polygonal or oval-shaped planes embedded in a rock matrix. Multi-scale fracture networks are usually gener- ated to map the fractures at a eld scale. Several techniques are used to generate the fracture networks (Kasiri and Bashiri, 2011): stochastic forward modeling, simulated annealing, fractal fracture networks, and fractal-wavelet-neural network model. The uid ow problem is solved at both fracture domain and matrix domain, while the ow between the fractures and matrix is controlled by transmissivity parameters. The eld-scale simulation with DFN is computationally dicult because of the complex- ity of the fracture network, and the heavy demand of detailed eld data for model calibration limited DFN's practical applications (Berkowitz, 2002). Dual continuum techniques are developed to implicitly simulate the fracture net- work in fractured reservoirs. The fractures are modeled with a separate continuum along with the rock matrix, and the ow between fracture and matrix is simulated with a transmissivity function. As one of the most famous dual continuum approaches, the dual porosity model was introduced in early 1960s by Warren and Root (1963). The fractured reservoir is characterized by two completely overlapping continua for the fracture and matrix blocks. The two continua possess separate state variables (e.g., pressure, saturations, etc.). Fluids can only move in the fractures, and matrix blocks act as source terms feeding the fractures. The mass transfer between the two media is modeled by a source-sink transfer function. Dual continuum methods have been widely used in petroleum reservoirs as because of good computational eciency. Facies models are widely used in geostatistics to construct reservoir models with complex geology. The properties (e.g., permeability, porosity, etc.) of the facies are populated into the facies maps, in order to generate parameter elds to be simulated by pixel-based reservoir simulators. Ping and Zhang (2013) used vector-based level set method to update the facies maps in fractured reservoirs, where the fracture facies and matrix facies are assigned with distinctive permeabilities and porosities. The multiple point statistics modeling (Strebelle, 2002; Caers, 2003) can also be used to generate facies maps for fractured reservoirs by training images with fractures. 19 1.2 History Matching and Reservoir Characteriza- tion for Fractured Reservoirs In addition to the prior knowledge of the geological characteristics of the reservoir and the near wellbore estimations from well logs, the relatively easily available production history of the wells is an invaluable source of knowledge to be used in characterizing the parameters in reservoir models. With the history matching techniques, the histor- ical production and pressures are compared with the simulated values from reservoir simulators to estimate petrophysical parameters and dynamic states. The increase of computational power and the adoption of geostatistics and Monte Carlo methods stimulated remarkable progress in the development of history matching algorithms in the past couple decades (Oliver and Chen, 2011). As one of the most popular advanced history matching methods applied in the petroleum industry, the Ensemble Kalman Filter (EnKF) is a robust and ecient technique to estimate model parameters and quantify the uncertainty of model predictions, while continuously assimilating time- dependent observations. Proposed by Evensen (2004), the EnKF is a Monte Carlo approach based on the Kalman lter, using an ensemble of model realizations to eval- uate the statistics of the states. It is capable to handle non-linear and large-scale problems. However, the analysis step in the EnKF is based on the rst and second order moments of the distributions, limiting the performance if the distributions are too far from Gaussian (Aanonsen et al., 2009). Especially in highly heterogeneous reservoirs, such as fractured reservoirs and channelized reservoirs, the petrophysical parameters and model predictions are usually of multimodal distribution. The rst and second order moments are inadequate to characterize such distributions, thus the EnKF can be incapable or inecient to provide correct parameter estimations and model predictions. In order to reduce the number of model parameters and to characterize complex geological features, where the reservoir models are usually represented by a distribu- tion of geological facies, various parameterization techniques have been developed in the eld of petroleum engineering (Aanonsen et al., 2009; Oliver and Chen, 2011). Liu and Oliver (2005a,b) used truncated pluri-Gaussian method for estimating reservoir 20 facies distribution with EnKF. Two Gaussian random elds Z1 and Z2 are dened on the reservoir domain, and the facies at any location is determined by a truncation map that assigns regions of the Z1-Z2 domain to a particular facies. Dovera and Della Rossa (2007) adopted the Gaussian mixture model in the EnKF. The pdfs of the parameters are described as weighted sums of Gaussian pdfs to be able to sample bimodal parameters and posteriors. Iglesias et al. (2014) used a small number of geometric parameters to characterize the simplied geometry of geological facies in layered and channelized reservoirs, and sampled the posterior distributions of both the geometric parameters and permeabilities within each facies by Markov Chain Monte Carlo methods. Nejadi and Leung (2012) constructed initial realizations of a dual-porosity model with discrete fracture network (DFN) models and upscaling techniques. The ensemble is updated by the EnKF with normal score transform. Jafarpour and McLaughlin (2008) proposed the parameterization with Discrete Cosine Transform (DCT) as a robust and computationally ecient approach to char- acterize facies structures in history matching. The DCT is a Fourier-based image compression technique which uses a set of predened basis functions that do not de- pend on the permeability covariance. Both model parameters and dynamic variables are parameterized by DCT. The coecients of the cosine basis functions are updated by the EnKF, and the facies map and permeability eld can be reconstructed via inverse DCT. Lorentzen et al. (2012) proposed a method in which distance functions are used to parameterize facies distributions. A distance function is dened as the shortest distance between a given position in the eld and the facies boundary. The distances are spatially smoother than the facies map. The distances are updated by the EnKF, and this method is able to estimate varying parameters within an arbitrary number of facies types. multiple point statistics (MPS) and facies modeling have also been used to char- acterize complex geology in history matching. The MPS can infer spatial patterns using many spatial locations of a given geometric template scanning a training im- age (Strebelle, 2002; Caers, 2003). Jafarpour and Khodabakhshi (2011) developed a probability conditioning method to condition the MPS simulation of facies maps with production data. Various conditional geological structures, such as uvial channels 21 and fractures, can be generated by MPS simulation constrained by the probability map, which is a probabilistic spatial description of facies distribution. The proba- bility map is statistically honored by the MPS simulation. The probability map is continuously updated by EnKF using ow measurements, and an ensemble of new facies maps is generated by MPS simulator at each step. Khodabakhshi and Ja- farpour (2011) quantied the uncertainty associated with the training image in the probability conditioning method. Several training images with dierent statistics are utilized to simulate facies maps, and the importance of dierent training images are quantied by history matching with production data. Level set method is a powerful tool to characterize the boundaries of the geologic facies (Iglesias and McLaughlin, 2011). Moreno and Aanonsen (2007) used level set formulation to model facies distribution. Facies boundaries are implicitly modeled through a level set function, and the boundary velocity represented by a Gaussian random eld is updated by the EnKF. Chang et al. (2010) assigned the values of level set function to a representing node system as Gaussian random numbers. The level set function values are linearly interpolated to the rest of the eld, and the facies type is determined by the sign of the level set function. Ping and Zhang (2013) developed a vector based level set method, in which fractures are represented by vectors dened by three parameters at each representing node: a level set function indicating the existence of a fracture starting from the node by its sign, the angle, and the length of the fracture. The representing nodes are deterministically positioned in the ow domain so that the most fractures can be represented. The parameter sets at the representing nodes are updated by the EnKF. Ping and Zhang (2014) also developed a vector-based level set parameterization for channelized reservoirs, where the channels are characterized by the paths connecting the representing nodes. The facies is determined by vectors consisting of level set function, real radius, and virtual radius at the representing nodes. The level set function indicates the facies type at a representing node, the real radius determines the size of the facies unit centered at the node, and the virtual radius controls the connectivity between the nodes. The vectors are updated by the EnKF in two-dimensional and three-dimensional synthetic examples. 22 The parameterization methods mentioned above improved the history matching for reservoirs with non-Gaussian parameters or complex facies distributions. Most of them are preferable for channelized reservoirs, and lack the constraints for the geometrical shapes and internal continuity of the geological features. In order to ex- plicitly represent the orientation and position of any arbitrary fracture in a reservoir, we develop a parameterization algorithm (Lu and Zhang, 2015) to enable the trans- form between the discontinuous facies map and parameter vectors on a continuous domain. The angle-radius formulation of straight lines in Hough transform is adopted to formulate the transform and inverse transform. By parameterizing the fractures into angle-radius parameters, with the help of additional parameters such as length and axial displacement, any arbitrary fracture in the ow domain can be represented by a set of parameters in the angle-radius domain. The proposed methodology will be described in Chapter 3, and some synthetic examples in Chapter 4 are provided to demonstrate the eectiveness of the method in estimating fracture distributions with production data. The assimilation of seismic data with the proposed param- eterization is discussed in Chapter 5, with synthetic examples to demonstrate the improved fracture characterization and model predictions. In Chapter 6, a few realis- tic examples based on the tracer tests at Soultz-sous-For^ ets geothermal reservoir are used to demonstrate the performance of proposed method in characterizing complex geological structure with fault network. 23 Chapter 2 Ensemble Based Methods for History Matching In this chapter, we review the ensemble Kalman lter and smoother, as well as their iterative formulations. 2.1 Ensemble Kalman Filter The Ensemble Kalman Filter (EnKF) (Evensen, 1994) is a Monte Carlo based sequen- tial data assimilation method that can continuously update the estimate of model parameters and the dynamic states with time-dependent observations. Based on the Kalman lter, the EnKF uses an ensemble of model realizations to evaluate the statis- tics of the states. It is capable to handle non-linear and large-scale problems. Since it avoids the explicit computation of covariance or the Jacobian of the model states, it is computationally ecient and can be non-intrusively implemented in commercial reservoir simulators. It is now one of the most widely used history matching methods for petroleum reservoirs. The EnKF simultaneously updates an ensemble of equally probable model state realizations, while the model state is dened as a vector containing parameters, state 24 variables, and observations, i.e.: y n = 2 6 6 6 4 m u n g n 3 7 7 7 5 2R Nm+Nu+Ng (2.1) where y n is the state vector at t = t n ;n = 1; 2;:::;N t ; m2 R Nm denotes the pa- rameters of the model, such as permeability and porosity; u2 R Nu represents the dynamic state variables simulated by the model, such as pressure and saturation; and g n 2R Ng denotes additional dynamic measurements, such as oil rate, water cut, etc. Let N e be the ensemble size, it is convenient to express the state vectors at time t n in matrix notation: Y n =fy 1 n ; y 2 n ;:::; y Ne n g2R NyNe (2.2) The EnKF uses a small ensemble to approximate the covariance between model pa- rameters and the dynamic observations. Thus the ensemble size N e should be small for eciency, but adequate for the complexity of the reservoir model. N e is usu- ally dependent on the number of model parameters, and the degree of uncertainty associated with the parameters and the observations. The initial ensemble Y n can be generated from prior information about the pa- rameters and model's initial conditions. The procedure of EnKF consists of two alternating steps: the forecast step and analysis step. In the forecast step, the reser- voir simulator is applied on each ensemble member independently to obtain the state variables and dynamic measurements at time t n with the model parameters at t n1 : y f;i n =F (y a;i n1 ); i = 1; 2;:::;N e (2.3) whereF () indicates the forward modeling of the reservoir simulator; and superscripts f and a indicate the forecast and assimilation procedure, respectively. 25 Then, in the analysis step, an updated ensemble of state vectors is computed as Y a n = Y f n + K n (d n HY f n ) (2.4) where K2R NyN d is the Kalman gain matrix, common for all the ensemble members; d n 2R N d Ne is the matrix containing observations perturbed by a white noise " n at time t n ; and H is a matrix used to transform the state matrix into the observation matrix, containing only 0s and 1s. The Kalman gain matrix is calculated using the following equations: K n = P f n H T (HP f n H T + R) 1 (2.5) hy f n i 1 N e Ne X i=1 y f;i n (2.6) P f n 1 N e 1 Ne X i=1 y f;i n hy f n i y f;i n hy f n i T (2.7) where the covariance matrix P f n of the state vectors is explicitly estimated from the ensemble; and R is the covariance matrix of the observation errors. The EnKF algorithm can be summarized as the following steps: 1. Generate N e independent realizations of state vector y a;i 0 ; i = 1; 2;:::;N e based on prior knowledge of model parameters and initial conditions. 2. For n = 1; 2;:::;N t , assimilate data by two steps: (a) Forecast step: compute the forecast state Y f n by running the reservoir simulation to obtain the most updated values for dynamic states and mea- surements. (b) Analysis step: update the states to Y a n using Eqs. 2.5-2.7. 2.2 Ensemble Smoother Proposed by van Leeuwen and Evensen (Van Leeuwen and Evensen, 1996), the En- semble Smoother (ES) is an alternative history matching method that computes a 26 single global update, rather than using recursive updates in time as the EnKF. EnKF and ES provide identical solutions for linear models and measurements (Evensen, 2004). Although then EnKF is superior to the ES in nonlinear models with chaotic dynamics (Van Leeuwen and Evensen, 1996), the ES is more ecient than the EnKF in complex reservoir applications because the dynamic variables (pressure, saturation, etc.) are not updated in the ES. Since the history-matched model output is obtained by a new reservoir simulation with the updated parameters, the inconsistent and non- physical updates of the dynamic variables and model restarts typical in the EnKF are avoided in the ES (Skjervheim et al., 2007). If the ensemble of the model state at time t n is dened as Eq. 2.2, the ensemble of the joint state from t 0 to t k in the ES is: f Y = 8 > > > < > > > : Y 0 . . . Y k 9 > > > = > > > ; (2.8) Then the mean removed ensemble perturbation matrix is: f Y 0 = f Y Y = f Y(I 1 Ny ) (2.9) The ensemble covariance e P2R NyNy is dened as: e P = f Y 0 ( f Y 0 ) T N e 1 (2.10) The linear ES analysis equation for f Y a is: f Y a = f Y + f H e P( f H e P f H T + f R) 1 ( e d f H f Y) (2.11) where f H is a matrix used to transform the state matrix into the observation matrix, 27 containing only 0s and 1s, and: e d = 8 > > > < > > > : d 1 . . . d m 9 > > > = > > > ; ; f R = 8 > > > < > > > : R 1 . . . R m 9 > > > = > > > ; (2.12) for measurements available at m measurement times. In the forecast step, the reservoir simulation is run from the initial condition to the last timestep with observation data, and the outputs from the simulation are used to construct the initial ensemble. All the data are assimilated in a single analysis step, when only the model parameters are updated. The reservoir simulation is often run again with the updated parameters to obtain the updated dynamic variables and measurements. 2.3 Iterative Methods The EnKF is ecient in nonlinear and large scale problems. However, for nonlinear problems, it may be impossible to update the dynamic variables, especially the uid saturations, without re-running the reservoir simulation. Usually the updated uid saturations can exceed the bounds and become unphysical. Moreover, in many prob- lems where the relationship between model parameters and observations are strongly nonlinear, especially in reservoirs with complex geology, the posterior pdf of the model parameters can be too far from Gaussian to be approximated with a Gaussian pdf, and the updates to model parameters are not as accurate as in linear problems. Iterative forms of the EnKF have been developed to overcome these diculties (Aanonsen et al., 2009). Nonlinear dynamic systems can be linearized by using a series of linear approximations to the nonlinear functions. Gu and Oliver (2007) pro- posed an ensemble approximation to the Randomized Maximum Likelihood (RML) Filter based on the Gauss-Newton formula to iteratively minimize the objective func- tion. The ensemble sensitivity matrix, which is the coecient matrix relating the changes in model parameters to the changes in measurements, is solved explicitly 28 by singular value decomposition. The sensitivity matrix is numerically unstable and causes low convergence rates for large-scale problems (Chen and Oliver, 2013). Sakov et al. (2012) proposed two formulations of iterative EnKF: The rst formulation, the iterative EnKF (IEnKF), is similar to the EnRML of Gu and Oliver (2007), but uses the ensemble square root lter instead of the traditional EnKF and rescaling of the ensemble anomalies to improve the stability of the sensitive matrix. The second for- mulation, the iterative extended Kalman Filter (IEKF), is an ensemble formulation of the extended Kalman Filter. Both the IEnKF and IEKF oer a considerable ad- vantage over the EnKF in strongly nonlinear problems. Bocquet and Sakov (2012) introduced a nite-size iterative EnKF (IEnKF-N) which is a new implementation of the IEnKF based on the Levenberg-Marquardt optimization algorithm. Iterative forms of the ES have also been developed to assist the history matching in nonlinear problems. Chen and Oliver (2012) used the ensemble randomized max- imum likelihood (EnRML) as an iterative ensemble smoother. The convergence was quite slow due to the poor approximation of the sensitivity matrix from a small ensem- ble. Chen and Oliver (2013) proposed an iterative ensemble smoother LM-EnRML based on a modied Levenberg-Marquardt algorithm. The explicit computation of the sensitivity matrix is avoided by singular value decomposition of the ensemble perturbations. The algorithm converged relatively quickly on multiple test cases, including the Brugge benchmark case. The main purpose of the LM-EnRML algorithm is to minimize the objective function: S j (m) = 1 2 (m m pr ) T C 1 M (m m pr ) + 1 2 (F (m) d) T C 1 D (g(m) d) (2.13) where m pr is the model parameters from the prior distribution with covariance C M , the perturbed observation d is sampled from a multivariant Gaussian pdf with mean d o bs and covariance C D , and the function F () denotes the model predictions from model parameters to data. The l th iteration of a Levenberg-Marquardt update for 29 the model parameter m is: m l = [(1 +! l )C 1 M + G T l C 1 D G l ] 1 [C 1 M (m l m pr ) + G T l C 1 D (F (m l ) d)] (2.14) where ! is the Levenberg-Marquardt regularization parameter, G is the sensitivity matrix of data to model parameters. The sensitivity matrix can be estimated by the ensemble: G l = C 1=2 D d l (m l ) 1 C 1=2 SC (2.15) where m l = C 1=2 SC (m l m l )= q N e 1 (2.16) d l = C 1=2 D (d l d l )= q N e 1 (2.17) Here C SC is a scaling matrix for model parameters, usually selected to be a diagonal matrix with diagonal elements equal to the variance of the prior pdf for each model parameter. Chen and Oliver (2013) replaced the C M in the Hessian term in Eq. 2.14 with P l , the ensemble estimation of model covariance from the updated ensemble to avoid the explicit computation of the sensitivity matrix G l . Eq. 2.14 becomes: m l = [(1 +! l )P 1 l + G T l C 1 D G l ] 1 [C 1 M (m l m pr ) + G T l C 1 D (F (m l ) d l )] (2.18) and P l is computed by P l = C 1=2 SC m l m T l C 1=2 SC (2.19) The Sherman-Woodbury-Morrison matrix inversion (Press et al., 1986) is applied to Eq. 2.18 to obtain: m l = [(1 +! l )P 1 l + G T l C 1 D G l ] 1 C 1 M (m l m pr ) P l G T l [(1 + l )C D + G l P l G T l ] 1 (F (m l ) d l ) (2.20) 30 Truncated singular value decomposition (Press et al., 1986) of d l and m pr gives d l = U d W d V T d (2.21) and m pr = U m W m V T m (2.22) Substituting Eqs. 2.15,2.16,2.19,2.21,and 2.22 to Eq. 2.20, the correction of the l th iterative update becomes m l C 1=2 SC m l V d [(1 + l )I d + W 2 d ] 1 V T d m T l U m W 2 m U T m C 1=2 SC (m l m pr ) C 1=2 SC m l V d W d [(1 +! l )I d + W 2 d ] 1 U T d C 1=2 D (F (m) d l ) (2.23) where I d is an identity matrix of dimension N d . Eq. 2.23 is referred as the LM- EnRML method. An approximation of the LM-EnRML algorithm was derived by neglecting the updates from the model mismatch term: m l = C 1=2 SC m l V d W d [(1 +! l )I d + W 2 d ] 1 U T d C 1=2 D (F (m) d l ) (2.24) The LM-EnRML methods showed good performance in tests where nonlinearity is large (Chen and Oliver, 2013). The incorporation of the Levenberg-Marquardt damp- ing parameter reduces the tendency to add model roughness at early iterations for nonlinear models, and the truncated singular value decomposition of the ensemble avoids the explicit computation of large matrices if the amount of data is large. 31 2.4 Covariance Localiztion in Ensemble Kalman Filters and Smoothers As an ensemble approximation to the Kalman lter, the standard EnKF estimates the covariance matrix with an ensemble of nite number of model realizations. The application of EnKF with small ensemble size suers from two critical problems: 1. Spurious correlation. Sampling error from the small size of ensemble can result in spurious correlations, that induces unnecessary updates to the state variables. Aanonsen et al. (2009) compared the covariance estimate from an ensemble size of 25 with that of 900 in a test example, showing that signicantly spurious correlations are observed in the covariance estimate from a small ensemble, especially for the distant observations to the state variable. 2. Rank-deciency. With the standard EnKF update, each updated ensemble member must be a linear combination of the corresponding initial ensemble, and the updates are restricted to the subspace spanned by the members of the forecast ensemble. If the size of ensemble is smaller than the amount of data to be assimilated, the rank of the approximation of covariance is less than or equal to the number of ensemble members. The limited degree of freedom prevents EnKF to assimilate a large number of independent data. Although the problems can be alleviated by using a large ensemble, the power of the EnKF is the ability to approximate the covariance with relatively small en- sembles. Consequently, Covariance localization is used to reduce the in uence of spurious correlations and allow larger amount of data to be assimilated. Covariance localization helps reduce the spurious correlations to the model parameters far from the observations. On the other hand, each parameter may be updated with a dier- ent linear combination of realizations, and the global update can be obtained from a much larger space than the subspace spanned by the ensemble predictions. In practice, covariance localization is usually applied by either a distance cut- o, or with a distance-dependent function. Houtekamer and Mitchell (1998) used a distance cuto to the Kalman gain so that only parameters within optimal distance of the observation were updated. Dong et al. (2006) updated the parameters at a grid 32 block by only the seismic data available at the same grid. However, discontinuities are induced to the analysis by the covariance cut-o, since much of the spatial continuity constraints in the prior ensemble is ignored. The distance-dependent localization was developed to provide a smoother updated state. The distance-dependent localization function is applied element-wise to the co- variance matrix using the Schur product, which is a element-wise product of matrices: Z P 2 6 6 6 6 6 6 6 4 1;1 1;2 ::: 1;n 2;1 2;2 ::: 2;n . . . . . . . . . . . . n;1 n;2 ::: n;n 3 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 4 p 1;1 p 1;2 ::: p 1;n p 2;1 p 2;2 ::: p 2;n . . . . . . . . . . . . p n;1 p n;2 ::: p n;n 3 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 4 1;1 p 1;1 1;2 p 1;2 ::: 1;n p 1;n 2;1 p 2;1 2;2 p 2;2 ::: 2;n p 2;n . . . . . . . . . . . . n;1 p n;1 n;2 p n;2 ::: n;npn;n 3 7 7 7 7 7 7 7 5 (2.25) where P is the covariance matrix of the state vectors, and Z is the matrix of distance- dependent localization functions in the form of P. The Kalman gain equation for the standard EnKF becomes: K n; = (Z P f n )H T [H(Z P f n )H T + R] 1 (2.26) In practical applications when the amount of data is large, the covariance matrix P is rarely explicitly computed and stored. Hence, the localization functions can be applied element-wise to the Kalman gain matrix in Eq. 2.5, instead of the covariance matrices (Hamill et al., 2001; Agbalaka and Oliver, 2008; Chen and Oliver, 2010): Y a n = Y f n + Z K n (d n HY f n ) (2.27) where the Kalman gain matrix can be solved by singular value decomposition to avoid the explicit inversion and storage of large matrices (Evensen, 2003). Similarly, Chen 33 and Oliver (2014) modied the LM-EnRML method in Eq. 2.24 to a formula with localization m l = ZfC 1=2 SC m l V d W d [(1 + l )I d + W 2 d ] 1 U T d g C 1=2 D (g(m) d l ) (2.28) by multiplying the coecient matrix element-wise by a localization function matrix Z. The algorithm was tested in a history matching for the Norne eld, oshore Norway (Chen and Oliver, 2014). A common choice of the localization function is the compactly supported fth- order piecewise rational function by Gaspari and Cohn (1999): (d) = 8 > > > < > > > : 1 4 d c 5 + 1 2 d c 4 + 5 8 d c 3 5 3 d c 2 + 1; 0dc 1 12 d c 5 1 2 d c 4 + 5 8 d c 3 + 5 3 d c 2 5 d c + 4 2 3 d c 1 ; c<d 2c 0; d> 2c (2.29) whered is the Euclidean distance between a grid point and the observation location, and c is the length scale of the function. The Gaspari-Cohn function has been used in petroleum engineering applications (Agbalaka and Oliver, 2008; Zhang and Oliver, 2009). However, the length scale c should be chosen wisely since it increases with the increase of ensemble size, as shown in the sensitivity studies (Hamill et al., 2001; Houtekamer and Mitchell, 2001). Chen and Oliver (2010) suggested the use of a localization function proposed by Furrer and Bengtsson (2007): (d) = 1 1 + [1 +c(0) 2 =c(d) 2 ]=N e (2.30) wherec() is the target covariance function, andN e is the ensemble size. The function tapers covariance functions with respect to distance and ensemble size, and its value approaches one when the ensemble size is very large. A spherical covariance function is selected to represent the correlation between model parameters and data because the localization function equals zero as the distance goes beyond the range (Chen and 34 Oliver, 2010). Nan and Wu (2011) used another localization function proposed by Furrer and Bengtsson (2007) for exponential covariance functions: c(d) = exp 3d ! (2.31) where is the practical range such that c() 0:05. The localization function is given as: (d) = exp " 3d 4 ( p 9 + 8N e 5) # (2.32) The function increases with ensemble size, and ensures correlation between model parameter and distant observation data, which is especially important in history matching of facies distributions to avoid the spatial discontinuity in the updated states (Aanonsen et al., 2009). 2.5 Assimilation of Time-Lapse Seismic Data The highly available production history provides temporal continuous observations, which is very valuable for the sequential data assimilation techniques, such as the EnKF. However, these observations are spatially sparse, that is, only available at lim- ited number of wells. The estimated reservoir parameters between wells are generally much more uncertain than those near wells. While in fractured reservoirs, the distri- bution of fractures between wells is highly uncertain, especially the fractures that do not connect to wells. Hence, observations with more areal coverage are necessary to estimate the fracture distributions away from wells. Seismic technologies have been used by petroleum geologists to map and interpret potential petroleum reservoirs. Seismic images are produced by generating, recording, and analyzing sound waves that travel through the Earth. Time-lapse seismic data is acquired by two or more seismic surveys shot at the same site over time, in order to image uid migration in a producing reservoir. The rst seismic survey, which is usually taken before the production, is referred as the reference survey, and the 35 others are called monitor surveys. If each of them are 3D seismic surveys, the process is also called 4D seismic. The time-lapse seismic data is of importance in reservoir monitoring, since the dierence between a monitor survey and the reference survey contains information about uid movement in the reservoir due to depletion and water injection. Increasing applications of assimilating seismic data in history matching for petroleum reservoirs had been observed in the past decade (Aanonsen et al., 2009; Oliver and Chen, 2011). Dong and Oliver (2005) used a quasi-Newton minimization method to incorporate time-lapse seismic data in history matching. Skjervheim et al. (2007) used subspace inversion to assimilate interpreted time-lapse seismic data with en- semble Kalman lter. The solution for the timestep with reference seismic survey is continuously updated in parallel with the current solution by ensemble Kalman smoother. The model parameters at a grid block are updated using only data at nearby grid blocks. Zhao et al. (2008) assimilated seismic data to assist in generating the facies maps with truncated pluri-Gaussian models. Chen and Oliver (2010) used seismic data in the ensemble-based closed-loop optimization for the Brugge eld, with distance-dependent covariance localization. In most of the applications, the seismic data is assimilated with the help of a rock physics model, combined with the reservoir simulator. The seismic impedance, which characterize the seismic response of the reservoir, is dened as: Z = b V p = s b K + 4 3 G = s b K + 4 3 2 b V 2 s (2.33) where b is the bulk density,V p is the P-wave velocity,K is the bulk modulus, andG is shear modulus. The seismic impedance can be modeled from reservoir parameters by Gassmann (1951) and Han et al. (1986) equations. The bulk density b can be dened as a volume average of the solid density and uid density: = ( o S o + w S w + g S g ) + (1) solid (2.34) where o , w , and g are densities of oil, water, and gas, respectively; S o , S w , and S g are the saturations of the corresponding phases; is porosity, and solid is the 36 density of solid. According to Gassmann (1951)'s theory, the bulk modulus K of the saturated rock can be estimated by the following equations: K =K grain K frame +Q (K grain +Q) (2.35) Q = K fluid (K grain K frame ) (K grain K fluid ) (2.36) whereK frame ,K grain , andK fluid are the moduli of rock frame, solid grain, and uid, respectively. The modulus of rock frame is computed by the grain modulus and porosity: logK frame = logK grain 4:25 (2.37) and the grain modulus is estimated by the moduli of sand, clay, and the shaliness, if the rock consists of shaly sand: K grain = 1 2 " K c + (1 )K s + K s K c K s +K c (1 ) # (2.38) where K s is sand modulus, K c is clay modulus, and is shaliness. The modulus of the mixed pore uid K fluid can be calculated by a weighted harmonic average of the individual uid bulk modulus: K fluid = 1= S w K w + S g K g + S o K o ! (2.39) where K w , K g , and K o are bulk modulus of water, gas, and oil. The shear wave velocity can be estimated by a correlation to porosity and shaliness (Han et al., 1986): V s = 3520:0 4910:0 1890:0 (2.40) Therefore, with Eqs. 2.33-2.40, the seismic impedance is a function of uid satura- tions if the other parameters are constant. In Fig. 2.1,, seismic impedance is plotted against water saturation for various porosity values, assuming gas saturation is zero 37 and the other parameters are constant. For any given porosity value, the impedance increase with water saturation, because of the low density oil is substituted by high density water. This relationship is very important in time-lapse seismic monitoring, where the substitution of oil by water can be clearly observed from the dierence between a monitor survey and the reference survey. Figure 2.1: Seismic impedance vs. water saturation for various poros- ity values, assuming gas saturation is zero and the other parameters are constant. In ensemble-based history matching, when seismic impedance data is assimilated, the rock physics model based on Eqs. 2.33-2.40, is solved with the reservoir simula- tion to provide predictions of the seismic surveys. The uid saturations updated by reservoir simulation, as well as the updated porosity, are inputs for the rock physics model. When time-lapse seismic data is used, the ensemble Kalman smoother is usually used to update the reference survey based on the updated reservoir static parameters (porosity, shaliness, etc.). 38 Chapter 3 Hough-Transform-Based Parameterization In this study, a Hough-transform-based parameterization for facies models in frac- tured reservoirs is developed to assist the history matching and uncertainty quanti- cation. The main motivation of this work is to improve the capability and eciency of the ensemble based history matching techniques for fractured reservoirs, while assuring the results from the history matching are geologically realistic. 3.1 Hough Transform The Hough transform is a feature extraction technique introduced by Hough (1962), and improved by Duda and Hart (1972). It is now widely used in image analysis, computer vision, and digital image processing. The technique identies instances of objects (straight lines, curves, ellipses, etc.) in a image by a voting procedure. The result of the voting procedure is stored in a so-called accumulator space, where objects can be identied by local maxima. Therefore, an image with discrete objects can be transformed into a relatively continuous accumulator map. Conventionally, a straight line in an image can be characterized by the slope- intercept equation: y =mx +b (3.1) 39 wherem is the slope of the line andb is the y-intercept. This denition of lines is quite straight forward, but is computationally inecient when being used in digital image processing, because the dimension of the parameter space can be quite large. The range of b can be from1 to1, and the vertical lines are dicult to characterize. In the Hough transform, lines are described with angle-radius rather than the slope- intercept parameters to simplify the computation. As shown in Fig. 3.1, when a xed origin is dened in an image (usually the bottom left corner), a straight line can be specied by the angle of its normal and its algebraic distance from the origin. The angle-radius equation of a straight line is: =xcos +ysin (3.2) If is restricted in the interval [=2;=2], every straight line in the x-y plane can be uniquely represented by a point (;) in the - plane, which is referred to as Hough space for the set of straight lines in two dimensions. X y θ ρ o Figure 3.1: Angle-radius parameterization of a straight line. is the angle of its normal, and is its algebraic distance from the origin o. Using this parameterization, an arbitrary pixel (x p ;y p ) on the image plane can be 40 transformed into a sinusoidal curve in the Hough space by: =x p cos +y p sin (3.3) while each point (;) on the curve denes a line passing through point (x p ;y p ). On the other hand, the points on the same line in the image plane will be transformed into sinusoidal curves that cross at the point ( ; ) which denes that line. For instance, in Fig. 3.2a, three pixels (P 1, P 2, and P 3) on the same straight line in an image can be transformed by Eq. 3.3 into three sinusoidal curves in the Hough space, as shown in Fig. 3.2b. The three curves cross at the same point Q(0:9; 20:0), indicating that the three pixels are co-linear points on that straight line. Therefore, the problem of detecting co-linear points can be converted to the problem of nding concurrent curves (Duda and Hart, 1972). Figure 3.2: (a) Pixels in a binary image are transformed into (b) sinusoidal curves in the Hough space by Eq. 3.3. The curves belonging to co-linear pixels will cross at a point in the Hough space which represents the straight line through the pixels. 41 Based on the parameterization of Hough transform in Eq. 3.3, the Hough trans- form is dened as: A(;) = Z +1 1 Z +1 1 I(x;y)(x cosy sin)dxdy (3.4) where A(;) is dened as the accumulator (Duda and Hart, 1972); I(x;y) is the original binary image; and is the Dirac delta function. Therefore, the accumulator A(;) gives the total number of sinusoidal curves that cross at point (;), and hence, the total number of points lying on the line dened by (;). If there is a line segment on an image, which can be represented by a point ( ; ) in the Hough space, all of the pixels on that line segment will contribute to the value of accumulator A( ; ). This line segment can be identied by the outstanding value of A( ; ). For instance, if the straight line in Fig. 3.2 is digitalized into an image with 40 40 pixels, an accumulator array can be obtained by Eq. 3.4, as illustrated in Table 3.1. The accumulator array is large for the entire Hough space, thus only part of the array is shown in the table. The largest entry in Table 3.1 occurs at (0:9; 20:0), corresponding to the straight line in Fig. 3.2a. An inverse Hough transform can be accomplished by nding peaks in the function A(;) in Hough space and drawing the straight lines represented by the peaks. 42 (Rad) 0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 0.0 2.5 1 1 5.0 3 4 7.5 2 1 4 10.0 3 3 5 2 12.5 3 4 3 5 4 15.0 2 4 5 7 5 4 17.5 3 3 5 9 8 5 4 20.0 3 4 5 10 31 10 5 2 22.5 3 3 5 9 6 5 4 25.0 2 3 5 3 1 2 27.5 3 4 3 30.0 3 2 32.5 35.0 Table 3.1: Accumulator array for the binary image in Fig. 1. Only part of the array within 2 [0:1; 1:5] and 2 [0:0; 35:0] is provided. The largest entry in the array occurs at (0:9; 20:0), corresponding to the straight line in Fig. 1. 3.2 Hough-Transform-Based Parameterization for Fractured Reservoirs In typical commercial reservoir simulators in petroleum industry, the reservoirs are discretized into regular grid blocks, and a set of parameters are assigned to each of the grid blocks. However, in reservoirs with more than one facies type, such as fractured reservoirs and channelized reservoirs, facies models are often used in the forecast and history matching. The facies maps can better represent the geological structures in 43 these reservoirs. However, it is dicult to estimate the facies maps by the ensemble based history matching methods, such as the EnKF. The primary problem is that the statistical distribution of facies maps is far from the Gaussian distribution, which is optimal for the performance of EnKF. Another problem, especially for fractured reservoirs, is that the geological structures in the facies maps are hard to be con- served after EnKF updates without proper constraints. Therefore, since the Hough transform can transform a binary image with arbitrary lines into a relatively con- tinuous parameter space, we consider using the parameters in the Hough space to represent the original facies maps. The angle-radius representation of straight lines in the Hough transform is utilized in this parameterization. The advantage of this pa- rameterization is its computational simplicity for arbitrary fractures in a facies map, and this parameterization acts as a constraint for the shape of fractures to ensure that the updated facies map is geologically realistic. For simplicity, the two-dimensional facies map of a fractured reservoir is treated as a binary image, where 0 means rock matrix and 1 is fracture. The petrophysical parameters, such as permeability and porosity, are assumed known and constant in each facies. Fractures are considered as line segments in the binary facies map. Each fracture can be determined by its Hough parameter and , as well as its length L and axial displacement D, as illustrated in Fig. 3.3. To accommodate for the common commercial reservoir simulators, in which the grid block at the upper left corner has the smallest index, the origin of a facies map is dened at the upper left corner, as shown in Fig. 3.3. The Hough parameter, the line segment lengthL, and the axial displacementD are dimensionless parameters, being the ratios of the actual length to the diagonal L d of the facies map. Thus, the maximum ranges of the rst three parameters are 2 [=2;=2], 2 [1; 1], and L2 (0; 1], respectively. The range of D is dependent on the other three parameters. 44 X y θ ρ o D L Figure 3.3: The Hough-transform-based parameterization for a fracture. is the angle of its normal; is its algebraic distance from the origin o; L is its length; and D is the axial displacement of its midpoint from its perpendicular foot. In Hough space, fractures are dened by a vector of parametersfH(;);L(;);D(;)g T . H(;) is an analog to the accumulator, but does not necessarily have the same si- nusoidal structure as dened in Eq. 3.4 because only the angle-radius representation of straight lines is adopted in this parameterization. So, H(;) is named as Hough function in this context to distinguish it from the accumulator in the original Hough transform. The Hough function acts like an indicator of the existence of fractures on a facies map. A fracture represented by (;) exists on a facies map only when the corresponding H(;)> 0. Thus, fractures can be added into or removed from a facies map by adjusting the values of the Hough function. To be computationally feasible for history matching techniques, the new parameter vectors are dened on a discrete Hough space. On a nm discrete Hough grid with indices i = 1; 2;:::;n and j = 1; 2;:::;m, the parameter vector at grid block (i;j) becomes: ' i;j = 8 > > > < > > > : H i;j L i;j D i;j 9 > > > = > > > ; (3.5) For any grid block (i;j), ifH i;j > 0, a fracture dened by i , j , and' i;j exists, where 45 i and j are the Hough coordinates at the center of grid block (i;j). Under this parameterization, any fracture on the facies map is uniquely dened on the Hough space by the parameter vector. The two ends (x 1 ;y 1 ) and (x 2 ;y 2 ) of the fracture on the facies map are computed by: x 1 = j L d cos i + (D(i;j) L(i;j) 2 )L d sin i (3.6) y 1 = j L d sin i (D(i;j) L(i;j) 2 )L d cos i (3.7) x 2 = j L d cos i + (D(i;j) + L(i;j) 2 )L d sin i (3.8) y 2 = j L d sin i (D(i;j) + L(i;j) 2 )L d cos i (3.9) With these two points, a fracture can be added into the facies map by setting the grid blocks intersected by the straight line between (x 1 ;y 1 ) and (x 2 ;y 2 ) as fracture facies. Therefore, a general inverse transform from a given ' eld on Hough space to a facies map consists of the following steps: 1. Initially, set all grid blocks in a facies map to matrix. 2. Identify' vectors in Hough space whose Hough function values are larger than 0. 3. Set corresponding grid blocks in the facies map to fracture facies given the vectors identied in Step 2, with the help of Eqs. 3.6-3.9. However, when the dimensions of the discrete Hough space are large, additional steps can be added into the inverse transform to constrain the abundance of fractures in facies maps, which will be introduced later in this dissertation. As discrete Hough space is used in this parameterization, there exists inconsis- tency between fracture characterization in continuous Hough space and that in low resolution discrete Hough space. For example, an arbitrary fracture in the ow do- main is characterized by the parameter vector rounded to the nearest discrete Hough 46 grid. Therefore, the fracture in the facies map generated by the inverse transform us- ing the rounded parameter vector can be slightly dierent from the original fracture, resulting in errors in the facies maps. This facies map error FE can be expressed as: FE = lim N c f !1 2 4 1 N c f N c f X k=1 N 1 X i=1 N 2 X j=1 i;j ( k ; k ) i;j ( ~ k ; ~ k ) 3 5 (3.10) where N c f is the number of fractures dened in continuous Hough space, i;j (;) is the facies index at grid block (i;j) in the facies map generated for the fracture (;), ( e ;e ) is the parameter vector (;) when rounded to the grid block in discrete Hough space, N 1 and N 2 are the numbers of grid blocks in ow domain along x and y direction, respectively. A discrete Hough space with very high resolution, instead of the continuous Hough space, is used to approximate this facies map error: g FE = 1 N hr f N hr f X k=1 N 1 X i=1 N 2 X j=1 i;j ( k ; k ) i;j ( ~ k ; ~ k ) (3.11) where N hr f is the total number of fractures dened in the high resolution discrete Hough space, ( e ;e ) is the parameter vector (;) in high resolution Hough space when rounded to the grid block in lower resolution discrete Hough space. This g FE is tested for N 1 N 2 ow domain with various Hough space resolution as N N . For instance, in Fig. 3.4, g FE is plotted for a 40 40 ow domain with dierent Hough space resolutions in whichN H =N =N . This trade-o curve reveals the relationship between the ow domain resolution and the Hough space resolution: with low Hough space resolutions, some of the fractures cannot be char- acterized correctly by the Hough space parameter sets by the rounding errors; while with high Hough space resolutions, it is dicult to distinguish in facies maps the fractures with similar characteristics. 47 Figure 3.4: The relationship between the approximation of facies map error g FE and Hough space resolution N H , for a 40 40 ow domain. N H is the resolution of a N H N H Hough space. Figure 3.5: The relationship between the approximation of facies map error g FE and Hough space resolution N H , for (a)60 60 ow domain and (b)80 80 ow domain. N H is the resolution of a N H N H Hough space. The trade-o curve remains similar in shape for dierent sizes of ow domains. As shown in Fig. 3.5, the g FE curves for a 60 60 ow domain (Fig. 3.5a) and a 80 80 ow domain (Fig. 3.5b) indicate the same g FEN H relationship as for the 40 40 ow domain. The shape of the curve is conrmed by extensive tests for 48 various ow domain resolutions. A recommended Hough space resolution is selected as N H = 100 in Fig. 3.4, yielding a 100 100 Hough space for a 40 40 ow domain. Such criteria ensures both the accuracy in characterizing the fractures and the eciency of the history matching with a relatively small parameter space. Another factor that controls the facies map error g FE is the aspect ratios of the ow domain grid and Hough space grid. Upon extensive tests, the trade-o curves are similar for dierent aspect ratios of the ow domain, as shown in Fig. 3.6, where two examples are demonstrated for a 20 40 ow domain (Fig. 3.6a) and a 40 20 ow domain (Fig. 3.6b). It can also be concluded from the tests that the facies map error g FE is least when N =N for any ow domain resolutions. Two examples are demonstrated in Fig. 3.7 for the relationship between g FE and the aspect ratio of Hough space grid. Figure 3.6: The relationship between the approximation of facies map error g FE and Hough space resolution N H , for (a)20 40 ow domain and (b)40 20 ow domain. N H is the resolution of a N H N H Hough space. The criteria is applied to various ow domain resolutions and the recommended Hough space resolutions are shown in Fig. 3.8. There is a clear linear relationship 49 Figure 3.7: The relationship between the approximation of facies map error g FE and Hough space resolution N H , for (a)20 40 ow domain and (b)40 20 ow domain. N H is the resolution of a N H N H Hough space. between the recommended Hough space resolution and ow domain resolution: N H = 2:5Max(N 1 ;N 2 ) (3.12) whereN H is the resolution of aN H N H Hough space, recommended for aN 1 N 2 ow domain. Although the recommended Hough space discretization yields more grid Figure 3.8: The linear relationship between the recommended Hough space resolution and ow domain resolution. N H is the resolution of a N H N H Hough space, recommended for a N 1 N 2 ow domain. 50 blocks than the ow model, the additional parameters are necessary to characterize the fractures as geometric objects in fractured reservoirs, and the number of param- eters in proposed parameterization will only grow linearly with the complexity of the reservoir model. It is also worth mentioning that Eq. 3.12 is developed by considering only the balance between the acceptable precision of fracture representation and sizes of the parameter space, independent of the dynamic data to be assimilated during history matching. Since the relationship between the facies map error and the Hough space resolution is a trade-o curve, as demonstrated in Fig. 3.4, it is still possible to slightly increase or decrease the Hough space resolution depending on the abundance of dynamic data to be assimilated. For instance, if large amount of production data with good spatial coverage is available in history matching, the coecient in Eq. 3.12 can be increased above the recommended 2:5 to further improve the accuracy of frac- ture characterization. However, rather than trading accuracy for smaller parameter space by using a coarser Hough space resolution when dynamic data is insucient, it is better to incorporate prior fracture statistics in order to signicantly reduce the parameter space without losing accuracy, which will be discussed later in this paper. In addition, not all grid blocks in the Hough space are active during the history matching, since the fractures represented by parameter vectors in some of the grid blocks always lie outside of the ow domain. Thus, a grid block (i;j) in Hough space is active if its coordinates ( i ; j ) satisfy: 8 < : j 2 h Ly L d sin i ; Lx L d cos i i ; i < 0 j 2 h 0; cos 4 i i ; i 0 (3.13) where L x and L y are the lengths of the ow domain along x and y direction, respec- tively; andL d is the diagonal of the ow domain. Instead of the 10,000 grid blocks in Hough space, only 4,504 grid blocks are active, as shown in Fig. 3.9. Only the pa- rameter vectors in these active grid blocks are updated in history matching, reducing the total number of parameters to be updated. 51 Figure 3.9: The Hough space is discretized into a 100 100 Cartesian grid. The black colored grid blocks are active during the history matching. 3.3 Implementation of Hough-Transform-Based Pa- rameterization in Ensemble-Based History Match- ing With the Hough-transform-based parameterization, the Hough-transform-based pa- rameter vectors in the Hough space, instead of the facies maps, are updated by the EnKF or ES. So the state vector in Eq. 2.1 becomes: y n = 8 > > > > > > > > > < > > > > > > > > > : [H] [L] [D] u n g n 9 > > > > > > > > > = > > > > > > > > > ; 2R N H +N L +N D +Nu+Ng (3.14) where the brackets [ ] indicate the collection of parameters assigned to all the grid blocks in the Hough space;u is the dynamic variables, such as pressure and saturation; 52 and g is the additional dynamic measurements, such as oil rate, water cut, etc. An inverse transform is performed to generate facies maps for the reservoir simulator whenever a simulation run is needed in the history matching process. In the context of EnKF/ES, the Hough space parametersH,L, andD become ran- dom variables. Performing a similar role to the levelset function (Chang et al., 2010; Ping and Zhang, 2013, 2014), the Hough function can be an independent Gaussian random variable. However, for a typical reservoir model, the Hough space resolu- tion can be quite large. If the Hough function is an independent standard Gaussian random variable, statistically the average number of fractures in each realization is (1 P (H 0))N H , where P (H 0) is the cumulative density function of the Gaussian variable H, and N H is the number of grid blocks in the Hough space. For instance, if there are 10,000 grid blocks in the Hough space, and the possibility that the Hough functions are larger than 0 is only 5%, then there will be 500 fractures in each realization. This can cause problems, such as poor covariance matrix in Eq. 2.7 or 2.10, and low eciency of the lter/smoother. The probability (1P (H 0)) cannot be too small, especially with a limited ensemble size in history matching; oth- erwise, that probability will have large uctuation due to sampling error. Therefore, the Hough function for ne Hough space grids is considered to be spatially correlated random elds, e.g., a Gaussian random eld with an exponential covariance function: C(x; y) = 2 exp 3jx y j 3jx y j ! (3.15) where C(x;y) is the autocovariance function between point x at (x ;x ) and y at (y ;y ) on the Hough space; 2 is the variance; and is the correlation length or practical range. The spatially correlated Hough function maps are relatively smooth elds with a much smaller number of peaks than the total number of grid blocks. A fracture dened by grid block (i;j) in the Hough space exists if its Hough function satises: H i;j > 0; and M(H i;j ) = 1 (3.16) whereM(H i;j ) = 1 whenH i;j is a local maximum comparing with the Hough functions 53 in the nearby grid blocks, otherwise M(H i;j ) = 0. Consequently, the average number of fractures can be easily constrained by adjusting the expectation of Hough function and its correlation length. A small expectation and large correlation length will result in fewer fractures in facies maps. The steps of inverse transform then become: 1. Initially, set all grid blocks in a facies map to matrix facies. 2. Identify ' vectors in Hough space that satisfy the conditions in Eq. 3.16. 3. Set corresponding grid blocks in the facies map to fracture facies given the vectors identied in Step 2, with the help of Eqs. 3.6-3.9. The work ow of history matching using the iterative smoother LM-EnRML and Hough-transform-based parameterization is shown in Fig. 3.10, as an example about how the parameterization is implemented in ensemble-based history matching meth- ods. The initial realizations of the Hough function, for instance, can be generated by Karhunen-Lo eve expansion (Zhang and Lu, 2004). The two dimensional covariance function in Eq. 3.15 can be decomposed into: C(x; y) = 1 X n=1 n f n (x)f n (y) (3.17) where n andf n are called eigenvalues and eigenfunctions, respectively, which can be solved analytically for a rectangular random eld on a regular Cartesian grid. The Hough function at point x (x ;x ) can be expressed in terms of the eigenvalues and eigenfunctions as: H(x) =m + 1 X n=1 n q n f n (x) (3.18) where m is the mean of the Hough function; and n are standard Gaussian random variables. Even though there are innite terms in Eq. 3.18, only the rst N terms are needed to generate the random eld, as discussed by Zhang and Lu (2004). Therefore, the initial Hough function random elds can be generated by independently sampling 54 Figure 3.10: The work ow of history matching using the iterative smoother LM-EnRML and Hough-transform-based parameterization. 55 a certain number of n from N(0; 1) and then populating into: H(x) =m + N X n=1 n q n f n (x) (3.19) The fracture length L and axial displacement D random elds are spatially in- dependent. Two fractures dened by two dierent parameter vectors on the Hough space should have irrelevant lengths and axial displacements. So, the initialL andD random elds are constructed by independently sampling Gaussian random numbers. 3.4 Integration of Prior Information Due to the strong non-linear relationship between the fracture distributions and the production history, as well as the limited observations at the sparsely distributed wells, multiple fracture distributions might yield equally good match to the observed data. Hence, prior information about the fracture distributions is considerably useful to reduce the uncertainty of the updated distribution and improve the convergence of the history matching algorithms. In this study, we consider two kinds of prior information: (1) Knowledge of the orientations of the fractures that pass through wells. The orientation of a fracture passing through a well can be estimated from core samples and well logs, e.g., downhole image logs. In such cases, that fracture is considered to exist in the facies map with certainty, and its Hough space coordinate 0 is known from the well logs. Its 0 is determined by: 0 =x w cos 0 +y w sin 0 (3.20) where x w and y w are the coordinates of the well in the facies map. However, the lengthL 0 and axial displacementD 0 of this fracture are unknown, except that it passes through the well. Thus, the initialL 0 in the corresponding grid block is still generated by sampling Gaussian random numbers, with the same statistics as the rest of the 56 Hough space. To ensure that the fracture passes through the well,D 0 is consequently dependent on the value of L 0 , satisfying D 0 2 N(m D ; 2 D ) and D 0 2 (D min ;D max ), where: m D = 8 < : xw 0 cos 0 sin 0 ; 0 < 4 or 0 > 4 0 sin 0 yw cos 0 ; 4 0 4 (3.21) D = L 0 6 (3.22) D min =m D L 0 2 (3.23) D max =m D + L 0 2 (3.24) Whenever a facies map is generated by the inverse transform, the known fractures are added independently from the others identied from the Hough function elds. (2) Known statistics of the fracture orientation in the ow domain. The gen- eral statistics of fracture orientation can be estimated from the regional earth stress regimes and information collected from the existing wells. Therefore, only part of the Hough space grid blocks will be used in the inverse transform, and updated during the history matching. For example, if the most possible direction of fractures in a reservoir is NE-SW with 2 [=4!;=4 +!], where ! is a small radius indicating the range of the orientations, only the grid blocks with 2 [=4!;=4 +!] in the Hough space are active during the history matching, and both the initial and updated fractures are in the NE-SW direction. This kind of prior information not only reduces the number of parameters in the Hough space, but also improves the history matching performance by adding extra constraints to the parameters. 3.5 Covariance Localization The ensemble-based history matching techniques, such as the EnKF and ES, ap- proximates the covariance matrix from an ensemble, which is often smaller than the number of state variables. As discussed in Chapter 2, this approximation suers from two weaknesses: 1) spurious correlations resulting from the sampling error with small 57 ensemble size; 2) the rank of the approximation of the covariance is less than or equal to the number of ensemble members, thus the limited degree of freedom prevents EnKF to assimilate a large number of independent data. Distance-dependent covariance localization is used to reduce the in uence of spu- rious correlations and allow larger amount of data to be assimilated. The localization function is dierent in Hough-transform-based parameterization since the distance between a variable in the Hough space and the observation point in the ow domain cannot be calculated directly. Therefore, the exponential localization function is de- ned between the Hough space and the ow domain by the Hough-transform-based parameterization: (d H ) = exp " 3d H 4 ( p 9 + 8N e 5) # (3.25) where d H is the distance in Hough-transform-based parameterization, and is the practical range. For a point (;) in the Hough space, as illustrated on Fig. 3.11, the distanced H between the parameter at (;) and the observation in the ow domain is dened as the perpendicular distance from the observation to the fracture represented by (;). The equation for d H is derived in terms of (;) and (x; y), which is the X y θ ρ o dH (x ,y ) o o Figure 3.11: Distance between a Hough space grid (;) and a data point (x o ;y o ) in ow domain. 58 location of the observation in the ow domain: d H =j (x cos + y sin)j (3.26) For instance, in a ow domain discretized into 41 41 regular grid, there are nine observation wells as shown in Fig. 3.12. As suggested by Eq. 3.12, a 100 100 discrete Hough space is used to parameterize the fractures in the ow domain. Fig. 3.13 illustrates the nine localization function map for the Hough space parameters to the observations at the nine wells, respectively, with = 0:3. The shapes of the localization function maps look like sinusoidal bands, whose magnitude is least for observations near the ow domain origin, and is largest for observations farthest from the origin. Figure 3.12: The ow domaine is discretized into 4141 regular grid. Nine observation wells marked as black dots on this map provide observation data to be assimilated to update the Hough space parameters. 59 Figure 3.13: The localization function maps for the Hough space parame- ters to the observations at wells (a)O1, (b)O2, (c)O3, (d)O4, (e)O5, (f)O6, (g)O7, (h)O8, and (i)O9, respectively. 60 Chapter 4 History Matching for a Synthetic Fractured Reservoir Model In this chapter, the proposed history matching techniques and the Hough-transform- based parameterization are tested on synthetic cases with dierent reference elds representing the fracture distributions in the reservoir. Since synthetic reference mod- els are used to provide production observations, the capability of the proposed history matching technique in both matching the production data and estimating the correct fracture distributions is evaluated and analyzed. The rst case is designed to demon- strate the eectiveness of this history matching technique to estimate the position of a single fracture in the reservoir. The improvements from covariance localization are investigated in the second case. The third case is similar to the rst one, but with a reference eld with two connected fractures. In the rst three cases all fractures in the reference facies maps pass through the wells; thus, the production history at the wells is strongly in uenced by, and thus sensitive to, the fractures in the ow domain. In the fourth and fth cases, in which the fractures in the reference facies maps are either not or only partially connected by wells, the parameters are less sensitive to the observed data, and the assimilation of prior information is demonstrated. 61 4.1 Model Description The reservoir model used in the cases in this chapter is a two-dimensional water ood- ing model with nine wells. The ow equations are solved with Eclipse (Schlumberger, 2007). The ow domain is a square of size 820ft 820ft, uniformly discretized into 41 41 square grids, as shown in Fig. 4.1. The thickness of the eld is 50ft, and the top of the eld is at the depth of 4,000 ft. All four sides of the eld are considered no ow boundaries. An injector (I in Fig. 4.1) is located in the center of the eld, injecting water at a xed rate of 1,000 B/D. Eight producers (P1 to P8 in Fig. 4.1) are located at the corners and the middle of the sides of the eld, producing at xed bottom hole pressure of 3,000 psi. Initially, the eld has a uniform pressure of 4,000 psi and a uniform oil saturation of 0.8. The total duration of the production is 30 days, divided evenly into 30 time steps, and observations of the oil rate and water cut at the eight producers and the Figure 4.1: The ow domain and well locations of the synthetic cases. Oil rate observations are available at all eight producers denoted by lled circles, and bottom hole pressure is available at the injector in the center of the domain. 62 bottom hole pressure at the injector are available at the end of each time step. The eld has only two facies: fracture (more specically, highly fractured strip consisting of intensive cluster of fractures) and rock matrix. Grid blocks inside the fracture facies have constant eective permeability and porosity of 30,000 md and 0.2, respectively. The permeability and porosity in grid blocks inside the rock matrix facies are 3 md and 0.1, respectively. The facies map of the eld is discretized in the same way as the ow domain, so the permeability and porosity of each grid block can be easily determined by assigning values with respect to the facies type of that grid block. The Hough space of the cases is a rectangular space with 2 [=2;=2] and 2 [1; 1]. As suggested by Eq. 3.12, the Hough space is discretized into a 100100 Cartesian grid, as shown in Fig. 4.2. However, not all grid blocks in the Hough space are active during the history matching, since the fractures represented by parameter vectors in some of the grid blocks always lie outside of the ow domain. Thus, a grid Figure 4.2: The Hough space is discretized into a 100 100 Cartesian grid. The black colored grid blocks are active during the history matching. 63 block (i;j) in Hough space is active if its coordinates ( i ; j ) satisfy: 8 < : j 2 h Ly L d sin i ; Lx L d cos i i ; i < 0 j 2 h 0; cos 4 i i ; i 0 (4.1) where L x and L y are the lengths of the ow domain along x and y direction, re- spectively; and L d is the diagonal of the ow domain. Out of the 10,000 grid blocks in Hough space, only 4,504 grid blocks are active, as shown in Fig. 4.2. Only the parameter vectors in these active grid blocks are updated in history matching. The initial Hough functionH(i;j) is generated by the Karhunen-Lo eve expansion with mean m =0:5, variance 2 = 0:5, and correlation length = = 0:3[L], where L is the length of the Hough space in the same direction of the correlation lengths. The average number of fractures is around ve in each initial realization. The initial fracture length L(i;j) and axial displacement D(i;j) on the Hough grid are generated by independently sampling Gaussian random numbers. The statistical mean of L(i;j) is 0.7 and the variance is 0.05. The initial distribution of D(i;j) is more complicated than that of L(i;j), due to the fact that the mean of the axial displacement of a fracture is dependent on the angle of its normal, in order to place the fractures inside of the ow domain. So, the mean of D(i;j) is determined by: D(i;j) = 0:5 sin( i 4 ) (4.2) and the variance of D(i;j) is 0.03. 4.2 History Matching Methods In fractured reservoir models, a slight change in the fracture distribution will cause signicant changes in the production history. The well performance can be completely dierent whether the well is connected by a fracture or not. The strong non-linear relationship between model parameters, dynamic variables, and production measure- ments makes it dicult to sequentially update the dynamic variables, i.e., pressure and saturations in EnKF. It is almost necessary to run the simulation from the initial 64 state to obtain the updated dynamic variables. As a result, the computational com- plexity in terms of the eort spent in reservoir simulation in EnKF will be O(N 2 ), which is expensive for large scale eld applications. Thus ES based methods are pre- ferred because all data are assimilated at a single step, with a complexity ofO(N). In addition, it is generally necessary to iterate the ES in order to match the production data (Chen and Oliver, 2012), especially for strongly non-linear problems. Therefore, the LM-EnRML (Chen and Oliver, 2013) is used in most of the illustrative cases for history matching. In some occasions when the observed data are inadequate to de- termine the fracture distributions, the iterative methods will suer from convergence issues. In such cases, the EnKF is still a viable method to sequentially update the model parameters. 4.3 Case 1: Eectiveness of The Hough-Transform- Based Parameterization in a Simple Reference Field In this case, a synthetic reference eld is designed to test the eectiveness of the proposed methodology to estimate fracture distributions in petroleum reservoirs. The facies map of the reference eld is demonstrated in Fig. 4.3. The black grid blocks indicate fracture facies. The producer P3 and injector I are connected by a fracture. Synthetic observed data are obtained by simulating the reservoir model with the reference facies map in Fig. 4.3a for 30 days. Only the observations in the rst ve timesteps are used in history matching, and they are assimilated at a single step in the LM-EnRML. No prior information is available in this case. The iterative ES converged after 17 steps of iteration. A total of 31 simulation runs (for the ve timesteps) were performed for each of the 100 realizations. Fig. 4.4 shows the reduction in the mean of data mismatch dened by: S d = 1 N e Ne X i=1 (F (e y i ) e d i ) T e R 1 (F (e y i ) e d i ) (4.3) 65 Figure 4.3: Case 1: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix; the producer P3 and the injector I are connected by a fracture. Figure 4.4: Case 1: Reduction in the mean of the data mismatch with iteration. 66 where F () denotes the reservoir simulation which compute the measurements from the model variables. The mean of data mismatch decreases substantially from 1:32 10 12 to 7:5410 6 . However, the reduction is not a linear process. The data mismatch declines slowly but steadily in the rst 15 steps because the initial guess is far from the solution. Since the ensemble updated after step 15 is close enough to the solution, the sensitivity matrix is less noisy and the objective function is more linear; thus, the iteration converges quickly to the solution in the last two steps. We computed and plotted the probability maps of fractures in Fig. 4.5 at the initial stage, after the 5 th iteration step, and after the 17 th iteration step when the iteration converges, with the equation: p(i;j) = 1 N e Ne X k=1 k (i;j); i = 1; 2;:::;N 1 ;j = 1; 2;:::;N 2 (4.4) where k is the facies map of the k th realization generated by the inverse transform from the ensemble. The probability maps are used to indicate the probability of a grid block in the ow domain to be inside of the fracture facies. In the initial state, the fractures from all realizations are uniformly distributed over the ow domain, rendering a at probability map with low values, as shown in Fig. 4.5a. After ve iteration steps, the probability map in Fig. 4.5b shows higher Figure 4.5: Probability maps of fractures for Case 1: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 17 th iteration step when the iteration converges. 67 probability of fractures in the area between P 3 and I, where the fracture in the reference facies map is located, than the rest of the eld. After the history matching, in Fig. 4.5c, the probability map exhibits a pattern of fracture distribution very similar to the reference facies map, with very high probability values (> 90%) at the grid blocks connecting P 3 and I. In addition to the orientation and position of the fracture, the updated probability map provides accurate estimation of its length, with the lower left tip of the fracture placed exactly on the injector. It can be noticed that there are other possible fractures with relatively high probability values on the probability map, such as those near P 1 and those near P 8. The production history, especially that of the oil rate and water cut, is not very sensitive to the existence of those short fractures. With additional observed data, it will be possible to rule out those articial fractures from the updated probability map. The initial facies maps of three randomly selected realizations are shown in the rst row in Fig. 4.6, and the updated facies maps are shown in the second and third rows in Fig. 4.6. Although the initial facies maps show randomly distributed fractures, all three updated facies maps contain a fracture similar to the one in the reference facies map. The short fractures in the lower right corner of the three updated facies maps are dierent, indicating relatively high uncertainty due to the low sensitivity of these fractures to the observed data. The mean Hough function eld at the initial stage, after the 5 th iteration step, and after the 17 th iteration step when the iteration converges, is plotted in Fig. 4.7. These mean elds are computed simply by taking the arithmetic averages of Hough functions over all the realizations in the ensemble. At the initial stage, the mean Hough function eld is nearly uniform at0:5, owing to the stationarity of the spatially correlated Gaussian random elds. According to the reference facies map, the fracture in Fig. 4.3 is represented by the parameter vector at point p (=4; 0:5) in the Hough space. After the 5 t h iteration step, the mean Hough function in Fig. 4.7b shows high values near pointp . The updated mean Hough function in Fig. 4.7c shows a obvious peak exactly at pointp , indicating the high probability of a fracture represented by the parameter vector at point p . The high values of Hough function below the peak at pointp in Fig. 4.7c represent the short fractures at the lower right 68 Figure 4.6: Facies maps of three randomly selected realizations for Case 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 69 Figure 4.7: Mean Hough function maps for Case 1: (a) at the initial stage, (b) after the 5 t h iteration step, and (c) after the 17 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. corner of the updated facies maps. The updated parameters were used to simulate the reservoir model for another 25 days with the same condition to investigate the predictability of the updated models. In Fig. 4.8, the model predictions for the production history before and after the updating are provided side by side to demonstrate the ability of the methodology to improve the model's predictability. The black colored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched parameters. The blue asterisks are observations. Figs. 4.8a and 4.8b are model predictions for the oil rate and water cut at P 2, respectively. There are no fractures passing through P 2; thus, the observed oil rates are very low and the observed water cut is 0 throughout the production period, as a result of the low matrix permeability. The gures demonstrate the randomness of the model predictions from the initial ensemble and the strong nonlinearity of the model in which the shapes of the oil rate curves are completely dierent whether or not the well is connected by a fracture. After updating, the match between the model predictions and observations is signicantly improved. The majority of the oil rate, water cut, and bottom hole pressure curves from the realizations match the observation curves, and the variations of the predictions are very small compared to those before updating. As for the updated predictions at P 3, where the fracture passes through, both the oil rate and 70 Figure 4.8: The predictions from all the ensemble members before and after the history matching for Case 1: oil rate at (a) P 2 and (c) P 3; water cut at (b)P 2 and (d)P 3; (e) the bottom hole pressure at the injector. The black colored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched pa- rameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. The vertical axes in the updated predictions in (a), (c), and (e) are adjusted to show graphs clearly. 71 water cut are also matched with the observed data. The match is also considerably improved for the bottom hole pressure at the injector, as is shown in Fig. 4.8e. Not only the production data in the rst 50 days match the observations when they are used in the history matching, the update model is also capable to provide accurate predictions for the future period when the observations are not assimilated. 4.4 Case 2: Improvements from covariance local- ization This case is designed to investigate the improvement to fracture estimations by co- variance localization. The reference facies map is the same as in Case 1. No prior information is available, and the initial ensemble is generated in the same way as Case 1. Covariance localization is used in the LM-EnRML to reduce the spurious correlation, as discussed in Section 2.4. Distance-dependent localization is adopt in the form of an exponential localization function: (d H ) = exp " 3d H 4 ( p 9 + 8N e 5) # (4.5) whered H , the distance between a Hough space grid (;) and an observation at (x; y) in ow domain, is dened by d H =j (x cos + y sin)j (4.6) The iteration converged at the 24 t h step, and the reduction in the mean of data mismatch is shown in Fig. 4.9. The data mismatch is reduced steadily from 1:33 10 12 to 1:01 10 7 . Similar to Case 1, the data mismatch drops rapidly after iteration step 21 because the parameters are close to the solution. The probability maps computed by Eq. 4.4 at the initial stage and after history matching are shown in Fig. 4.10. Compared with Fig. 4.5c in Case 1, the updated probability map in Fig. 4.10c has a better match with the reference facies map: there are no artifacts with signicant probability on the updated probability map, when 72 covariance localization is used. It is easier to see the improvement by covariance localization on the mean Hough function maps (Fig. 4.11). In the updated mean Hough function map (Fig. 4.11c), there is only one signicant peak at (=4; 0:5), which represents the fracture in the ref- erence facies map. The rest of the mean Hough function map has much smaller values than this peak, while in Case 1 a peak appears at around (=4; 0:9) after updating due Figure 4.9: Case 2: Reduction in the mean of the data mismatch with iteration. Figure 4.10: Probability maps of fractures for Case 2: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 24 th iteration step when the iteration converges. 73 to false correlation (Fig. 4.7c). Therefore, covariance localization considerably allevi- ates the spurious correlation resulting from the relatively small ensemble, providing a better estimation of the fracture distribution with minimal artifacts. 4.5 Case 3: History Matching for Reference Field with Two Fractures In this case, the reference facies map contains two connected fractures to test the performance of the proposed method in more complicated geology, as illustrated in Fig. 4.12. The initial ensemble is generated in the same way as Case 1. Synthetic observed data are obtained by simulating the reservoir model with the reference facies map in Fig. 4.12 for 30 timesteps. The observations in the rst ve timesteps are used to update model parameters in the LM-EnRML. No prior information is available in this case. The iteration converged at the 20 t h step, and the reduction in the mean of data mismatch is shown in Fig. 4.13. The data mismatch is reduced steadily from 1:12 10 11 to 3:63 10 7 . Similar to Case 1, the data mismatch drops rapidly after iteration step 16 because the parameters are close to the solution. Although the geology is more complicated than Case 1, it did not take many more steps to converge. Figure 4.11: Mean Hough function maps for Case 2: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 24 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 74 The probability maps computed by Eq. 4.4 at the initial stage and after history Figure 4.12: Case 3: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix; the producer P 2 is connected to P 4 and P 5 by two fractures, respectively. Figure 4.13: Case 3: Reduction in the mean of the data mismatch with iteration. 75 matching are plotted in Fig. 4.14. The probability map after the 5 th iteration step (Fig. 4.14b) already shows large probabilities of fractures between P 2 and P 4, as well as those between P 2 and P 5. The updated probability map clearly shows the two fractures similar to those in the reference facies map. The estimated fracture connectingP 2 andP 5 is slightly dierent in the angle from the reference, due to the discretization of the Hough space and the ow domain. It is noticeable that there are some grid blocks with relatively high probability values, indicating the possibility of the existence of fractures, near the upper corners of the probability map although the production history of the wells is much less sensitive to these fractures. In Fig. 4.15a, the mean Hough function eld is uniform at the initial stage. The updated mean Hough function elds (Figs. 4.15b and 4.15c) show obvious peaks at p 1 (=4; 0:25) and p 2 (=4; 0:25), that represent the two fractures in the reference facies map (Fig. 4.12). The facies maps of three randomly selected realizations before and after updating in Fig. 4.16 also demonstrated that the history matching with the proposed param- eterization can successfully identify the main fractures conditional on the production history, although the initial fractures are randomly distributed in the facies maps. The model predictions with the updated parameters are shown in Fig. 4.17. Not only the oil rate and water cut at theP 2 (Figs. 4.17a and 4.17b) andP 5 (Figs. 4.17c Figure 4.14: Probability maps of fractures for Case 3: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 20 th iteration step when the iteration converges. 76 and 4.17d) that are connected by fractures, but also those at the other producers, e.g.,P 7 (Figs. 4.17e and 4.17f), are matched with the observed data available in the rst 5 days after the history matching. The ensemble is also proven to be able to predict the reservoir behavior for another 25 days under current conditions without the assimilation of additional observed data. 4.6 Case 4: Integration of Prior Information In this case, the reference facies map contains two connected fractures, while one of them is connected with P3, as illustrated in Fig. 4.18. Since the sensitivity of the fractures' characteristics to the production history is not as strong as the previous cases, additional prior information is introduced into the history matching process: (1) assuming that theP 3 is known to be connected by a fracture with the orientation of ==4; and (2) the fractures in this ow domain are known to be in the NE-SW and NW-SE directions, thus the active grid blocks in the Hough space is limited to those that satisfy 2f[0:9;0:6][ [0:6; 0:9]g. As a result, there are 997 active grid blocks in the Hough space. The observations in the rst 50 days are used in LM-EnRML to iteratively update the model parameters. The LM-EnRML converges after the 16 th iteration step, and the reduction in the mean of data mismatch is shown in Fig. 4.19. The mean of data mismatch is Figure 4.15: Mean Hough function maps for Case 3: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 20 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 77 Figure 4.16: Facies maps of three randomly selected realizations for Case 3: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 78 Figure 4.17: The predictions from all the ensemble members before and after the history matching for Case 3: oil rate at (a) P 2; (c) P 5, and (e) P 7, as well as water cut at (b) P 2, (d) P 5, and (f) P 7. The black colored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. The vertical axes in the updated predictions in (a), (c), and (e) are adjusted to show graphs clearly. 79 reduced to as low as 1:54 10 7 at step 11, indicating a faster convergence rate than Case 1 to Case 3, due to the smaller number of active Hough space grid blocks and the integration of prior information. The probability maps are shown in Fig. 4.20. Dierent from the last two cases, the probability map at the initial stage (Fig. 4.20a) is not uniform over the ow domain, in the presence of the prior information about the fracture passing through P 3. That fracture's orientation is known, but its length and axial position are un- certain, thus the probabilities of the fracture become smaller away from P 3. After the 5 th iteration step, the fracture that passes through P 3 is already identied with its correct length (Fig. 4.20b), although there is still some variation at the SW tip. However, there are no high probability grid blocks along the other fracture at this stage, because the sensitivity of the NE-SW fracture is much stronger. After the iteration converged, the NE-SW fracture is clearly identied on the probability map Figure 4.18: Case 4: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix. There is a fracture in NE- SW direction that passes through P 3, and another fracture in NW-SE direction that does not pass through any wells but connects the former one. 80 (Fig. 4.20c), although one end of this fracture is not connected to a well. The slightly lower probabilities in the grid blocks in the direction perpendicular to the former one indicate the existence of the second fracture, with larger uncertainties in both its orientation and length. The mean Hough function elds are plotted in Fig. 4.21. Although there are several peaks in the mean Hough function eld after the 5 th iteration step (Fig. 4.21b), there is only one signicant peak atp 1 (=4; 0:3) when updated (Fig. 4.21c). Figure 4.19: Case 4: Reduction in the mean of the data mismatch with iteration. Figure 4.20: Probability maps of fractures for Case 4: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 16 th iteration step when the iteration converges. 81 This peak represents the NW-SE fracture in the reference facies map. There is no noticeable peak at the point representing the NE-SW fracture that passes P 3, since that fracture is added as prior information with certainty, and the Hough function at point p 2 (=4; 0:25) is no longer sensitive to the production history, while the length and axial displacement still are. From the facies maps of three randomly selected realizations in Fig. 4.22, it is obvious that initially the fractures are in the two directions suggested by the prior information (Figs. 4.22a, 4.22b and 4.22c), while the positions of the fractures are still random on the facies maps. When the iteration converged, the three randomly selected facies maps (Figs. 4.22g, 4.22h and 4.22i) show the identical NE-SW fracture that passes throughP 3. However, variations in orientation, length, and position of the NW-SE fracture are observed on the three facies maps. Although the characteristics of the NW-SE fracture are not as sensitive to the production history as those that pass through wells, they are successfully identied by the history matching process with the estimation of their uncertainty due to the lack of decisive observed data. The model predictions at P 2, P 3, and P 7 before and after updating are plotted in Fig. 4.23. The predictions from the updated model match both the observed data in the rst 5 days that are used in history matching, and the observations in the other 25 days. As the NW-SE fracture in the reference facies map does not pass through any wells, the uncertainty in the updated model predictions is slightly larger Figure 4.21: Mean Hough function maps for Case 4: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 16 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 82 Figure 4.22: Facies maps of three randomly selected realizations for Case 4: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 83 than the last two cases. 4.7 Case 5: History Matching for a Reference Field with Fractures Away from Wells In this case, as seen in the reference facies map (Fig. 4.24), none of the two frac- tures pass through any wells. The sensitivity of the fractures' characteristics to the production history is even weaker than Case 4. The fractures are known to be in the NE-SW and NW-SE directions; thus, the active grid blocks in the Hough space are in the range of 2f[0:9;0:6][ [0:6; 0:9]g. With only the production history at the eight producers and the injector, this inverse problem with thousands of parameters is strongly ill-posed, rendering it dicult to obtain a unique solution to the fracture distribution. Consequently, the objective function can have multiple local minima (Aanonsen et al., 2009). The iterative methods encounter diculties in such prob- lems, suering from slow convergence rate, or even not being able to converge at all. In such situations, the EnKF is preferred as a history matching method to se- quentially update model parameters, when a smaller number of observations are as- similated at a time. Observations in the rst 20 timesteps are assimilated by EnKF, and the probability maps after the 20 th timestep of three separate EnKF assimila- tion runs are plotted in Fig. 4.25. The three updated probability maps indicate the relative positions of the two fractures similar to that in the reference facies map, and the fractures are generally within the upper right quadrant of the ow domain. However, the three probability maps suggested three dierent sets of fractures, while the predictions of the all three updated models match the production history equally well. In such cases, more timesteps of observed data and additional indirect measure- ments, e.g., tracer and seismic data, will help to reduce the uncertainty of the model parameters and improve the estimation of the fracture distributions. Our ongoing research is addressing this issue. 84 Figure 4.23: The predictions from all the ensemble members before and after the history matching for Case 4: oil rate at (a) P 2, (c) P 3, and (e) P 7, as well as water cut at (b) P 2, (d) P 3, and (f) P 7. The black colored sections of the lines are the period (rst 5 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 25 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. The vertical axes in the updated predictions in (a), (c), and (e) are adjusted to show graphs clearly. 85 4.8 Summary The Hough-transform-based parameterization is tested in two dimensional synthetic water ooding examples to evaluate its eectiveness to estimate the fracture distribu- tions. The impacts of the complexity of the reference fracture distribution, the sensi- tivity of the parameters to the observations, and the prior information, are examined in terms of the quality of the estimated fracture distribution and the predictability of the updated models. The examples indicate that the proposed method is an ecient and robust way to estimate fracture distribution in fractured reservoirs. The fracture distribution, in- cluding the number of fractures, the orientation, position and length of the fractures, are successfully estimated after the assimilation of production history at the wells. The estimated fracture distributions can be expressed as probability maps, which provide probabilities of the estimated fractures and are valuable for future reservoir Figure 4.24: Case 5: The reference facies map. Black grid blocks are fracture facies, and the rest are rock matrix. 86 Figure 4.25: Probability maps after updating for Case 5 from three sepa- rate EnKF assimilation runs. management and optimization. The sensitivity between the parameters and the pro- duction history depends on the distance between the fractures and the wells. High sensitivity is observed when the fractures pass through wells, while the sensitivity is lower when the fractures are away from the wells, rendering larger uncertainties in the estimated probability maps. Additional observations other than the production history at the wells may be essential to accurately estimate the fractures in the area away from the wells. Prior information about the fracture distribution plays an important role in the es- timation of fracture distributions, especially when the observed data are inadequate for an accurate estimation. Two kinds of prior information are assimilated in this chapter: (1) known fractures that pass through wells; and (2) statistics of regional fracture orientations in the ow domain. Both prior information considerably im- prove the updated estimation of the fractures, as well as the eciency of the history matching. The second type of prior information also signicantly reduces the number of parameters to be updated. The predictability of the model improves after the history matching. The updated model outputs match the production history assimilated during the history match- ing. With the help of the accurate characterization of the fractures by the proposed parameterization, the calibrated models can provide good reservoir responses and improved predictions for the periods beyond the available production history. 87 Chapter 5 Assimilation of Time-lapse Seismic Data with Hough-Transform-Based Parameterization In this chapter, the proposed methodology is tested on synthetic cases in which time- lapse seismic data is assimilated along with production data. With the help of the good spatial coverage of seismic data, the synthetic cases are designed to investigate the ability of the proposed parameterization to characterize the fractures that are not connected by wells. The rst case is designed to demonstrate the eectiveness of this history matching technique to estimate the position of a single fracture in the reservoir. The second case is more complicated with two connected fractures. In both cases all fractures in the reference facies maps do not pass through the wells; thus, the production history at the wells is not sensitive to the fractures in the ow domain, yielding uncertain interpretation of the position of the fractures. The time-lapse seismic data provides essential information about the uid migration in the entire reservoir during the production, which improves the characterization of the fractures. 88 5.1 Model Description The reservoir model used in the cases in this section is the same two-dimensional water ooding model as in Chapter 5. The ow equations are solved with Eclipse (Schlumberger, 2007). The ow domain is a square of size 820ft 820ft, uniformly discretized into 41 41 square grids, as shown in Fig. 5.1. The thickness of the eld is 50ft, and the top of the eld is at the depth of 4,000 ft. All four sides of the eld are considered no ow boundaries. An injector (I in Fig. 5.1) is located in the center of the eld, injecting water at a xed rate of 1,000 B/D. Eight producers (P1 to P8 in Fig. 5.1) are located at the corners and the middle of the sides of the eld, producing at xed bottom hole pressure of 3,000 psi. Initially, the eld has a uniform pressure of 4,000 psi and a uniform oil saturation of 0.8. The total duration of the production is 300 days, divided evenly into 30 time steps, and observations of the oil rate and water cut at the eight producers and the Figure 5.1: The ow domain and well locations of the synthetic cases. Oil rate observations are available at all eight producers denoted by lled circles, and bottom hole pressure is available at the injector in the center of the domain. 89 bottom hole pressure at the injector are available at the end of each time step. The eld has only two facies: fracture and rock matrix. Grid blocks inside the fracture facies have constant permeability and porosity of 30,000 md and 0.2, respectively. The permeability and porosity in grid blocks inside the rock matrix facies are 3 md and 0.1, respectively. The facies map of the eld is discretized in the same way as the ow domain, so the permeability and porosity of each grid block can be easily determined by assigning values with respect to the facies type of that grid block. The Hough space of the cases is a rectangular space with 2 [=2;=2] and 2 [1; 1]. As suggested by Eq. 3.12, the Hough space is discretized into a 100100 Cartesian grid. The initial Hough space parameter elds are generated in the same way as in Chapter 4. 5.2 History Matching Methods In addition to the nonlinearity of the models in fractured reservoirs, the assimila- tion of large amount of data is a challenge to the ensemble based history matching methods, when time-lapse seismic data is included in the observations. As discussed in Chapter 2, the matrix inversion in the Kalman gain is not computationally fea- sible if the amount of data is large. The iterative ensemble smoother LM-EnRML (Chen and Oliver, 2013, 2014) with covariance localization is used to assimilate both production data and time-lapse seismic data. The formulations with truncated sin- gular value decomposition avoid the explicit computation of the matrix inversion, and the localization discussed in Section 3.5 helps to reduce spurious correlation and rank-deciency. In order to assimilate seismic data, the rock physics model discussed in Section 2.5 is incorporated in the LM-EnRML to update the seismic predictions from petrophysics parameters. The parameters for the rock physics model are given in Table 5.1. The seismic data is assimilated as the dierence of acoustic impedance between the reference survey at initial state and the monitor survey, since the impedance dierence is better correlated with the saturation changes in the reservoir due to production. 90 Parameter Value Shaliness 0.2 Bulk modulus of sand (Pa) 3:8 10 10 Bulk modulus of clay (Pa) 2:12 10 10 Bulk modulus of oil (Pa) 6:71 10 8 Bulk modulus of water (Pa) 2:39 10 9 Density of solids (kg=m 3 ) 2650 Density of oil (kg=m 3 ) 707 Density of water (kg=m 3 ) 986 Table 5.1: Parameters used in the rock physics model for seismic impedance 5.3 Case 1: Assimilation of Time-Lapse Seismic Data in a Simple Reference Field In this case, a synthetic reference eld is designed to test the eectiveness of the proposed methodology to estimate fracture distributions with both production data and seismic data. The facies map of the reference eld is demonstrated in Fig. 5.2. The black grid blocks indicate fracture facies. There is one fracture between the producer P3 and injector I, but not connected to any wells. Synthetic observed data are obtained by simulating the reservoir model with the reference facies map in Fig. 5.2a for 50 days, with seismic impedance data available before the production and at 20 th day. Only the observations in the rst twenty days, as well as the seismic dierence data, is assimilated at a single step in the LM-EnRML. Prior information is introduced into the history matching process: the fractures in this ow domain are known to be in the NE-SW direction, thus the active grid blocks in the Hough space is limited to those that satisfy 2 [0:6; 0:9]. As a result, there are 499 active grid blocks in the Hough space. After 13 steps of iteration, the iterative ES converged. Fig. 5.3 shows the reduction in the mean of data mismatch. The mean of data mismatch decreases 91 substantially from 4:55 10 10 to 2:73 10 9 . The data mismatch declines slowly but steadily in the rst 9 steps and drops quickly in the last few steps. We computed and plotted the probability maps of fractures in Fig. 5.4 at the Figure 5.2: Case 1: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix. Figure 5.3: Case 1: Reduction in the mean of the data mismatch with iteration. 92 initial stage, after the 5 th iteration step, and after the 13 th iteration step when the iteration converges. Initially, the fractures from all realizations are uniformly distributed over the ow domain in the NE-SW directions as indicated by the prior information, as shown in Fig. 5.4a. After ve iteration steps, the probability map in Fig. 5.4b shows higher probability of fractures in the area between P 3 and I, where the fracture in the reference facies map is located, than the rest of the eld. When the iteration converged, the probability map in Fig. 5.4c exhibits a pattern of fracture distribution very similar to the reference facies map, with very high probability values (> 90%) at the grid blocks between P 3 and I. In addition to the orientation and position of the fracture, the updated probability map provides accurate estimation of its length, and the fractures are not connected to either the producer or the injector. It can be noticed that there are other possible fractures with relatively low probability values on the probability map, such as those near P 8. The production history and seismic data is not very sensitive to the existence of those fractures that are perpendicular to the ow path between the producer and injector. With additional observed data, it will be possible to rule out those articial fractures from the updated probability map. The initial facies maps of three randomly selected realizations are shown in the Figure 5.4: Probability maps of fractures for Case 1: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 13 th iteration step when the iteration converges. 93 rst row in Fig. 5.5, and the updated facies maps are shown in the second and third rows in Fig. 5.5. Although the initial facies maps show randomly distributed fractures, all three updated facies maps contain a fracture similar to the one in the reference facies map. The short fracture in the lower right corner only appears in one of the realizations, indicating relatively high uncertainty due to the low sensitivity of these fractures to the observed data. The mean Hough function eld at the initial stage, after the 5 th iteration step, and after the 13 th iteration step when the iteration converges, is plotted in Fig. 5.6. These mean elds are computed simply by taking the arithmetic averages of Hough functions over all the realizations in the ensemble. At the initial stage, the mean Hough function eld is nearly uniform at0:5, owing to the stationarity of the spatially correlated Gaussian random elds. According to the reference facies map, the fracture in Fig. 5.2 is represented by the parameter vector at point p (=4; 0:5) in the Hough space. After the 5 t h iteration step, the mean Hough function in Fig. 5.6b shows high values near pointp . The updated mean Hough function in Fig. 5.6c shows a obvious peak exactly at pointp , indicating the high probability of a fracture represented by the parameter vector at point p . There are no other Hough function peaks in Fig. 5.6c, thanks to the covariance localization in the history matching method. The updated parameters were used to simulate the reservoir model for another 30 days with the same condition to investigate the predictability of the updated models. In Fig. 5.7, the model predictions for the production history before and after the updating are provided side by side to demonstrate the ability of the methodology to improve the model's predictability. The black colored sections of the lines are the period (rst 20 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 30 days with the history matched parameters. The blue asterisks are observations. The gures demonstrate the randomness of the model predictions from the initial ensemble and the strong nonlinearity of the model in which the shapes of the oil rate curves are completely dierent whether or not the well is connected by a fracture. After updating, the match between the model predictions and observations is signicantly improved. The 94 Figure 5.5: Facies maps of three randomly selected realizations for Case 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 95 majority of the oil rate, water cut, and injector bottomhole pressure curves from the realizations match the observation curves, and the variations of the predictions are very small compared to those before updating. Not only the production data in the rst 20 days match the observations when they are used in the history matching, the update model is also capable to provide accurate predictions for the future period when the observations are not assimilated. The observations of the mean seismic impedance dierence data at the end of the 20 th day is shown in Fig. 5.8a, and the mean impedance dierence eld at the initial stage and after the 13 th iteration step when the iteration converges, is plotted in Fig. 5.8b and Fig. 5.8c, respectively. These mean elds are computed simply by taking the arithmetic averages of impedance dierence over all the realizations in the ensemble, while the impedance dierence is the arithmetic dierence between the reference seismic survey at the initial state and the monitor survey at the end of the 20 th day. At the initial stage, the mean impedance dierence eld is well correlated to the water saturation distribution in the reservoir, with high water saturation near the injector and low water saturation near the boundaries. There is no indications to the existence of fractures in Fig. 5.8b. The updated mean impedance dierence in Fig. 5.8c matches the observations as shown in Fig. 5.8a, where grid blocks with high impedance dierence values between P 3 and the injector indicate the fracture, which provides the preferential ow path between the two wells and increases the Figure 5.6: Mean Hough function maps for Case 1: (a) at the initial stage, (b) after the 5 t h iteration step, and (c) after the 13 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 96 Figure 5.7: The predictions from all the ensemble members before and after the history matching for Case 1: oil rate at (a) P 2 and (c) P 3; wa- ter cut at (b) P 2 and (d) P 3; (e) the bottom hole pressure at the injec- tor. The black colored sections of the lines are the period (rst 20 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 30 days with the his- tory matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. 97 water saturation in those grid blocks. Fig. 5.9 gives the seismic impedance elds from the monitor survey at the end of the 20 th day. Compared with the impedance dierence maps in Fig. 5.8, the seismic monitor surveys are poorly correlated with the fracture distribution. With high porosity and high water saturation in the grid blocks in the fracture facies, the impedance is similar to the matrix blocks with low porosity and low water saturation, Figure 5.8: Mean impedance dierence maps at the end of the 20 th day for Case 1: (a) the observations, (b) at the initial stage, and (c) after the 13 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. Figure 5.9: Mean impedance maps at the end of the 20 th day for Case 1: (a) the observations, (b) at the initial stage, and (c) after the 13 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 98 as discussed in Section 2.5. Therefore, the impedance dierence data is preferred for characterizing the fractures. 5.4 Case 2: Assimilation of Time-Lapse Seismic Data for Reference Field with Two Fractures In this case, the reference facies map contains two connected fractures that are not connected to any wells, as illustrated in Fig. 5.10. Prior information is introduced into the history matching process: the fractures in this ow domain are known to be in the NE-SW and NW-SE directions, thus the active grid blocks in the Hough space is limited to those that satisfy 2f[0:9;0:6][ [0:6; 0:9]g. As a result, there are 997 active grid blocks in the Hough space. Moreover, the fractures in the NE-SW direction are assumed to be known in the northeast quadrant of the reservoir, to improve the eciency of the history matching. Figure 5.10: Case 2: The reference facies map: black grid blocks are fracture facies, and the rest are rock matrix. 99 The initial ensemble is generated in the same way as case 1. Synthetic observed data are obtained by simulating the reservoir model with the reference facies map in Fig. 5.10a for 50 days, with seismic impedance data available before the production and at 20 th day. Only the observations in the rst twenty days, as well as the seismic dierence data, is assimilated at a single step in the LM-EnRML. The iteration converged at the 15 t h step, and the reduction in the mean of data mismatch is shown in Fig. 5.11. The data mismatch is reduced steadily from 5:7610 8 to 5:0510 7 . The convergence is considerably improved by the assimilation of seismic data, compared with Case 4 in Chapter 4. The probability maps computed by Eq. 4.4 and the mean Hough function elds at the initial stage and after history matching are plotted in Fig. 5.12. Initially, the majority of the fractures in NE-SW direction are found in the northeast quadrant of the reservoir in Fig. 5.12a, as stated in the prior information used to generate the initial ensemble. The probability map after the 5 th iteration step (Fig. 5.12b) already shows large probabilities of the two fractures in the reference facies map, though the NE-SW fracture shows larger uncertainty than the other one. There are still artifacts with low probabilities in the southwest quadrant. The updated probability map clearly shows the two fractures similar to those in the reference facies map, and the artifacts are eliminated by the iteration and the localization. The estimated fractures Figure 5.11: Case 2: Reduction in the mean of the data mismatch with iteration. 100 are slightly dierent in the angle from the reference, due to the discretization of the Hough space and the ow domain. Compared with Case 4 in Chapter 4, the estimated fracture distribution is much more accurate. In Fig. 5.13a, the mean Hough function eld is uniform at the initial stage. The updated mean Hough function elds (Figs. 5.13b and 5.13c) show obvious peaks at p 1 (=4; 0:23) and p 2 (=4; 0:49), that represent the two fractures in the reference facies map (Fig. 5.10). The facies maps of three randomly selected realizations before and after updating Figure 5.12: Probability maps of fractures for Case 2: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 15 th iteration step when the iteration converges. Figure 5.13: Mean Hough function maps for Case 2: (a) at the initial stage; (b) after the 5 t h iteration step; and (c) after the 15 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 101 in Fig. 5.14 demonstrate that the history matching with the proposed parameter- ization can successfully identify the fractures conditional on the production history and seismic data, although the initial fractures are randomly distributed in the facies maps and the fractures are not connected to wells. The model predictions with the updated parameters are shown in Fig. 5.15. The updated oil rate and water cut at the P 2 (Figs. 5.15a and 5.15b) and P 3 (Figs. 5.15c and 5.15d), as well as the bottomhole pressure at the injector(Figs. 5.15e), are matched with the observed data available in the rst 20 days after the history matching. The ensemble is also proven to be able to predict the reservoir behavior for another 30 days under current conditions without the assimilation of additional observed data. The observations of the mean seismic impedance dierence data at the end of the 20 th day is shown in Fig. 5.16a, and the mean impedance dierence eld at the initial stage and after the 15 th iteration step when the iteration converges, is plotted in Fig. 5.16b and Fig. 5.16c, respectively. At the initial stage, the mean impedance dierence eld is well correlated to the water saturation distribution in the reservoir, with high water saturation near the injector and low water saturation near the boundaries. There is no indications to the existence of fractures in Fig. 5.8b as a result of the randomly distributed fractures in the initial ensemble. The updated mean impedance dierence in Fig. 5.16c matches the observations as shown in Fig. 5.16a, where grid blocks with high impedance dierence values in the northeast quadrant of the reservoir indicate the two connected fractures, which provide the preferential ow path between the two wells and increase the water saturation in those grid blocks. 5.5 Summary Time-lapse seismic data is assimilated with production data to improve the char- acterization of fractures that are not connected to wells. With the help of the distance-dependent localization designed for the Hough-transform-based parameteri- zation, large amount of seismic data is eciently assimilated. The seismic impedance dierence data between the reference survey and monitor survey is used because of 102 Figure 5.14: Facies maps of three randomly selected realizations for Case 2: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps when the iteration converges. 103 Figure 5.15: The predictions from all the ensemble members before and after the history matching for Case 2: oil rate at (a) P 2 and (c) P 3, as well as water cut at (b) P 2 and (d) P 3; (e) the bottom hole pressure at the injector. The black colored sections of the lines are the period (rst 20 days) within which the model is updated with observed data, while the green colored sections are model predictions for another 30 days with the history matched parameters. The blue asterisks represent observations. The dashed vertical line denotes the end of the history matching period. 104 its good correlation with the change of uid saturations. The methodology is tested in two dimensional synthetic water ooding examples to evaluate its eectiveness to estimate the fracture distributions. The impacts of the complexity of the reference fracture distribution, the sensitivity of the parameters to the observations, and the prior information, are examined in terms of the quality of the estimated fracture distribution and the predictability of the updated models. The examples indicate that the assimilation of seismic data enables the proposed method to characterize fractures that are not connected to wells. Compared with Case 4 in Chapter 4, the accuracy of the estimated fracture distribution is considerably improved, and the associated uncertainty is signicantly reduced. The sensitivity between the parameters and the production history depends on the distance between the fractures and the wells. But the sensitivity between the fractures away from the wells and the seismic impedance dierence data is relatively strong, as the seismic dierence data provides valuable observations about the uid migration with good areal coverage. Prior information about the fracture distribution plays an important role in the estimation of fracture distributions. It helps to improve the eciency of the history matching and to reduce the parameters to be updated. The statistics of regional Figure 5.16: Mean impedance dierence maps at the end of the 20 th day for Case 2: (a) the observations, (b) at the initial stage, and (c) after the 15 t h iteration step when the iteration converges. The three maps are plotted with the same color scale. 105 fracture orientations acts as an additional constraint to the estimation of the fracture distribution, which is important in reducing the uncertainty associated with the non- unique interpretations of the observations. The predictability of the model improves after being updated with production history and seismic data. The updated model outputs match both the production history and seismic surveys assimilated during the history matching. With the help of the accurate characterization of the fractures by the proposed parameterization, the calibrated models can provide good reservoir responses and improved predictions for the periods beyond the available production history. 106 Chapter 6 Characterization of Structural Control in Geothermal Reservoirs In this chapter, more realistic case studies are developed to investigate the eec- tiveness of the proposed Hough-transform-based parameterization in characterizing complex geology. A good application is in the geothermal reservoirs, where the suc- cess of the thermal production relies on the fault or fracture network to provide the ow capacity for the hot uids. The rst case is designed to demonstrate the eectiveness of this history matching technique to estimate the position of a single fault in a geothermal reservoir based on a tracer test at the Soultz-sous-For^ ets enhanced geothermal system (EGS) site. The other two cases are more complicated with a fault network from the depth of 5000ft to 17000ft at Soultz-sous-For^ ets. Although the measured data about the dip, strike, and depth of the fracture zones intercepted by wells are available, the extensions and connectivity of the faults are uncertain. The production data and tracer test data are used to estimate the hydraulic connections between the wells, and the extensions of the faults. 107 6.1 Summary of Soultz-sous-For^ ets EGS site The Soultz geothermal area is located in the Upper Rhine valley, about 70 km north of Strasbourg, France (G erard et al., 2006). The crystalline basement is overlain by a 1.5 km layer of sediments, and consists of Paleozoic granites with fractures ranging from micro-cracks to large normal faults. As an EGS reservoir, the rock matrix is used as a heat exchanger and the ow capacity is provided by the natural or hydraulically induced fractures. The stimulation of the wells can also result in opening of preexisting fractures (Vogt et al., 2012). The temperature of 200 C was encountered at the depth of 5000 m. Four wells (GPK1, GPK2, GPK3, and GPK4) had been drilled in the area since 1987 (G erard et al., 2006). Wells GPK2, GPK3, and GPK4 were directionally drilled to be 600-700 m apart at the target depth, in order to avoid the hydraulic short-circuiting between the wells. 6.2 Case 1: History Matching for a Tracer Test at Soultz-sous-For^ ets This case study is based on a tracer test performed among the wells GPK2, GPK3, and GPK4 (Fig. 6.1) at the depth of about 5000 m (Sanjuan et al., 2006). The wells GPK2 and GPK3 intersect a fault (80 in dip and striking to 255 N), while GPK4 misses the fault. Hence, the continuation of the fault toward GPK4 is unknown. The reservoir model is considered as a two-dimensional single phase geothermal model with three wells. The ow domain is a rectangle of size 6000ft 9000ft 1500ft, centered at the depth of 5000 m. The initial temperature of the entire reservoir is 392 F (200 C), and the reservoir pressure is at the hydraulic static pressure. The reservoir is discretized into 4060 regular grids, as shown in Fig. 6.2. The white grid blocks in Fig. 6.2 indicate rock matrix of the granite, the gray grid blocks represent the hydraulically stimulated zone around the wellbores, and the fault is indicated by the black grid blocks. The eective permeabilities in the individual facies are given in Table 6.1., as well as the other hydraulic and thermal properties of the reservoir. In this reference facies map, GPK2 and GPK3 are connected by the fault, but GPK4 108 Parameter Value Eective permeability of rock matrix (md) 0.001 Eective permeability of fracture (md) 100 Eective permeability of stimulated zone (md) 10 Eective porosity 0.0005 Thermal conductivity (Btu D 1 ft 1 F 1 ) 34.67 Heat capacity (Btu lb 1 F 1 ) 492 Diusion coecient (ft 2 D 1 ) 465 Table 6.1: Parameters used in the geothermal simulations is isolated. The ow equations are solved with Eclipse (Schlumberger, 2007). During the tracer test, GPK2 produces 392 F hot water at 6470 B/D, and GPK4 Figure 6.1: Deep borehole conguration at Soultz EGS site, based on the sketch from Gessner et al. (2009). The fault is shown as the blue ellipse, and the trajectories of the wells are shown as black curves. 109 produces at 1680 B/D, while the produced water (with tracer in the solution) is reinjected by GPK3 at 158 F . The uorescein tracer is injected at GPK3 for the rst day with a concentration of 0.0618 Mol/B, and the circulation is carried on for another 150 days without adding more tracer. Adsorption and reaction of the tracer to the reservoir rock is assumed to be trivial. The water rate, bottomhole pressure, and produced tracer concentration are available at GPK2 and GPK4 at the end of each day. Water injection rate and bottomhole pressure is available at the injector GPK3. The synthetic observations of production data and tracer concentration are generated by simulating the reservoir model with the reference facie map in Fig. 6.2. In addition to the nonlinearity of the models in fractured reservoirs, the assimila- tion of tracer test data induces extra nonlinearity between the fracture distribution and tracer concentration. Therefore, iterative history matching method is necessary to estimate the fracture distribution and to match the observations. The iterative en- semble smoother LM-EnRML (Chen and Oliver, 2013, 2014) with distance-dependent Figure 6.2: The ow domain and well locations of Case 1. Color of the grid blocks indicates the facies distribution in the reservoir. 110 covariance localization is used to assimilate both production data and tracer concen- tration data from the entire duration of the tracer test. Prior information is introduced into the history matching process: (1) assuming that GPK3 is known to be connected by a fault with the orientation of ==2; and (2) the fractures in this ow domain are known to be in the N-S directions, thus the active grid blocks in the Hough space is limited to those that satisfy 2 [0:2; 0:2]. After 9 steps of iteration, the LM-EnRML converged. Fig. 6.3 shows the reduc- tion in the mean of data mismatch. The iteration converged relatively fast, due to the incorporation of prior information. The probability maps of fractures are plotted in Fig. 6.4 at the initial stage, after the 5 th iteration step, and after the 9 th iteration step when the iteration converges. Initially, the fractures from all realizations are uniformly distributed over the ow domain, except the fault at GPK3 as indicated by the prior information, as shown in Fig. 6.4a. The position of the prior fault is largely uncertain as randomly generated in the initial ensemble. After ve iteration steps, the probability map in Fig. 6.4b shows higher probability of fractures that connect GPK2 and GPK3, where the fault in the reference facies map is located, than the rest of the eld. When the iteration converged, the probability map in Fig. 6.4c exhibits a pattern of fracture distribution very similar to the reference facies map. In addition to the orientation and position Figure 6.3: Case 1: Reduction in the mean of the data mismatch with iteration. 111 of the fault, the updated probability map provides accurate estimation of its length, and the fault is not connected to GPK4. The initial facies maps of three randomly selected realizations are shown in the rst row in Fig. 6.5, and the updated facies maps are shown in the second and third rows in Fig. 6.5. All three updated facies maps contain a fault similar to the one in the reference facies map. The randomly distributed fractures in the initial facies maps disappear after the updating. The updated parameters were used to simulate the reservoir model with the same condition to obtain the predictions from the updated models. In Fig. 6.6, the model predictions for the tracer concentration history before and after the updating are provided side by side for both the producers. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. At GPK2, which is connected to the injector by the fault, the updated tracer con- centration predictions match the observations, and the uncertainty of the predictions is signicantly reduced. On the other hand, although the match between the model predictions and observations is improved at GPK4, which is not connected by the fault, the uncertainty level of the updated predictions is larger than that at GPK2. The model predictions for water rate and bottomhole pressure at GPK2, GPK4, Figure 6.4: Probability maps of fractures for Case 1: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 9 th iteration step when the iteration converges. 112 Figure 6.5: Facies maps of three randomly selected realizations for Case 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 113 Figure 6.6: The predictions of tracer concentration from all the ensemble members before and after the history matching for Case 1: tracer concen- tration at (a)(b) GPK2 and at (c)(d) GPK4. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 114 and GPK3 before and after updating are plotted in Fig. 6.7. The predictions from the updated model match the observed data. The water rates at the wells remains constant at the rate target before and after updating. The updated predictions of bottomhole pressure match the observations, and the variance of the predictions from the ensemble decreases signicantly. 6.3 Case 2: Characterization of Large-Scale Fault Network at Soultz-sous-For^ ets In this case study, the reservoir model is based on a large scale structural study of the granite basement at Soultz-sous-For^ ets (Place et al., 2011). The study covered a vast area from almost the top of the basement to the target depth of the wells. A major fault (GPK3-FZ4770) was identied with a vertical extension of 4km, intersecting the smaller faults at various depth. The conceptual connection plan by Place et al. (2011) is developed in this case study to construct a inclined two-dimensional reservoir, in order to investigate the hydraulic connections among the faults. Although the wells intersect multiple faults at dierent depths, the intersects are treated as dierent wells in this case. 6.3.1 Model Description The reservoir model is considered as a inclined two-dimensional single phase geother- mal model with ve wells. The ow domain is a rectangle of size 12,000ft 18,000ft 1000ft. The top of the reservoir is at the depth of 4600ft in the north end, and at the depth of 16,400ft in the south end. The initial temperature of the entire reservoir is 392 F (200 C), and the reservoir pressure is at the hydraulic static pressure. The reservoir is discretized into 4060 regular grids, as shown in Fig. 6.8. The white grid blocks in Fig. 6.8 indicate rock matrix of the granite, the gray grid blocks represent the hydraulically stimulated zone around the wellbores, and the fault is indicated by the black grid blocks. The eective permeabilities in the individual facies are given in Table 6.1., as well as the other hydraulic and thermal properties of the reservoir. In 115 Figure 6.7: The predictions from all the ensemble members before and after the history matching for Case 1: water rate at (a) GPK2, (c) GPK4, and (e) GPK3, as well as bottomhole pressure at (b) GPK2, (d) GPK4, and (f) GPK3. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 116 this reference facies map, I1 and P 4 are connected by the major fault F 4. P 1, P 2, andP 3 are connected by three smaller faultsF 1,F 2, andF 3, respectively. There are two other faultsF 5 andF 6, that are not connected to wells, butF 5 is communicating with the major fault F 4. The ow equations for this geothermal reservoir are solved with Eclipse (Schlumberger, 2007). During the tracer test, P 1 produces 392 F hot water at 1000 B/D, P 2 produces at 1900 B/D, P 3 produces at 2250 B/D, and P 4 produces at 3000 B/D, while the produced water (with tracer in the solution) is reinjected by I1 at 158 F . The uo- rescein tracer is injected atI1 for the rst day with a concentration of 0.0618 Mol/B, and the circulation is carried on for another 150 days without adding more tracer. Adsorption and reaction of the tracer to the reservoir rock is assumed to be trivial. The water rate, bottomhole pressure, and produced tracer concentration are available at the producers at the end of each day. Water injection rate and bottomhole pressure is available at the injector. The synthetic observations of production data and tracer Figure 6.8: The ow domain and well locations of Case 2. Color of the grid blocks indicates the facies distribution in the reservoir. Names fo the wells are marked as red, and names of the faults are marked as blue. 117 concentration are generated by simulating the reservoir model with the reference fa- cies map in Fig. 6.8. The iterative ensemble smoother LM-EnRML (Chen and Oliver, 2013, 2014) with distance-dependent covariance localization is used to assimilate both production data and tracer concentration data from the entire duration of the tracer test. 6.3.2 Scenario 1: History Matching with Prior Information at Producing Wells In this scenario, only the prior information at the producing wells are applied in the history matching process:: (1) the major faultF 4 is added as a deterministic fracture to all ensemble members throughout the history matching; (2) producerP 1 is known to be connected by a fault with the orientation of = 0:436; (3) producer P 2 is known to be connected by a fault with the orientation of = 0:942; (4) producer P 3 is known to be connected by a fault with the orientation of = 1:030;and (5) the fractures in this ow domain are known to be generally in the NNE and NEE directions, thus the active grid blocks in the Hough space is limited to those that satisfy 2 [0; 1:57]. After 5 steps of iteration, the LM-EnRML converged. Fig. 6.9 shows the reduc- tion in the mean of data mismatch. The iteration converged relatively fast, due to the incorporation of prior information. In Fig. 6.10, The probability maps of fractures are plotted at the initial stage, af- ter the 1 st iteration step, and after the 5 th iteration step when the iteration converges. Initially, the fractures from all realizations are uniformly distributed over the ow domain, except the faultsF 1,F 2,F 3, andF 4, as indicated by the prior information (Fig. 6.10a). The position of the prior faults F 1, F 2, and F 3 are largely uncertain as randomly generated in the initial ensemble. When the iteration converged, the probability map in Fig. 6.10c exhibits a pattern of fracture distribution very similar to the reference facies map in Fig. 6.8, where F 1, F 2, and F 3 are connected to F 4. The uncertainty of the position and length of the faults are reduced in the updated 118 probability map. However, the faults F 5 and F 6, that are not directly connected to wells, are not identied by the history matching. These two faults are not in the main ow paths among the wells, and the observations of production data and tracer concentration at the wells are not very sensitive to the existence of the faults. The initial facies maps of three randomly selected realizations are shown in the rst row in Fig. 6.11, and the updated facies maps are shown in the second and third rows in Fig. 6.11. All three updated facies maps contain the faults F 1, F 2, Figure 6.9: Scenario 1: Reduction in the mean of the data mismatch with iteration. Figure 6.10: Probability maps of fractures for Scenario 1: (a) at the initial stage; (b) after the 1 st iteration step; and (c) after the 5 th iteration step when the iteration converges. 119 F 3, and F 4, similar to those in the reference facies map. The randomly distributed fractures in the initial facies maps disappear after the updating. The updated parameters were used to simulate the reservoir model with the same condition to obtain the predictions from the updated models. In Fig. 6.12, the model predictions for the tracer concentration history before and after the updating are provided side by side for producersP 1 andP 2. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. The updated tracer concentration predictions match the observations, and the uncertainty of the predictions is signicantly reduced. Since P 1 and P 2 are relatively far from the injector, there are still certain degree of uncertainty associated with the updated predictions. On the other hand, in Fig. 6.13, where the model predictions for the tracer concentration history before and after the updating are provided side by side for producersP 3 andP 4, the model predictions for the tracer concentration are considerable improved after updating, with much smaller uncertainties than that at P 1 and P 2, since P 3 and P 4 are closer to the injector. The model predictions for water rate and bottomhole pressure at P 1,P 4, andI1 before and after updating are plotted in Fig. 6.14. The predictions from the updated model match the observed data. The water rates at the wells show great variations before updating, because some of the faults in the initial ensemble are connected to the major fault F 4, while others are not. The updated water rates match the observations almost perfectly. The updated predictions of bottomhole pressure match the observations, and the variance of the predictions from the ensemble decreases signicantly. 6.3.3 Scenario 2: History Matching with Prior Information for All Faults In this scenario, the prior information for all the faults in Fig. 6.8 are applied in the history matching process:: (1) the major faultF 4 is added as a deterministic fracture to all ensemble members throughout the history matching; (2) producerP 1 is known to be connected by a fault with the orientation of = 0:436; (3) producer P 2 is 120 Figure 6.11: Facies maps of three randomly selected realizations for Sce- nario 1: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 1 st iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 121 Figure 6.12: The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 1: tracer concentration at (a)(b) P 1 and at (c)(d) P 2. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 122 Figure 6.13: The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 1: tracer concentration at (a)(b) P 3 and at (c)(d) P 4. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 123 Figure 6.14: The predictions from all the ensemble members before and after the history matching for Scenario 1: water rate at (a) P 1, (c)P 4, and (e) I1, as well as bottomhole pressure at (b) P 1, (d) P 4, and (f) I1. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 124 known to be connected by a fault with the orientation of = 0:942; (4) producer P 3 is known to be connected by a fault with the orientation of = 1:030; (5) a fault at = 1:047 and = 0:225 with uncertain length and axial position, is known to exist to represent F 5; (6) a fault at = 0:471 and = 0:588 with uncertain length and axial position, is known to exist to represent F 6; and (7) the fractures in this ow domain are known to be generally in the NNE and NEE directions, thus the active grid blocks in the Hough space is limited to those that satisfy 2 [0; 1:57]. After 12 steps of iteration, the LM-EnRML converged. Fig. 6.9 shows the reduction in the mean of data mismatch. The iteration converged relatively fast, due to the incorporation of prior information. The probability maps of fractures are computed and plotted in Fig. 6.10. Ini- tially, the fractures from all realizations are uniformly distributed over the ow do- main, except the faults introduced by the prior information (Fig. 6.10a). The position of the prior faults F 1, F 2, F 3, F 4, F 5, and F 6 are largely uncertain as randomly generated in the initial ensemble. When the iteration converged, the probability map in Fig. 6.10c exhibits a pattern of fracture distribution very similar to the reference facies map in Fig. 6.8, where F 1, F 2, and F 3 are connected to F 4. The uncertainty of the position and length of the faults connected by wells are reduced in the updated probability map. However, the position and length of the faults F 5 andF 6, that are Figure 6.15: Scenario 2: Reduction in the mean of the data mismatch with iteration. 125 not directly connected to wells, are still uncertain after updating. These two faults are not in the main ow paths among the wells, and the observations of production data and tracer concentration at the wells are not very sensitive to the existence of the faults. With the help of the prior information, which provides extra constrains to these faults, the faultF 6 is characterized to be disconnected fromF 4, showing an improvement to the initial ensemble. The initial facies maps of three randomly selected realizations are shown in the rst row in Fig. 6.17, and the updated facies maps are shown in the second and third rows in Fig. 6.17. All three updated facies maps contain the faults F 1, F 2, F 3, F 4, F 5, and F 6, similar to those in the reference facies map. However, the length and axial position ofF 5 is still uncertain after updating, as shown in the three realizations. The randomly distributed fractures in the initial facies maps disappear after the updating. The updated parameters were used to simulate the reservoir model with the same condition to obtain the predictions from the updated models. In Fig. 6.18, the model predictions for the tracer concentration history before and after the updating are provided side by side for producersP 1 andP 2. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. The Figure 6.16: Probability maps of fractures for Scenario 2: (a) at the initial stage; (b) after the 5 th iteration step; and (c) after the 12 th iteration step when the iteration converges. 126 Figure 6.17: Facies maps of three randomly selected realizations for Sce- nario 2: (a), (b), and (c) the facies maps of three realizations at the initial stage; (d), (e), and (f) the corresponding updated facies maps after the 5 th iteration step; (g), (h), and (i) the corresponding updated facies maps when the iteration converges. 127 updated tracer concentration predictions match the observations, and the uncertainty of the predictions is signicantly reduced. Since P 1 and P 2 are relatively far from the injector, there are still certain degree of uncertainty associated with the updated predictions. On the other hand, in Fig. 6.19, where the model predictions for the tracer concentration history before and after the updating are provided side by side for producersP 3 andP 4, the model predictions for the tracer concentration are considerable improved after updating, with much smaller uncertainties than that at P 1 and P 2, since P 3 and P 4 are closer to the injector. The model predictions for water rate and bottomhole pressure at P 1,P 4, andI1 before and after updating are plotted in Fig. 6.20. The predictions from the updated model match the observed data. The water rates at the wells show great variations before updating, because some of the faults in the initial ensemble are connected to the major fault F 4, while others are not. The updated water rates match the observations almost perfectly. The updated predictions of bottomhole pressure match the observations, and the variance of the predictions from the ensemble decreases signicantly. 6.4 Summary The Hough-transform-based parameterization is used in the history matching for tracer tests in geothermal reservoirs based on the Soultz-sous-For^ ets EGS project, where the uid ow between the injector and producers are controlled by fault net- works with uncertain extension. Production history and tracer concentration data is assimilated by the LM-EnRML with distance-dependent covariance localization. The prior information helps to constrain the positions of the faults and to improve the convergence rate of the iterative smoother. 128 Figure 6.18: The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 2: tracer concentration at (a)(b) P 1 and at (c)(d) P 2. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 129 Figure 6.19: The predictions of tracer concentration from all the ensemble members before and after the history matching for Scenario 2: tracer concentration at (a)(b) P 3 and at (c)(d) P 4. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 130 Figure 6.20: The predictions from all the ensemble members before and after the history matching for Scenario 2: water rate at (a) P 1, (c)P 4, and (e) I1, as well as bottomhole pressure at (b) P 1, (d) P 4, and (f) I1. The gray colored lines are the model predictions from the ensemble members, and the blue lines are observations. 131 Chapter 7 Conclusions and Future Work 7.1 Conclusions It is a challenge to estimate fracture distribution in fractured reservoirs by conven- tional history matching methods because of the complex geology and strong hetero- geneity in such reservoirs. Without proper constraints to the geometry and continuity of the fractures, it is dicult to obtain a geologically realistic solution even when the model predictions are matched with the observations, since the history matching techniques can give non-unique solutions for the ill-posed problems. In this study, we propose a parameterization based on the Hough transform, allow- ing the assisted history matching methods to directly update the fracture distributions conditional on observed data. The fractures in binary facies maps are characterized by the angle-radius representation of straight lines, the same technique used in the Hough transform. Each fracture is dened by the angle of its normal, its algebraic distance from the origin in the facies map, and its length and axial displacement dened on a continuous Hough space. The existence of the fractures in a facies map is determined by a Hough function also dened on the Hough space. Hence the bi- nary facies maps can be replaced by new parameter vectors that can be Gaussian distributed on the Hough space. Instead of the discontinuous facies maps, the new parameter vectors are updated by the assisted history matching methods. Whenever a simulation run is needed during the updating process, a new set of facies maps is 132 generated from the parameter vectors by the inverse transform. When coupled with history matching methods, such as EnKF and ES, this parameterization satises the Gaussian assumption of these methods, and the facies maps used in model simulation or prediction are always geologically realistic with the help of the inverse transform. Moreover, the proposed method can work seamlessly with the commercial pixel-based reservoir simulators, such as the Eclipse, without the needs to modify the simulators. The methodology is tested in two dimensional synthetic water ooding examples to evaluate its eectiveness to estimate the fracture distributions with production data. The impacts of the complexity of the reference fracture distribution, the sensitivity of the parameters to the observations, and the prior information, are examined in terms of the quality of the estimated fracture distribution and the predictability of the updated models. The examples indicate that the proposed method is an ecient and exible way to estimate fracture distribution in fractured reservoirs. The estima- tion of the fracture distribution, including the number of fractures, the orientation, position and length of the fractures, improves signicantly when production history at the wells is assimilated. The Gaussian distributed Hough space parameters ensure good eciency in the history matching using the EnKF or ES. The estimated fracture distributions can be expressed as probability maps, which provide probabilities of the estimated fractures and are valuable for future reservoir management and optimiza- tion. The sensitivity between the parameters and the production history depends on the distance between the fractures and the wells. High sensitivity is observed when the fractures pass through wells, while the sensitivity is lower when the fractures are away from the wells, rendering larger uncertainties in the estimated probability maps. Additional observations other than the production history at the wells may be essen- tial to accurately estimate the fractures in the area away from the wells. Assimilation of the additional observations is being developed by our ongoing research. Prior information about the fracture distribution plays an important role in the estimation of fracture distributions, especially when the observed data are inade- quate for an accurate estimation. Two kinds of prior information are assimilated in this study: (1) known fractures that pass through wells; and (2) statistics of regional 133 fracture orientations in the ow domain. Both prior information considerably im- prove the updated estimation of the fractures, as well as the eciency of the history matching. The second type of prior information also signicantly reduces the number of parameters to be updated. The predictability of the model improves after being updated with production history. The updated model outputs match the production history assimilated during the history matching. With the help of the accurate characterization of the fractures by the proposed parameterization, the calibrated models can provide good reservoir responses and improved predictions for the periods beyond the available production history. Furthermore, compared with the applications of `big loop' history matching in fractured reservoirs using DFN models, the proposed method characterizes and updates the fracture distribution directly by production history, thus prevents the costly return to DFN models where fracture objects are typically generated. The assimilation of time-lapse seismic data with the Hough-transform-based pa- rameterization is investigated with a few two dimensional synthetic water ooding examples. The impacts of the complexity of the reference fracture distribution, the sensitivity of the parameters to the observations, and the prior information, are exam- ined in terms of the quality of the estimated fracture distribution and the predictabil- ity of the updated models. The examples indicate that the assimilation of seismic data enables the proposed method to characterize fractures that are away from wells. The characterization of those fractures and the reduction of the associated uncer- tainty benet from the time-lapse seismic data, which provides valuable information about the uid migration in the reservoir, especially in the vast area away from wells. The method is also tested in some realistic examples based on the Soultz-sous- For^ ets EGS project, where the uid ow between the injector and producers are con- trolled by fault networks with uncertain extension. The iterative ensemble smoother LM-EnRML with Hough-transform-based parameterization is coupled with a geother- mal reservoir model to assimilate tracer test data and production data. The fault network is successfully characterized, though the posterior uncertainty of some faults are relatively high because they are away from the main ow path. 134 7.2 Future Work In the future, this research work can be furthered along the following directions: 1) The work might be extended to three dimensional reservoir models. In order to characterize three dimensional fractures, the spherical coordinates can be used to dene the orientation of fractures, resulting in a three dimensional Hough space. Hence, a fracture is represented by a point in the Hough space with three coordinates: strike of its normal, dip, and distance from the origin. Assumed to be ellipsoid, the size of a fractures can be characterized by the length along the long and short axes. The longitudinal and lateral displacement, as well as the axial rotation of the fracture on the fracture plane, can be determined by three additional parameters. Therefore, a three-dimensional fracture will be represented by a point in the three dimensional Hough space with six parameters (including the Hough function). 2) In this study, the properties (permeability, porosity, etc.) in the facies models are constant in each facies type. In the future work, the methodology can be developed to update the homogeneous permeability and porosity in each facies, together withe the facies maps. And furthermore, it is of interest to investigate the possibility of updating spatial-variant properties within each facies. 3) The iterative ensemble smoother is necessary for the characterization of fracture distributions and the match of observation data. The process can be time-consuming if applied to complex fracture networks where the convergence of the iteration can be slow. It is worth to investigate the possibilities of using polynomial expansion methods or surrogates to approximate the model predictions, since the Hough function eld is Gaussian distributed and generated by Karhunen-Lo eve expansion. 135 Nomenclature A accumulator C autocovariance function C D covariance matrix of observations C M covariance matrix of model parameters C SC scaling matrix for model parameters D axial displacement of a fracture d vector of noisy observations f eigenfunction F () reservoir model G sensitivity matrix g vector of dynamic measurements G shear modulus, m=Lt 2 , Pa H operator matrix H Hough function of a fracture I identity matrix I binary image 136 K Kalman gain matrix K bulk modulus, m=Lt 2 , Pa K c modulus of clay, m=Lt 2 , Pa K fluid modulus of uid-in-place, m=Lt 2 , Pa K frame modulus of rock frame, m=Lt 2 , Pa K g modulus of gas, m=Lt 2 , Pa K grain modulus of rock grain, m=Lt 2 , Pa K o modulus of oil, m=Lt 2 , Pa K s modulus of sand, m=Lt 2 , Pa K w modulus of water, m=Lt 2 , Pa L length of a fracture L d diagonal of ow domain, L, ft L x length of ow domain along x direction, L, ft L y length of ow domain along y direction, L, ft m vector of model parameters M function to nd local maxima N 1 number of grid blocks in ow domain along x direction N 2 number of grid blocks in ow domain along y direction N d number of observations N e number of ensemble members N g number of dynamic measurements 137 N m number of model parameters N t number of timesteps N u number of dynamic variables N y number of elements in state vector p probability map of fracture distributions P covariance matrix of state vectors R covariance matrix of observation errors S d mean of data mismatch S g gas saturation S o oil saturation S w water saturation u vector of dynamic variables V p P-wave velocity, L=t, m=s V s S-wave velocity, L=t, m=s y state vector Y ensemble of state vectors Z localization function in matrix form Z seismic impedance, m=L 2 t, Pas=m practical range shaliness Dirac delta function 138 standard Gaussian random variable localization function the angle of a straight line's normal eigenvalue ! Levenberg-Marquardt regularization parameter the algebraic distance from a straight line to the origin b bulk density, m=L 3 , kg=m 3 g density of gas, m=L 3 , kg=m 3 o density of oil, m=L 3 , kg=m 3 solid density of solid, m=L 3 , kg=m 3 w density of water, m=L 3 , kg=m 3 standard deviation porosity ' parameter vector in Hough-transform-based parameterization facies map 139 Bibliography Aanonsen, S. I., Nvdal, G., Oliver, D. S., Reynolds, A. C., and Vall es, B. (2009). \The ensemble kalman lter in reservoir engineering{a review." SPE Journal, 14(03), 393{412 SPE-117274-PA. Agbalaka, C. C. and Oliver, D. S. (2008). \Application of the enkf and localiza- tion to automatic history matching of facies distribution and production data." Mathematical Geosciences, 40(4), 353{374. Berkowitz, B. (2002). \Characterizing ow and transport in fractured geological media: A review." Advances in water resources, 25(8), 861{884. Bocquet, M. and Sakov, P. (2012). \Combining in ation-free and iterative ensemble kalman lters for strongly nonlinear systems." Nonlinear Processes in Geophysics, 19(3), 383{399. Caers, J. (2003). \History matching under training-image-based geological model constraints." SPE Journal, 8(03), 218{226 SPE-74716-PA. Chang, H., Zhang, D., and Lu, Z. (2010). \History matching of facies distribution with the enkf and level set parameterization." Journal of Computational Physics, 229(20), 8011{8030. Chen, Y. and Oliver, D. S. (2010). \Ensemble-based closed-loop optimization applied to brugge eld." SPE Reservoir Evaluation & Engineering, 13(01), 56{71. Chen, Y. and Oliver, D. S. (2012). \Ensemble randomized maximum likelihood method as an iterative ensemble smoother." Mathematical Geosciences, 44(1), 1{26. 140 Chen, Y. and Oliver, D. S. (2013). \Levenberg{marquardt forms of the iterative ensemble smoother for ecient history matching and uncertainty quantication." Computational Geosciences, 17(4), 689{703. Chen, Y. and Oliver, D. S. (2014). \History matching of the norne full-eld model with an iterative ensemble smoother." SPE Reservoir Evaluation & Engineering, 17(02), 244{256. Dong, Y., Gu, Y., and Oliver, D. S. (2006). \Sequential assimilation of 4d seismic data for reservoir description using the ensemble kalman lter." Journal of Petroleum Science and Engineering, 53(1), 83{99. Dong, Y. and Oliver, D. S. (2005). \Quantitative use of 4d seismic data for reservoir description." SPE Journal, 10(01), 91{99 SPE-84571-PA. Dovera, L. and Della Rossa, E. (2007). \Ensemble kalman lter for gaussian mix- ture models." EAGE Petroleum Geostatistics 2007, 10-14 September 2007, Cascais, Portugal, Utrecht, The Netherlands: Extended Abstracts Book, EAGE Publica- tions BV, A16. Duda, R. O. and Hart, P. E. (1972). \Use of the hough transformation to detect lines and curves in pictures." Communications of the ACM, 15(1), 11{15. Evensen, G. (1994). \Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics." Journal of Geophys- ical Research: Oceans (1978{2012), 99(C5), 10143{10162. Evensen, G. (2003). \The ensemble kalman lter: Theoretical formulation and prac- tical implementation." Ocean dynamics, 53(4), 343{367. Evensen, G. (2004). \Sampling strategies and square root analysis schemes for the enkf." Ocean Dynamics, 54(6), 539{560. Furrer, R. and Bengtsson, T. (2007). \Estimation of high-dimensional prior and posterior covariance matrices in kalman lter variants." Journal of Multivariate Analysis, 98(2), 227{255. 141 Gaspari, G. and Cohn, S. E. (1999). \Construction of correlation functions in two and three dimensions." Quarterly Journal of the Royal Meteorological Society, 125(554), 723{757. Gassmann, F. (1951). \Elastic waves through a packing of spheres." Geophysics, 16(4), 673{685. G erard, A., Genter, A., Kohl, T., Lutz, P., Rose, P., and Rummel, F. (2006). \The deep egs (enhanced geothermal system) project at soultz-sous-for^ ets (alsace, france)." Geothermics, 35(5), 473{483. Gessner, K., K uhn, M., Rath, V., Kosack, C., Blumenthal, M., and Clauser, C. (2009). \Coupled process models as a tool for analysing hydrothermal systems." Surveys in geophysics, 30(3), 133{162. Gu, Y. and Oliver, D. S. (2007). \An iterative ensemble kalman lter for multiphase uid ow data assimilation." SPE Journal, 12(04), 438{446 SPE-108438-PA. Hamill, T. M., Whitaker, J. S., and Snyder, C. (2001). \Distance-dependent ltering of background error covariance estimates in an ensemble kalman lter." Monthly Weather Review, 129(11), 2776{2790. Han, D.-h., Nur, A., and Morgan, D. (1986). \Eects of porosity and clay content on wave velocities in sandstones." Geophysics, 51(11), 2093{2107. Hough, P. V. (1962). \Method and means for recognizing complex patterns (Decem- ber). Houtekamer, P. L. and Mitchell, H. L. (1998). \Data assimilation using an ensemble kalman lter technique." Monthly Weather Review, 126(3), 796{811. Houtekamer, P. L. and Mitchell, H. L. (2001). \A sequential ensemble kalman lter for atmospheric data assimilation." Monthly Weather Review, 129(1), 123{137. Iglesias, M., Lin, K., and Stuart, A. (2014). \Well-posed bayesian geometric inverse problems arising in subsurface ow." Inverse Problems, 30, 114001. 142 Iglesias, M. A. and McLaughlin, D. (2011). \Level-set techniques for facies identi- cation in reservoir modeling." Inverse Problems, 27(3), 035008. Jafarpour, B. and Khodabakhshi, M. (2011). \A probability conditioning method (pcm) for nonlinear ow data integration into multipoint statistical facies simula- tion." Mathematical Geosciences, 43(2), 133{164. Jafarpour, B. and McLaughlin, D. B. (2008). \History matching with an ensemble kalman lter and discrete cosine parameterization." Computational Geosciences, 12(2), 227{244. Kasiri, N. and Bashiri, A. (2011). \Status of dual-continuum models for naturally fractured reservoir simulation." Petroleum Science and Technology, 29(12), 1236{ 1248. Khodabakhshi, M. and Jafarpour, B. (2011). \Multipoint statistical characterization of geologic facies from dynamic data and uncertain training images." SPE Reservoir Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE, 9- 11 October (9-11 October). SPE-146935-MS. Liu, N. and Oliver, D. S. (2005a). \Critical evaluation of the ensemble kalman lter on history matching of geologic facies." SPE Reservoir Evaluation & Engineering, 8(06), 470{477 SPE-92867-PA. Liu, N. and Oliver, D. S. (2005b). \Ensemble kalman lter for automatic history matching of geologic facies." Journal of Petroleum Science and Engineering, 47(3), 147{161. Lorentzen, R. J., Nvdal, G., and Shaeirad, A. (2012). \Estimating facies elds by use of the ensemble kalman lter and distance functions{applied to shallow-marine environments." SPE Journal, 3(01), 146{158 SPE-143031-PA. Lu, L. and Zhang, D. (2015). \Assisted history matching for fractured reservoirs by use of hough-transform-based parameterization." SPE Journal (published online). 143 Moreno, D. and Aanonsen, S. (2007). \Stochastic facies modelling using the level set method." EAGE Petroleum Geostatistics 2007, 10-14 September 2007, Cascais, Portugal, Utrecht, The Netherlands: Extended Abstracts Book, EAGE Publica- tions BV, A18. Nan, T. and Wu, J. (2011). \Groundwater parameter estimation using the ensemble kalman lter with localization." Hydrogeology Journal, 19(3), 547{561. Nejadi, S. and Leung, J. (2012). \Integration of production data for estimation of natural fracture properties in tight gas reservoirs using ensemble kalman lter." SPE Canadian Unconventional Resources Conference, Calgary, Alberta, Canada (30 October-1 November). SPE-162783-MS. Oliver, D. S. and Chen, Y. (2011). \Recent progress on reservoir history matching: a review." Computational Geosciences, 15(1), 185{221. Ping, J. and Zhang, D. (2013). \History matching of fracture distributions by ensem- ble kalman lter combined with vector based level set parameterization." Journal of Petroleum Science and Engineering, 108, 288{303. Ping, J. and Zhang, D. (2014). \History matching of channelized reservoirs with vector-based level-set parameterization." SPE Journal, 19(03), 514 { 529 SPE- 169898-PA. Place, J., Garzic, E., Geraud, Y., Diraison, M., and Sausse, J. (2011). \Characteriza- tion of the structural control on uid ow paths in fractured granites." Thirty-Sixth Workshop on Geothermal Reservoir Engineering. Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T. (1986). Numer- ical recipes: the art of scientic computing, Vol. 1. Cambridge University Press. Sakov, P., Oliver, D. S., and Bertino, L. (2012). \An iterative enkf for strongly nonlinear systems." Monthly Weather Review, 140(6), 1988{2004. Sanjuan, B., Pinault, J.-L., Rose, P., G erard, A., Brach, M., Braibant, G., Crouzet, C., Foucher, J.-C., Gautier, A., and Touzelet, S. (2006). \Tracer testing of the 144 geothermal heat exchanger at soultz-sous-for^ ets (france) between 2000 and 2005." Geothermics, 35(5), 622{653. Schlumberger (2007). Eclipse reference manual 2007.2. Schlumberger Limited. Skjervheim, J. A., Evensen, G., Aanonsen, S. I., et al. (2007). \Incorporating 4d seismic data in reservoir simulation models using ensemble kalman lter." SPE Journal, 12(03), 282{292 SPE-95789-PA. Strebelle, S. (2002). \Conditional simulation of complex geological structures using multiple-point statistics." Mathematical Geology, 34(1), 1{21. Van Leeuwen, P. J. and Evensen, G. (1996). \Data assimilation and inverse methods in terms of a probabilistic formulation." Monthly Weather Review, 124(12), 2898{ 2913. Vogt, C., Marquart, G., Kosack, C., Wolf, A., and Clauser, C. (2012). \Estimating the permeability distribution and its uncertainty at the egs demonstration reservoir soultz-sous-for^ ets using the ensemble kalman lter." Water Resources Research, 48(8). Warren, J. and Root, P. J. (1963). \The behavior of naturally fractured reservoirs." Society of Petroleum Engineers Journal, 3(03), 245{255. Zhang, D. and Lu, Z. (2004). \An ecient, high-order perturbation approach for ow in random porous media via karhunen{loeve and polynomial expansions." Journal of Computational Physics, 194(2), 773{794. Zhang, Y. and Oliver, D. S. (2009). \History matching using a hierarchical stochastic model with the ensemble kalman lter: a eld case study." SPE Reservoir Simula- tion Symposium, Society of Petroleum Engineers. Zhao, Y., Reynolds, A. C., and Li, G. (2008). \Generating facies maps by assimilating production data and seismic data with the ensemble kalman lter." SPE Symposium on Improved Oil Recovery, Society of Petroleum Engineers. 145

Abstract (if available)

Abstract Natural or hydraulically induced fractures provide the primary flow capacity in many reservoir, such as tight sand and shale gas reservoirs, where the permeability of the rock matrix is very low compared with the conventional reservoirs. Hence, the success of production in fractured reservoirs is significantly dependent on the knowledge of the spatial distribution of fractures. However, it is a challenge to estimate fracture distribution by conventional history matching methods because of the complex geology of such reservoirs. Although there have been great progress in assisted history matching techniques during the last couple of decades, estimating fracture distributions in fractured reservoirs is still inefficient due to the strong heterogeneity and spatial discontinuity of model parameters. It is well known that the ensemble-based history matching methods, such as the Ensemble Kalman Filter (EnKF), performs best for models with Gaussian distributed parameters, since two-point statistics is used to represent the relationship between the model parameters and the observations. The performance of the EnKF can be far from optimum for the non-Gaussian distributed parameters, such as facies indicator, effective permeability, and porosity. Although the geometric shapes of fractures may be properly generated at the initial step, the geometrical shape and spatial continuity of the features are difficult to preserve after updating, resulting in geologically unrealistic fracture distribution maps. ❧ In this research, we develop an assisted history matching method for fractured reservoirs with a Hough-transform-based parameterization. The facies maps of fractured reservoirs are parameterized into Hough function fields in a discrete Hough space, while each grid block in the Hough space represents a fracture defined by its two Cartesian coordinates: angle θ of its normal and its algebraic distance ρ from the origin in the flow domain. The length and axial position of the fractures are defined by two additional parameters nested on the same grid. The Hough function value of each grid block in the Hough space is used as the indicator of the existence of the fracture in the facies map. Facies maps can be uniquely generated from the Hough space parameters by an inverse transform, which guarantees that no matter how the Hough space parameters are modified, the fractures in facies maps are always geologically realistic. ❧ In assisted history matching with Hough-transform-based parameterization, the parameter fields in the Hough space, instead of the facies maps, are updated with production history and seismic data. Before each fluid flow simulation during the history matching, an inverse transform is performed to generate facies maps for the reservoir simulator. Point-wise prior information such as known fractures discovered from well log data, as well as the statistics of regional fracture orientations, can be honored by the inverse transform. The effectiveness of this method is investigated with two dimensional synthetic waterflooding examples. Seismic data is also assimilated in some of the examples to improve the characterization of fractures away from wells. The fracture distributions in the reference maps are identified by this method, and updated models are capable to provide improved predictions for prolonged periods of production. A few realistic examples based on the tracer tests at Soultz-sous-Forêts geothermal reservoir are used to demonstrate the performance of proposed method in characterizing complex geological structure with fault network.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

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

PDF

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

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

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

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

Subsurface model calibration for complex facies models

$Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data$

PDF

Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Modeling and simulation of complex recovery processes

PDF

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

PDF

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

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

PDF

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

PDF

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

PDF

Deep learning for subsurface characterization and forecasting

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

Asset Metadata

Creator Lu, Le (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 09/03/2015

Defense Date 08/06/2015

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

Tag ensemble Kalman filter,fractured reservoirs,history matching,Hough transform,OAI-PMH Harvest

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Zhang, Dongxiao (committee chair), Ershaghi, Iraj (committee member), Ghanem, Roger (committee member)

Creator Email lelu.usc@gmail.com,lelu@usc.edu

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

Unique identifier UC11273108

Identifier etd-LuLe-3870.pdf (filename),usctheses-c40-174032 (legacy record id)

Legacy Identifier etd-LuLe-3870.pdf

Dmrecord 174032

Document Type Dissertation

Format application/pdf (imt)

Rights Lu, Le

Type texts

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

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

Repository Name University of Southern California Digital Library

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