History Matching

INSTITUTE OF PHYSICS PUBLISHING Inverse Problems 21 (2005) 565–590

INVERSE PROBLEMS

doi:10.1088/0266-5611/21/2/009

History matching problem in reservoir engineering using the propagation–backpropagation method Pedro González-Rodr´ıguez, Manuel Kindelan, Miguel Moscoso and Oliver Dorn Mathematics Department, Universidad Carlos III de Madrid, Avenida de la Universidad 30, Leganes 28911, Spain

Received 14 June 2004, in final form 29 November 2004 Published 23 February 2005 Online at stacks.iop.org/IP/21/565 Abstract This paper describes the application of the adjoint method to the history matching problem in reservoir engineering. The history matching problem consists in adjusting a set of parameters, in this case the permeability distribution, in order to match the data obtained with the simulator to the actual production data in the reservoir. Several numerical experiments are presented which show that our method is able to reconstruct the unknown permeability distribution in a reliable and efficient way from synthetic data provided by an independent numerical forward modelling code. An efficient and flexible regularization scheme is introduced as well, which stabilizes the inversion and enables the reservoir engineer to incorporate certain types of prior information into the final result. (Some figures in this article are in colour only in the electronic version)

1. Introduction History matching techniques are used in reservoir modelling for estimating unknown properties of a reservoir, such as porosity and permeability, from measured data. The data are usually taken at the production wells and might consist of pressure or flow data. Proper characterization of the reservoir heterogeneity is a crucial aspect of any optimal reservoir management strategy. It helps to better understand the reservoir behaviour so that its performance can be predicted and controlled with higher reliability. History matching can be carried out either manually (by a trial-and-error process), or automatically by computing a set of parameter values so as to minimize a prescribed cost function. Generally, the reservoir simulator uses a spatial grid, and the history matching process is used to determine the permeability and/or porosity in each grid cell. This can be done manually when a small number of parameters has to be retrieved from the observed production data. However, for large scale models, where more than 50 000 parameter values are unknown, this procedure is infeasible. Automatic strategies for history matching were initially based on the perturbation method. This method computes the single grid block sensitivity 0266-5611/05/020565+26$30.00 © 2005 IOP Publishing Ltd Printed in the UK

565

566

P González-Rodr´ıguez et al

coefficients in order to evaluate the change of an objective function under small perturbations of the permeability of each cell [27]. A direct method of obtaining these sensitivity coefficients requires a number of forward solutions per iteration equal to the number of reservoir parameters to be determined. In consequence, this approach is very time consuming for large geophysical models. A major breakthrough was achieved by Jacquard and Jain [26]. They used an electric circuit analogue of a reservoir, to compute the sensitivity coefficients for pressure to changes in the permeability map. The method was computationally efficient since it required only one simulation per observation point. It was based on a linear, single phase, reservoir model and it was later rigorously derived in terms of the diffusion equation [9]. The resulting optimization problem was solved by a linear programming procedure which computed the reservoir parameters including constraints. The work by Chavent et al represented another major breakthrough in automatic history matching [11]. They applied an adjoint technique from optimal-control theory to this application. For the same linear, single phase, reservoir model as used in [26] they derived a corresponding adjoint equation, and computed the sensitivity coefficients by just solving two PDEs (the equation for the reservoir model, and the adjoint equation) and some integrals involving the pressure and the adjoint variable. Secondary oil recovery techniques involve the simultaneous flow of up to three fluid phases. The water, oil and gas flows are governed by a set of coupled nonlinear partial differential equations and therefore the derivations in [9, 26] had to be extended. Several researchers have paid attention to this and other coupled inverse problems [8, 19, 32, 50]. The combination of spatially dependent parameters and nonlinearity results in a very challenging inverse problem. Wasserman et al [52] modified the original work by Chavent et al [11] to treat multiphase effects approximately with a ‘pseudo’ single-phase model. They applied the adjoint method directly to the set of ordinary differential equations resulting from the finite-difference discretized model. The adjoint method for the case of a two-phase, nonlinear, reservoir model was derived in [53] starting from the discretized equations. Based on a set of adjoint operation rules, Sun and Yeh [45] introduced a general procedure for deriving the adjoint state equations in a straightforward way for general coupled inverse problems. As examples, they applied their method to problems of saltwater intrusion and two-phase flow. There are other more recent techniques which have been applied successfully in reservoir engineering. We only want to mention two approaches in the following. The first one is an alternative streamline approach which has been proposed in [28, 47, 51]. A method is introduced in that work for estimating the sensitivity coefficients semi-analytically with one single forward simulation, which significantly reduces the numerical cost. This gives rise to a large and sparse, but ill-conditioned, linear system. With appropriate regularization this system can be solved efficiently [2, 25]. As a second approach we mention the use of geostatistical methods [7, 21, 29] as well as genetic algorithms [24, 40]. These methods generate different realizations of a geological model which is derived from well-log and seismic data, and select the realization that best matches the pressure or production history. Although accurate and fast simulators are available, the large number of realizations to be tested usually makes this procedure computationally quite demanding. Adjoint techniques, which we are focusing on in this paper, are particularly useful in large scale inverse problems where relatively few independent experiments can be performed for gathering data but many parameters need to be reconstructed. Since typically only one experiment is performed in history matching due to the simultaneous production process, the adjoint technique is therefore much faster in this application. We mention that adjoint techniques have been applied recently with great success also in other applications of medical

History matching problem in reservoir engineering

567

and geophysical imaging. See for example the articles [3, 15, 33, 49] dealing with examples from ultrasound tomography, optical tomography and microwave imaging with medical applications, or [16, 17, 23] describing applications in geophysical monitoring or prospecting. More general information regarding adjoint techniques for solving large scale inverse problems can be found for example in [35, 48]. One important difference between the history matching problem and these applications is that in the history matching problem the experimental setup cannot easily be changed in order to obtain independent data. As already mentioned, typically only one field experiment is available due to the production process. As a consequence, the information content of the data is quite low, and the history matching problem is highly ill-posed. Therefore, when developing reconstruction schemes for this application, an important component will be the incorporation of efficient and flexible regularization tools, as well as the capability of incorporating as much ‘prior information’ as possible into the search for a good candidate for the parameter distribution. Prior information can consist of results which have been obtained with different physical experiments such as seismic imaging, core analysis, well-logs or general geological information. Therefore, we will put a strong emphasis on the regularization tools which we want to use in our reconstruction method. Good overviews of general regularization techniques for inverse problems can be found, for example, in [18, 22, 34, 48]. In this paper, we derive the adjoint formulation as needed here for a two-phase reservoir, starting directly from the partial differential equations of the nonlinear reservoir model. The resulting method is a propagation–backpropagation algorithm for history matching which is used to estimate the permeability distribution inside the reservoir based on the water flow extracted at a small number of production wells. Starting from a homogeneous permeability distribution, we calculate the difference between the computed and the given production data. These residuals are numerically backpropagated into the reservoir by solving the corresponding adjoint equation, and the correction to the guess is directly calculated from the result. The process is repeated iteratively until a convergence criterion is satisfied. As mentioned, a particular emphasis is put on the derivation and discussion of flexible and useful ‘regularization schemes’ for this method, which can be used for stabilizing the reconstruction process and, in addition, for incorporating certain types of prior information into the reconstruction. In contrast to the more standard approach of Tikhonov–Phillips regularization, we will not change our cost functional for the derivation, but instead will restrict the search for a minimum to a smaller subspace of functions with certain properties. This has the advantage that, upon convergence, only the data misfit has been minimized, which is our primary objective. For describing the flow dynamics in the reservoir (our basic ‘forward model’), we use a simplified Black-oil model [46]. In order to avoid the so-called inverse crime, the data are generated with a streamline method, while during the reconstruction process we use an independent IMPES method to solve the equations (see appendix B for details of these two different forward modelling codes). Here our approach differs slightly from the more standard approach of simply adding statistically random noise to the data which have been created with the same simulator. We believe that, by using a different simulator for creating data, we can incorporate some component of ‘systematic noise’ which might render the data more realistic. Purely random noise often can easily be removed (at least partially) by a simple filtering of the data, whereas systematic noise might be harder to cope with. The noise content of our data using this strategy is typically about 3%. This paper is organized as follows. In section 2 we give a short introduction to our setup of the reservoir characterization problem. In section 3 a mathematical description of the underlying flow equations is given. Section 4 introduces the inverse problem and gives a

