University of Southern California Dissertations and Theses

Control of displacement fronts in porous media by flow rate partitioning

(USC Thesis Other)

Control of displacement fronts in porous media by flow rate partitioning

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content CONTROL OF DISPLACEMENT FRONTS IN POROUS MEDIA BY FLOW RATE PARTITIONING by Mohsen Heidary-Fyrozjaee A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) May 2008 Copyright 2008 Mohsen Heidary-Fyrozjaee Dedication This dissertation is lovingly dedicated to my parents, who in spite of their limitations and obstacles in life always continued to support and inspire me as I pursued my dreams. ii Acknowledgments My deepest thanks and respect first go to my thesis advisor, professor Yannis C. Yortsos whose support, guidance, and friendship were essential in my progress as a student. I am also grateful to the members of my committee, professor Ershaghi and professor Sadhal for their time, encouragement and expertise. I would like to thank my friends and the individuals who helped me to get where I am today. Finally I would like to thank Kija who believed in me and gave me confidence and support to finish my school. iii Table of Contents Dedication ii Acknowledgments iii List of Figures vii Abstract xi 1 Chapter 1: Introduction 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Chapter 2: Control of Displacement Fronts in Homogeneous Media 12 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1 The Governing Potential Equation . . . . . . . . . . . . . . . . 13 2.2.2 Front Kinematic Equation . . . . . . . . . . . . . . . . . . . . 15 2.2.3 Fredholm Integral Equation of the First Kind . . . . . . . . . . 15 2.3 Homogeneous and Isotropic Porous Media . . . . . . . . . . . . . . . . 16 2.3.1 Pressure Solution . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.2 Front Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.1 Numerical Procedure . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.2 Front/Particle Tracking . . . . . . . . . . . . . . . . . . . . . . 25 2.4.3 Example I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4.4 Example II . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4.5 Example III . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.5.1 Degree of Ill-conditioning . . . . . . . . . . . . . . . . . . . . 35 2.5.2 Optimum Regularization Parameter . . . . . . . . . . . . . . . 36 iv 2.5.3 Admissible Fronts . . . . . . . . . . . . . . . . . . . . . . . . 40 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3 Chapter 3: Control of Displacement Fronts in Heterogeneous Media 50 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2.1 Two-Well Problem . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3 Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.1 Two-Well Problem . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.2 Interpolation Methods . . . . . . . . . . . . . . . . . . . . . . 56 3.3.3 Front/Particle Tracking . . . . . . . . . . . . . . . . . . . . . . 58 3.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.1 Anisotropic Porous Media . . . . . . . . . . . . . . . . . . . . 59 3.4.2 Heterogeneous Porous Media . . . . . . . . . . . . . . . . . . 72 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4 Chapter 4: The Inverse Problem 91 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.2 Fredholm Integral Equation of the First Kind . . . . . . . . . . . . . . 92 4.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.2.2 Discrete Ill-posed Problems/Ill-conditioning . . . . . . . . . . . 93 4.2.3 Singular Value Expansion . . . . . . . . . . . . . . . . . . . . 96 4.2.4 Singular Value Decomposition . . . . . . . . . . . . . . . . . . 97 4.3 Solution Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.3.1 Least Square Method . . . . . . . . . . . . . . . . . . . . . . . 98 4.3.2 The Picard Condition . . . . . . . . . . . . . . . . . . . . . . . 99 4.4 Regularization Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.4.1 Tikhonov Regularization Technique . . . . . . . . . . . . . . . 103 4.4.2 Truncated SVD . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.4.3 Semi-Iterative Conjugate Gradients . . . . . . . . . . . . . . . 106 4.4.4 The Regularization Parameter . . . . . . . . . . . . . . . . . . 106 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5 Chapter 5: Analysis of Forward Problem 108 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.2 Semi-Infinite Reservoir . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.3 Method of Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.4 Special Case: One Harmonic Front Solution . . . . . . . . . . . . . . . 114 5.4.1 With Constant Coefficient . . . . . . . . . . . . . . . . . . . . 114 5.4.2 With Variable Coefficient . . . . . . . . . . . . . . . . . . . . . 116 5.5 Special Case: Two Harmonics Front Solution . . . . . . . . . . . . . . 118 5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 v 6 Chapter 6: Conclusions and Recommendations 120 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.2 Future Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . 122 References 124 vi List of Figures 2.1 A simplified schematic of the front control problem . . . . . . . . . . . 14 2.2 A simple schematic of the front location and its correspondent system of equations relating the velocity at the front to the injection rates along the injection well . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3 Desired front location as a function of time for example I . . . . . . . . 26 2.4 Solution to Example I: rate solutions obtained by standard method ver- sus rate solutions obtained by regularization method . . . . . . . . . . . 28 2.5 Validation of the solution to Example I . . . . . . . . . . . . . . . . . . 29 2.6 Solution to Example II . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.7 Validation of the solution to Example II . . . . . . . . . . . . . . . . . 32 2.8 Solution to Example III . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.9 Validation of the solution to Example III . . . . . . . . . . . . . . . . . 34 2.10 How fast singular values decay, determines the degree of ill-conditioning of the system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.11 The L-curve of front in example I at the10 th time step . . . . . . . . . 39 2.12 Three choices of the optimum regularization parameters for example I . 40 2.13 The effect of different regularization parameters on the solution . . . . . 41 2.14 Three choices of the optimum regularization parameter for example III . 42 2.15 The effect of different regularization parameters on the solutions . . . . 43 2.16 Different values of¾ and² defines fronts with different sharpness. The front for example I used previously is also shown on both figures., (a) effect of varying².,(b) effect of varying¾ . . . . . . . . . . . . . . . . 44 vii 2.17 The desired front dynamics of example I with four different character- istic forms corresponding to different values of² and¾ . . . . . . . . . 46 2.18 The rate solutions of example I with four different characteristic forms corresponding to different values of² and¾ . . . . . . . . . . . . . . . 47 2.19 The comparison between desired and achieved front profiles of exam- ple I with four different characteristic forms corresponding to different values of² and¾ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.1 Super-imposition of multiple sources in potential flow . . . . . . . . . . 53 3.2 Schematic of interpolation process . . . . . . . . . . . . . . . . . . . . 56 3.3 Validation of the numerical method vs. analytical one . . . . . . . . . . 60 3.4 The desired front progression in time in examples I and III at larger time intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.5 Homogeneous porous media for example I (a) the corresponding rate profile solution. , (b) the desired front progression versus the achieved one 63 3.6 Anisotropic system with K y Kx = 0:1 for example I (a) the correspond- ing rate profile solution , (b) the desired front progression versus the achieved one . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.7 Anisotropic system with Ky K x = 0:01 for example I (a) the correspond- ing rate profile solution , (b) the desired front progression versus the achieved one . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.8 Anisotropic system with K y Kx =0:0001 for example I (a) the correspond- ing rate profile solution., (b) the desired front progression versus the achieved one . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.9 Homogeneous porous media (a) the corresponding rate profile solution. , (b) the desired front progression versus the achieved one . . . . . . . . 68 3.10 Anisotropic system with K y Kx = 0:1 for example III (a) the correspond- ing rate profile solution. , (b) the desired front progression versus the achieved one . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.11 Anisotropic system with Ky K x = 0:01 for example III (a) the correspond- ing rate profile solution . , (b) the desired front progression versus the achieved one . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 viii 3.12 Anisotropic system with Ky K x = 0:0001 for example III (a) the corre- sponding rate profile solution . , (b) the desired front progression versus the achieved one . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.13 Permeability field increasing iny direction and the desired front dynamics 75 3.14 Front propagation corresponds to the permeability field increasing in y direction showed in Figure 3.13(a)., (a) for the equal injection rate policy .,(b) and for rates inversely proportional to the permeability at the injection points . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.15 The obtained rate and the comparison between desired and achieved front 77 3.16 Permeability field and the desired front dynamics . . . . . . . . . . . . 79 3.17 Front propagation corresponds to the permeability field increasing in y direction showed in Figure 3.16(a)., (a) for the equal injection rate policy .,(b) and for rates inversely proportional to the permeability at the injection points . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.18 The obtained rate solution and the comparison between desired and achieved front . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.19 Permeability field and the desired front dynamics . . . . . . . . . . . . 83 3.20 Front propagation corresponds to the permeability field increasing in y direction and showed in Figure 3.19(a)., (a) for the equal injection rate policy .,(b) for rates inversely proportional to the permeability at the injection points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.21 The obtained rate solution and the comparison between desired and achieved front . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.22 Permeability field and the desired front dynamics . . . . . . . . . . . . 87 3.23 Front propagation corresponds to the permeability field varies in bothx andy directions shown in Figure 3.22(a)., (a) for the equal injection rate policy .,(b) and for rates inversely proportional to the permeability at the injection points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.24 The obtained rate solution and the comparison between desired and achieved front . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.1 Discrete Picard plot for the unperturbed (top) and perturbed (bottom) data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 ix 4.2 Discrete Picard plot is not satisfied for the integral equation 4.16 . . . . 102 4.3 The generic form of the L-curve plot . . . . . . . . . . . . . . . . . . . 105 5.1 The semi-infinite reservoir . . . . . . . . . . . . . . . . . . . . . . . . 110 5.2 The amplitude of the harmonic decreases inversely proportional to its frequency. The results for first six harmonics are presented here. As time increases, the sensitivity of the front to the rate fluctuations decreases inversely proportional to the frequency . . . . . . . . . . . . . . . . . . 117 x Abstract In modern reservoir management, optimizing the recovery efficiency of displacement processes is a fundamental problem with important economic ramifications. Given a well configuration, the optimum displacement efficiency may be achievable through controlling injection rates. Smart wells and intelligent completion have gained signif- icant attention in recent years in the field of dynamic optimization of Enhanced Oil Recovery processes. Implementing separate downhole control valves enables the appli- cation of optimal rate control in real time, through the appropriate valve settings in injection and production wells. In this thesis, we consider a fundamental question in rate-control flow problems, namely the control of a displacement front via flow-rate partition in a horizontal well. The specific question we address is how to partition the flow rate within the well, so that the displacement front can be steered according to pre-determined dynamics, or as is necessary. The ability to steer a front at will is of obvious significance. Given the reservoir geology, one can for example, steer the front away from flow obstacles, affect a piston-like displacement in heterogeneous formations or otherwise control the displacement process as needed. xi The objective of this thesis is to investigate the steering of the front by partitioning the injection rates along the well. First, the problem is addressed in a simple rectangu- lar geometry based on potential flow. It is assumed that the flow is potential and that the displacement is at a unit mobility ratio. These assumptions are for mathematical convenience only and can be relaxed. They allow, however, significant insight into the problem. When the reservoir is homogeneous and isotropic, an integral equation in an analytical form is derived, the solution of which determines the desired injection rate profile. A similar approach applies for an anisotropic or a heterogeneous system, except that the kernel in the integral equation must be determined numerically. This is done either by repeated calculations of the Green’s function in a heterogeneous system or for a modified two well system. For the solution of the integral (Fredholm) equation, a reg- ularization technique is necessary. However, it is found that numerical instabilities do develop, even with the use of regularization, for later times, when transverse cross-flow is large. Conversely, the instabilities diminish with a more stratified structure. The appli- cation of the technique is highlighted and extended to heterogeneous systems where its complexity and the solution were studied. Also in a forward displacement problem, the effect of rate fluctuations on the front dynamics were analyzed in a semi-infinite porous media. We considered the special case of injection rate with one harmonic term. It was found that the front is less sensitive to the fluctuations in rate for larger harmonic terms. The results of this study find applications to the rapidly emerging field of smart wells and the optimization of displacement problems in oil reservoirs using flow rate control. xii Chapter 1 Introduction 1.1 Introduction Oil production from subsurface reservoirs is achieved by a number of processes, includ- ing oil and gas expansion or the injection of a displacing fluid. Both methods may be used simultaneously. When the displacement of oil by another fluid is the main control- ling process, the problem arises of the advancement of the interface between the fluid in place and the displacing fluid. Typically, the displacement fluid is injected through a number of injection wells and the displaced fluid is produced through a number of production wells. Water flood, gas flood and steam injection are examples of such dis- placement processes. Injection and production wells arrangement, varying the injection and production rates of an individual well, and varying the rates of individual segments along the injection and production wells, all affect the flood front movement through the porous medium. Flood front control is a key to achieve an optimum displacement efficiency. The common well configuration includes the well know five-spot and line drive patterns for vertical wells, and toe-to-heel water flooding, air injection and thermal flooding for horizontal wells. Usually geologic, economic and other field practical limitations deter- mine the well configuration. With a pre-determined well configuration, the basic control parameter for controlling the interface progress between fluids is the injection rate allocation of various injection 1 wells as well as the rate partitioning within the individual wells. Once the ability to steer a front is accomplished, one will be able to achieve the optimum recovery efficiency of the displacement process. At present, the problem of steering the flood front at will by manipulating injection rates along the well has not been systematically investigated. To achieve a high recov- ery efficiency, there exist two conventional approaches: one is to allocate the injection rates equally to all injection wells, the other is to apply variable injection rates. The former practice assumes that the reservoir is homogeneous and the geometric pattern is symmetric, while the later takes into account the inherent heterogeneity of the porous medium. As a well trajectory cuts through the reservoir, it is intuitive to partition the injection rate along the well trajectory, so as to control and steer the flood front accord- ingly. Smart wells and intelligent completion allow for implementing separate downhole control valves along the well trajectory, thereby making possible a rate partitioning. Injection at constant and equally allocated rates as well as variable rates are not the only policies used in the displacement process to achieve high recovery. An intuitive and efficient way of achieving high recovery would be to determine the best rate partitioning along the well trajectory. This may allow steering the flood front as desired in order to achieve an optimum displacement efficiency. 1.2 Background The subject of optimization for reservoir development has been addressed by numer- ous scholars and professionals within the petroleum industry as early as 1919 [Ros19]. In the early stages, determining the optimum number of wells and well spacing were the objectives of the optimization [Ros19, Phe29, BW30, Cor30, Has30]. Muskat 2 [WBM34, MW34] was the first to explore mathematically the problem of wells spacing and injecton/production well configuration. With the assumption of Darcy’s law and its analogy to the law of electrical conduction, he modeled the performance of five-spot, seven-spot and line drive flood patterns and determined the optimum recovery at the time of breakthrough. In the problem of well spacing, Muskat [Mus40] concluded that the effect of well spacing on flood efficiency is insignificant. Rather, the appropriate well spacing must be relegated to the economics of the project. Following that work, and for the next two decades, focus shifted toward develop- ing techniques for determining the water intake zones in injection wells [Pfi48, CBJ46, Geo59, Jos57] and water entry zones within production wells [SB48, Dal49]. The objec- tive was to achieve higher oil production by avoiding injection through thief zones which was causing early water breakthrough in production wells leading early well abundant. Optimization of field operations was introduced by Aronofsky [Aro62] and Lee and Aronofsky [LA58], where a time discretized optimization process applied to a set of single and double well reservoirs. With the net profit over eight years as the objec- tive function, the authors applied a linear programming approach to obtain the optimum schedule of production rates for each time period and for all reservoirs. In other similar studies, the optimal scheduling of production and compression in the development of gas field has been considered [CHD69, ME78, RG74, Hup74]. In the optimization phase, Cooksey [CHD69] and Murray [ME78] used linear programming, while Rosenwald et al. [RG74] adopted mixed-integer programming and Huppler [Hup74] applied nonlin- ear programming to solve the optimum reservoir development plan. Bernard and Horne [SH83] used a commercial reservoir simulator to model the reservoir dynamics instead of a simple linear representation of the relationship between well rates and pressure. 3 Again linear programming was used for the optimization phase of the study. Bitten- court et al. [BH97] and Hardy et al. [HRK98] employed a hybrid algorithm based on direct methods such as Genetic Algorithms (GA) and Polytope search, to the production scheduling of an oil and gas field. Dynamic optimization of the displacement process has gained significant attention in the last two decades. Optimal control methods are the backbone of the dynamic opti- mization. Its application in displacement processes has been introduced by Ramirez and his colleagues [RFL84, FR84, LRQ93]. The authors applied optimal control the- ory of distributed systems to determine the best injection policy of a surfactant slug for a tertiary oil recovery chemical flood and the steam quality of a steam flooding. A steepest descent gradient method was used to achieve the optimum solution of the prob- lem. The control variable to be determined was the size or the quality of the injected agents (the slug size of a miscible flood or the steam quality of steam injection). The authors coupled a numerical simulation with the optimal control strategy and achieved an improvement of 20% in recovery metrics defined over the initial policy for a pair of injection and production wells in a symmetric and homogeneous reservoir. The process of water displacing oil gained more attention for purposes of dynamic optimization, since it is a less complex process and the injection fluid composition is fixed. Here, the only control parameters are the rate allocation of the injection fluid to the injectors and of produced fluids to the producers [SY01]. In addition, recovery efficiency in terms of an objective function can be defined in such a way that will lead to corresponding optimal rate profiles. Asheim [Ash88] applied a general optimiza- tion technique to different production scenarios, such as natural water drive and water drive by injection. The objective function was the NPV (Net Present Value), with corre- sponding improvements found in the range of 2-11%. His method combined reservoir 4 simulation with the control theory practice of implicit differentiation. In the reservoir optimization approach, a 2-D, two-phase simulator was used. Asheim used the IMPES technique, instead of a linear form of the pressure-production relationship, which is only applicable for the single phase reservoir and was common at the time [Aro62, ME78]. In the IMPES technique, at each time step, the pressure is solved implicitly, and then the saturation is calculated explicitly using the known pressure. The hypothetical reser- voir models used by Asheim included two cases: the first case is a reservoir with two injection and one production wells, the second case has only production wells, and a natural water drive. In the optimized production scheme, injection and production rates are changing in time to achieve higher recovery compared to the constant injection and production rates policy. Sudaryanto and Yortsos [Sud98, SY00] provided a systematic approach for the dynamic optimization of displacement problems, in which the objective function was the displacement efficiency at breakthrough. They formulated the control problem for mul- tiple injection wells and one production well. Using potential flow in 2-D, in the absence of dispersion and gravity, Sudaryanto and Yortsos applied optimal control methodology and showed that the optimal injection policy is the ”bang-bang” type. This means that the control parameters, namely the injection rates at the individual source points, must take their extreme values (maximum or minimum admissible rate). Then the problem becomes finding optimal switch times for each well. For the case of three wells there is only one switch time. On the other hand, depending on the problem type, multiple switch times are possible. The authors achieved higher displacement efficiency for the optimized case compared to the conventional equal rate injection policy. For the particles that define the front, the problem is formulated as a non-linear dynamic. For homogeneous media, this formulation is analytic. In the heterogeneous 5 media the authors have used a superposition approach to separate time and spatial depen- dences, which reduces extensively the computational expenses. By conducting sensitiv- ity studies in the heterogeneous media, the authors showed that the range of improve- ment in displacement efficiency depends on the particular realization of the permeabil- ity field, as well as on the presence or absence of spatial correlation. Nevertheless, the advantage of the bang-bang control policy compared to constant rate injection is still valid, even in the absence of information on the heterogeneity. Smart wells and intelligent completion have gained significant attention in recent years in the field of dynamic optimization of enhanced oil recovery processes. The abil- ity to implement separate downhole control valves enables the application of optimal rate control in real time through the appropriate valve settings in injection and produc- tion wells. Advancement in the field of intelligent completion has opened new opportu- nities in the field of optimization. Yetan [YDA02] applied an optimization technique based on conjugate gradients to optimize the operation of non-conventional wells. Non-conventional wells (wells with an arbitrary trajectory and a number of branches or laterals) provide great potential in maximizing oil recovery. The objective of his study was to develop a general proce- dure for the optimal deployment of non-conventional wells. Because of the complexity of the problem and a large number of variables, gradient based methods are not suit- able. Genetic Algorithm (GA), a stochastic search engine based on the general principle of Darwinian evolution, was used as an optimization algorithm along with a number of acceleration or “helper“ routines such as Artificial Neural Network. Yetan linked the optimization algorithm with a commercial simulator to model the inflow control devices. The main constraints of the optimization problem were the following; the proposed well must lie completely within the reservoir, no trajectory can intersect any 6 other wells or lateral trajectories and the proposed trajectories must honor any speci- fied anisotropy/stress directions. The optimal well type and recovery improvement vary depending on reservoir model and objective function. An improvement of up to 30% was achieved. In recent works, Brouwer et al. [BJvdSS + 01, BJ02] published a series of papers investigating the optimization of water flooding in heterogeneous reservoirs having mul- tiple segments along both injection and production wells. The model consisted of two smart horizontal wells, in both injection and production, with multiple Inflow Control Valves, (ICVs). In the static approach, optimization of the water flooding was based on heuristic algorithms resulting in fixed valve settings, which were kept fixed during the displacement process. Based on the productivity index (PI) of each segment, to delay water breakthrough, the algorithms were selected such that the flow paths were moved away from the higher permeability zones. In one specific scenario, where there is a high permeability strip perpendicular to the well trajectory, a higher improvement of sweep efficiency was achieved by applying the static optimization approach. The approach was then extended to dynamic conditions, using optimal control theory to maximize ultimate recovery by continuously controlling the valve settings. As in all previous studies, significant improvements in ultimate recovery were obtained compared to both the non-optimized scenario and the static optimized scenario. There are few studies on the subject of dynamic optimization of the displacement process using streamline technology. Grinestaff et al. [Gri99, GC00] used streamline simulation to highlight the waterflood efficiency in Northwest Fault Block of the Prude- hole Bay. For each producer the fraction of support from each injector was defined as a bundle of streamlines connecting injectors to the producers. Because the well rates 7 change during the production, the streamline allocation information is dynamic. The process must be done on a history matched reservoir. Alhutalhi et al. [AODG06] provided a streamline-based approach to determine opti- mal injection and production rates, maximizing the sweep efficiency and delaying the arrival time of the water front. The main goal was to equalize the arrival times of the water front at all producers in a selected sub-region of the waterflood project. The arrival time to the producer is defined as the time required for the front to reach the producer from its current position. It is calculated by the average of the fastest (top 20%) stream- lines. The assumption is that the streamlines do not shift significantly because of a small perturbation in rates. In addition, the water front arrival time at produceri will be mainly sensitive to the production at the produceri and relatively insensitive to the rates of the other wells. Streamlines were used to calculate the sensitivity of arrival time to injection rates analytically in one single run. For the field application, Alhutali et al. modified the arrival time as the time required to travel from the injector to the producer in order to continue the optimization even after the breakthrough. Thiele and Batycky [TB06] applied streamlines to define dynamic well allocation factors (WAF) between injector and producer wells. The WAF information was used to determine each phase rates at either end of each injector-producer pair. With this information at hand, they proposed the idea of injection efficiency for injector as well as each injector-producer pair. The underline idea is very similar to Grinestaff [Gri99, GC00]. The major difference with Alhutali’s work and other gradient based methods is that there is no formal attempt to optimize the field performance based on a minimization of an objective function. The method is solely based on a heuristic reservoir engineering approach. For every injector, the efficiency is defined as the ratio of offset oil production to water injection. After all injector-producer pairs injection efficiencies are known, then 8 with a target update rate method, water is reallocated from inefficient well connectors to the efficient ones. The incremental change in rate is a function of injection efficiency of the well. The reference efficiency is arbitrary and is selected to be the field average injection efficiency. The new rate is updated at certain number of times (time steps), where the wells with an injection efficiency higher than the reference field injection efficiency will see an increase in rate, while the wells with injection efficiency lower than the injection efficiency of the field will see a decrease in rate. The authors have not addressed an automated way of forecasting and allocating the injection and production rates. 1.3 Motivation Regardless of the different methodologies being used, one can state that in the above studies, the main objective is to control the interface propagation. For instance, by maximizing Net Present Value in a waterflood, one intends to achieve injection and production rate profiles that delay water production and enforce a uniform displacement. Therefore, the ability of having control over the front and its propagation has significant value in achieving optimum displacement efficiency. The objective of this thesis is to investigate the control of a displacement front through porous medium by partitioning the injection rate along the horizontal well. We address the fundamental question in rate flow problem. The specific question is what is the injection rate policy at various injection points along the well in order for a particular front to follow a pre-described dynamics in a possibly heterogeneous porous medium. The injection points are independent and can be controlled separately. Clearly, partition- ing injection rates along the well trajectory directly influences the dynamics of the fluids 9 in the porous medium. To understand the dynamics of the displacement process and the extend in which one can control the front, the problem is considered for a potential flow (mobility ratio of one) for a simple two dimensional homogeneous rectilinear sys- tem. With the insight offered by the potential flow solution in the homogeneous porous media, the more general problem involving heterogeneous reservoirs are explored. 1.4 Thesis Outline The thesis is organized as follows: In chapter 2, the general form of the front control problem is presented. The front control problem is considered in a two-dimensional rectangular porous medium with injection well at the left side and production at the right. The porous medium is homogeneous and isotropic. We utilized the method of separation of variables to describe the pressure field and coupled it with the front kine- matics to arrive at a non-linear relation in the form of Fredholm Integral Equation of the First Kind. Through several numerical experiments, we describe the front steering tech- nique and its ability to control fronts by partitioning injection rates along the well. In addition the degree of ill-conditioning, source of errors and the optimum regularization parameters are studied. In chapter 3, we consider the same front control of a first-contact miscible displacement at unit mobility ratio in a similar geometry except that the porous medium is heterogeneous. We illustrate the ability to control the front through different heterogeneous porous media. In the chapter we utilize the equal injection policy and the injection rate inversely proportional to the permeability at injection points and com- pare their front propagation with the controlled front propagation scenario. Through the comparison, we show the potential applicability of the front control techniques in maximizing efficient sweeps. In chapter 4, we present a review of inverse problems and 10 their solution techniques that provides the background for the front control problem in obtaining the rate solution. The inverse problem, the Fredholm Integral Equation of the First Kind, its discretized form, its analysis, its solution and the necessary regulariza- tion techniques are all presented. In chapter 5, we present an asymptotic analysis of the displacement front in the forward problem. The injection rate is assumed to have single harmonic term in a semi-infinite homogeneous porous medium. We present an analysis that can be used to explain the diminishing front sensitivity to the injection rate at larger times where front moves away from injection side. We conclude the study in chapter 6, where future recommendations on the control of the displacement front through porous media by partitioning injection rates along the well are discussed. 11 Chapter 2 Control of Displacement Fronts in Homogeneous Media 2.1 Introduction In this chapter, we consider the control of the displacement front through porous media by partitioning injection rates accordingly along the well. We assume a first contact mis- cible displacement at unit mobility ratio and a two-dimensional, areal and homogeneous reservoir with two horizontal wells, one injector at the left side and one producer at the right. In the absence of gravity and dispersion, potential flow governs the dynamics of the displacement process, which simplifies the mathematical formulation considerably. The objective of the problem is to provide an answer to the question of how to con- trol or steer a displacement front through porous media by variable injection rates along the well. The well can be either vertical or horizontal and is equipped with individu- ally controllable segments (or perforation intervals). The coupling of the pressure field solution and the front kinematics equation leads to an integral form, a Fredholm Integral Equation of the First Kind. Solutions to this equation are presented. The control parameters to be determined, are the injection rates along the injection well trajectory. These rates are to be determined based on the desired trajectory of the front dynamics, which is known in advance. In the absence of gravity and dispersion, and for a homogeneous porous medium, an analytical formulation is derived that relates 12 the control parameter to the front dynamics. As the front moves away from the injection well, however, an instability in the solution develops and a regularization technique is necessary. Several examples are used to assess the robustness of the solution technique and to determine the sensitivity of the control parameter to the sharpness of the front dynamics. The extent to which the regularization method is effective and the source of errors causing instability in solution are discussed as well. In the next section and its subsections, the general problem statement and its mathe- matical formulation will be presented, along with the Fredholm Integral Equation of the First Kind as it corresponds to the general problem statement. 2.2 Problem Formulation In this section a general methodology and formulation in a simple rectangular geometry is presented. A simplified schematic of the problem considered is shown in Figure 2.1 with the appropriate boundary conditions and a typical front location. The injection points along the well at the left boundary are independent of each other and therefore can be controlled individually. The producing well is at the right boundary and is at constant pressure zero. 2.2.1 The Governing Potential Equation Under the assumption of incompressible miscible fluids of equal mobility and in the absence of dispersion and gravity effects, the governing equation in the reservoir is @ @x (k(x;y) @P @x )+ @ @y (k(x;y) @P @y )=0 (2.1) 13 Figure 2.1: A simplified schematic of the front control problem where the P is the normalized flow potential and k(x;y) is the permeability at point (x;y). For the geometry of Figure 2.1, the appropriate boundary conditions in dimen- sionless notations are ¡ @P @x j x=0 = q(y;t) (2.2) P j x=1 = 0 (2.3) ¡ @P @y j y=0 = 0 (2.4) ¡ @P @y j y=1 = 0 (2.5) The injection points are located along the injection well located at the left side of the two-dimensional reservoir considered. The unknown rate profile along the well trajec- tory isq(y;t). The production well at the right side of the reservoir is at zero pressure. In the following sections the homogeneous and isotropic reservoirs will be addressed 14 first, where we obtain analytical solutions, and then investigate an extension to the more general and heterogeneous reservoir. 2.2.2 Front Kinematic Equation In displacement process, the interface of the displacing and displaced fluid is the locus of all the displacing agent particles that form the interface. This interface can be written as F(x;y;t)=0 (2.6) wheret is time. The front can be also recast as F(x;y;t)=x¡f(y;t) (2.7) where f is the front location. Because the front moves with the fluid, the material derivative is zero dF(x;y;t) dt = @F @t +v x @F @x +v y @F @y =0 (2.8) where we rescaled the time appropriately. Equivalently, using equation 2.7 we have @f @t =v x ¡v y @f @y (2.9) 2.2.3 Fredholm Integral Equation of the First Kind The solution of the equation 2.9 describes the evolution of the front in time. However, the pressure field must be either previously determined or determined simultaneously with the front. Here by substituting an appropriate form of the pressure solution of 15 equation 2.1, as described in the next section, into equation 2.9, we obtain a generic Fredholm Integral Equation of the First Kind @f(y;t) @t = Z 1 0 K(0;´;f(y;t);y;t)q(´;t)d´ (2.10) where K(0;´;f(y;t);y;t) is the kernel and q(´;t) is the injection rate at the injection point (0;´). The above equation is a coupling of the pressure field and the evolving front as a function of time. The kernel of the integral equation is a function of the front dynamics and the reservoir properties. Equation 2.10 is a Fredholm Integral Equation of the First Kind, the solution of which will allow for the control profile, namely the injection rates along the well, to be determined for proper time intervals. 2.3 Homogeneous and Isotropic Porous Media In the case of a homogeneous medium, the transverse and flow directions are scaled by H andL respectively, in which case the dimensionless parameterR L = L H q Kv K L is equal to one. The dimensionless parameterR L is a key parameter that controls the degree of vertical flow equilibrium and suggested by Lake [Lak89] and analytically shown to be so by Yortsos et al. [Yor91, YYS02]. H and L denote the width and length of the porous medium and K v and K L are the values of the transverse and the flow direction permeability. For a general anisotropic porous medium where R L 6= 1, one should replace equations 2.1 and 2.9 respectively by the following 1 R 2 L @ 2 P @ 2 x + @ 2 P @ 2 y =0 (2.11) 16 @f @t =¡ 1 R 2 L @P @x + @P @y @f @y (2.12) For a two-dimensional rectilinear and homogeneous porous media, an analytical expres- sion is possible for the pressure by the method of Separation of Variables. 2.3.1 Pressure Solution In a homogeneous two-dimensional porous media, the governing pressure equation is the Laplace’s equation. Its solution can be uniquely determined if the value of the func- tion (Dirichlet boundary condition), its normal derivative (Newmann boundary condi- tion), or a combination of both are specified. As depicted in the schematic of the prob- lem in Figure 2.1, there are one Dirichlet and three Newmann boundary conditions. We restate the problem @ 2 P @ 2 x + @ 2 P @ 2 y =0 (2.13) ¡ @P @x j x=0 = q(y;t) (2.14) P j x=1 = 0:0 (2.15) ¡ @P @y j y=0 = 0:0 (2.16) ¡ @P @y j y=1 = 0:0 (2.17) Because the problem is linear, the solution can be obtained by a superposition [Iva79, Ozi80] P(x;y)= Z 1 0 G(0;´;x;y)q(´;t)d´ (2.18) where the kernel G(0;´;x;y) can be obtained as follows: The Separation of Variables assumes that the solution to the partial differential equation 2.13 can be written as a 17 product of the solution of two separate ordinary differential equations with respect to each of the two independent variables P(x;y)= 1 X n=0 X n (x)Y n (y) (2.19) We substitute the above into the Laplace equation 2.13 to find ¡ X " n (x) X n (x) = Y " n (y) Y n (y) =¡¸ 2 n (2.20) where the eigenvalues¸ n need to be determined. The two homogeneous boundary con- ditions aty =0 andy =1 allow us to fix the eigenvalues. The general solution is Y n (y)=B n cos(¸ n y) (2.21) where¸ n =n¼. Likewise ¡ X " n (x) X n (x) =¡¸ 2 n (2.22) with solution n=0; X 0 (x)=A 0 (x)+B 0 (2.23) n6=0; X n (x)=C n sinh(¸ n x)+D n cosh(¸ n x) (2.24) Applying the boundary conditions at x = 0 and x = 1 allows us to determine the unknown coefficients X 0 (x)=A 0 (x¡1) (2.25) 18 A similar procedure can be applied forX n (x), to arrive at the following expression: X n (x)=A n sinh(¸ n (1¡x)) (2.26) Combining the two we have P(x;y)=A 0 (x¡1)+ 1 X n=1 A n sinh(¸ n (1¡x))cos(¸ n y) (2.27) Expression 2.27 is the general pressure solution that satisfies the boundary conditions at y = 0,y = 1 andx = 1. There is one more boundary condition that must be satisfied. Applying the Fourier cosine series technique to the boundary condition functionq(y;t) allow us to obtain an expression for the coefficientA n . ¡ @P @x j x=0 =q(y;t) =¡A 0 + 1 X n=1 A n ¸ n cosh(¸ n )cos(¸ n y) = 1 2 a 0 + 1 X n=1 a n cos(¸ n y) (2.28) Where a 0 = 2 Z 1 0 q(y 0 ;t)dy 0 (2.29) a n = 2 Z 1 0 q(y 0 ;t)cos(¸ n y 0 )dy 0 (2.30) 19 Comparing the cosine series representation of the injection rate profile with the har- monic terms, then the coefficientsA n can be determined A 0 = ¡ 1 2 a 0 (2.31) A n = a n ¸ n cosh(¸ n ) (2.32) The pressure solution can be expressed in the following form: P(x;y)=¡ 1 2 a 0 (x¡1)+ 1 X n=1 a n ¸ n cosh(¸ n ) sinh(¸ n (1¡x))cos(¸ n y) (2.33) or as the following form: P(x;y)=(1¡x) Z 1 0 q(y 0 ;t)dy 0 + 1 X n=1 2 R 1 0 q(y 0 ;t)cos(¸ n y 0 )dy 0 ¸ n cosh(¸ n ) sinh(¸ n (1¡x))cos(¸ n y) (2.34) By further simplification and rearrangement of the summation and the integral, the final form of the pressure solution is P(x;y)= Z 1 0 G(0;´;x;y)q(´)d´ = Z 1 0 [(1¡x)+ 1 X n=1 2sinh(¸ n (1¡x))cos(¸ n y) ¸ n cosh(¸ n ) cos(¸ n ´)]q(´)d´ (2.35) 20 2.3.2 Front Kinematics From the pressure solution we can express the velocity in thex andy directions analyt- ically and substitute in the front kinematics equation. For a homogeneous medium, the velocity is the pressure gradient along the spatial coordinate @f @t =v x ¡v y @f @y =¡ @P @x + @P @y @f @y (2.36) Using 2.35, we reach the following expressions forv x andv y ; v x = Z 1 0 [1+ 1 X n=1 2cosh(¸ n (1¡x))cos(¸ n y) cosh(¸ n ) cos(¸ n ´)]q(´)d´ (2.37) v y = Z 1 0 [ 1 X n=1 2sinh(¸ n (1¡x))sin(¸ n y) cosh(¸ n ) cos(¸ n ´)]q(´)d´ (2.38) Substituting into equation 2.36 we finally find @f @t = Z 1 0 [1+ 1 X n=1 2cosh(¸ n (1¡x))cos(¸ n y) cosh(¸ n ) cos(¸ n ´)]q(´)d´ ¡ @f @t Z 1 0 [ 1 X n=1 2sinh(¸ n (1¡x))sin(¸ n y) cosh(¸ n ) cos(¸ n ´)]q(´)d´ (2.39) and if we re-arrange we obtain the final equation @f @t = Z 1 0 f1+ 1 X n=1 f 2 cosh(¸ n ) (cosh(¸ n (1¡x))cos(¸ n y) ¡f y sinh(¸ n (1¡x))sin(¸ n y))cos(¸ n ´)ggq(´)d´ (2.40) 21 Integral equation 2.40 is an inverse problem relating the injection rates to the front veloc- ity at the front location. It is a Fredholm Integral Equation of the First Kind and its solution enables us to determine the appropriate injection rate profile according to the prescribed front evolution. In its compact form it is the same as equation 2.10 @f(y;t) @t = Z 1 0 K(0;´;f(y;t);y;t)q(´;t)d´ (2.41) wheref(y;t) is the front location and the kernelK(0;´;f(y;t);y;t) is a function of the front dynamics K(0;´;f(y;t);y;t)=1+ 1 X n=1 f 2 cosh(¸ n ) (cosh(¸ n (1¡x))cos(¸ n y) ¡f y sinh(¸ n (1¡x))sin(¸ n y))cos(¸ n ´)g (2.42) Equation 2.41 must be transformed into a system of equations to make them suitable for numerical evaluation. The discretized form of equation 2.41 is f t (y i ;t)= n X j=1 K(0;´ j ;f(y i ;t);y i ;t)q(´ j ;t)±´ j (2.43) Equations 2.43 describe the non-linear relationship between the front velocity at each point to the rates of all the injection points along the injection well trajectory. The discrete form of this relationship is shown in Figure 2.2. The blue line represents a typical front location at time t where f t (´;t) is the velocity at the point (x;y) on the front. We are interested in obtaining the injection rate profile along the well trajectory, namelyq(´;t) when the desired front location and its velocities are given as a function of time. 22 q(η 1 ,t) q(η 2 ,t) q(η n ,t) K(η 1 ,x,y,t) K(η 2 ,x,y,t) K(η n ,x,y,t) f t (x,y) = * f t (x,y) y−axis: q(η,t) o o o o o o o o o o o o o o No Flow x-axis: No Flow P=0 The non-linear relation between the front velocity and injection rate can be written as a system of equations. o o Figure 2.2: A simple schematic of the front location and its correspondent system of equations relating the velocity at the front to the injection rates along the injection well The next section illustrates the problem and its solution through multiple examples. It also presents the discussion on the necessary regularization to obtain stability. Admis- sible fronts and how far one can control the front effectively through the porous media are discussed as well. 23 2.4 Numerical Results Next, we illustrate the front steering problem by injection rate partitioning through three examples. For each example, we define specific front dynamics. The porous media is homogeneous and the flow is potential in the two-dimensional domain as depicted in Figure 2.1. In all the examples, the injection well is at the left side and flow rate at any point along the injection well is independent of the other injection points. The right hand side is at constant pressure of zero. Top and bottom sides of the medium are no-flow boundaries. 2.4.1 Numerical Procedure Because the front dynamics is known as a function of time, the system of equations 2.43 can be solved at each time instant. We consider a finite set of time intervals and solve the system of equations for each time interval ¢t. In addition, the kernel is the sum of infinite number of harmonic and hyperbolic terms. Each term is calculated analytically and based on our investigation, a fifty or higher number of terms will converge to the solution. In all examples we consider a hundred and twenty terms for the summation purposes. The simple collocation was utilized in all examples and the domain is divided into hundred and twenty cells in each direction. Because the front is exact then there is no error in front dynamics term in the equation. The only error associated with the system of equations is the discretization error associated with the numerical calculations. In obtaining the rate solutions, the Tikhonov regularization technique is applied at each time step. The optimum regularization parameter at each time step was determined as explained in chapter 4. 24 2.4.2 Front/Particle Tracking Given rate solutions corresponding to the controlled front as a function of time, a for- ward solving algorithm is used to track the position of the front. This is done to confirm that the rate solutions can produce the front dynamics. This allows us to validate our solution. We start at time zero, apply the given rate at each time step and solve for the pressure. Then, based on the pressure field, we calculate the velocity components in the x andy directions. The front particles are moved along the normal velocity at the front. In two-dimensional porous media, the front location at timet+±t is x f (t+±t)=x f (t)+v x ±t (2.44) y f (t+±t)=y f (t)+v y ±t (2.45) where ±t is a time interval. For the forward problem, we simulate the dynamics of the fronts by simply tracking the fronts as a function of time. The assumption made is that the fronts are represented by a finite number of particles. This approach is used both in homogeneous and heterogeneous porous media. Given the solution methods and the validation procedure, three front examples are introduced in the next sections, illustrating front steering problems and the related dis- cussion. 25 2.4.3 Example I The first example corresponds to a front evolving according to f(y;t) = 1 2 t(¾¡1)tanh( y¡ 1 2 ² )+ 1 2 t(1+¾) (2.46) f t (y;t) = 1 2 (¾¡1)tanh( y¡ 1 2 ² )+ 1 2 (1+¾) (2.47) f y (y;t) = 1 2² t(¾¡1)f1¡tanh 2 ( y¡ 1 2 ² )g (2.48) where the parameters² and¾ dictate the sharpness of the front. A characteristic example is shown in Figure 2.3, corresponding to²=0:35 and¾ =0:5, respectively. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases Figure 2.3: Desired front location as a function of time for example I The question is to achieve an injection rate profile that will result into the front dynamics shown in Figure 2.3. However, unregulated rate profile solutions applying 26 standard methods are unstable (and meaningless) as shown in Figure 2.4(a). On the other hand, the regularized rate profile solutions are stable as shown in Figure 2.4(b). In the regularized rate solution case, the front is well controlled and with good accu- racy. Indeed, using the obtained rate solution into a forward solving algorithm with front tracking, produces the results shown in Figure 2.5(a). In this example, the comparison between desired and achieved is good. The regularization allowed us to obtain a stable solution for a substantially larger values of time than would otherwise been possible. 2.4.4 Example II The second example is a rather simpler front compared to the first. The mathematical expression for the front and its velocity inx andy directions are f(y;t) = t+0:215tsin(¼y) (2.49) f t (y;t) = 1+0:215sin(¼y) (2.50) f y (y;t) = +0:215tcos(¼y) (2.51) The middle part of the front moves faster with respect to both upper and lower sections. Figure 2.6(a) describes the front movement as time increases. In order to steer the front the regularization method was used to obtain corresponding rate profile along the injection well. The corresponding regularized rate solution is shown in Figure 2.6(b). As the front moves and changes shape, the rate profile has the same variation but with a larger amplitude. Optimum regularization parameters are shown in Figure 2.7(b). At each time step, we applied an approach similar to L-curve until no further improvements in solution were possible (see below). As shown in Figure 2.7(a), the front is well controlled until it reaches half way through the domain. As time increases and the front 27 −3 −2 −1 0 1 2 x 10 14 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves: Standard solution x-axis Rate along injection well (a) Rate solution using Least Squared method 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (b) Rate solution achieved by Tikhonov regularization method Figure 2.4: Solution to Example I: rate solutions obtained by standard method versus rate solutions obtained by regularization method 28 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x-axis y-axis, Injection side time increases -Desired front location: Blue -Front by appling injection rate: Red (a) The front dynamics calculated by the injection rate matches well the desired front dynamics 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 time Optimal regularization parameter (b) The optimum regularization parameter values as time increases Figure 2.5: Validation of the solution to Example I 29 moves further, one notes an increasingly weakening in the control of the front. This is an important result and will be discussed in later sections. 2.4.5 Example III The third front profile has the sinusoidal form f(y;t) = t+®tcos(¯¼y) (2.52) f t (y;t) = 1+®cos(¯¼y) (2.53) f y (y;t) = ¡®¯sin(¯¼y) (2.54) corresponding to a front translating at a y¡dependent velocity. The amplitude ® and the frequency ¯ dictate the form of the front. A characteristic example is shown in Figure 2.8(a), corresponding to® = 0:15 and¯ = 4, respectively. In comparison to the previous example, this front is more complex in its dynamics. In this example, for the time shown, the front is well controlled and with good accu- racy. In fact, using the obtained injection rate solution into a forward solving algorithm with front tracking, produces the result shown in Figure 2.9(a). The comparison between desired and achieved is good although one notes that there is a progressive deterioration as time increases. Again, regularization enabled us to to achieve a stable solution for a substantially larger values of time that would have otherwise been possible. Similar to the previous example, we observe that as the front propagates through the domain, the sensitivity of the front dynamics to the rate changes diminishes. The deterioration of the solution as time increases is a reflection of the extensive cross-flow that develops and will be explained in later sections. 30 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases (a) Desired front location as a function of time 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (b) Rate solution achieved by Tikhonov regularization method Figure 2.6: Solution to Example II 31 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x-axis y-axis, Injection side time increases - Desired front: Blue - Front location by injection rate : Red (a) The front dynamics calculated by the injection rate matches well the desired front dynamics 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 time Optimal regularization parameter (b) The optimum regularization parameter values as time increases. Figure 2.7: Validation of the solution to Example II 32 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases (a) Desired front location as a function of time 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (b) Rate solution achieved by Tikhonov regularization method Figure 2.8: Solution to Example III 33 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x-axis y-axis, Injection side time increases - Desired front : Blue - Front location by injection rate : Red (a) The front progression calculated by the injection rate matches the desired front progression 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 time Optimal regularization parameter (b) The optimum regularization parameter values as time increases Figure 2.9: Validation of the solution to Example III 34 2.5 Discussion 2.5.1 Degree of Ill-conditioning Because the Fredholm integral equation of the first kind is an ill-posed problem, its dis- retized form is a discrete ill-conditioned problem. The solution is sensitive to the small perturbations of the kernel, the non-linear relationship between velocity at the front and rates at the injection points. The ill-conditioning property of the inverse problem pre- vents us to obtain solutions by standard methods such as LU and QR factorization. This does not mean, that a solution can not be obtained. Depending on the degree of ill- conditioning, one can achieve an approximate solution to the underline exact solution. As explained in chapter 4.1, singular value decomposition is a powerful tool to examine the characteristic elements of the coefficient matrix. One such characteristic element is the singular value of the coefficient matrix [Wah80]. The diminishing of the singular values can be used to characterize the degree of ill-conditioning. The behavior of the singular values and their corresponding singular functions are directly connected to the degree of ill-conditioning. The faster the singular value decays to zero, the higher the degree of ill-conditioning. For instance, in example I as time increases the singular value of the kernel dimin- ishes faster as depicted in Figure 2.10(a) for the first twelve time steps. Similar behav- ior is observed for example III shown in Figure 2.10(b), where we observe an even faster decay. This explains the difficulty in steering front in example III, as it has more complex front dynamics in comparison to the dynamics of the front in example I. As the front moves away from the injection side, the singular value diminishes faster and that directly translates into a higher degree of ill-conditioning of the control problem. 35 Because of this behavior prior to any solution, one must have an approximate of the degree of ill-conditioning. It was shown previously that how fast the singular values decay can be utilized to determine the degree of ill-conditioning. However, one may ask for the source of instabilities that adversely affects the degree of ill-conditioning. Math- ematically, the integral operator is the source of the instability as the ill-conditioning is directly connected with the compactness of the operator [Han92c] in inverse problems. The theory of linear compact operators is a generalization of the Fredholm’s the- ory of integral equations in formulating boundary value problems ([HY02] chapter 3) and are very well behaved when applied in forward directions. In the forward problem, the kernel has a smoothing effect on the injection rates in a sense that high frequency components are smoothed out by the integration. In inverse problems the opposite hap- pens, where high frequency components get amplified and cause instabilities during the inversion (differentiation) process. This translates as the diminishing of the sensitivity of the front to the injection rate for the problem, as implied in the Figure 2.10(a) and Figure 2.10(b), with an increasing slope of decay in singular value as time increases. In addition, in discretized form, the system of equations is affected by discretization errors, measurement errors and rounding errors. 2.5.2 Optimum Regularization Parameter A key element in obtaining a proper solution to an ill-conditioned system of equations is the choice of an optimized regularization parameter. It enables us to recover information associated with the larger singular values and filter out the high frequency components 36 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time increaases (a) Decay in singular values of the coefficient matrix of Example I 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time increases (b) Decay in singular values of the coefficient matrix of Example III Figure 2.10: How fast singular values decay, determines the degree of ill-conditioning of the system 37 associated with the smaller singular values. Here, we examined the L-curve and quasi- optimality criterion for choosing the regularization parameter at each time step. In addi- tion, we determined manually the optimum regularization parameter for each time step for the front of example I and III. We observed that, one can use one parameter value and achieve very similar solutions for rate profiles. The L-curve is a plot of the side con- straintlogkL(x¡x o )k versus the error normlogkAx¡bk for the Tikhonov and similar regularization methods. With its definition, the value of the parameter corresponding to the maximum curvature represents the near optimum value for the regularization param- eter [Han92b]. Figure 2.11 depicts the L-curve for choosing the optimum regularization parameter at the 10 th time step corresponding to the front of example I. An interesting aspect of the L-curve is that, one can visually see the balance between regularization error and perturbation errors. However, for the front steering problem the L-curve method and its corner underestimate the optimum regularization parameter value resulting in less filtering effect in respect to the quasi-optimality criterion and manually determined optimum values. Another parameter choice method that looks for a balance between the regulariza- tion error and the perturbation error is the quasi-optimality criterion. Morozov [Mor84] derived the quasi-optimality criterion as the minimum distance between two successive regularized solutions Q(¸)=k¸ dx ¸ d ¸ k 2 (2.55) where one can reach a minimum where the solutionx ¸ changes from being dominated by regularization errors to be dominated by perturbation errors. One problem with quasi- optimality criterion method is that the function to be optimized may have several local 38 10 −6 10 −4 10 −2 10 0 10 2 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 Residual norm || A x − b || 2 Side constraint λ|| x−x o || 2 L-curve- Example one time step 10 The λ of 0.0035 corresponding to the corner is the optimum reg. parameter Figure 2.11: The L-curve of front in example I at the10 th time step minimum [Flo02], therefore it might be difficult to reach the global minimum. The quasi-optimality criterion was used to determine the optimum regularization parame- ter in example I. For the front in example I, Figure 2.12 depicts the desired front and three different regularization parameters used to obtain the solution. As shown in Figure 2.13, the rate solutions obtained are very similar. Still one notes that the case where¸ is chosen by quasi-optimality criterion, the front control is accurate. Also there is a minor difference in rate profiles of the three regularization approaches. All three approxima- tions of the underline solution are acceptable as the solution is not unique. Similarly, Figure 2.14 shows the desired front in example III and the three regularization parame- ter choices. The front control is the same for each of regularization choices even though the rate solutions are different. We see that in example III, shown in Figure 2.15 the rate solutions can be different but with the same accuracy in front control. There are two reasons for that: one is the complexity of the third front compared to the front in 39 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases (a) Desired front for example I 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 time step Optimum regularization parameter λ λ by quasi−optimal criterion λ obtained manually λ of 0.1 (b) Three Optimum¸ being used Figure 2.12: Three choices of the optimum regularization parameters for example I example I and the other is the faster increase in the degree of ill-conditioning (to be shown in the next section) related to the vertical transverse flow effect. 2.5.3 Admissible Fronts The issue of controllability and admissibility of the front steering problem considered in this study is discussed through the front of example I. We have shown that the front control problem is described by f t = Z 1 0 K(0;´;f(y;t);y;t)q(´;t)d´ =Lq(´;t) (2.56) whereL is the integral operator operating on the injection rate profileq(´;t). To get the rate solution of the inverse problem 2.56, a necessary condition is that the operatorL is invertible q(´;t)=L ¡1 (f t ) (2.57) 40 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (a) Rate solution corresponding to ¸ by the quasi- optimum criterion 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x-axis y-axis, Injection side time increases -Desired front location: Blue -Front by appling injection rate: Red (b) Front control performance corresponding to ¸ by the quasi-optimum criterion 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (c) Rate solution corresponding to ¸ determined manually by trial and error 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front (blue) vs front obtained by applying calculated rate(red) x-axis y-axis, Injection side time increases (d) Front control performance corresponding to ¸ determined manually by trial and error 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (e) Rate solution corresponding to¸ of 0.1 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front(blue) vs front obtained by applying calculated rate(red) x-axis y-axis, Injection side time increases (f) Front control performance corresponding to¸ of 0.1 Figure 2.13: The effect of different regularization parameters on the solution 41 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases (a) Desired front for example III 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 time step Optimum regularization parameter λ λ by quasi−optimal criterion λ obtained manually λ of 0.1 (b) Three Optimum¸ being used Figure 2.14: Three choices of the optimum regularization parameter for example III Admissible fronts are ones for which L ¡1 (f t ) > 0 and R 1 0 L ¡1 (f t ) = constant. The ability to control the front also depends on the complexity of the front. By complexity, we mean the sharpness of the front velocity. Let us consider the front of example I again where its front velocity is f t (y;t)= 1 2 (¾¡1)tanh( y¡ 1 2 ² )+ 1 2 (1+¾) (2.58) Figure 2.16(a) and Figure 2.16(b) show front velocity profiles of example I with different values of ² and ¾. Figure 2.16(a) shows the effect of varying ² on the sharpness of the front. Figure 2.16(b) shows the effect of varying¾ on the front shape and its sharpness. Here, smaller values of¾ lead to velocity values at the lower part which are faster than in the upper part. A combination of these two parameters determines the specific front profile and its evolution. In both figures, the front of example I controlled previously for ² = 0:35 and ¾ = 0:5 is shown for comparison purposes. Figure 2.17 depicts the four desired front evolution scenarios with²=0:1 and¾ =0:1,²=0:1 and¾ =0:5,²=0:5 42 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (a) Rate solution corresponding to ¸ by the quasi- optimum criterion 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front(blue) vs front obtained by applying calculated rates(red) x-axis y-axis, Injection side time increases (b) Front control performance corresponding to ¸ by the quasi-optimum criterion 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (c) Rate solution corresponding to ¸ determined manually by trial and error 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x-axis y-axis, Injection side time increases - Desired front : Blue - Front location by injection rate : Red (d) Front control performance corresponding to ¸ determined manually by trial and error 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (e) Rate solution corresponding to¸ of 0.1 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front(blue) vs front obtained by applying calculated rates(red) x-axis y-axis, Injection side time increases (f) Front control performance corresponding to¸ of 0.1 Figure 2.15: The effect of different regularization parameters on the solutions 43 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X-axis Y-axis; Injection side The front velocity profiles in example one with different ² and σ =0.1 ε = 0.1 ε = 0.2 ε = 0.3 ε = 0.4 ε = 0.5 ε = 0.35, σ = 0.5 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X-axis Y-axis; Injection side The front velocity profiles in example one with different σ and ²=0.1 σ = 0.1 σ = 0.2 σ = 0.3 σ = 0.4 σ = 0.5 ε = 0.35, σ = 0.5 (b) Figure 2.16: Different values of ¾ and ² defines fronts with different sharpness. The front for example I used previously is also shown on both figures., (a) effect of varying ².,(b) effect of varying¾ 44 and ¾ = 0:1, and ² = 0:5 and ¾ = 0:5 that represents different level of complexity in front dynamics as time increases. Figure 2.18 shows the corresponding rate profiles obtained after applying regularization. For the times shown, we observe that the front is well controlled and with accuracy for the front with ² = 0:5 and ¾ = 0:5, while the front is not well controlled for the front defined by ² = 0:1 and ¾ = 0:1. Using the obtained rate solution in a forward front tracking algorithm produces the results in Figure 2.19. The comparison between desired and achieved for the four fronts shows that the controllability of the front depends on the complexity of the front dynamics. 45 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics: ²=0.1 and σ =0.1 X-axis Y-axis, Injection well time increases (a) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics: ²=0.1 and σ =0.5 X-axis Y-axis, Injection well time increases (b) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics: ²=0.5 and σ =0.1 X-axis Y-axis, Injection well time increases (c) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics: ²=0.5 and σ =0.5 X-axis Y-axis, Injection well time increases (d) Figure 2.17: The desired front dynamics of example I with four different characteristic forms corresponding to different values of² and¾ 46 −2 −1 0 1 2 3 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves: ²= 0.1 and σ = 0.1 X-axis Rate along injection well (a) 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile : ²=0.1 and σ =0.5 X-axis Rate along injection well (b) 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile : ²=0.5 and σ =0.1 X-axis Rate along injection well (c) 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile : ²=0.5 and σ =0.5 X-axis Rate along injection well (d) Figure 2.18: The rate solutions of example I with four different characteristic forms corresponding to different values of² and¾ 47 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red) : ² = 0.1 and σ = 0.1 x-axis y-axis, Injection side time increases (a) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red): ² = 0.1 and σ = 0.5 X-axis Y-axis, Injection side time increases (b) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red): ²= 0.5 and σ = 0.1 X-axis Y-axis, Injection side time increases (c) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red): ² = 0.5 and σ = 0.5 X-axis Y-axis, Injection side time increases (d) Figure 2.19: The comparison between desired and achieved front profiles of example I with four different characteristic forms corresponding to different values of² and¾ 48 2.6 Summary We considered the problem of control of a displacement front in porous medium, via flow-rate partition in a well. With the assumption of potential flow and that the displace- ment is at unit mobility ratio, we derived a formalism that allows to control displacement fronts of pre-determined dynamics using rate control. The formulation derived in this section results into the solution of an integral equation, which is in an analytical form, in a rectangular geometry and a homogeneous porous medium. For the solution of the inte- gral equation, a regularization technique is necessary. However it is found that even with the use of the regularization numerical instabilities do develop at larger times. We illus- trated the front control approach through multiple examples and showed that the front control becomes more difficult as the dynamics of the front gets more complex. It was found that regularization parameter is an important element of the solution technique. Multiple approaches of choosing optimum regularization parameters were examined. In addition, as the front moves away from injection side, there is an increase in instabilities, and that is translated into increasing transverse cross-flow effect. 49 Chapter 3 Control of Displacement Fronts in Heterogeneous Media 3.1 Introduction In this chapter we consider extension of the front control problem in a two-dimensional, heterogeneous porous medium. It is assumed that gravity and dispersion effects are absent and that we have independent injection points at the left boundary. The problem considered is the same to that discussed in the previous chapter except that the displace- ment takes place in a heterogeneous porous medium. Similar to the homogeneous case, the problem will be formulated as a Fredholm Integral Equation of the First Kind with injection rates as control parameters. The rate solution will be obtained following the same procedure as described previously. Two objectives will be pursued. First, in connection to the results from the homo- geneous case, the controllability of the front for larger times and transverse flow effects will be investigated for anisotropic porous media. Then we will apply the front control to the general heterogeneous porous media. 50 3.2 Mathematical Model The formulation approach is similar to potential flow in homogeneous medium case. The velocity in the x and y directions at the front can be represented by superimpos- ing the velocity effect of each injection point along the injection side of the boundary. This was achieved analytically in homogeneous case, however they must be determined numerically for a heterogeneous porous medium. Because of the linearity of the problem, the velocity terms v x and v y can be repre- sented as the superposition of an infinite number of velocity corresponding to the indi- vidual injection point along the injection side. One can apply a similar superposition technique, as elaborated in Sudaryanto and Yortsos [SY01] to determine the velocities at front numerically. The idea is to express the displacement as the superposition of the response of individual well segments of point source strength, along the injection trajectory. Therefore, one can replace the velocity termsv x andv y in equation 3.2 with v x = N wi X i=1 ® i v i x (3.1) v y = N wi X i=1 ® i v i y (3.2) Here, ® i is the correspondent injection rate for each doublet (along each injection seg- ment), and N wi is the number of well segments along the injection well (theoretically infinite, but in practice finite). At the production well, segments are represented as one well by imposing the same Dritchlet boundary condition for all segments. The velocities v i x andv i y must be determined numerically. Because the fluids are equal viscosity and are 51 incompressible, velocities v i x and v i y need to be determined only once. By substituting velocities into for each time step, we arrive at @f @t = N wi X i=1 ® i v i x ¡ N wi X i=1 ® i v i y @f @y = N wi X i=1 (v i x ¡ @f @y v i y )® i = N wi X i=1 K j;i ® i (3.3) where K j;i is the kernel relating the injection rate profile to the front dynamics. The solution of 3.3 gives the rate partition along the injection well as time increases corre- sponding to the desired front progression. Because it is the discrete form of the Fred- holm Integral Equation of the First Kind, regularization is necessary in obtaining the rate solution. 3.2.1 Two-Well Problem The displacement can be described in terms of the dynamics of the displacement fronts where the effect of each individual injection source can be superimposed as graphi- cally illustrated in Figure 3.1. The velocity response of each individual injection point must be determined numerically for heterogeneous porous medium. We present a brief description of the method used by Sudaryanto and Yortsos [SY01]. 52 Figure 3.1: Super-imposition of multiple sources in potential flow The velocity termsv i x andv i y must be determined by the Darcy’s law v i x (x;y) = ¡k(x;y) @p i @x (3.4) v i y (x;y) = ¡k(x;y) @p i @y (3.5) where k is the permeability of the porous medium and p i is the pressure field that cor- responds to the two well problem, with i segment of the injection side injects at +1 and the production side fixed at constant pressure zero (which will produce¡1 for the incompressible problem at hand). For the velocities, one must calculate the pressure fieldp i by solving @ @x (k(x;y) @p i @x )+ @ @y (k(x;y) @p i @y )=0 (3.6) subject to the appropriate boundary conditions. 53 3.3 Numerical Solution 3.3.1 Two-Well Problem In general the governing equation 3.6 must be solved numerically, using finite difference method. The basis for finite difference method is the construction of a discrete grid where the continuous derivative terms in the partial differential equation are replaced by equivalent finite difference terms [Fle91]. In this study, we apply a block-centered grid system with constant and equal cell sizes in x and y directions. The pressure solution are at the center of the grid blocks. The spatial discretization of the equation 3.6 is ( kA x 4x ) i+ 1 2 ;j (p i i+1;j ¡p i i;j )¡( kA x 4x ) i¡ 1 2 ;j (p i i;j ¡p i i¡1;j )+ ( kA y 4y ) i;j+ 1 2 (p i i;j+1 ¡p i i;j )¡( kA y 4y ) i;j¡ 1 2 (p i i;j ¡p i i;j¡1 )=¡q i i;j (3.7) whereA x andA y are the cross sectional area of grid interface in thex andy directions. Here,q i i;j is the injection rate at grid block(i;j) and is equal to+1 if the block includes the injection segmenti and zero otherwise. The point(i;j) is the center of the grid block. And points(i§ 1 2 ;j) and(i;j§ 1 2 ) represents the interface boundary of two adjacent grid blocks ini andj directions. The above discrete system of equations allow us to obtain the pressure values at the center of each cells. The values of some properties need to be evaluated at the boundaries between adjacent blocks. To do so, we use harmonic averaging to evaluate( kA x 4x ) i§ 1 2 ;j and( kAy 4y ) i;j§ 1 2 , as shown ( kA x 4x ) i§ 1 2 ;j = 2A x i;j k i;j A x i§1;j k i§1;j 4x(A x i;j k i;j +A x i§1;j k i§1;j ) (3.8) 54 ( kA y 4y ) i;j§ 1 2 = 2A y i;j k i;j A y i;j§1 k i;j§1 4y(A y i;j k i;j +A y i;j§1 k i;j§1 ) (3.9) With some rearrangement, equation 3.7 can be written in the following form L i;j p i i¡1;j +R i;j p i i+1;j +C i;j p i i;j +U i;j p i i;j+1 +D i;j p i i;j¡1 =¡q i i;j (3.10) where L i;j = ( kA x 4x ) i¡ 1 2 ;j ; R i;j =( kA x 4x ) i+ 1 2 ;j ; U i;j =( kA y 4y ) i;j+ 1 2 ; D i;j =( kA y 4y ) i;j¡ 1 2 C i;j = ¡(L i;j +R i;j +U i;j +D i;j ) Equation 3.10 is written for every grid block. This produces N number of equations withN unknowns for pressure and can be presented in terms of a pentadiagonal matrix. Standard methods in linear algebra such as LU factorization can be used to solve for the unknown pressure field. After calculating the pressure field, Darcy’s velocity terms can be obtained at discrete points v i x i§ 1 2 ;j = k i§ 1 2 ;j ( p i i;j ¡p i i§1;j 4x ) (3.11) v i y i;j§ 1 2 = k i;j§ 1 2 ( p i i;j ¡p i i;j§1 4y ) (3.12) Since the pressure field of the system does not vary with time, then the corresponding velocity field also does not vary and only needs to be determined once for the injection segmenti. 55 3.3.2 Interpolation Methods By the finite difference method, the pressure and velocity solutions are only available at discrete points. The solutions are known at the center of the cells and the velocities are calculated on the interfaces between the cells. Therefore, the velocities at any points (x;y) in the equation 3.3 must be interpolated from the known values at discrete points. The multi-linear interpolation method presented below will be used to determine the Darcy’s velocities at a point(x;y) on the desired front by the four nearest known values of the velocities. A simple schematic of the multi-linear interpolation at unknown points, where its surrounding values are given, is shown in Figure 3.2. Given values at discrete x y f(x,y) Known values of f(x,y) Interpolated point Figure 3.2: Schematic of interpolation process points, the following multi-linear interpolation can be used to calculate the unknown values of the function at any point(x;y) 56 f(x;y)=(1¡h x )(1¡h y )f(x o ;y o )+h x (1¡h y )f(x 1 ;y o ) +(1¡h x )h y f(x o ;y 1 )+h x h y f(x 1 ;y 1 ) (3.13) where h x = x¡x o x 1 ¡xo , h y = y¡y o y 1 ¡yo , x o · x · x 1 and y o · y · y 1 . For the present problem, the block-centered grid system with equal grid sizes are used and the pressures are being determined at grid block centers. Then having the pressure field, the velocity values at the cell boundaries are calculated through Darcy’s relation. The velocity v i x andv i y at an arbitrary point (x,y) on the front or within the domain are calculated from the following relations: v i x (x;y)=(1¡h 1x )(1¡h 1y )v i x i¡ 1 2 ;j¡1 +h 1x (1¡h 1y )v i x i+ 1 2 ;j¡1 +(1¡h 1x )h 1y v i x i¡ 1 2 ;j +h 1x h 1y v i x i+ 1 2 ;j (3.14) v i x (x;y)=(1¡h 2x )(1¡h 2y )v i y i¡1;j¡ 1 2 +h 2x (1¡h 2y )v i y i;j¡ 1 2 +(1¡h 2x )h 2y v i y i¡1;j+ 1 2 +h 2x h 2y v i y i;j+ 1 2 (3.15) where h 1x = x¡(i¡1)4x 4x ; (i¡0:5)4x·x·(i+0:5); i=1;2;3;::: h 1y = y¡(j¡0:5)4y 4y ; (j¡1)4y·y·(j+1); j =1;2;3;::: h 2x = x¡(i¡0:5)4x 4x ; (i¡1)4x·x·(i+1); i=1;2;3;::: h 2y = y¡(j¡1)4y 4y ; (j¡0:5)4y·y·(j+0:5); j =1;2;3;::: 57 3.3.3 Front/Particle Tracking The computational procedure of front tracking is similar to the method described pre- viously for front control in the homogeneous system. To validate the accuracy of the solution, the obtained rate solutions corresponding to the controlled front as a function of time, will be used in the forward solving algorithm to track the position of the front. Since the velocities v i x and v i y in x and y directions are known from two well solution, we use the equation 2.44 and 2.45 to advance the front (particles) as time progresses. 3.4 Numerical Results The general problem considered is the same as in the previous chapter. Various scenarios of heterogeneous reservoirs will be used in the numerical experiments. These include the special cases of anisotropic reservoir, layered reservoir and reservoirs with specified heterogeneity, such as, decreasing permeability asy increases and a permeability streak at the middle of the reservoir parallel to thex axis. In addition we apply front steering techniques to more general heterogeneous medium generated by Sequential Gaussian Simulation methods. By conducting numerical experiments for different types of heterogeneity, we may investigate the effect of heterogeneity on the front control. In addition we can study the controllability of the front and intend to provide results regarding how far one can control the front away from injection well. This is a fundamental question in displace- ment processes where injection rate partitioning is the only available factor in achieving optimum displacement efficiency. 58 First, we will reproduce the front control problem of example I. This will assure that the numerical approach is accurate and can lead us to the same result as the analytical approach. In the numerical approach, we discretized the rectangular and homogeneous domain into71£71 grid blocks. The upper and lower boundaries are no-flow, and the right side boundary is at pressure zero. The injection wells are located at the first block to the right side of the left boundary as graphically shown in Figure 3.1. The results of the numerical approach for example I is depicted in Figure 3.3. The results obtained by the numerical approach are as accurate as those by the analytical method. In both methods, the quasi-optimal criterion approach was used in determining the optimum regularization parameter for each time step. 3.4.1 Anisotropic Porous Media The anisotropic porous medium is a special case where the permeability is direction dependent. For our numerical experiments, we consider examples I and III where we conduct the displacement for a larger period of time. We introduce an anisotropy with transverse-to flow-direction permeability ratios of 0:1, 0:01 and 0:0001 and examine front control in both examples. The transverse direction is they direction, and the flow direction is the x direction. Figure 3.4 shows the desired front progression in time for the two examples. The previous chapter showed that as the front moves away from the injection well, the ability to control the front diminishes as time increases. It was postulated, that one reason for the instability was the development of extensive cross-flow effect as the front moves further away from the injection points. 59 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regulariza- tion method by analytical method 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile variation as front evolves: By numerical method x-axis Rate along injection well (b) Rate solution achieved by Tikhonov regulariza- tion method by numerical method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x-axis y-axis, Injection side time increases -Desired front location: Blue -Front by appling injection rate: Red (c) Achieved front vs desired one by analytical method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved (red) : Example one by numerical method x-axis y-axis, Injection side time increases (d) Achieved front vs desired one by numerical method 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 time Optimal regularization parameter (e) The optimum regularization parameter as time increases by analytical method. 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 The optimized parameter at each time step time Optimal regularization parameter (f) The optimum regularization parameter as time increases by numerical method. Figure 3.3: Validation of the numerical method vs. analytical one 60 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases (a) Desired front progression in time: example I 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics x-axis y-axis, Injection well time increases (b) Desired front progression in time: example III Figure 3.4: The desired front progression in time in examples I and III at larger time intervals In early times (when the front is very close to injection points), cross-flow is limited and the problem behaves as in a non-communicating layered systems. In fact the early- time solution for the rate shown in Figure 3.5(a) is identical to that corresponding to a layered system, which for the front in example I has the simple form q(y;0)= 1 2 (¾¡1)tanh( y¡ 1 2 ² )+ 1 2 (1+¾) (3.16) The Figure shows that as time increases, the rate profile has the same variation but with a larger amplitude. At the same time, at the larger times the intensity of the cross-flow increases and the sensitivity of the front to the rate diminishes, causing instability. Introducing anisotropy where the permeability in the y direction is lower that the permeability in the x direction by a magnitude of 0:1 or 0:01, reduces the effect of the extensive cross-flow and diminishes the instability as shown in Figure 3.6(b) and Fig- ure 3.7(b) respectively. As the system approaches conditions of a non-communicating layer, the flow rate profile becomes closer to equation 3.16 and control can be extended 61 at much later times. This behavior is more pronounced in the more extreme case shown in Figure 3.8(a) and Figure 3.8(b), corresponding to the anisotropic medium with trans- verse to flow direction permeability ratio of0:0001. The rate profile becomes indepen- dent of time and it is identical to equation 3.16. The second example corresponds to the front for example III corresponding to a front translating at a y-dependent velocity. The desired evolution of the front at different times is shown on the right side of Figure 3.4. The corresponding rate solution obtained after applying regularization is shown in Figure 3.9(a) for the homogeneous porous medium. Again, at early times the front is accurately controlled. As time increases, one notes that there is a progressive deterioration in front control. The comparison between desired and achieved front is shown in Figure 3.9(b). The deterioration of the solution as time increases is a reflection of the cross-flow that develops. In early times (when the front is very close to the injection points), the cross-flow is insignificant and the problem behaves as in a non-communicating layered systems. In fact the early-time solution for the rate shown in Figure 3.9(a) is identical to that corresponding to a layered system, which for the front in example two has the simple form q(y;0)=1+®cos(¯¼y) (3.17) As time increases, the rate profile has the same shape but with a larger variation in its values as shown in the Figure 3.9(a). At the same time, at the larger times the intensity of the cross-flow increases and the sensitivity of the front to the rate diminishes, causing instability. 62 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Injection rate profile: Homogeneous medium x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red) : Homogeneous medium x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.5: Homogeneous porous media for example I (a) the corresponding rate profile solution. , (b) the desired front progression versus the achieved one 63 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile : Anisotropic medium ky kx =0.1 x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red): Anisotropic medium ky kx = 0.1 x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.6: Anisotropic system with Ky K x = 0:1 for example I (a) the corresponding rate profile solution , (b) the desired front progression versus the achieved one 64 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile : Anisotropic medium ky kx =0.01 x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red) : Anisotropic medium ky kx = 0.01 x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.7: Anisotropic system with Ky K x =0:01 for example I (a) the corresponding rate profile solution , (b) the desired front progression versus the achieved one 65 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile: Anisotropic medium ky kx =0.0001 x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red): Anisotropic medium ky kx = 0.0001 x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.8: Anisotropic system with Ky K x = 0:0001 for example I (a) the corresponding rate profile solution., (b) the desired front progression versus the achieved one 66 By decreasing the transverse to flow direction permeability ratio, we can limit the cross-flow effect and the front can be controlled for a larger values of time. As the sys- tem approaches conditions of a non-communicating layer, the flow rate profile becomes closer to equation 3.17 and control can be extended to a larger time interval. This is shown in Figure 3.10(b) and Figure 3.11(b) for transverse to flow direction permeabil- ity ratios of 0:1 and 0:01. This behavior is more pronounced in the more extreme case shown in Figure 3.12(a) and Figure 3.12(b), corresponding to the anisotropic medium with transverse to flow direction permeability ratio of0:0001. The rate profile becomes independent of time and is identical to equation 3.17. 67 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile as front evolves: Homogeneous medium x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red) : Homogeneous medium x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.9: Homogeneous porous media (a) the corresponding rate profile solution. , (b) the desired front progression versus the achieved one 68 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile : Anisotropic medium Ky Kx =0.1 x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red) : Anisotropic medium Ky Kx = 0.1 x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.10: Anisotropic system with Ky K x = 0:1 for example III (a) the corresponding rate profile solution. , (b) the desired front progression versus the achieved one 69 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile: Anisotropic medum ky kx =0.01 x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red) : Anisotropic medium ky kx = 0.01 x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.11: Anisotropic system with Ky K x = 0:01 for example III (a) the corresponding rate profile solution . , (b) the desired front progression versus the achieved one 70 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile: Anisotropic medium ky kx =0.0001 x-axis Rate along injection well (a) Rate solution achieved by Tikhonov regularization method 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs achieved front (red) : Anisotropic medium ky kx = 0.0001 x-axis y-axis, Injection side time increases (b) The front progression achieved by using obtained rate solution Figure 3.12: Anisotropic system with Ky K x =0:0001 for example III (a) the corresponding rate profile solution . , (b) the desired front progression versus the achieved one 71 3.4.2 Heterogeneous Porous Media In this section, by conducting numerical experiments in various heterogeneous reser- voirs, the effect of heterogeneity on the control of front by injection rate partitioning is presented. It is assumed that the heterogeneity is known and the vertical permeability is one tenth of the horizontal permeability. The later assumption delays the cross-flow effect as the front moves away from injection side and allows us to investigate the front control for larger values of times. This assumption is fairly correct and is accepted in the representation of heterogeneous reservoirs in real field applications. In general, reservoir heterogeneity in particular high permeability zones adversely affect the displacement efficiency causing poor sweep and early water breakthrough. We first conduct numerical experiments where heterogeneity is smooth but represents channels as well as cases with increasing permeability in the y direction. We apply the front control technique and study the applicability of the technique in steering the front as well as its applicability in achieving optimized efficiency. Then we test the front steering technique to fields with less smooth heterogeneity generated by Sequential Gaussian Simulations. To assess the effectiveness of our technique, we compare results obtained by front control with two other injection rate policies, an equal injection rate policy, and one in which rates are inversely proportional to the permeability at the injection points along the injection well. For the forward displacement problem, we applied front tracking and moved the particles defining the front in time. The amount of injection at each time step is the same for all three injection scenarios. 72 Heterogeneity in Vertical Direction At first we consider two numerical experiments where the permeability variation in the y direction is smooth and there is no change in permeability values in the x direction. The examples include a decreasing permeability by increasing depth and the existence of a channel, respectively. For the first example, let us consider the permeability field shown in Figure 3.13(a) where permeability is increasing in they direction according to K x (x;y)= 1 2 +exp(µy) (3.18) with µ = 0:06 and K y (x;y) = 0:1K x (x;y). The desired front is the same as example I but with values of ² = 0:65 and ¾ = 0:65, as shown in Figure 3.13(b). The front is moving faster at the bottom and slower at the top although it is less sharper compared to the ones used previously. First, we consider the equal injection rate scenario and inject equally at all injection points along the injection well and track the front in a forward simulation. One expects to observe early water breakthrough at the top of the production side as the permeability values are larger. This is shown in Figure 3.14(a). As expected, the front reaches the production side at the top while at the lower permeability section the front is still less than half the distance away from the injection points. In order to avoid early breakthrough and enforce a more uniform displacement, one common practice is to inject at a rate inversely proportional to the permeability at the injection side, namely injecting at higher rate in lower permeability sections and at a lower rate at higher per- meability areas of the perforations along the injection trajectory. Figure 3.14(b) shows 73 front propagation in a permeability field described in Figure 3.13(a) corresponding to injection rates inversely proportional to the permeability. Then, we applied the front control technique and obtained rate profiles correspond- ing to the desired front propagation as depicted in Figure 3.13(b). The obtained rate solution is used in a front tracking algorithm and produced the results shown in Fig- ure 3.15(b). Again, regularization is found to be necessary for a stable rate solution to be obtained in this and all remaining front control numerical experiments. The compar- ison between desired and achieved is good although there is progressive deterioration as time increases. One notes that the front is steered opposite to the front propagation corresponds to the equal injection rate scenario. The ability to steer the front gives us the advantage of both delay in breakthrough and control over the direction one desires to direct the front. For the second experiment, let us consider the permeability field shown in Fig- ure 3.16(a), where permeability is defined by K x (x;y)= 1 2 +ysin(¼y) (3.19) with permeability in the y direction K y (x;y) = 0:1K x (x;y). The desired front is the same as above with ² = 0:65 and ¾ = 0:65 and shown in Figure 3.16(b). We first consider the forward problem with equal injection rate policy. One expects to see a faster front propagation in the areas of higher permeability and an early breakthrough. This is indeed shown in Figure 3.17(a). As expected, the front arrives at the production side through the higher permeability streak perpendicular to the well trajectories, while at the lower permeability section the front is behind and at the time of breakthrough the fluid in-place is bypassed. 74 1 11 21 31 41 51 61 1 11 21 31 41 51 61 Permeability field K x (x,y)= 1 2 +(exp(0.06∗y)) and K y =0.1K x X-axis Y-axis 5 10 15 20 25 30 35 (a) Permeability increases iny direction 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics X-axis Y-axis, Injection well time increases (b) The Desired front propagation Figure 3.13: Permeability field increasing iny direction and the desired front dynamics 75 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X-axis Y-axis, Injection side Front propagation – Equal injection scenario time increases (a) The front propagation - equal rate injection policy 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X-axis Y-axis, Injection side Front propagation – Rates inversely proportional to permeability time increases (b) The front propagation - permeability proportional rate injection policy Figure 3.14: Front propagation corresponds to the permeability field increasing in y direction showed in Figure 3.13(a)., (a) for the equal injection rate policy .,(b) and for rates inversely proportional to the permeability at the injection points 76 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile X-axis Rate along injection well (a) Obtained rate solution 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red) X-axis Y-axis, Injection side time increases (b) Desired front versus achieved front Figure 3.15: The obtained rate and the comparison between desired and achieved front 77 To avoid early breakthrough and achieve a more uniform displacement, one common practice is to inject at a rate inversely proportional to the permeability at each injection point. By injecting at higher rate in lower permeability sections and lower rate at higher permeability points along the injection trajectory, the front will be propagate at the same speed at all locations, in the absence of cross-flow Figure 3.17(b) shows front propa- gation corresponding to the injection rates inversely proportional to the permeability. The policy of injection rate being inversely proportional to the permeability field at the perforation improves the displacement efficiency compared to the equal injection rate policy, although one notes a tendency of the front to deviate from uniform propagation at larger times. As the front moves away from injection, extensive cross-flow effects cause the displacement front to move toward higher permeability regions even though it is not significant in this example. In parallel, we applied the front control technique and determined rate profiles cor- responding to the desired front propagation shown in Figure 3.16(b). The obtained rate solution is used in a front tracking algorithm to produce the results shown in Fig- ure 3.18(b). The comparison between desired and achieved is good although there is a progressive deterioration as time increases . One notes that the front is steered oppo- site to the front propagation corresponds to the equal injection rate scenario. The front steering technique applied here gives us the advantage of both delay in breakthrough and control over the direction where one desires to steer the front. Heterogeneity in Both Directions In the above two numerical experiments, the permeability varies only in they direction, and is constant in thex direction (layered system). 78 1 11 21 31 41 51 61 1 11 21 31 41 51 61 Permeability field K x (x,y)= 1 2 +y ∗sin(π∗y) X−axis Y−axis 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 (a) Permeability increases iny direction 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics X-axis Y-axis, Injection well time increases (b) The Desired front dynamics Figure 3.16: Permeability field and the desired front dynamics 79 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X-axis Y-axis, Injection side Front propagation – Equal injection scenario time increaes (a) The front propagation - equal rate injection policy 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X-axis Y-axis, Injection side Front propagation – Rates inversely proportional to permeability time increases (b) The front propagation - permeability proportional rate injection policy Figure 3.17: Front propagation corresponds to the permeability field increasing in y direction showed in Figure 3.16(a)., (a) for the equal injection rate policy .,(b) and for rates inversely proportional to the permeability at the injection points 80 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile X-axis Rate along injection well (a) Obtained rate solution 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red) X−axis Y−axis, Injection side time increases (b) Desired front versus achieved front Figure 3.18: The obtained rate solution and the comparison between desired and achieved front 81 The following two examples will consider permeability fields representing hetero- geneity in both directions. In both examples the permeability fields are generated by Sequential Gaussian Simulation. The injection points are located at the left side and production points at the right side with pressure at zero. The desired front shown in Fig- ure 3.19(b) is the same as the one used above with ² = 0:65 and¾ = 0:65. Again, we first consider the forward displacement problem where we allocate equal injection rate at all injection points. One expects to see a faster front propagation in the areas of higher permeability and an early breakthrough. This is shown in Figure 3.20(a). As expected, the front progresses toward the production side through the higher permeability regions, while at the lower permeability regions, the front is moving at a slower pace. One common practice of avoiding early breakthrough is to inject at a rate inversely proportional to the permeability at each injection point. This injection rate policy is applied here, but as shown in Figure 3.20(b), the front propagates very similar to the case of the equal injection rate policy. In fact, the displacement front for the times shown reveals that the rates inversely proportional to permeability is not beneficiary over the equal injection rate policy in delaying the early breakthrough. Because the permeability varies in bothx andy directions as shown in Figure 3.19(a), rates inversely proportional to the permeability values at the sand face do not represent the rates corresponding a uniform front propagation. We applied the previous front control technique to obtain rate profiles correspond- ing to the desired front propagation shown in Figure 3.19(b). The obtained rate solution is used in a forward front tracking algorithm and produced the results shown in Fig- ure 3.21(b). The comparison between desired and achieved is good for the times shown. 82 10 20 30 40 50 60 10 20 30 40 50 60 Permeability field X-axis Y-axis, Injection side 10 20 30 40 50 60 70 80 90 100 (a) Permeability varies in bothx andy directions 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics X-axis Y-axis, Injection well time increases (b) The Desired front dynamics Figure 3.19: Permeability field and the desired front dynamics 83 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front propagation – Equal injection rate scenario X-axis Y-axis, Injection side time increases (a) The front propagation - equal rate injection policy 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front propagation – Rates inversely proportional to the permeability X-axis Y-axis, Injection side time increases (b) The front propagation - permeability proportional rate injection policy Figure 3.20: Front propagation corresponds to the permeability field increasing in y direction and showed in Figure 3.19(a)., (a) for the equal injection rate policy .,(b) for rates inversely proportional to the permeability at the injection points 84 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile X-axis Rate along injection well (a) Obtained rate solution 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red) X-axis Y-axis, Injection side time increases (b) Desired front versus achieved front Figure 3.21: The obtained rate solution and the comparison between desired and achieved front 85 The second numerical experiment corresponds to the permeability field shown in Figure 3.22(a). Injecting equal rates along the injection points produces results shown in Figure 3.23(a). The front propagates faster in the regions with higher permeability and will cause an early breakthrough if simulation continues. A different injection pol- icy must be applied if one intends to prevent poor displacement efficiency. We also apply the injection scenario where the rates are inversely proportional to the permeabil- ity at the injection points. This rate policy produces results shown in Figure 3.23(b). A quick comparison between the front propagation of the two injection policies shows no advantage of either methods. The front control technique was used to steer the front according to the pre-determined front dynamics shown in Figure 3.19(b). The obtained rate profiles shown in Figure 3.23(a). The obtained rate solution is used in a forward front tracking algorithm and produced the results shown in Figure 3.24(b). The compar- ison between desired and achieved is good for the times shown. 86 10 20 30 40 50 60 10 20 30 40 50 60 Permeability field X-axis Y-axis, Injection side 10 15 20 25 30 35 40 45 50 55 60 (a) Permeability varies in bothx andy directions 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired front dynamics X-axis Y-axis, Injection well time increases (b) The Desired front dynamics Figure 3.22: Permeability field and the desired front dynamics 87 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front propagation – Equal injection rate scenario X-axis Y-axis, Injection side time increases (a) The front propagation - equal rate injection policy 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front propagation– Rates are inversely proportional to permeability X-axis Y-axis, Injection side time increases (b) The front propagation - permeability proportional rate injection policy. Figure 3.23: Front propagation corresponds to the permeability field varies in both x andy directions shown in Figure 3.22(a)., (a) for the equal injection rate policy .,(b) and for rates inversely proportional to the permeability at the injection points 88 0 0.5 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Obtained injection rate profile X-axis Rate along injection well (a) Obtained rate solution 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Desired (blue) vs. achieved front (red) X-axis Y-axis, Injection side time increases (b) Desired front versus achieved front Figure 3.24: The obtained rate solution and the comparison between desired and achieved front 89 3.5 Summary We considered the problem of control of a displacement front in porous media, through flow-rate partition in a well. With the assumption of potential flow and in the absence of gravity and dispersion, we presented a formalism similar to the one previously described, that allows to control displacement fronts of a pre-determined dynamics using rate control. The problem considered is the same except that the displacement takes place in a heterogeneous porous media. The problem is formulated as a Fredholm Inte- gral Equation of the First Kind with the injection rates as the control parameters. The velocity inx andy directions are determined numerically and multi-linear interpolation methods are used to determine the velocities at each point(x;y) in the domain. For the solution of the integral equation, a regularization technique was used. In the previous chapter, it was postulated that the deterioration of the solution as time increases is a reflection of the cross-flow that develops. Through numerical experi- ments, we showed that introducing an anisotropy (with transverse-to flow-direction per- meability ratios of 0:1, 0:01 and 0:0001) decreases the cross-flow effect and limits the instabilities and the front control can be extended at much larger times. Also, we stud- ied the effect of heterogeneity on the control of front by injection rate partitioning. We demonstrated the potential application of the front control technique in displacement process through numerical examples. To assess the effectiveness of the method, we compared the result obtained by front control technique with the front propagation of two other injection policies, the equal injection rate policy and injection rates inversely proportional to the permeability at the injection points. We showed that, for a general heterogeneous media, the front control approach leads to achieve a controlled displace- ment as the others lead to an early breakthrough through higher permeability regions. 90 Chapter 4 The Inverse Problem 4.1 Introduction The purpose of this chapter is to give a brief overview of inverse problems, their charac- teristics and special treatments required in their solution. We introduce a set of mathe- matical tools currently available to analyze and solve the inverse problems. The Inverse Problems belong to a broader class of parameter estimation problems. These type of problems are common in many science and engineering applications where the goal is to estimate unknown parameters (either the input to the physical system or the phys- ical systems’s structure itself), given measurements of the outputs that only indirectly related to the unknown parameters. Some science and engineering examples include Earthquake location problem [NST85], Medical Imaging [Nat02], Satellite Naviga- tion, Remote sensing [Mil99], Geophysical imaging [Men89], Ground water modeling [Sun94] and Reservoir Engineering Applications [Oli92, Pen97]. For instance, in Darcy’s experiment, for a column of porous medium with unknown permeability, given the applied pressure difference at both inlet and outlet and by mea- suring the injection rate, one can obtain an estimation to the average permeability of the medium. The above example is the simplest in its form and more complicated exam- ples include reservoir characterization by means of history matching [MCWS93] and permeability estimation through pressure transient analysis [Oli90]. 91 In this chapter, the focus is mainly on a specific class of inverse problems, namely Fredholm Integral Equation of the First Kind [Bit95]. In comparison to the forward problems, the inverse problems usually do not have a straight forward solution, there- fore certain analytical tools are required to analyze the behavior of the system and its parameters. In the following sections, relevant analytical tools and the theoretical back- ground of the inverse problems are presented within the framework of Fredholm Integral Equation of the First Kind. 4.2 Fredholm Integral Equation of the First Kind 4.2.1 Introduction Fredholm Integral Equations of the First Kind arises in many science and engineering applications involving inverse problems. They belong to a broader class of Integral Equations [Kon91] where a kind of functional equation that includes the integral of an unknown function. Some examples include geophysics, antenna design, medical imaging, mathematical biology and remote sensing. It is a classical example of a linear ill-posed problem. The generic form of the equation is Z b a k(s;t)f(t)dt=g(s) c·s·d (4.1) where the kernelk(s;t), represents the mathematical model of the underlying physical process. The right hand side g(s) (the measured function) is known, while f(t) is the unknown function to be determined. There is a difficulty associated with the equation 4.1 when solving for unknown f(t). In the forward problem the integration has a smoothing effect on the function 92 f(t) especially for the high frequency components. It smoothes out the edges of the input function and damps the effect of high frequency components. In contrast, the high frequency components are amplified during the solution of the inverse problem and cause instability in the solution. The ill-posedness of the Fredholm integral equation of the first kind is illustrated below through a simple example with the following analysis. Let us consider the forward integral of equation 4.1. If one perturbs the solutionf by a very small value¢f(t) ¢f(t)=²sin(2¼pt); p=1;2:::; ²=constant and substitutes the perturbed solution into the equation 4.1 then the corresponding per- turbation of the data functiong is given by ¢g(s)=² Z b a k(s;t)sin(2¼pt)dt; p =1;2:::; Based on the Riemann-Lebesgue lemma , when p!1 then ¢g ! 0. Therefore, the ratio of k¢fk k¢gk can become extremely large by having the integer p large enough. The above example shows the sensitivity of the Fredholm Integral equation of the first kind to high frequency perturbations. For finite dimensional case which will be elaborate in next section the above mentioned ratio will be bounded but it may become very large. 4.2.2 Discrete Ill-posed Problems/Ill-conditioning The Fredholm equation 4.1 in its original form is a continuous ill-posed problem, where both kernel and data functions are defined within an interval. However, in real engineer- ing applications, instead of continuous form of equation 4.1, an equation with a discrete 93 right hand side is being used where there is a finite set of kernels k i and a finite set of measuredg i available. Therefore, the equation 4.1 is only continuous on one variablet Z b a k i (t)f(t)dt=g i i=1;¢¢¢ ;m (4.2) To solve such system of equations, one need to desceretize the integral into a finite set of discrete points. One can write the integral in a summation form with a finite number of terms applying appropriate collocation techniques. n X j k i j(t)f j (t)±t=g i i=1;¢¢¢ ;m j =1;¢¢¢ ;n (4.3) The goal is to solve for unknown functionx of the following systems of equations. Ax=b A2IR m£n (4.4) The above linear system of equations, can be solved for the unknownx through standard methods if the equation is well-posed. In the case of inverse problems and especially for the solution of the integral equation, the standard solution methods fail to compute the solution to the equation 4.4 and special treatments are necessary in obtaining accurate estimation of the functionx through regularization techniques. Hadamard [Had23], studied inverse problems and considers the equation 4.4 ill- posed if one or more of the following criterion do not apply: ² There exists a solution for any datab in data space ² The solution is unique ² The inverse mapping from datab to the unknownx is continuous 94 The first two conditions are related to the existence of the inverse of the kernel or the coefficient matrix A and will be explored in the subsequent sections. The third condi- tion requires a strong continuity between data space b and the unknown space x and is necessary but not sufficient condition for the stability of the solution [TC]. The third condition emphasizes the fact that the solution must depend continuously on the data b. In practice the data measurement (e.g measurement of the front location) has bias errors that in principle can destabilize the true solution. For a well-conditioned prob- lem, the relative error propagation from data spaceb to the solution spacex is controlled by the condition number of the coefficient matrix A. If ¢b is a change in data and ¢x the corresponding change of x, then the following relation explains the relative error propagation k¢xk kxk ·cond(A) k¢bk kbk (4.5) Thecond(A) is called the condition number of the coefficient matrix A and is defined as the ratio of the largest singular value of matrix A to the smallest. The process of obtain- ing the singular values are presented in the next section. A large condition number is not in favor of the solution and the solution becomes very sensitive to small perturba- tions in the measured data. Ill-conditioned coefficient matrices have very large condition numbers. There is a very close relationship between the ill-posedness of the Fredholm equation 4.1 and the ill-conditioning of the coefficient matrix A discretized by means of collocation methods.The degree of difficulty associated to inverse problems is not well defined and is problem dependent. In the next section we describes analytical and numerical tools that allow one to better analyze the inverse problem and determine the extend of the ill-posedness. 95 4.2.3 Singular Value Expansion Singular Value Expansion [Han92c], SVE is a powerful analytical tool for analysis of Fredholm integral equation of the first kind. With the definition of SVE, any square integrable kernel K can be presented as an infinite sum shown below: K(s;t)= 1 X i=1 ¹ i u i (s)v i (t) (4.6) The singular functions ofK areu i andv i , which are orthogonal with respect to the usual inner product. (u i ;u j ) = (v i ;v j ) = ¾ ij , where(;) is defined as(Á;!) = R Á(t)!(t)dt. The quantities¹ i are the singular values ofK and are always positive and can be written in a non-increasing order¹ 1 ¸¹ 2 ¸¹ 3 ¢¢¢¸ 0. One of the important relations between singular values and singular functions is the following: Z b a K(s;t)v i (t)dt=¹ i u i (s) (4.7) namely that any singular functionv i is mapped onto the correspondingu i with the sin- gular values¹ i as the amplification of this mapping. If we insert equation 4.7 and 4.6 into the original form of Fredholm integral equation 4.1, Z b a K(s;t)f(t)dt= Z b a 1 X i=1 ¹ i u i (s)v i (t)f(t)dt = 1 X i=1 ¹ i Z b a v i (t)f(t)dtu i (s) = 1 X i=1 ¹ i (v i ;f)u i (s) =g(s) (4.8) 96 Then it is not difficult to reach the following relationship 1 X i=1 ¹ i (v i ;f)u i (s)= 1 X i=1 (u i ;g)u i (s) (4.9) This leads to the following expression for the solution to the Fredholm integral equation of the first kind. f(t)= 1 X i=1 (u i ;g) ¹ i v i (4.10) Equation 4.10 shows that the functionf can be expressed in terms of singular functions v i and the corresponding expansion coefficients (u i ;g) ¹ i of the coefficient matrix. The SVE provides us a powerful tool to get insight on the kernel and its structure and characterize the solutionf by analyzing the coefficients (u i ;g) ¹ i and the functionv i . The properties of singular values and singular functions are not arbitrary and are strongly connected to the properties of the kernel of the integral. 4.2.4 Singular Value Decomposition Singular Value Decomposition, SVD is a discrete analog of the singular value expansion and has the same properties. We use SVD to analyze the discretized form of equation 4.1 as a system of equationsAx = b. In contrast to the equation 4.6 for SVE, the SVD of a m£n coefficient matrixA can be written as A= min(m;n) X i=1 ¾ i u i v T i (4.11) 97 and similarly the solution to theAx=b can be represented as below. x k ´ k X i=1 (u T i ;b) ¾ i v i (4.12) Similarly, the other relations derived for singular value expansion are also applicable for the singular value decomposition except for the summation, where the number of terms is finite. This is important because depending on the number of terms being used, we differentiate between different regularization methods based on singular values such as Truncated SVD and Generalized SVD. 4.3 Solution Techniques There is a numerical difficulty associated with the solution of the linear system of equa- tionAx=b, the discretized form of the Fredholm integral equation. The solution is very sensitive to the perturbations ofA andb because of the ill-posedness associated with the integral illustrated in previous section. In this section we present the least square solu- tion method and explore the difficulty in obtaining a solution for the discrete ill-posed problem. Then certain analysis and regularization methods are introduced in order to obtain acceptable solutions. 4.3.1 Least Square Method Considering the system of equations 4.4 in the form of a lest square problem where the solution of the following optimization problem min x kAx¡bk 2 A2IR m£n (4.13) 98 is to be found. Using SVD, similar to equation 4.10 a lest-square solution can be calcu- lated using singular value and singular functions x LSQ = n X i=1 (u i ;b) ¹ i v i (4.14) Having a good estimate solution for x, equation 4.14 can not be used directly to solve when the coefficinet matrix is coming from a discrete ill-posed problem. There are two criteria that are the direct result of the SVD and reveals the difficulty associated with the x LSQ solution and the degree of the ill-conditioning of the matrixA. ² the singular value ofA decays gradually to zero ² the ratio between the largest and the smallest nonzero singular values is large The presence of small singular values leads to more oscillations in the singular functions u i andv i and therefore the direct use of equation 4.14 leads to unstable solutions. One of the characteristics of the solution is that its norm (the size of the solution vector) is significantly higher than the norm of the exact solution. With the same reasoning, the common standard techniques such as Cholesky, LU decomposition and QR factorization will fail to compute a solution to the integral equation. One possible way to improve the solution is to add constraint to the size of the solution as well as filter out the high frequency components associated with the small singular values. Then it is feasible to reach a meaningful solution reasonably close to the underlying exact solution. 4.3.2 The Picard Condition Equation 4.10 produced by singular value decomposition has two important elements. In order to have a smooth solution f(t), the right hand side of equation 4.1, g(s) is 99 required to be smoother than the functionf(t). Therefore, the following condition, the Picard condition must satisfy. 1 X i=1 ( (u i ;g) ¹ i ) 2 <1 (4.15) The absolute value of the coefficients (u i ;g) must decay faster than the corresponding singular values ¹ i in order that a solution exists. The importance of the Picard condi- tion in its connection with the Fredholm integral equation of the first kind requires this condition to be checked every time one tries to solve the integral equation. It shows the severity of ill-posedness associated with the smoothing effect of the kernel under the integral. The discretized form of the Picard condition is called the discrete Picard condition. The smaller singular values and their corresponding singular functions must be filtered out by the solution technique. This prevents the oscillatory behavior of the solution. In fact, the quantities ¹ i and (u i ;g) can be used in a plot that visually allow us to inspect wether the discrete Picard condition is satisfied for the system of equations produced from a collocation of integral equations. The importance of Picard condition can be understood by the following example taken from [Han92c]. A one-dimensional image restoration can be written as a Fredholm Integral Equation of the First Kind. Once the kernel and the blurring functions are given we can produce a system of equations by implementing the integration in a forward process and arrive at the final blurred pic- ture. Now if we desire to restore the original picture, we must solve the inverse problem which is in the form of an integral equation and is an ill-posed problem. In Figure 4.1, the top graph shows the behavior of an unperturbed image restoration where the(u i ;b) decays faster than the singular values, while the bottom one shows the behavior of the same problem with additional noise added to the blurred picture. In order to regularize 100 0 5 10 15 20 25 30 35 40 10 −20 10 −10 10 0 10 10 i Discrete Picard plot μ i |u i T b| |u i T b|/μ i 0 5 10 15 20 25 30 35 40 10 −20 10 0 10 20 i Discrete Picard plot μ i |u i T b| |u i T b|/μ i Figure 4.1: Discrete Picard plot for the unperturbed (top) and perturbed (bottom) data the solution, we should dampen the components that destabilize the solution and use the first few eigenvalues for the perturbed image restoration instead. The Picard condi- tion in its discrete or continuous form is a diagnostic tool to be applied prior to apply any solution techniques. If the Picard condition does not satisfy even a few numbers of singular values, then it is not possible to reach a meaningful solution. This analysis allow us to see if the solution does exist but will not determine if the instability comes form the integral equation itself, the decsritization or from data sources. The following integral equation adopted from reference [Han92c] has no solution as it does not satisfy the discrete Picard condition depicted in Figure 4.2. Z 1 i=1 1 s+t+1 f(t)dt=1 (4.16) 101 0 5 10 15 20 25 30 35 40 10 −15 10 −10 10 −5 10 0 10 5 10 10 i Discrete Picard plot μ i |u i T b| |u i T b|/μ i Figure 4.2: Discrete Picard plot is not satisfied for the integral equation 4.16 4.4 Regularization Methods The purpose of using regularization methods is to diminish or dampen the contributions to the solution by the smaller singular values. It is helpful to incorporate any prior information and knowledge about the desired solution. The most common approach is to introduce a side constraint that limits the size of the solution. This is the goal of applying any regularization methods. By introducing the side constraint, one must seek an estimation to the exact solution rather than the exact solution itself. The regularized solution must satisfy a balance between minimizing the side constraint and minimizing the residual normkAx¡bk 2 . In the following sections we consider the numerical treatment of discrete ill-posed problems and briefly review the relevant regularization methods. These methods include 102 direct and iterative procedures. Here, we focus on direct methods where the solution is described by a direct computation (e.g. Tikhonov method, Truncated SVD) and explain their solution procedures. The iterative Conjugate Gradients methods will be presented as well. 4.4.1 Tikhonov Regularization Technique The most common and direct form of regularization is the Tikhonov regularization technique. Tikhonov [TA77] introduced a regularized form of the solution to the Fredholm integral equation, namely a side constraint term that prevents unstable solutions. In theory, regularization is applicable to the integral equation proposed by Tikhonov in his original work. However, in practice it is simpler to implement it to the linear system produced by descretization of the integral equation. The generalized form of the Tikhonov regularization [TA77, TG91] is to find a solu- tion that minimizes a weighted combination of the residual norm and the introduced side constrained of x ¸ =argminfkAx¡bk 2 2 +¸ 2 kL(x¡x o )k 2 2 g (4.17) where the matrix L is either the identity matrix I n or a banded matrix with full rank. The regularization parameter, ¸ controls the amount of importance given to the side constraint compared to the minimization of the residual norm. It is apparent that larger ¸ correspond to a large amount of weight given to the side constraint with the cost of larger residual norm , while smaller¸ have the opposite effect. The role of¸ is to control the sensitivity of the solution x ¸ to the small perturbations in A and b. With singular 103 value decomposition concept, the regularized solutionx ¸ withx o = 0 can be written in the following from. x ¸ = n X i=1 f i (u i ;g) ¹ i v i ifL=I n (4.18) x ¸ = n X i=1 f i (u i ;g) ¹ i v i + n X i=1 (u i ;g)v i ifL6=I n (4.19) Equation 4.18, whereL=I n is called the standard Tikhonov method and equation 4.19 is called the general Tikhonove method. The filter factors, f i must go to zero as the singular values¹ i decreases. This property of filter factors filters out the effect of small singular values. The major difference among the various regularization techniques is actually the way that these methods determine the filter factors. In applying Tikhonov regularization methods, the filter factors are defined as eitherf i = ¹ 2 i (¹ 2 i +¸ 2 ) for the stan- dard Tikhonov regularized method orf i = ° 2 i (° 2 i +¸ 2 ) for the generalized Tkhonov regular- ization method. The idea behind Tikhonov regularization is to introduce a way to control the tradeoff between the residual norm and the size constraint. The optimum value of ¸ is the one that balances the two. This is the concept behind the L-curve technique that visually helps to choose the best¸. The L-curve is a plot of the sidelogkL(x¡x o )k versus the error normlogkAx¡bk for the Tikhonov and similar regularization methods. Figure 4.3 is taken from Hasnen [Han92c] shows the generic form of the L-curve and the corner point corresponding to the optimum regularization parameter. 104 Figure 4.3: The generic form of the L-curve plot 4.4.2 Truncated SVD The second regularization method reviewed here is the Truncated Singular Value Decomposition [Han90, Var79] where simply we truncate the SVD-based expres- sion 4.12 and consider only k number of terms for summation purposes. We remove the smaller singular values before they dominate and cause instability in the solution x k ´ k X i=1 (u T i ;b) ¾ i v i (4.20) Here, integerk acts as a regularization parameter. Hansen [Han92b] showed that when the matrix L is the identity matrix then the truncated SVD is basically the same as the Tikhonov regularization method. 105 4.4.3 Semi-Iterative Conjugate Gradients The application of iterative methods in solving discrete inverse problems are not com- parable to the application of direct methods [AP81]. The Conjugate Gradient method is robust and popular for solving large sparse system of equations when the prob- lem is well-posed [Han92c]. The only interesting fact that attracts this method for solving ill-posed problems is that when applied to the unregulated normal equations Ax=b, the low-frequency components of the solution tend to converge faster the high- frequency components. Therefore, the Conjugate Gradients method has some regular- ization effect implicitly. Here the number of iterations act as the regularization parame- ter [Han92a, Bra86]. 4.4.4 The Regularization Parameter A key element that controls the solution of the regularized inverse problem is the regu- larization parameter. It is important to have an optimum regularization parameter as it balances the regularization and residual errors. If the regularization parameter is smaller than the optimal value, then the solution will be masked by the error contribution from smaller singular values. When the regularization parameter is larger than the optimal value, over filtering occurs that leads to the loss of information. The regularization parameter in Truncated SVD is the k number of singular values considered to reach the solution. In Semi-Iterative Conjugate Gradient method, the regu- larization parameter is the number of iteration being used to obtain the regularized solu- tion. There are different techniques such as discrepancy principal, the quasi-optimality criterion [Mor84], and the generalized cross-validation available to select the optimal 106 regularization parameter. The concept of the optimum regularization parameter for the Tikhonov method is shown in Figure 4.3. 4.5 Summary In this chapter, we presented solution techniques for a class of inverse problems, reviewed analytical tools and numerical methods to analyze them. As described in the previous chapters, the mathematical formulation of the control of front through porous medium falls in the category of the Fredhlom Integral Equation of the First Kind. To obtain the solution of the integral equation, an inverse problem must be solved, which is ill-posed and requires regularization. We introduced analytical tools that can be used for diagnostic purposes. In addition, numerical tools are described that allow us to arrive at a solution to the discrete form of the inverse problems. 107 Chapter 5 Analysis of Forward Problem 5.1 Introduction In the previous chapters, the front steering problem was addressed in a simple rectan- gular geometry based on potential flow. The objective was to answer the fundamental question in rate-control flow problems, namely the control of a displacement front via flow rate partition in a well. We arrived at a front steering technique to partition the rate within the well, so that the displacement front can be steered accordingly. When the reservoir is homogeneous and isotropic, an integral equation in an analytical form was derived, the solution of which prescribes the desired injection rate profiles. A sim- ilar approach applied to an anisotropic or a heterogeneous porous medium except the integral equation is determined numerically. For the solution of the integral equation, a regularization technique is necessary. However, it was found that numerical instabilities do develop, even with the use of regularization, for later times when the front moves away from the injection well. This behavior is seen in all homogeneous and heterogeneous porous medium examples. We hypothesized that the numerical instabilities develops at larger times when transverse cross-flow is large. It was shown that the instabilities diminish with a more stratified structure. In fact, the deterioration of the solution at the larger times (as the front moves away from injection point) can be damped and diminished by decreasing the vertical 108 cross-flow. This is done by choosing small values of R L (which is the same as small values of Ky K x ), where the problem behaves as a non-communication layered system. The fact is, as the front evolves and moves away from injection points, the ability to control its dynamics by varying injection rates diminishes. Therefore, it is important to quantify and explain the sensitivity of the front to the rate fluctuation at larger times. In the following sections the forward problem is considered and an analytical analysis of the effect of the injection rate fluctuation on the front at larger times is presented. For convenience only, a homogeneous porous medium is considered, where an analytical expression can be obtained for direct analysis purposes. 5.2 Semi-Infinite Reservoir In order to analyze the forward problem analytically, a semi-infinite reservoir with boundary conditions similar to the one used in the previous chapters is considered. With the assumption of incompressible miscible fluids of equal mobility and in the absence of dispersion and gravity effects, the governing equation in the semi-infinite homogeneous system is the Laplacian equation. @ 2 P @ 2 x + @ 2 P @ 2 y =0:0 (5.1) The following boundary conditions in dimensionless notations are applied to the system: ¡ @P @x j x=0 =q(y;t) P j x!1 =0 109 Figure 5.1: The semi-infinite reservoir ¡ @P @y j (y=0) =0 ¡ @P @y j (y=1) =0 The pressure field can be expressed in an integral form containing the injection rate q(y), either by Green’s function or the Separation of Variables method. Similar to the procedure described in section 2.3.1, the pressure value at any point (x;y) within the domain and on the boundary can be analytically expressed as P(x;y)=¡x+ X f 2cos(¸ n y)exp(¡¸ n x) ¸ n Z 1 0 q(´)cos(¸ n y)d´g (5.2) where ¸ n = n¼ is the eigenvalue fixed by the homogeneous boundary conditions at y =0 andy =1. In addition,n is the number of terms andn=1;¢¢¢ ;1. 110 Assume that the front is given by the general equation below that relates the front progression and the velocity field F(x;y;t)=x(t)¡f(y;t)=0 (5.3) where as a function of time, the front (the moving surface) evolves. Taking the substan- tial derivative of the moving surface leads to the following relationship between the front velocity and the local pressure change at the front. This relationship can be written: @f @t =¡ @P @x + @P @y @f @y (5.4) where @P @x =¡1¡ X f2cos(¸ n y)exp(¡¸ n x) Z 1 0 q(´)cos(¸ n ´)d´g (5.5) @P @y =¡ X f2sin(¸ n y)exp(¡¸ n x) Z 1 0 q(´)cos(¸ n ´)d´g (5.6) 5.3 Method of Characteristics For the forward problem at hand, with a given injection rate profile along the injection well, the pressure field need to be calculated only once. Then front tracking algorithm can be used to move the particles defining the front forward at each time step. In con- trast to the inverse problem, the forward problem is well-posed and the front solution is unique so there is no need for regularization technique. Here, the objective is to arrive at an analytical solution for the front in the form ofx(t) = f(y;t) corresponding to an injection rate profile. 111 The Method of Characteristics is used to produce a solution of the kinematics equa- tion 5.4. The method of characteristics is a technique that simplifies the first order partial differential equation into ordinary differential equations along the characteristic lines. Equation 5.4 can be written as: f t +G(f;y)f y +H(f;y)=0:0 (5.7) whereG(f;y) andH(f;y) are G(f;y)=¡ X fa n sin(¸ n y)exp(¡¸ n f)g H(f;y)=1+ X fa n cos(¸ n y)exp(¡¸ n f)g where a n =2 Z 1 0 q(´)cos(¸ n ´)d´ The two characteristic equations are dy dt =G(f;y) (5.8) df dt =H(f;y) (5.9) and one must solve them first using the initial conditionf(y;t) = f(y;0). The relation df dy can be written as: df dy = df dt dy dt = H(f;y) G(f;y) = 1+ P fa n cos(¸ n y)exp(¡¸ n f)g ¡ P fa n sin(¸ n y)exp(¡¸ n f)g (5.10) 112 Consider©(y;f) defined as ©(y;f)= X f a n sin(¸ n y)exp(¡¸ n f) ¸ n g (5.11) Since d©= @© @y dy+ @© @f df = X fa n cos(¸ n y)exp(¡¸ n f)gdy¡ X fa n sin(¸ n y)exp(¡¸ n f)gdf (5.12) one can arrive at d© dy = 1. This leads to the solution of the new equation Á = y + constant, where the constant is the proper initial condition for the front. ©(y;f)=y¡´¡ X f a n sin(¸ n ´)exp(¡¸ n f o (´)) ¸ n g (5.13) where 0 · ´ · 1. To solve for the front f(y;t), we use the equation dy dt = G(f;y), which can be interpreted to give Z y ´ dy P fa n sin(¸ n y)exp(¡¸ n f)g =t (5.14) It is difficult to find an analytical expression for equation 5.14. To analyze the effect of injection rate fluctuations on front dynamics, for an special case, the case of one harmonic is examined below. The result then can be extended to the more general case. 113 5.4 Special Case: One Harmonic Front Solution 5.4.1 With Constant Coefficient The boundary condition imposed on the left can be represented as a fourier series. q(y)=1+ 1 X 1 fa n cos(¸ n y)g (5.15) wherea n is a n =2 Z 1 0 q(y)cos(¸ n y)dy Each individual term of the summation of the injection rate and its corresponding ampli- tude can be selected to solve for the front propagation. For an analytical result this procedure allows the injection rate to be a constant plus one harmonic term, for example q(y)=1+a n cos(¸ n y) a n =2 Z 1 0 q(y)cos(¸ n y)dy With this assumptions, the two relations 5.13 and 5.14 become a n sin(¸ n y)exp(¡¸ n f) ¸ n =¡y+´+ a n sin(¸ n ´)exp(¡¸ n f o (´)) ¸ n (5.16) Z y ´ dy a n sin(¸ n y)exp(¡¸ n f) =t (5.17) 114 Substituting the right-hand-side of 5.16 into the denominator of the integral of the equa- tion 5.17, we can integrate to find Z y ´ dy a n sin(¸ n y)exp(¡¸ n f) = Z y ´ dy ¡y+´+ a n sin(¸ n ´)exp(¡¸ n f o (´)) ¸n =¡lnj a n sin(¸ n ´) ¸ n y) +´¡yj y ´ =ln a n sin(¸ n ´)exp(¡¸ n f o (´)) ¸n ansin(¸n´)exp(¡¸nfo(´)) ¸n ¡y+´ =t (5.18) Assuming the initial front location f o (´) = 0. This leads us to the following set of equations: a n sin(¸ n y)exp(¡¸ n f) ¸ n =¡y+´+ a n sin(¸ n ´) ¸ n (5.19) a n sin(¸ n ´)exp(¡¸ n t) ¸ n =¡y+´+ a n sin(¸ n ´) ¸ n (5.20) from which f(y;t)= 1 ¸ n logf a n sin(¸ n y) a n sin(¸ n ´)exp(¸ n t) g = 1 ¸ n logf sin(¸ n y) sin(¸ n ´) g+t (5.21) This can also be represented as ¸ n [f(y;t)¡t]=log sin(¸ n y) sin(¸ n ´) (5.22) Analysis of this equation shows a decrease in the amplitude of harmonic in propor- tion inverse to its frequency. The front reaches a traveling wave solution and for higher harmonic terms the fluctuation in injection rate will diminish faster as the wavelength is 115 smaller. This is shown graphically in Figure 5.2 for six harmonics. As time increases, the sensitivity of the front to the fluctuations in the rate for larger values of harmonics is diminishing. This phenomenon explains why as the front progresses we loose control over its steering by varying injection rates. 5.4.2 With Variable Coefficient Consider now a problem similar to the previous with the additional assumption that the injection rate has variable coefficient a n (t), namely q(y;t) = 1+a n (t)cos(¸ n y). Through additional manipulation and simplification, one can reformulate equations 5.13 and 5.14: a n (t)sin(¸ n y)exp(¡¸ n f) ¸ n =¡y+´+ a n (t)sin(¸ n ´)exp(¡¸ n f o (´)) ¸ n (5.23) a n (t) Z y ´ dy a n (t)sin(¸ n y)exp(¡¸ n f) = Z t 0 a n (t)t (5.24) Then we arrive at the following solution for the front f(y;t)= 1 ¸ n logf sin(¸ n y) sin(¸ n ´) g+ 1 a n (t) ¿ (5.25) where¿ = R t 0 a n (t)dt. The front solution with a time-varying rate differs from the front solution with constant rate in the oscillating behavior of the front fluctuation. 116 −0.1 −0.05 0 0.05 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front location at time t = 0.025 Front location [f(y,t)−t] Y-axis, Injection side n=1 n=2 n=3 n=4 n=5 n=6 (a) −0.1 −0.05 0 0.05 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front location at time t = 0.075 Front location [f(y,t)−t] Y-axis, Injection side n=1 n=2 n=3 n=4 n=5 n=6 (b) −0.1 −0.05 0 0.05 0.1 0.15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front location at time t = 0.25 Front location [f(y,t)−t] Y-axis, Injection side n=1 n=2 n=3 n=4 n=5 n=6 (c) −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front location at time t = 0.75 Front location [f(y,t)−t] Y-axis, Injection side n=1 n=2 n=3 n=4 n=5 n=6 (d) −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front location at time t = 1.5 Front location [f(y,t)−t] Y-axis, Injection side n=1 n=2 n=3 n=4 n=5 n=6 (e) −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Front location at time t = 3.0 Front location [f(y,t)−t] Y-axis, Injection side n=1 n=2 n=3 n=4 n=5 n=6 (f) Figure 5.2: The amplitude of the harmonic decreases inversely proportional to its fre- quency. The results for first six harmonics are presented here. As time increases, the sensitivity of the front to the rate fluctuations decreases inversely proportional to the frequency 117 5.5 Special Case: Two Harmonics Front Solution Consider now the case where the case where the injection rate consists of more than one harmonics,q(y)=1+a n cos(¸ n y)+a m cos(¸ m y) wheren<m. One can reformulate equations 5.13 and 5.14 as shown below a n sin(¸ n y)exp(¡¸ n f) ¸ n + a m sin(¸ m y)exp(¡¸ m f) ¸ m =¡y+´+ a n sin(¸ n ´)exp(¡¸ n f o (´)) ¸ n + a m sin(¸ m ´)exp(¡¸ m f o (´)) ¸ m (5.26) Z y ´ dy a n sin(¸ n y)exp(¡¸ n f)+a m sin(¸ m y)exp(¡¸ m f) =t (5.27) At large times, the following is applicable a n sin(¸ n y)exp(¡¸ n f) ¸ n ¿ a m sin(¸ m y)exp(¡¸ m f) ¸ m namely, we can neglect the m th order harmonic compared to the n th . Thus, a similar can be found f(y;t)= 1 ¸ n logf sin(¸ n y) sin(¸ n ´)+ na m man sin(¸ n ´) g+t (5.28) showing that the high-frequency behavior (m th harmonic) is suppressed compared to the lower frequencies. 118 5.6 Summary In this chapter, we considered a semi-infinite homogeneous porous media, where we solved for the displacement front dynamics given the injection rates. This is a for- ward problem and is well-posed. We obtained a general expression describing the front dynamics. However, in order to quantify and analyze the effect of rate fluctuations on the front dynamics, we considered the special case of injection rate with one harmonic term. The results show that at larger times (when the front is far from the injection side) the front is less sensitive to the fluctuations in rate for larger harmonic terms. Similarly we obtained the solution to the front dynamics for variable rate of one harmonic. The front solution with a time-varying rate differs from the front solution with constant rate in the oscillating behavior of the front fluctuation. 119 Chapter 6 Conclusions and Recommendations 6.1 Conclusions In this thesis, we studied the control of a displacement front through porous media by partitioning injection rates along the well for a two-dimensional rectilinear geometry. It was assumed that the injection points along the injection well are independent of each other and can be controlled separately. The production side of the system is at constant pressure. We considered a first contact miscible displacement at mobility ratio of one and it was assumed that the effects of gravity, dispersion and capillarity are absent and the flow is incompressible. We studied a fundamental question in rate-control flow problems, namely the con- trol of a displacement front via flow-rate partition in a well. The specific question we addressed was how to partition the flow rate within the well, so that the displacement front can be steered according to pre-determined dynamics. The problem was addressed in a simple rectangular geometry based on potential flow. When the reservoir is homo- geneous and isotropic, an integral equation in an analytical form was derived, where its solution determines the desired injection rate profile. We showed that the injection rates and the front kinematics can be related through the integral equation where the ker- nel is computed analytically. The solution to the integral equation is indeed an inverse problem and regularization methods were applied to obtain rate solutions. 120 Numerical experiments of the front control in homogeneous reservoirs for various fronts dynamics were conducted to test the ability of the front steering technique. The results of the numerical experiments showed that we can control the front by partition- ing the rates along the injection well. In addition, we applied a front tracking algorithm and used the obtained rate solution to produce the front dynamics. Comparison between achieved and desired fronts were accurate for early times. However, as the front moves away from the injection points, there is a deterioration in the solution. The decreas- ing ability to control the front continues as the front progresses in time and that was believed to be due to the development of cross-flow effects. We showed that by intro- ducing anisotropy in the vertical direction, one can limit the cross-flow and control the front for larger times. Through applying different methods in determining the optimum regularization parameters, it was found that regularization is necessary to obtain the rate solution at each time step. However, in all numerical experiments, one single value of the parameter can produce an accurate estimation of the rate solution of all time steps. In addition the degree of ill-conditioning is found to be related to the complexity of the kernel function which includes both the degree of heterogeneity of the media and the sharpness of the front. In the homogeneous porous media, the extent to which one can control the front depends on the sharpness of the front. This was shown through numerical experiments of fronts that had different sharpness. We showed that for a heterogeneous porous medium, a similar technique can be applied to relate the injection rates to the front dynamics. The potential application of the front control in heterogeneous porous media was presented through numerical examples. For different heterogeneity generated by Sequential Gaussian Simulation, we compared the front propagation of three injection rate policies. For the front control approach, the desired front selected to propagate resembling a uniform displacement. 121 The result compared to the equal injection policy at the injection well and to the injection rate inversely proportional to the permeability at the injection points. Results showed that, for a general heterogeneous media, the front control approach leads the ability to obtain a uniform displacement as the two other scenarios lead to early breakthrough at higher permeability regions. As the front evolves and moves away from the injection points (at larger times), the ability to control the front is weaken. We considered the forward problem in a semi- infinite and homogeneous porous media and showed that the sensitivity of the front to the rate fluctuations diminishes inversely proportional to the frequency. As a result, the high-frequency components of the rate has no effect on the front. 6.2 Future Recommendations The study of the front control in displacement process by injection rate partition along the well presented in this thesis is not yet fully developed. There is room for further analysis of the front control problem. Below is the list of recommended studies that are relevant to the work presented in this thesis. We considered the front control problem in a two-dimensional rectangular porous media, where the displacement is at unit-mobility ratio, where we neglected the effect of gravity and dispersion. This allowed us to achieve an integral equation in an analytical form for the homogeneous porous media. In addition, a similar approach applies for a heterogeneous system, except that the kernel in the integral equation is determined numerically. We recommend a future study to address the front control problem of fluids with different mobility ratios. We also recommend the study of the displacement front control in the presence of gravity. 122 We considered the problem of steering the front at will by manipulating injection rates only for the case with one injection well equipped with independent injection inter- vals and in a two-dimension porous media. We recommend further study of applying this concept for three-dimensional porous media with multiple wells. 123 References [AODG06] A.H. Alhuthali, D. Oyerinde, and A. Datta-Gupta. Optimal waterflood management using rate control. San Antonio, TX, 24-27 September 2006. SPE Annual Technical Conference and Exhibition. [AP81] R.S. Anderssen and P.M. Prenter. A formal comparison of methods pro- posed. for the numerical solution of first kind integral equations. J. Aust. Math. Soc. Series B, 22, 1981. [Aro62] J.S. Aronofsky. Linear programming a problem-solving tool for petroleum industry management. Journal of Petroleum Technology, 14(7):729–736, July 1962. [Ash88] H. Asheim. Maximization of water sweep efficiency by controlling pro- duction and injection rates. London, United Kingdom, 16-19 October 1988. European Petroleum Conference, Society of Petroleum Engineers. [BH97] A. Bittencourt and Roland N. Horne. Reservoir development and design optimization. San Antonio, TX, October 1997. SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [Bit95] A.V . Bitsadze. Integral equations of First Kind. World Scientific Pub- lishing Co. Pte. Ltd., Steklov Institute for Mathematics Moscow, Russua, 1995. [BJ02] D.R. Brouwer and J.D. Jansen. Dynamic optimization of water flooding with smart wells using optimal control theory. Aberdeen, United King- dom, 29-31 October 2002. European Petroleum Conference, Society of Petroleum Engineers. [BJvdSS + 01] D.R. Brouwer, J.D. Jansen, van der S. Starre, C. Berentsen, and van C. Kruijsdijk. Recovery increase through water flooding with smart well technology. The Hague, Netherlands, 21-22 May 2001. SPE European Formation Damage Conference, Society of Petroleum Engineers. 124 [Bra86] H. Brakhage. On ill-posed problems and the method of conjugate gradi- ents. Proc. of Inverse and Ill-posed Problems, pages 165–175, 1986. [BW30] A.H. Bell and E.W. Webb. Repressuring the selover zone at seal beach and the effect of proration spaces. AIME Transactions, 86:240–245, 1930. [CBJ46] W. Cerini, R. Willis Battles, and P.H. Jones. Some factors influencing the plugging characteristics of an oil-well injection water. AIME Trans- actions, 165:52–63, 1946. [CHD69] A.R. Cooksey, J.H. Henderson, and John Dempsey. Total computer sim- ulation of a gas producing complex. Journal of Petroleum Technology, 21(8):942–948, August 1969. [Cor30] C.S. Corbett. Equilateral trianular system of well spacing Banach spaces. AIME Transactions, 86:168–173, 1930. [Dal49] C.R. Dale. Bottom hole flow surveys for determination of fluid and gas movements in wells. Petroleum Transactions, AIME, 186:205–210, 1949. [Fle91] C. Fletcher. Computational Techniques for Fluid Dynamics, volume 1. Springer, 1991. [Flo02] R. Floberghagen. Lunar Gravimetry: Revealing the Far-Side. Springer, 1 edition edition, May 31 2002. [FR84] Z. Fathi and Fred Ramirez. Optimal injection policies for enhanced oil recovery: Part 2- surfactant flooding. SPE Journal, 24(3):333–341, June 1984. [GC00] G.H. Grinestaff and D.J. Caffrey. Water management: A case study of the northwest fault block area of prudehule bay , alaska, using streamline simulation and traditional waterflood analysis. Dallas, TX, 1-4 October 2000. SPE Annual Technical Conference and Exhibition. [Geo59] H. George. Gas injectivity logging. Midland, TX, May 1959. Permian Basin Oil Recovery Conference, Society of Petroleum Engineers. [Gri99] G.H. Grinestaff. Waterflood pattern allocations: Quantifying the injector to producer relationship with streamline simulation. Anchorage, Alaska, 22-24 May 1999. SPE Western Regional Meeting. [Had23] J. Hadamard. Lectures on Cauchy’s Problem in Linear Partial Differen- tial Equations. Yale University Press, New Heaven, USA, 1923. 125 [Han90] P.C. Hansen. Truncated svd solutions to discrete ill-posed problems withill-determined numerical rank. SIAM Journal on Scientific and Sta- tistical Computing, 11(3), May 1990. [Han92a] M. Hanke. Regularization with differential operators. an iterative approach. J. Numer. Funct. Anal. Optim., 13, 1992. [Han92b] P.C. Hansen. Analysis of discrete ill-posed problems by means of the l-curve. SIAM Press, 34(4):561 – 580, December 1992. [Han92c] P.C. Hansen. Numerical tools for analysis and solution of fredholm inte- gral equations of the first kind. Inverse Problems, (6):849–872, Decem- ber 1992. [Has30] W.P. Haseman. A theory of well spacing (banach) spaces. AIME Trans- actions, 86:146–155, 1930. [HRK98] T.J. Harding, N.J. Radcliffe, and P.R. King. Hydrocarbon production scheduling with genetic algorithms. SPE Journal, 3(2):99–107, June 1998. [Hup74] J.D. Huppler. Scheduling gas field production for maximum profit. SPE Journal, 14(3):279–294, June 1974. [HY02] G.W. Hanson and Alexander B. Yakovlev. Operator Theory for Electro- magnetics: An Introduction. Springer, Printed in New York, 2002. [Iva79] S. Ivar. Green’s functions and boundary value problems. Wiley, New York, 1979. [Jos57] B. Joseph. Selective plugging of waterflood input wells theory, methods and results. Journal of Petroleum Technology, 9(3):26–31, March 1957. [Kon91] J. Kondo. Integral equations. Oxford applied mathematics and comput- ing series, Science Council of Japan, 1991. [LA58] A.S. Lee and J.S. Aronofsky. A linear programming model for scheduling crude oil production. Journal of Petroleum Technology, 10:51–54, July 1958. [Lak89] L.W. Lake. Enhanced oil recovery. Prentice Hall, New York, 1989. [LRQ93] W. Liu, F. Ramirez, and Yu Feng Qi. Optimal control of steamflooding. SPE Advanced Technology Series, 1(2):73–82, July 1993. 126 [MCWS93] E.M. Makhlouf, Wen Chen, Mel L. Wasserman, and John H. Seinfeld. A general history matching algorithm for three-phase, three-dimensional petroleum reservoirs. SPE Advanced Technology Series, 1(2):83–92, July 1993. [ME78] J.E. Murray and T.F. Edgar. Optimal scheduling of production and com- pression in gas fields. Journal of Petroleum Technology, 30(1):109–116, January 1978. [Men89] W. Menke. Geophysical Data Analysis: Discrete Inverse Theory, vol- ume 45 of International Geophysics series. Elsevier, revised edition edi- tion, 1989. [Mil99] A.S. Milman. Mathematical Principles of Remote Sensing: Making Inferences from Noisy Data. Sleeping Bear Press, London, March 1999. [Mor84] V .A Morozov. Methods for Solving Incorrectly Posed Problems. Springer, Printed in Berlin, 1984. [Mus40] M. Muskat. Principles of well spacing Banach spaces. AIME Transac- tions, 136:37–56, 1940. [MW34] M. Muskat and R.D. Wyckoff. A theoretical analysis of water-flooding networks Banach spaces. AIME Transactions, 107:62–76, 1934. [Nat02] F. Natterer. The mathematics of computerized tomography. Medical Physics, 29(1):107–108, June 2002. [NST85] G.D. Natale, R. Scarpa, and E. Terzini. Earthquake location inverse prob- lem: a review and an application to the aftershocks of the southern italy 23 november 1980 earthquake. Il Nuovo Cimento C, 8(1):1–25, January 1985. [Oli90] D.S. Oliver. The averaging process in permeability estimation from well- test data. Soc. Pet. Eng. Formation Evaluation, 5(3):319–324, September 1990. [Oli92] D. Oliver. Estimation of radial permeability distribution from well test data. Soc. Pet. Eng. Formation Evaluation, 7(4):290–295, 1992. [Ozi80] M.N. Ozisik. Heat Conduction. Wiley-Interscience, Printed in the United States of America, 1980. [Pen97] Yue-Jun Peng. An inverse problem in petroleum exploitation. Inverse Problems, 13(6), December 1997. 127 [Pfi48] R.J. Pfister. An improved water-input profile instrument. AIME Transac- tions, 174:269–285, 1948. [Phe29] R. Phelps. Analytical principles of the spacing of oil and gas wells Banach spaces. AIME Transactions, 82:90–103, 1929. [RFL84] F. Ramirez, Zohreh Fathi, and Jean Luc. Optimal injection policies for enhanced oil recovery: Part 1- theory and computational strategies. SPE Journal, 24(3):328–332, June 1984. [RG74] G. Rosenwald and D.W. Green. A method for determining the optimum location of wells in a reservoir using mixed-integer programming. SPE Journal, 14(1):44–54, February 1974. [Ros19] J.H. Roswell. Staggering locations for oil wells Banach spaces. AIME Transactions, 61:611–616, 1919. [SB48] D. Silverman and A.R. Brown. Location of points of water entry in oil wells. Petroleum Transactions, AIME, 174:286–304, 1948. [SH83] A.B. See and Roland Horne. Optimal reservoir production scheduling by using reservoir simulation. SPE Journal, 23(5):717–726, October 1983. [Sud98] B. Sudaryanto. Optimization of Displacement Efficiency of Oil Recovery in Porous Meida Using Optimal Control Theory. Petroleum engineering, University of Southern California, The Graduate School, University Park, Los Angeles,CA 90007, December 1998. [Sun94] Ne-Zheng Sun. Inverse Problems in Groundwater Modeling, volume 6 of Theory and Applications of Transport in Porous Media. Springer, 1994. [SY00] B. Sudaryanto and Y .C. Yortsos. Optimization of fluid front dynamics in porous media using rate control i. equal mobility fluids. Physics of Fluids, 12(7):1656–1670, 2000. [SY01] B. Sudaryanto and Y .C. Yortsos. Optimization of displacements in porous media using rate control. New Orleans, Louisiana, 30 September-3 Octo- ber 2001. SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [TA77] A.N. Tikhonov and V .Y . Arsenin. Solutions of ill-posed problems. Hal- sted Press, 1977. [TB06] M. Thiele and R.P. Batycky. Using streamline-derived injection effi- ciencies for improved waterflood management. SPEREE, 9(2):187–196, April 2006. 128 [TC] S.M. Tan and Fox Colin. Introduction to inverse problems. chapter 1, pages 1–14. Physics 707 Inverse Problems. [TG91] A.N. Tikhonov and A.V . Goncharsky. Ill-posed problems in natural sci- ences. Brill Academic Publishers, December 1991. [Var79] J. Varah. A practical examination of some numerical numerical methods for linear discrete ill-posed problems. SIAM Review, 21(3), 1979. [Wah80] G. Wahba. Ill-posed problems: numerical and statistical methods for mildy,moderately and severely ill-posed problems with noisy data. Tech- nical Report 595, Department of Statistics, University of Wisconsin, 1980. [WBM34] R.D. Wyckoff, H.G. Botset, and M. Muskat. The mechanics of porous flow applied to water-flooding problems Banach spaces. AIME Transac- tions, 103:219–260, 1934. [YDA02] B. Yetan, Louis J. Durlofsky, and Khalid Aziz. Optimization of smart well control. London, United Kingdom, 4-7 November 2002. SPE Inter- national Thermal Operations and Heavy Oil Symposium and Interna- tional Horizontal Well Technology Conference, Society of Petroleum Engineers. [Yor91] Y .C. Yortsos. A theoretical analysis of vertical equilibrium. University of Southern California, 6-9 October 1991. SPE Annual Technical Con- ference and Exhibition, Dallas, TX. [YYS02] Z.M. Yang, Yannis C. Yortsos, and Dominique Salin. Asymptotic regimes in unstable miscible displacements in random porous media. Advances in Water Resources, (25):885–898, 2002. 129