568


theoretical derivation of the basic algorithm which we use for the inversion. The algorithm itself is outlined in section 5. In section 6 we propose an efficient regularization scheme for the method. Various numerical experiments are presented in section 7 which demonstrate the performance of the reconstruction method in different situations. In section 8 we draw some conclusions and indicate some directions for future research. In appendix A a proof of theorem 1 of the paper is given. In appendix B we briefly describe the two basic forward modelling codes (the IMPES method and the streamline method).

2. Description of the problem Our basic flow model (the ‘Black-oil model’) consists of partial differential equations which govern the unsteady flow of different fluid phases in the reservoir medium. The petroleum engineer often uses this model for understanding the dynamics of petroleum reservoirs and petroleum production in order to design an efficient operational strategy. We consider here the case of ‘secondary recovery’ where water is injected through several injection wells conveniently located in order to enhance oil production. The behaviour of the reservoir is usually influenced by many factors (permeability, porosity, relative permeability, etc) which are never known precisely. Therefore, the engineer uses the best values available and compares the results from the simulator with the field-recorded histories at the wells. Generally, this history matching will show discrepancies which the engineer tries to minimize by modifying one or several of the parameters which define the reservoir (permeability, porosity, etc). Once the simulator matches the production data, it is used to predict its future behaviour and to design alternative plans of operation of the reservoir. It should be pointed out that there is not a unique set of parameters to match production data, so that a perfect match does not mean that the reservoir is correctly characterized. In fact, if after a perfect match the simulator is used to predict future behaviour, the actual performance may differ from the predicted one so that it is necessary to monitor periodically the predicted versus the actual performance in order to update reservoir characterization. We use inverse problem techniques to optimize this history matching process. For simplicity, we consider the two-dimensional case and we select the distribution of rock permeability throughout the reservoir as the parameter that we try to adjust. Therefore, we assume that all the other parameters needed to define the reservoir behaviour are known. The direct problem refers to the resolution of the equations describing the flow within the reservoir assuming that the permeability distribution is known. For instance, the left side of figure 1 shows a five-spot layout with an injector well (o) at the centre (location x(i) 1 ) and four production wells (x) at the corners of a two-dimensional reservoir (locations (p) xj , j = 1, . . . , 4, j = 1 being the well at the upper left corner and numbered in the clockwise direction). Also shown is the real permeability distribution in milli-Darcys (mD). The water injected at the injection well displaces the oil in the reservoir towards the production wells. Time resolution of the flow equations provides the time evolution of pressure and flow at each point of the reservoir. Of particular interest is the oil and water flow rate at each production well. The right-hand side of figure 1 shows the time history of water flow rate (Qw ) at each well obtained by solving the direct problem. Note that water arrival occurs first at well four since it is surrounded by a region of high permeability. In the inverse problem we assume that the water flow rate at each well is known but the permeability distribution is unknown. We will start with an initial permeability guess (typically some constant distribution) and will iteratively modify the permeability distribution until the actual water production rate at each well is matched by the simulator.

History matching problem in reservoir engineering mD 1800

600 x(p) 2

x(p) 1

(i) x1

300

1200 1000

200

−1 6

1400

400

Qw × 10 (s )

1600

500

800

100 x3

200

400

1.2 0.6

Producer 1

1.2 0.6

Producer 2

1.2 0.6

Producer 3

2.4 1.2

Producer 4

15 10 5

(p)

x(p) 4

0 0

569

600 600

Injector 1 0

100 200 Time (Days)

(a)

300

(b)

Figure 1. (a) Permeability distribution: 5-spot example; (b) extracted water flows at the producers (shown in (a) with crosses) and injected water flow at the injector (shown in (a) with a circle).

3. The mathematical model In order to model the time evolution of the flow in a reservoir we use a simplification of the Black-oil model [46]. We consider two incompressible phases (water and oil) in a porous medium where the effects of gravity and capillary pressure are neglected. Then, the governing equations for the multiphase incompressible flow in a reservoir ⊂ Rn (n = 2, 3) can be written as −∇ · [T ∇p] = Q φ

∂Sw − ∇ · [Tw ∇p] = Qw ∂t

in

× [0, tf ]

(1)

in × [0, tf ],

(2)

where p(x, t) and Sw (x, t) are the unknowns of the problem which represent the pressure and the water saturation at position x and time t respectively. The water saturation Sw measures the volume fraction of water. φ(x) is the porosity, and T and Tw are the transmissibilities, which are known functions which depend linearly on the permeability K, the parameter to be reconstructed, and nonlinearly on Sw , Tw = K(x)

Krw (Sw ) , µw

To = K(x)

Kro (Sw ) , µo

T = Tw + T o .

(3)

In (3), Krw (Sw ), Kro (Sw ), µw and µo denote the relative permeabilities and the viscosities of each phase, respectively. Hereafter, the subindex w stands for ‘water’, while the subindex o t) define the total flow and the water flow at the wells, stands for ‘oil’. Q(x, t) and Qw (x, respectively. They are given by Np Ni (i) (p) (p) (i) pwbj − p δ x − xj + cT pwbj − p δ x − xj Q = cT j =1

Np Ni (i) (p) (p) (i) Qw = cT pwbj − p δ x − xj + cTw pwbj − p δ x − xj , j =1

(4)

j =1

j =1

(5)

570

P González-Rodr´ıguez et al (p)

where x(i) j , j = 1, . . . , Np , j , j = 1, . . . , Ni , denote the locations of the Ni injector wells, x (p)

(i) , pwbj are the imposed well bore denote the locations of the Np production wells, and pwb j pressures at the Ni injector wells and at the Np production wells, respectively. Here, c is a (i) (p) pwbj are larger (smaller) than the constant that depends on the well model [12]. Since pwb j reservoir pressure at the injector (production) wells, Q and Qw are positive (negative) at the injector (production) wells. Equation (2) is the conservation law for water in a porous medium and equation (1) is obtained by combining the conservation laws for water and oil in order to eliminate the time derivative term. It is assumed that the flow obeys Darcy’s law K(x)K rl (Sw ) ul (x, t) = − ∇p(x, t), l = w, o µl which defines the velocity of each phase in the medium. Equations (1) and (2) are solved with the following initial and boundary conditions:

x , 0) = Sw0 ( x) Sw (

in ,

(6)

x) p( x , 0) = p (

in ,

(7)

∇p · ν = 0

on

(8)

0

∂,

where ν is the outward unit normal to ∂. The boundary condition (8) implies no flux across the boundary. Equations (1)–(8) define the direct problem for the dynamic production history at the extraction wells. The properties of the porous media are given by K(x) and φ(x). The properties of the fluids are defined by µw , µo , Krw (Sw ) and Kro (Sw ). The well bore pressures (i,p) pwbj are known functions of time at the wells’ positions. 4. The inverse problem 4.1. A propagation–backpropagation inversion method In reservoir characterization, typically, one tries to estimate the permeability distribution by matching production data. The porosity distribution and the relative permeabilities are usually assumed known from core analysis [47, 51]. For the mathematical analysis of the above-described problem, we need to specify some function spaces that will be used in the formulation of the problem. These function spaces will also become important when deriving our regularization tools. We will use suitably chosen Hilbert spaces throughout our derivation. We will denote the space of permeability distributions K by P, which is defined by P = L2 (), equipped with the usual L2 inner product. In section 6 we will introduce an additional space Pˆ for the permeability distributions. The set (p) (p) (p) of measurement locations (‘well-locations’) will be denoted by + := x1 , x2 , . . . , xNp . At each of these positions, the water flow is measured during a time 0 t tf , such that the data space D is given by D = (L2 ([0, tf ]))Np . Our forward operator M is given as M : P −→ D

M[K] = Qw [K]|+ ×[0,tf ] ,

(9)

where Qw is obtained by solving the direct problem for a given permeability distribution K ˜ (equations (1)–(8)). For some guess K of the permeability, and given the measured data G (water flow rate) at the production wells, we can furthermore define the residual operator R : P −→ D by ˜ R[K] = M[K] − G. (10)


571

Equation (10) describes the mismatch between these physically measured data and the data corresponding to a guess K. In the inverse problem, we ideally want to find a permeability distribution Kˆ in P such that ˆ = 0. R[K] (11) ˜ are in the range of M. Using real This equation has a solution in the situation where the data G data, we cannot be sure whether this is the case. Therefore, we generalize our criterion for a solution. Defining the least-squares cost functional J (K) = 12 R(K)2L2 ,

(12)

we are searching for a minimizer of this cost functional, which can be zero in the situation ˜ ∈ range(M). Otherwise, it will be a positive number. A standard method for finding where G a minimizer of the cost functional (12) is to start a search with some initial guess K (0) , and to find descent directions of (12) in each step of an iterative scheme. Popular choices for descent directions are for example the gradient direction, conjugate gradient directions, Newton- or quasi-Newton directions (see, for example, [13, 36, 48] for details). We will, in the following, derive update directions for minimizing (12) using formulation (11). These update directions will have the useful property that they can easily be generalized in order to incorporate efficient regularization schemes in our algorithm. In order to find an ‘update’ (or ‘correction’) δK for our permeability K we linearize the nonlinear operator R (assuming that this linearized operator R [K] exists and is well-defined) and write R[K + δK] = R[K] + R [K]δK + O(δK2 ).

(13)

The linearized operator R [K] is often called the Frechet derivative of R at K. (See for example [14] and references therein for some formal derivations of Frechet derivatives in different applications.) It is also closely related to the ‘sensitivity functions’ of the parameter profile with respect to the data. Using (13) we want to look for a correction δK such that R[K + δK] = 0. Neglecting terms of order O(δK2 ) in (13), this amounts to solving R [K]δK = −R[K].

(14)

Certainly, due to the ill-posedness of our problem, this equation needs to be handled with care. Treated as an ill-posed linear inverse problem, a classical solution of (14) will be the minimum-norm solution δKMN = −R [K]∗ (R [K]R [K]∗ )−1 R[K],

(15)

where R [K]∗ is the adjoint operator of R [K] with respect to our chosen spaces P and D [34]. In applications with very few data, this form has the useful property that it avoids contributions in the solution which are in the (often non-empty) null-space of the (linearized) forward operator R [K]. Using (13) it can be verified by direct calculation that (16) J (K + ωδKMN ) = J (K) − ωR(K)2D + O δKMN 2P such that (15) is also a descent direction of the least-squares cost functional (12). In our application the operator C = (R [K]R [K]∗ )−1 is very ill-conditioned, such that a regularized version needs to be used. This can be, for example, Cˆ = (R [K]R [K]∗ + λI )−1 , where λ is some regularization parameter and I is the identity operator. Unfortunately, in practice, both C and Cˆ are very expensive to calculate and to apply to the residuals R. Typically, a direct calculation of the operator Cˆ would require us to solve as many forward and adjoint problems as we have independent data values. Iterative schemes (such as the gradient

572


or conjugate gradient method) for applying this operator to a given vector are possible as well, but usually they converge only slowly. In addition, we also have the possibility of just using a rough approximation of it in each step of the inversion. This is much easier and faster to do and might yield good results. We will investigate these possibilities in our future research. When using a very large regularization parameter λ, the contribution of R [K]R [K]∗ can be neglected and we end up with essentially (i.e. up to the scaling factor λ−1 ) calculating δK = −R [K]∗ R(K).

(17)

For this update direction we have

J (K + ωδK) = J (K) − ωR [K]∗ R(K)2P + O δK2P

(18)

such that it is also a descent direction for (12). We will use this update direction throughout the paper, with some important modifications described below. Our goal will be to derive and test efficient schemes for applying the adjoint linearized residual operator to given data (the basic propagation–backpropagation scheme), and moreover to derive and evaluate a new regularization scheme for this backpropagation technique. We will emphasize, however, that the propagation–backpropagation scheme as well as the regularization scheme will also be applicable directly to the update directions given in (15), which will be treated in our future work. A standard method for deriving regularization schemes is to explicitly try to minimize a cost functional which incorporates, in addition to the usual least-squares data misfit, a Tikhonov–Phillips regularization term: 1 η R(K)2D + K2α , (19) 2 2 where η > 0 is the regularization parameter and .α indicates some norm or semi-norm, e.g. Kα = ∇KL2 [20, 48]. Using this approach, the cost functional is changed significantly with the goal of obtaining in a stable way a global minimizer. We do not want to take this route, but prefer instead to keep working with the original least-squares cost functional (12) which only involves the data fit. We will minimize this cost functional by restricting the search to elements of a smaller function space, which is an alternative form of regularization. The regularization scheme will be derived and discussed in detail in section 6. In the following, we will present the basic structure of our inversion method, and will derive practical ways of applying the adjoint linearized residual operator R [K]∗ to vectors R in the data space D. This will lead us to the propagation–backpropagation technique which is applied in this paper. JT P (K) =

4.2. Computation of the operator R [K]: the linearized problem Let us consider a small perturbation δK in the permeability distribution K that leads to small perturbations W and q in the saturation and the pressure, respectively. Here we assume that the pressure remains nearly unchanged so that ∇q is negligible. This is so because the pressure is a smooth function compared to the saturation. Using a heuristic approach to derive an expression for R , we introduce K + δK and Sw + W in (2) and we neglect second-order terms. Then, W solves the initial value problem ∂Tw ∂Qw δK δK ∂W −∇ · Qw + ∇ · Tw ∇p W ∇p − W = in (20) φ ∂t ∂Sw ∂Sw K K W (x, 0) = 0

in ,

(21)


573

where Sw and p are the solutions of (1)–(8). From the value of W we derive the linearized response of the data to a perturbation δK in the permeability distribution, which is given by ∂Qw W . (22) R [K]δK = ∂Sw + ×[0,tf ] 4.3. Computation of the operator R [K]∗ : the adjoint problem Here, we derive an expression for the adjoint operator R [K]∗ applied to a function ρ in the data space. The operator R [K]∗ is defined by R [K]δK, ρ D = δK, R [K]∗ ρ P .

(23)

We assume that the inner products in the parameter space P and in the data space D are given by f, g D =

Np j =1

tf

fj gj dt,

A, B P =

0

AB dx,

(24)

where fj = f xpj , t and gj = g xpj , t , j = 1, . . . , Np , are time functions defined at the production well positions xpj . We formulate the basic result of this section in the following theorem, which is derived in appendix A. Theorem 1. Let ρ ∈ D be an arbitrary function in the data space. Then R [K]∗ ρ is given by tf

Tw 1 ∗ (25) ∇p∇z − z Qw dt, R [K] ρ = K K 0 where z is the solution of the adjoint equation   Np ∂Qw ∂z ∂Tw (p) −φ + ∇p∇z − z − ρδ x − xj  =0 ∂t ∂Sw ∂Sw j =1 z(x, tf ) = 0

in ,

in

(26) (27)

and Sw and p are the solutions of (1)–(8). Note that Qw is nonzero only at the well locations. Therefore, when we assume in the mathematical derivation of the theorem that the permeability is known directly at the wells (a realistic assumption), the second term in (25) disappears and we only have to evaluate the first term in order to calculate the update in the rest of the domain . This will be the approach we use in our numerical reconstructions. Note that, as typical for the adjoint scheme, the system (26), (27) physically models some kind of backpropagation with respect to the linearized forward model. The residuals are applied at the production wells as artificial injectors, and backpropagated backward in time (note the minus sign in front of the time derivative in (26) and the prescribed final value conditions in (27), compared to a plus sign in (20) and initial values in (21)) and in space by the system (26), (27). Equation (25) uses these backpropagated fields to extract an update direction by combining forward and adjoint fields at each location.

574


5. The basic algorithm The basic inversion algorithm (without regularization) can be summarized as follows. Assume that the nth approximation K (n) to the true permeability distribution has been obtained. Then: ˜ on + × [0, tf ] by solving the direct (i) Compute the residual R[K (n) ] = M[K (n) ] − G problem (1)–(8). (ii) Apply ρ = R[K (n) ], as computed in (i), at the production wells and ‘backpropagate’ by solving the adjoint problem (26), (27) for this ρ. (iii) Combine the results of (i) and (ii) for calculating the update δK (n) from (25). Let p be the pressure as calculated in (i), and let z be the solution of (26), (27) as calculated in (ii). Since we assume that K is known at the well locations, and Qw is zero in the rest of the domain, (25) simplifies to tf Tw δK (n) = ∇p∇z dt. (28) K 0 (iv) The new approximation to the true permeability distribution is K (n+1) = K (n) + ωδK (n) , where ω is some step-length to be chosen properly. An efficient scheme for practically choosing ω (a variant of a line-search technique) will be explained in section 7.1.2. (v) Compute the residual R[K (n+1) ] and proceed with step (i) until some convergence criterion is achieved. 6. Regularization and smoothing 6.1. Smoothing with function spaces We have presented above the basic algorithm which recovers L2 functions of permeability from given data such that the misfit in the data is minimized. This procedure does not incorporate ˆ explicit regularization (except of the stabilizing procedure incorporated in the operator C). In some situations, it might be necessary or desirable to restrict the search for permeability functions to a smaller subset of L2 , for example of smoothly varying functions. This might be so in order to regularize the reconstruction algorithm, or in order to take into account some prior information or assumptions on the solution we are looking for. For example, the reservoir engineer might know or assume that the permeability distribution in some region is fairly smoothly varying. Or, he might only have very few data available for the inversion, so that he wants to select a smoothly varying profile as a regularized form of the reconstructed permeability distribution. This can easily be done in our framework. Instead of looking for permeability distributions in L2 (), let us assume now that we require the permeability to be an element of the smaller subspace H1 () := {m ∈ L2 (), ∂i m ∈ L2 ()

for

i = 1, 2, 3}.

This Sobolev space is usually equipped with the standard norm 1/2 m1,1 := m2L2 + ∇m2L2 and the standard inner product m1 , m2 1,1 := m1 , m2 L2 + ∇m1 , ∇m2 L2 . For reasons explained below, we will instead prefer to work with the equivalent norm 1/2 , α, β > 0 mα,β := αm2L2 + β∇m2L2


575

and its associated inner product m1 , m2 α,β := α m1 , m2 L2 + β ∇m1 , ∇m2 L2 . A proper choice of the weighting parameters α and β will allow us to steer the regularization properties of our algorithm in an efficient and predictable way. Let us denote the new parameter space H1 (), when equipped with the weighted norm .α,β , by Pˆ . When using this modified space in our algorithm, we also have to adjust the operators acting on it, in particular the adjoint of the linearized residual operator. This operator is now required to map from the data space D into Pˆ . Moreover, the minimum norm solution of (14) is now taken with respect to the weighted norm .α,β , which clearly gives us a different candidate. The necessary adjustments for our algorithm can be done as follows. Denote as before by R [K]∗ ζ the image of ζ ∈ D under application of the adjoint linearized residual operator as calculated in section 4.3, considered as an operator mapping from D into P = L2 (). Denote furthermore by R [K]◦ ζ its image under the adjoint linearized residual operator with respect to the newly defined weighted inner product, mapping into the smaller space Pˆ . With a straightforward calculation, using the definitions of the two adjoint operators R [K]x, ζ D = x, R [K]∗ ζ P = x, R [K]◦ ζ Pˆ , (29) it follows that R [K]◦ ζ = (αI − β )−1 R [K]∗ ζ, (30) where we supplement the inverted differential operator (αI − β )−1 by the boundary condition ∇(R [K]◦ ζ ) · n = 0 on ∂. The symbol I stands for the identity, and stands for the Laplacian operator. Equation (30) can easily be derived by applying Green’s formula to the right-hand side equality in (29). In practice, the ratio γ = β/α (which can be considered as a ‘regularization parameter’) is an indicator for the ‘smoothing properties’ of our scheme. The larger this ratio, the more weight is put on minimizing the derivatives of our solution. Therefore, by properly choosing this ratio, we can steer the smoothness properties of our final reconstruction to a certain degree. In our numerical experiments, we will choose this ratio once, when starting the algorithm, and keep it fixed during the iterations. The other free parameter, say α, will be chosen in each individual step to scale the update properly. In our numerical experiments, we choose α such that R [K]◦ ζ L2 = R [K]∗ ζ L2 is satisfied for the current update. This possibility of scaling the updates is the main reason for keeping the parameter α throughout the calculations instead of simply putting it to 1 right at the beginning. When testing and comparing the performance of different regularization parameters γ it is practically useful (in particular for the line-search method) that the order of magnitude of the calculated values of R [K]◦ ζ does not depend too much on γ . Note also that the new search directions using this modified adjoint operator are still descent directions for the least-squares cost functional (12), as can be verified easily by replacing P by Pˆ in (16) and (18). Practically, the scheme is implemented as follows: γ is fixed regularization parameter Define = R [K]∗ ζ. Solve (I − γ )ϕ = , ∇ϕ · n = 0 on ∂. (31) ϕ ϕ Define = α with α = (such that = ) Then we have (αI − β ) = , ∇ · n = 0 on ∂, with β = αγ . Put R [K]◦ ζ = .

576


We mention that applying this regularization scheme amounts to applying the postprocessing operator (αI − β )−1 to the updates calculated in the previous ‘unregularized’ scheme. Therefore, the effect of the regularization is similar to filtering the updates with a carefully designed (iteration-dependent) filtering operator. In the following, we want to give an interesting additional interpretation of this regularization scheme. Define the cost functional a b c Jˆ () = 2L2 + ∇2L2 + − 2L2 (32) 2 2 2 with = R [K]∗ ζ . Here, the third term penalizes the misfit between the unregularized update direction δK = R [K]∗ ζ and the new candidate , whereas the first two terms penalize roughness of . The gradient direction for this cost functional is [(a + c)I − b ] − c (where the Laplace operator is again understood to be accompanied by the boundary condition ∇ · n = 0 on ∂). Therefore, a necessary condition for the minimum can be stated as [(a + c)I − b ] = c.

(33)

Choosing c = 1, b = β 0 and a = α − 1 0 this amounts to calculating = (αI − β )−1 ,

(34)

which is equivalent to (30). Therefore, applying function space regularization as described above can be interpreted as minimizing the cost functional (32) with specifically chosen parameters a, b and c. 6.2. Smoothing with the heat kernel An alternative (and slightly more ad hoc) approach to regularization and smoothing is to use well-known concepts from image processing. Denote = R [K]∗ ζ for some residual vector ζ ∈ D. Then, we can convolve these (unregularized) updates calculated by (15) with a Gaussian kernel of variance σ > 0

|x|2 1 exp − (35) fσ (x) = 4π σ 4σ which produces the smoothed update ˆ = fσ ∗ = fσ (x − y)(y) dy.

(36)

Practically, this can be done by solving the initial value problem for the heat equation vt − v = 0 v(0) =

for t ∈ [0, τ ]

(37)

on with τ = σ and with suitably chosen boundary conditions, and putting ˆ = v(τ ).

(38)

Here, the smoothing time τ can be considered as a regularization parameter: for τ = 0 no regularization takes place, whereas with increasing τ the updates become increasingly smoothed. In image processing, this procedure (in a generalized form) is sometimes referred to as ‘defining a scale space’, with τ > 0 being the scale. Although the scheme described above is only ad hoc, it turns out that a similar scheme can be derived in a more rigorous way from the cost functional (32). This will be demonstrated in the following.


577

Let us assume that we again want to minimize (32) for the above given choice of the parameters c = 1, b = β 0 and a = α − 1 0. Now we want to use a gradient method for finding the minimum, starting with the initial guess (0) = . Using the gradient direction for (32) derived above, we get the iteration rule (n+1) = (n) − θ [(αI − β )(n) − ]

(39)

where θ is the (fixed) step-size in a given step n of the iteration. This can be written as (n+1) − (n) (40) = β (n) + ( − α(n) ) θ which is just one step of a finite-difference time discretization of a modified heat equation (37) vt − β v = ( − αv) v(0) =

for t ∈ [0, τ ]

(41)

with additional time-dependent heating source − αv and fixed time-step δt = θ . The choice of the free parameter α in this iteration can be done according to the same lines as described in the algorithm in section 6.1. This strategy leads to the practical algorithm: γ Define for ϕ (n+1) end for Define Then we have Put

is fixed regularization parameter = R [K]∗ ζ. n = 1, . . . , N : = ϕ (n) + θ [γ ϕ (n) + ( − ϕ (n) )] (Upon convergence we have ϕ (N) = (I − γ )−1 ) (N)

= ϕα with α = (αI − β ) = , R [K]◦ ζ = .

(42)

ϕ (N)

(such that = ) ∇ · n = 0 on ∂, with β = αγ .

It turns out that this procedure has good regularization properties even if we choose N fairly small (e.g. N = 5 or N = 10), which amounts to choosing the regularization time τ in (41) small. In that case, we will not exactly calculate R [K]◦ ζ = (αI − β )−1 , but we will have only a mildly smoothed form of = R [K]∗ ζ . Moreover, for small values of N the additional heating term − (n) is typically small (recall that we start the iteration with (0) = ), such that we can safely neglect it and arrive at the regularization scheme by the heat kernel (37). Although this yields only a very crude approximation to R [K]◦ ζ , the results achieved with this scheme are usually quite satisfactory when used in each step of our iterative scheme for solving the inverse problem. We complete this section by mentioning that, throughout this paper, we will always use the scheme described in section 6.1 (regularization by function spaces) for the regularization. However, most of our calculations have also been tested with the alternative scheme described above using a small value of N, which also gave us good results. 7. Numerical experiments In our numerical experiments we use two different numerical schemes for the forward modelling: the IMPES method and the streamline method. Both are described in appendix B. For the numerical solution of the adjoint problem (26), (27) we have implemented an explicit upwind finite difference scheme. We model a reservoir of 600 × 600 m2 which is discretized by a 25 × 25 uniform spatial grid. The typical time-step in the discretization is between 2 h and one day, and the reservoir is monitored for a duration between 50 and 100 days. When using the streamline method, each cell of the grid is intersected by at least one

578

P González-Rodr´ıguez et al 1 K rw K

0.9

ro

0.8 0.7

Kr

0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

Sw

Figure 2. Relative permeabilities versus saturation as used in the numerical experiments.

streamline. The parameters used to describe the reservoir are φ = 0.213, µw = 8.2 × 10−4 Pa s, µo = 7.9 × 10−4 Pa s. As boundary and initial conditions for the forward (p) (i) = 3500 psi and pwb = 2000 psi ∀p at modelling we use Sw0 (x) = 0.187, p0 = 3000 psi, pwb all wells. The relative permeabilities are nonlinear functions of saturation which are shown in figure 2. For the value c in formulae (4) and (5) we use c = dx dy2πln( rd ) (see also [12, 37]) 1

rw

where rd = 0.14[dx 2 + dy 2 ] 2 and rw = 0.108 m is the well radius. Here dx and dy denote the discretization lengths of the grid in the x and y directions, respectively. In the following we show permeability estimations for several examples with two different reference permeability distributions. In the first case we estimate a smooth permeability distribution, whereas in the second example we estimate a permeability distribution with sharp discontinuities. In both cases we use four injector wells and nine production wells arranged as an array of so-called ‘five-spot patterns’. 7.1. First example: reconstruction of a smoothly varying function We first investigate the reconstruction of the smoothly varying permeability distribution which ˜ =Q ˜ w at the producers is plotted in was already shown in figure 1. The extracted water flow G figure 3(b) with solid lines. In our initial experiment (sections 7.1.1 and 7.1.2), we compute the synthetic data by using the IMPES method and perform the reconstruction by using the same forward modelling code (the IMPES method). In the later section 7.1.3, and in all further numerical experiments presented in section 7.2, we will use instead data which have been generated by the streamline method, whereas the reconstruction is done by using the IMPES method. By doing this, we can avoid the so-called ‘inverse crime’, and also can investigate the performance of the algorithm in the case of very noisy data. 7.1.1. Basic reconstruction without explicit regularization. The well configuration and the estimated permeability distribution at the final iteration are shown in figure 3(a). Our initial model, a uniform permeability distribution corresponding to the mean permeability value (1400 mD) serves to generate the initial data M[K (0) ], plotted in figure 3(b) with dot-dashed


600

mD 1800

400

1400

C 200

1000

(b)

(12,12) (12,300) (12,588)

Qw × 106 (s −1)

(a)

579

(300,12) (300,300) (300,588) (588,12) (588,300)

A 0 0 B

2.4 0

600 200 D

400

600

(588,588) 20

40 Time (Days)

60

(d) 8

(c) 2000

Flow residual (a.u.)

A 2000

mD

B

2000 0 2000

C

D 0

200

400

600

6

4

2

0

20

40 Iterations

60

80

Figure 3. Synthetic example with four injectors and nine producers: (a) estimated permeability; (b) extracted water flows at the producers (solid lines, dot-dashed lines and dashed lines represent the extracted water flow for the reference model, for the initial model and for the estimated model, respectively); (c) cross sections (solid lines and dashed lines represent the reference permeability model and the estimated model, respectively); (d) residuals at each iteration. The reference permeability is shown in figure 1(a).

lines. The calculated extracted flows at the producers after 100 iterations are shown with dashed lines. The agreement is very good, indicating the validity of our approach. Note that the calculated responses match the ‘true’ production data over almost the whole time history. The small differences might be due to pressure effects that are not taken into account. Figure 3(d) shows the root mean square error in the measurements at each iteration. It decreases monotonically until becoming stationary (up to small fluctuations) at a small value, after about 80 iterations. For comparison purposes, we plot in figure 3(c) several cross sections of the real permeability (solid lines) and of the estimated permeability (dashed lines). Note that no explicit regularization has been applied in this numerical experiment. Accordingly, the reconstruction does not look very smooth. Moreover, a block structure can be observed in the reconstruction which is due to the use of five-spot patterns in the experimental setup. Each block reflects the flow from one injector to the corresponding neighbouring producer. In this example we have chosen the relaxation parameter ω to be constant. Usually, its value is determined by trial and error. We have chosen it such that the maximum of the relative update δK (0) /K (0) in the first iteration is 1% of the constant permeability value of the initial guess.

580

P González-Rodr´ıguez et al mD

600

1400

400

1000

200

0 0

(b) 8 1800

600 200

400

600


(a)

6

4

2

0

20

40 Iterations

60

80

Figure 4. Same as figure 3 but with variable ω and with regularization: (a) estimated permeability; (b) residuals at each iteration with variable ω (solid line) and with constant ω (dashed line).

7.1.2. Reconstruction with regularization and with a line-search variant. In the next example, we allow ω to be variable in order to speed up the convergence. We use a variant of a line search technique for finding a suitable update. Our goal is to take a maximally possible step which still decreases the mismatch between calculated and measured data. Therefore, step (iv) in the algorithm is replaced by K (n+1) = K (n) + ω(n) δK (n) , n = 0, 1, 2, . . .. For starting this scheme, ω(0) is chosen as before. We introduce a suitable factor γ > 1 for manipulating the relaxation parameter ω. If the residual R[K (1) ] decreases when applying the chosen value for ω(0) , we accept this step, and set ω(1) = γ ω(0) for the next iteration. Otherwise, we reject the update and try again with ω(0) → ω(0) /γ . In this case, we continue reducing ω(0) until we find an ω(0) for which the newly calculated residual R[K (1) ] is reduced. This step is then accepted, and we continue with the next iteration searching for a suitable ω(1) , assuming as initial value for the search the final value for ω(0) . The algorithm stops if at least 15 consecutive trials are rejected, i.e. if choosing ω(n+1) = γ −15 ω(n) does not lead to a decrease of the residual in some iteration step n. In our numerical experiments this scheme has shown to reduce the total computational time for the reconstruction significantly. We note here that Bulishev et al [6] have used a similar procedure for choosing the length of the iterative step in a gradient method. Alternative line search techniques are possible, and will be investigated in our future research. Figure 4(a) shows the estimated permeability after 20 of these modified iterations when using the value γ = 1.5. Observe in figure 4(b) that now the residual decreases much faster. In addition to the line search technique, we have also applied to this reconstruction the regularization scheme explained in section 6.1. Using here the same well configuration and reference permeability distribution as in figure 3, it is apparent that the estimated permeability distribution is now smoother than the one obtained in figure 3(a). 7.1.3. Reconstruction from data created by an independent STREAMLINE method. In the numerical experiments shown so far we have used the same forward modelling code for creating the data as we have used for the reconstruction task, without any additional noise added to these data. Certainly, in order to really evaluate the stability and practical usefulness of a reconstruction algorithm, it is necessary to apply the algorithm also to data corrupted with different types of noise. Noisy data with only random Gaussian noise can be created by simply adding random numbers of a certain magnitude to the previously considered noiseless data.

History matching problem in reservoir engineering (a)

mD

600

400

8


6

600

1400

400

600 0

200

400

0

600

1800

1400

1000

600 0

200

400

600

(d) 20

4

2

0

mD

200

Error in permeability (mD)

(c)

(b) 1800

1000

200

0

581

5

10 15 Iterations

20

25

15

10

5

0

5

10 15 Iterations

20

25

Figure 5. Reconstruction of the smooth permeability distribution shown in figure 1 from data created with an independent streamline method. Regularization is applied with a value β/α = c0 . Top left: reconstruction after 10 iterations; top right: reconstruction after 25 iterations. Bottom left: evolution of residuals and actual noise level of data. Bottom right: evolution of L2 -error of reconstruction.

In order to include some possible systematic noise component into our investigation, we have decided not to take this route, but instead to use a completely independent simulation technique for calculating our data. For this purpose, we use a streamline method that is not designed to produce data with the highest accuracy. In fact, we have measured that the disagreement with data obtained with the IMPES method is approximately 3%. In summary, the streamline created data are used as the input for our reconstruction scheme. Then, we use the IMPES method throughout the iterations for carrying out the reconstruction task. We believe that this is also an excellent test for evaluating the expected performance of the code when applied to real data. In order to regularize the scheme, we select a fixed ratio of β/α = c0 throughout the iterations and apply the postprocessing operator (αI − β )−1 in each update as described before. The results after 10 and 25 iteration steps are shown on the top left and top right images of figure 5. In the bottom left image of this figure we have plotted the evolution of the norm of the residuals during the iteration. It can be observed that the residuals are continuously decreasing, as expected from the algorithm. In order to get an indication of the actual level of noise contained in the data, we have also plotted the L2 -misfit between these data and the data we would get using the IMPES method (the horizontal line in the figure). For the value c0 chosen here for our reconstruction we observe some form of semiconvergence of the algorithm, as often occurs in iterative methods applied to noisy data.

582


In order to illustrate this, we show in the bottom right image of figure 5 the evolution of the L2 -error of the reconstruction compared to the reference permeability distribution. During the early iterations of the method this error obviously decreases. Approximately at a point where the residuals hit the noise level, this error starts to increase again, although the residuals are still decreasing. Our interpretation of this phenomenon is that at this point the algorithm starts fitting the noise rather than the data. The top left image of figure 5 shows the reconstruction at the point where the L2 -norm of the image error turns from decreasing to increasing (which can be considered as the ‘optimal stopping point for the algorithm’), and the top right image shows the final reconstruction which we get when keeping the algorithm running until the residuals do not decrease anymore. We see that the image obtained at the ‘optimal stopping point’ is a very good reconstruction of the reference permeability distribution, whereas in the image obtained at the final iteration, noise-related artefacts have degraded the reconstruction. Certainly, in real applications the noise level of the data can only be estimated, and the optimal stopping point must be chosen according to some appropriate criterion (e.g. the ‘Morozov discrepancy principle’). Moreover, it must be taken into account that the evolution of the L2 -error of the reconstruction itself needs to be considered with care. Since the problem at hand is vastly underdetermined, it does not really make sense at all to speak about ‘the correct solution’ without any further prior information available. We want to point out the role of prior information in our scheme in the following numerical experiment. In order to investigate this semi-convergence property further, we have increased in this experiment the ratio β/α to 10 times its previous value, i.e. β/α = 10c0 . The results of the reconstructions are displayed in figure 6. Now the residuals decrease about as much as in the previous situation. However, the L2 -norm of the error in the reconstruction (bottom right image of figure 6) does not show the semi-convergence behaviour here. The reconstruction error becomes practically stationary at some stage of the algorithm and we get a good final reconstruction in a stable way. The reason for this improvement is that now we have increased the smoothing property of the inversion algorithm by increasing the ratio β/α. Since the reference permeability distribution is fairly smooth, this amounts to putting more weight on correct prior information, such that we arrive at a better reconstruction than before. 7.2. Second example: reconstruction of a function with sharp discontinuities In contrast to the previous situation where we tried to reconstruct a very smooth unknown permeability distribution from the data, we will now investigate the performance of the numerical scheme when the true permeability distribution is highly discontinuous with large sharp edges of high contrast. The experimental setup is the same as before. Figure 7 shows the true permeability distribution and the well locations. The data have been calculated with the streamline method, and the numerical reconstruction is done by using the IMPES method for forward modelling. In this reconstruction, our initial guess is a constant permeability distribution of 1400 mD. Figure 8 shows the results of the experiment when using the ratio β/α = c0 (the same value as before) for the regularization at two different iteration steps of the reconstruction. The lower left image in figure 8 shows the evolution of the residual during the reconstruction, together with the theoretical ‘noise level’ of the data. The semiconvergence behaviour of the reconstruction is hardly noticeable here. One reason might be that the original profile is not as smooth as it was in the previous example. Certainly, also here we can display the results after an intermediate stage of the reconstruction, here after 19 iterations (top left image of figure 8), and at the final stage after about 70 iterations (top right image) when the iteration

History matching problem in reservoir engineering (a)

mD

600

400

600

1400

400

600 200

400

1000

600 200

400

10

20 Iterations

600

(d) 20

6

4

2

0

1800

1400

0 0

600

8

mD

200

Error in permeability (mD)


(c)

(b) 1800

1000

200

0 0

583

10

20 Iterations

15

10

5

0

30

30

Figure 6. Reconstruction of the smooth profile shown in figure 1 from data created with an independent streamline method. Regularization is applied with a value β/α = 10c0 . Top left: reconstruction after 15 iterations; top right: after 37 iterations. Bottom left: evolution of residuals and actual noise level of data. Bottom right: evolution of L2 -error of reconstruction.

mD

600

1400

1000

400

600 200 200 0 0

200

400

600

Figure 7. Reference permeability distribution for the second example. The figure also shows the location of the four injectors (◦) and the nine producers (×).

stops according to our stopping criterion. It is difficult to say which one could be selected as the ‘better’ reconstruction of the original profile.

584


(a) 600

mD

2000

(b) 600

1500

400

mD

1500

400

1000

1000

200

200

500

500 0

0

200

400

0

600

(c) 15

10

5

0

0

200

400

600

(d)

Error in permeabilty (mD)


2000

10

20

30 40 Iterations

50

60

25 20 15 10 5 0

10

20

30 40 Iterations

50

60

70

Figure 8. Reconstruction of the non-smooth profile shown in figure 7 from data created with an independent streamline method. Regularization is applied with a value β/α = c0 . Top left: reconstruction after 19 iterations; top right: after 70 iterations. Bottom left: evolution of residuals and actual noise level of data. Bottom right: evolution of L2 -error of reconstruction.

In the same way as in the previous example we want to compare the behaviour of the code when we change the ‘regularization parameter’ to β/α = 10c0 . As mentioned before, this amounts to putting more weight on the smoothing properties of the algorithm. Figure 9 shows the results in this situation. Again, the lower left image shows the evolution of the residual norm during the iterations, together with the theoretical ‘noise level’ of the data. The bottom right image shows the evolution of the L2 -error in the reconstruction. In contrast to the previous situation where we tried to recover the smooth profile, now the L2 -error of the image does not stabilize when using the higher value for β/α. Instead, we see in this situation a semi-convergence behaviour for this larger value for β/α. Therefore, in a way, the behaviour of the algorithm seems to be reversed compared to the previous smooth example. This makes sense since now our reference permeability distribution is highly discontinuous, and enforcing a smooth reconstruction is not the right way to proceed if this prior knowledge is available. Using the smaller value for β/α = c0 could be the better choice in this situation. We have displayed in the top left image of figure 9 the reconstruction at the point where the L2 -error in the image changes from decreasing to increasing (after 17 iterations), and in the top right image the final reconstruction of the algorithm. We want to mention at this point that, when it is known a priori that the permeability distribution is highly discontinuous, or if it is desired for different reasons to reconstruct a discontinuous profile, different regularization methods might be particularly useful here.

History matching problem in reservoir engineering (a) 600

mD

585 (b) 2000

1500

400

mD

600

2000

1500

400

1000

1000

200

200

500 0

0

200

400

0

200

400

600

(d) 28 Error in permeability (mD)


0

600

(c) 15

10

5

0

500

10

20 Iterations

30

26 24 22 20 18

5

10

15 20 Iterations

25

30

35

Figure 9. Reconstruction of the non-smooth profile shown in figure 7 from data created with an independent streamline method. Regularization is applied with a value β/α = 10c0 . Top left: reconstruction after 17 iterations; top right: after 36 iterations. Bottom left: evolution of residuals and actual noise level of data. Bottom right: evolution of L2 -error of reconstruction.

In particular the use of total variation regularization [1, 10], Huber norm regularization [4], the use of a Mumford Shah functional [41] or a level set based shape reconstruction scheme [17, 31, 43] would be appropriate tools to use. See also [39] for a description (and alternative treatment) of such a situation in reservoir engineering. We plan to investigate some of these regularization schemes in our future research in order to compare their performances with the reconstructions shown here. 8. Conclusions We have introduced in this paper a new reconstruction scheme for the history matching problem which is based on an adjoint method. Starting from some initial guess, this iterative scheme calculates successive corrections for the permeability distribution by numerically backpropagating the residuals corresponding to the latest best guess. This so-called ‘adjoint scheme’ avoids the usually expensive calculation of large sensitivity matrices in each step of the inversion. A line-search variant is used to find optimal steps for the updates. In particular, new efficient regularization schemes have been derived and tested which enable the reservoir engineer to incorporate prior information into the inversion which stabilizes the reconstruction process.

586


The numerical experiments presented in this paper show that this algorithm is able to reconstruct unknown permeability distributions in a fast and efficient way from relatively few data which have been created by an independent forward modelling scheme. In the case of very noisy data we observe typical semi-convergence behaviour. However, it can be controlled to a certain degree by introducing correct prior information into the reconstruction. In particular, for a smooth profile this leads to a significant improvement of the results. The presented scheme can easily be extended to a more realistic 3D situation. Since it is based on repeated forward and adjoint calculations of the reservoir equations, all that is needed are reliable forward and adjoint solvers in 3D. These solvers can be integrated into the inversion code in a black-box fashion, such that the reservoir engineer can easily switch back and forth between different forward modelling codes if necessary. The same holds true for the two regularization schemes presented in this paper. Since they have the form of easy to implement post-processing tools (applied in each step of the iterative inversion), they can be used as black-box filters which can be modified and exchanged independently from the rest of the algorithm in order to incorporate correct or assumed prior information into the inversion. In our future work we plan to investigate alternative regularization schemes, such as level-set techniques, which are better suited for the presence of sharp discontinuities. Acknowledgments We thank the anonymous referees for useful comments which improved the paper. Funding for this work was provided by the Dirección de Tecnolog´ıa y Soporte Técnico, Repsol-YPF. OD acknowledges support by the Ministerio de Ciencia y Tecnologia, Spain, through a Ramón y Cajal grant. Appendix A. Proof of theorem 1 In this appendix we prove theorem 1. Let z be a function defined in × [0, tf ]. Then, using (20) and (23) we can obtain δK, R [K]∗ ρ P = R [K]δK, ρ D + C1 + C2 ,

(A.1)

where ∂Tw ∂W ∂Qw −∇ · C1 = z φ W ∇p − W dt dx ∂t ∂Sw ∂Sw 0 tf δK δK Tw ∇p + Qw dt dx. C2 = − z ∇· K K 0

tf

(A.2) (A.3)

Note that C1 + C2 = 0. Integrating by parts (A.3) and applying the divergence theorem tf

Tw 1 (A.4) ∇p∇z − z Qw dt dx δK C2 = K K 0 where (8) has been used. Thus, defining tf

Tw 1 ∗ ∇p∇z − z Qw dt R [K] ρ = K K 0

(A.5)

we have C2 = δK, R [K]∗ ρ P , if R [K]δK, ρ D + C1 = 0.

(A.6)


587

Note that z is not determined yet. We want to select z so that (A.6) is satisfied. In order to find a formula for z we apply partial integration and the divergence theorem to (A.2). Using (8) and (21) tf ∂z ∂Tw ∂Qw + z(x, tf )φW (x, tf ) dx + W −φ ∇p∇z − z dt dx. (A.7) C1 = ∂t ∂Sw ∂Sw 0 We impose that z(x, tf ) = 0 in (A.7) so the first integral vanishes and therefore (A.6) reads tf ∂z ∂Tw ∂Qw dt dx = 0. (A.8) W −φ ∇p∇z − z R [K]δK, ρ D + + ∂t ∂Sw ∂Sw 0 On the other hand, using (22) and (24) Np t Np tf f ∂Qw ∂Qw (p) R [K]δK, ρ D = Wj ρj dt = Wρ δ x − xj dt dx, (A.9) ∂Sw ∂Sw 0 j =1 0 j =1 so (A.8) yields     Np tf   ∂Q ∂z ∂Tw w (p) + W −φ ∇p∇z − z − ρδ x − xj  dt dx = 0.  ∂t ∂Sw ∂Sw  0 j =1 Therefore, by choosing z to be the solution of the following adjoint problem:   Np ∂Qw ∂z ∂Tw (p) ∇p∇z − z − ρδ x − xj  =0 + −φ ∂t ∂Sw ∂Sw j =1 z(x, tf ) = 0,

(A.10)

(A.11) (A.12)

where Sw and p are the solutions of (1), (2) with boundary and initial conditions (6)–(8), we have found a z which satisfies our requirements. With this, we have proven the theorem. Appendix B. Numerical methods In this paper we have used two alternative methods to solve equations (1) and (2) describing the direct problem. The first alternative is the IMPES (implicit pressure explicit saturation) method which is based on finite differences and was first proposed by Sheldon et al [44, 46]. It is widely used in commercial reservoir simulators and consists of solving at each time step equation (1) implicitly to derive the pressure distribution, and then computing the saturation distribution at the new time step by solving equation (2) with an explicit finite difference scheme (such as the upwind method) where the values of transmissibilities, T, and pressure, p, are computed at the previous time step. The second alternative is the streamline method (see, for example, [5]), which instead of solving equation (2) in the 2D or 3D space, solves it along one-dimensional streamlines. The method became an attractive alternative for reservoir simulation after the work of Pollock [38], which proposed a very efficient technique to compute the streamlines. The main advantage of the method lies in its speed, since instead of computing the evolution of saturation in a twodimensional or three-dimensional grid, it is only necessary to solve a set of one-dimensional problems to derive the evolution of saturation along streamlines. Thus, it is a method very well suited for inverse problems, where many direct problems have to be solved for convergence and, therefore, very efficient methods are needed in order to solve realistic problems in three dimensions.

588


To derive the method it is useful to write equation (2) in terms of the total flow velocity t) = −T (x, t)∇p(x, t). Introducing this variable in (2) and using (1) leads to ut (x,

∂Sw Tw Tw +u t ∇ = Qw + Q. (B.1) φ ∂t T T This equation can be written along the streamlines using as the independent variable the time of flight, s φ τ= dσ (B.2) 0 |ut | where s represents distance along the streamline. Since the streamline has the direction of the velocity ut , the derivative along a streamline can be obtained from the gradient by ut ∂ = · ∇. ∂s |ut |

(B.3)

Also, from (B.2), ∂ φ ∂ = . ∂s |ut | ∂τ

(B.4)

Thus, ∂ = ut · ∇ (B.5) ∂τ and equation (B.1) can be written as a one-dimensional conservation law in terms of the variable time, t, and time of flight, τ , along the streamlines,

1 Tw ∂Sw ∂ TTw + = Qw + Q . (B.6) ∂t ∂τ φ T φ

The right-hand side is zero except at injection and production wells. Therefore, the procedure to solve the direct problem is the following: (i) Compute the pressure at time tp by solving equation (1) with finite differences in a 2D or 3D grid. (ii) From the pressure distribution compute the velocity distribution ut . (iii) Use Pollock’s method to derive a set of streamlines (at least one streamline must go through each grid cell). (iv) Project the values of the saturation on each grid cell at time tp , in order to derive the saturation distribution along each streamline. (v) Along each streamline solve the conservation law (B.6) with time step ts in order to advance the saturation distribution along each streamline from time tp to time tp + tp . (vi) Project the values of saturation along the streamlines at time tp + tp in order to compute the values of saturation at each grid cell. (vii) Use the saturation values at each cell to compute the values of T at time tp + tp . (viii) If tp = tp + tp is smaller than the total simulation time go to (i). The method takes advantage of the fact that the pressure varies more slowly than the saturation and, therefore, the time step for computation of the pressure, tp , can be taken much larger than the time step for saturation, ts . When other effects, such as gravity or capillarity forces, which cannot be written in terms of τ are included, the mapping between streamlines and grid cells has to be carried out at each time step, and the method loses its efficiency.


589

References [1] Acar R and Vogel C R 1994 Analysis of total variation penalty methods Inverse Problems 10 1217–29 [2] Agarwal B and Blunt M 2003 A streamline-based method for assisted history matching applied to an Arabian Gulf field SPE ATCE (Denver, CO, 5–8 Oct.) SPE paper 84462 [3] Arridge S R 1999 Optical tomography in medical imaging Inverse Problems 15 R41–93 [4] Ascher U and Haber E 2004 Computational methods for large distributed parameter estimation problems with possible discontinuities Inverse Problems, Design and Optimization Symposium 2004 (Rio de Janeiro, Brazil) [5] Batycky R P, Blunt M J and Thiele M R 1997 A 3-D field scale streamline-based reservoir simulator SPERE 12 246–54 [6] Bulishev A E, Souvorov A E, Semenov S Y, Svenson R H, Nazarov A G, Sizov Y E and Tatsi G P 2000 Threedimensional microwave tomography. Theory and computer experiments in scalar approximation Inverse Problems 16 863–75 [7] Caers J 2002 SPE Annual Technical Conference and Exhibition paper 77429 [8] Carrera J 1987 State of the art of the inverse problem applied to the flow and solute transport equations Groundwater Flow and Quality Modeling (NATO ASI Series) [9] Carter R D, Kemp L F Jr, Pierce A C and Williams D L 1974 Performance matching with constraints Soc. Petrol. Eng. J. 14 187–96 [10] Chan T F and Tai X 2003 Identification of discontinuous coefficients in elliptic problems using total variation regularization SIAM J. Sci. Comput. 25 881–904 [11] Chavent G, Dupuy M and Lemmonier C 1975 History matching by use of optimal control theory Soc. Petrol. Eng. J. 15 74–86 [12] Crichlow H B 1976 Modern Reservoir Engineering—A Simulation Approach (Englewood Cliffs, NJ: Prentice Hall) [13] Dennis J E and Schnabel R B 1996 Numerical Methods for Unconstrained Optimization and Nonlinear Equations (SIAM Classics in Applied Mathematics series no 16) (reprint) [14] Dierkes T, Dorn O, Natterer F, Palamodov V and Sielschott 2002 Frechet derivatives for some bilinear inverse problems SIAM J. Appl. Math. 62 2029–113 [15] Dorn O 1997 A transport–backtransport method for optical tomography Inverse Problems 14 1107–30 [16] Dorn O, Bertete-Aguirre H, Berryman J G and Papanicolaou G C 1999 A nonlinear inversion method for 3D electromagnetic imaging using adjoint fields Inverse Problems 15 1523–58 [17] Dorn O, Miller E L and Rappaport C 2000 A shape reconstruction method for electromagnetic tomography using adoint fields and level sets Inverse Problems 16 1119–56 [18] Engl H W, Hanke M and Neubauer A 1996 Regularization of Inverse Problems (Mathematics and Its Applications series no 375) (Dordrecht: Kluwer) [19] Ewing R, Lin T and Falk R 1987 Inverse and ill-posed problems in reservoir simulation Inverse and IllPosed Problem (Notes and Reports on Mathematics in Science and Engineering) (New York: Academic) pp 483–97 [20] Ewing R E, Pilant M S, Wade J G and Watson A T 1995 Identification and control problems in petroleum and groundwater modeling Control Problems in Industry (Progress in Systems and Control Theory) ed I Lasciecka and B Morton (Basle: Birkhauser) pp 119–49 [21] Gomez-Hernandez J J, Sahuquillo A and Capilla J E 1997 Stochastic simulation of transmissivity fields conditional to both transmissivity and piezometric data—I: theory; II: demonstration on a synthetic aquifer J. Hydrol. 203 162–88 [22] Groetsch C W 1993 Inverse Problems in the Mathematical Sciences (Braunschweig: Vieweg) [23] Haber E, Ascher U and Oldenburg D 2000 On optimization techniques for solving nonlinear inverse problems Inverse Problems 16 1263–80 [24] Harding T J, Radcliffe N J and King P R 1996 Optimization of production strategies using stochastic search methods SPE European 3-D Reservoir Modelling Conference paper 35518 [25] He Z, Datta-Gupta A and Yoon S 2002 Streamline-based production data integration with gravity and changing fields conditions Soc. Petrol. Eng. J. 423–36 [26] Jacquard P and Jain C 1965 Permeability distribution from field pressure data Soc. Petrol. Eng. J. 5 281–94 [27] Jahns H O 1966 A rapid method for obtaining a two-dimensional reservoir description from well pressure response data Soc. Petrol. Eng. J. 237 315–27 [28] Kulkarni K N and Datta-Gupta A 2000 Estimating relative permeability from production data: a streamline approach Soc. Petrol. Eng. J. 5 402–11 [29] Landa J L, Kamal M M and Jenkins C D 1996 SPE Annual Technical Conference and Exhibition paper 36511 [30] Landa J L and Home R N 1997 SPE Annual Conference paper 38653

590


[31] Litman A, Lesselier D and Santosa F 1998 Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set Inverse Problems 14 685–706 [32] Mishra S and Parker S C 1989 Parameter estimation for couple unsaturated flow and transport Water Resour. Res. 25 385–96 [33] Natterer F and Wübbeling F 1995 A propagation–backpropagation method for ultrasound tomography Inverse Problems 11 1225–32 [34] Natterer F 2001 The Mathematics of Computerized Tomography (SIAM Classics in Applied Mathematics series no 32) (reprint) [35] Natterer F and Wübbeling F 2001 Mathematical Methods in Image Reconstruction (SIAM Monographs on Mathematical Modeling and Computation) (Philadelphia, PA: SIAM) [36] Nocedal J and Wright S J 1999 Numerical Optimization (New York: Springer) [37] Peaceman D W 1978 Interpretation of well-block pressure in numerical reservoir simulation Soc. Petrol. Eng. J. 3 183–94 Peaceman D W 1978 Trans. AIME 253 [38] Pollock D W 1988 Semianalytical computation of path lines for finite difference models Ground Water 26 743–50 [39] Rahon D, Edoa P F and Masmoudi M 1997 Inversion of geological shapes in reservoir engineering using well tests and history matching Proc. SPE 38656 (San Antonio, Oct 5–6) [40] Romero C and Carter J N 2002 Develop. Petrol. Sci. 51 323–63 [41] Rondi L and Santosa F 2001 Enhanced electrical impedance tomography via the Mumford–Shah functional ESAIM: Control Optim. Calculus Variations 6 517–38 [42] Sagar R K, Kelkar M G and Thomson L G 1995 Reservoir Description by Integrating Well-Test Data and Spatial Statistics SPEFE (Dec.) p 267 [43] Santosa F 1996 A level-set approach for inverse problems involving obstacles ESAIM: Control Optim. Calculus Variations 1 17–22 [44] Sheldon J W, Harris C D and Bavly D 1960 A method for general reservoir behavior simulation on digital computers Presented at SPE 35th Annual Meeting (Denver, Oct. 2–5) paper SPE 1521-G [45] Sun N and Yeh W G 1990 Coupled inverse problems in groundwater modelling 1. Sensitivity analysis and parameter identification Water Resour. Res. 26 2507–25 [46] Thomas G W 1982 Principles of Hydrocarbon Reservoir Simulation (Englewood Cliffs, NJ: Prentice-Hall) [47] Vasco D M, Yoon S and Datta-Gupta A 1999 Integrating dynamic data into high-resolution reservoir models using streamline-based analytic sensitivity coefficients SPE J. 4 389–99 [48] Vogel C R 2002 Computational Methods for Inverse Problems (SIAM Series ’Frontiers in Applied Mathematics’) [49] Vogeler M 2003 Reconstruction of the three-dimensional refractive index in electromagnetic scattering by using a propagation–backpropagation method Inverse Problems 19 739–53 [50] Wagner B J and Gorelick S M 1987 Optimal groundwater quality management under parameter uncertainty Water Resour. Res. 23 1162–74 [51] Wang Y and Kovscek A R 2000 Streamline approach for history matching production data Soc. Petrol. Eng. J. 5 353–62 [52] Wasserman M L, Emanuel A S and Seinfeld J H 1975 Practical applications of optimal-control theory to history-matching multiphase simulator models Soc. Petrol. Eng. J. 15 347–55 [53] Wu Z, Reynolds A C and Oliver D S 1999 Conditioning geostatical models to two-phase production data Soc. Petrol. Eng. J. 4 142–55

History Matching

Recommend Documents