Abstract (if available)

Abstract In modern reservoir management, optimizing the recovery efficiency of displacement processes is a fundamental problem with important economic ramifications. Given a well configuration, the optimum displacement efficiency may be achievable through controlling injection rates. Smart wells and intelligent completion have gained significant attention in recent years in the field of dynamic optimization of Enhanced Oil Recovery processes.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

A study of dispersive mixing and flow based lumping/delumping in gas injection processes

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

The effects of gravity and thermal gradient on the drying processes in porous media

PDF

Dynamics of combustion fronts in porous media

$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

PDF

Pore Network Simulations Of Immiscible Displacements In Heterogeneous And Anisotropic Porous Media

PDF

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

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

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

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

PDF

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

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

Efficient control optimization in subsurface flow systems with machine learning surrogate models

PDF

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

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

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

PDF

Application of data-driven modeling in basin-wide analysis of unconventional resources, including domain expertise

PDF

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

Asset Metadata

Creator Heidary-Fyrozjaee, Mohsen (author)

Core Title Control of displacement fronts in porous media by flow rate partitioning

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 04/10/2008

Defense Date 03/13/2008

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

Tag fluid fronts,front control,intelligent wells,OAI-PMH Harvest,porous media,rate partitioning

Language English

Advisor Yortsos, Yannis C. (committee chair), Ershaghi, Iraj (committee member), Sadhal, Satwindar S. (committee member)

Creator Email heidarif@usc.edu

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

Unique identifier UC1134424

Identifier etd-HeidaryFyrozjaee-20080410 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-53344 (legacy record id),usctheses-m1095 (legacy record id)

Legacy Identifier etd-HeidaryFyrozjaee-20080410.pdf

Dmrecord 53344

Document Type Dissertation

Rights Heidary-Fyrozjaee, Mohsen

Type texts

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

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu