Analysis of Geometric Selection of the Data-Error Covariance Inflation for ES-MDA

12/03/2018
by   Alexandre A. Emerick, et al.
Petrobras
0

The ensemble smoother with multiple data assimilation (ES-MDA) is becoming a popular assisted history matching method. In its standard form, the method requires the specification of the number of iterations in advance. If the selected number of iterations is not enough, the entire data assimilation must be restarted. Moreover, ES-MDA also requires the selection of data-error covariance inflations. The typical choice is to select constant values. However, previous works indicate that starting with large inflation and gradually decreasing during the data assimilation steps may improve the quality of the final models. This paper presents an analysis of the use of geometrically decreasing sequences of the data-error covariance inflations. In particular, the paper investigates a recently introduced procedure based on the singular values of a sensitivity matrix computed from the prior ensemble. The paper also introduces a novel procedure to select the inflation factors. The performance of the data assimilation schemes is evaluated in three reservoir history-matching problems with increasing level of complexity. The first problem is a small synthetic case which illustrates that the standard ES-MDA scheme with constant inflation may result in overcorrection of the permeability field and that a geometric sequence can alleviate this problem. The second problem is a recently published benchmark and the third one is a field case with real production data. The data assimilation schemes are compared in terms of a data-mismatch and a model-change norm. The first norm evaluates the ability of the models to reproduce the observed data. The second norm evaluates the amount of changes in the prior model. The results indicate that geometric inflations can generate solutions with good balance between the two norms.

READ FULL TEXT VIEW PDF

Authors

page 16

page 17

page 21

page 22

page 26

page 27

06/02/2022

Hybrid iterative ensemble smoother for history matching of hierarchical models

The choice of the prior model can have a large impact on the ability to ...
05/16/2019

4D Seismic History Matching Incorporating Unsupervised Learning

The work discussed and presented in this paper focuses on the history ma...
10/10/2019

On k-error linear complexity of binary sequences derived from Euler quotients modulo 2p

We consider the k-error linear complexity of binary sequences derived fr...
03/31/2018

Improving Portfolios Global Performance with Robust Covariance Matrix Estimation: Application to the Maximum Variety Portfolio

This paper presents how the most recent improvements made on covariance ...
05/21/2014

Sparse Precision Matrix Selection for Fitting Gaussian Random Field Models to Large Data Sets

Iterative methods for fitting a Gaussian Random Field (GRF) model to spa...
08/06/2021

Adjusting PageRank parameters and Comparing results

The effect of adjusting damping factor α and tolerance τ on iterations n...
07/29/2009

Cooperative Training for Attribute-Distributed Data: Trade-off Between Data Transmission and Performance

This paper introduces a modeling framework for distributed regression wi...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

Abstract

The ensemble smoother with multiple data assimilation (ES-MDA) is becoming a popular assisted history matching method. In its standard form, the method requires the specification of the number of iterations in advance. If the selected number of iterations is not enough, the entire data assimilation must be restarted. Moreover, ES-MDA also requires the selection of data-error covariance inflations. The typical choice is to select constant values. However, previous works indicate that starting with large inflation and gradually decreasing during the data assimilation steps may improve the quality of the final models.

This paper presents an analysis of the use of geometrically decreasing sequences of the data-error covariance inflations. In particular, the paper investigates a recently introduced procedure based on the singular values of a sensitivity matrix computed from the prior ensemble. The paper also introduces a novel procedure to select the inflation factors. The performance of the data assimilation schemes is evaluated in three reservoir history-matching problems with increasing level of complexity. The first problem is a small synthetic case which illustrates that the standard ES-MDA scheme with constant inflation may result in overcorrection of the permeability field and that a geometric sequence can alleviate this problem. The second problem is a recently published benchmark and the third one is a field case with real production data. The data assimilation schemes are compared in terms of a data-mismatch and a model-change norm. The first norm evaluates the ability of the models to reproduce the observed data. The second norm evaluates the amount of changes in the prior model. The results indicate that geometric inflations can generate solutions with good balance between the two norms.

1 Introduction

Iterative forms of the ensemble smoother (ES) are rapidly becoming popular assisted history matching (AHM) techniques with a growing number of field applications reported in the literature in the last few years. An evidence of the increasing interest of the industry is the recent appearance of commercial AHM softwares based on iterative smoothers. Among these methods, the ES-MDA (Emerick and Reynolds, 2013a) is a popular choice because of the simplicity of its formulation. In a standard implementation of ES-MDA, it is necessary to specify in advance the number of data assimilations, , and the set of data-error covariance inflation factors,

. In fact, classifying ES-MDA as an iterative ES is debatable because its termination criterion does not rely on a convergence check

(Mannseth, 2015). A very common choice for ES-MDA is to set for (Emerick and Reynolds, 2013a, b; Emerick, 2016a). Unfortunately, if the selected value for is not sufficient to obtain a reasonable data match, it is necessary to restart the entire data assimilation with a larger value. In practice, however, it is desirable to keep as small as possible to reduce computational cost. Moreover, it has been suggested in the literature (Emerick and Reynolds, 2013a; Le et al., 2016; Emerick, 2016a) that starting with large value for and gradually decreasing this value during the data assimilation steps improves the quality of the final models. The rationality for this procedure is based on the conjecture that at the beginning of the data assimilation the predicted data are far from the observations, so it is beneficial to restrict the changes in the models to avoid overcorrection. This conjecture is supported by previous experiences with the application of Gauss-Newton and Levenberg-Marquardt methods for history matching; see, e.g., (Oliver et al., 2008, Chap. 8).

The problem of selecting the inflation factors for ES-MDA has been investigated in previous works. Emerick (2016a) presented a ES-MDA algorithm where and are selected according to the progress of the average data-mismatch objective function. This algorithm was tested in a field case requiring data assimilations. However, the final results were not significantly superior to the standard method with . Inspired by the work of Iglesias and Dawson (2013), Le et al. (2016) proposed an ES-MDA scheme that adaptively select and using the regularization condition proposed in (Hanke, 1997, 2010). They also investigated an alternative procedure based on limiting the maximum change in the model over a data assimilation step. They tested both schemes in a modified version of the PUNQ-S3 case and concluded that the standard ES-MDA with 8 and 16 data assimilations and constant inflation resulted in over and undershooting of the permeability field (very high and low values of permeability). On the other hand, the proposed schemes obtained smoother permeability realizations. Unfortunately, these adaptive schemes tend to be computationally too demanding when applied to large-scale problems because they often require several repetitions of the analysis (inversion) step, which becomes problematic for big models with large number of data points. Moreover, these schemes usually require several data assimilation steps, often more than ten. Recently, Rafiee and Reynolds (2017) proposed to select the inflation factors in a geometrically decreasing order where the first inflation is based on the regularization condition of (Hanke, 1997). For a simple test problem, they showed that the geometric scheme resulted in smoother permeability fields compared to constant inflations.

Here, the objective is to investigate the use of geometric inflations. This paper also introduces a novel scheme to select the inflations. The goal is to develop a robust procedure to improve the data assimilation results without compromising the computational efficiency and ease of implementation of the method.

The remaining of the paper is organized as follows: the next two sections reviews basic concepts of regularization for linear and nonlinear inverse problems. The section after that presents the ES-MDA equations emphasizing that the inflation factors plays the same role of the regularization parameter described in the inverse problems literature. This section also reviews the method proposed in (Rafiee and Reynolds, 2017) and presents the new scheme for selecting a geometric sequence of . Both methods are tested in three reservoir history-matching problems. This first problem is similar to the case used in (Rafiee and Reynolds, 2017) and it was designed to show overcorrection in the permeability field using constant inflation factors. The second one is a more realistic history-matching exercise proposed in (Maschio et al., 2013; Avansi and Schiozer, 2015). The third problem is the same field case used in (Emerick, 2016a). The last section of the paper summarizes the main conclusions of the work.

2 Regularization for Linear Inverse Problems

Consider the problem of finding such that

(1)

where

is a vector with given observations and

is a matrix. A unique solution exists only if and , where and denote the range and the null space of , respectively. From a practical perspective, however, it is necessary to find only an approximate solution of (1) because is typically corrupted with noise. An approximate solution of (1) can be obtained by solving the least-squares problem

(2)

All norms in this paper refer to the -norm, i.e., . The linear least-squares problem (2) is ill-posed in the sense that small errors in may cause errors of arbitrary size in . Instead, one can seek for the solution of the regularized least-squares problem

(3)

where is the regularization parameter. Eq. 3 is known as Tikhonov regularization. The solution of (3) exists; see, e.g., (Kaipio and Somersalo, 2005), and can be written as

(4)

The choice of the regularization parameter controls the quality of the estimates. There are several methods for choosing

presented in the literature; see, e.g., (Hamarik et al., 2012) for a detailed discussion. Perhaps the most well-known methods are the ones based on the Morosov’s Discrepancy Principle (MDP) (Morozov, 1984), which states that it is not reasonable to expect the approximate solution to yield an error smaller than the noise level in the data vector. The MDP can be formulated as selecting such that

(5)

In the above equation, is the noise level in and is an extra parameter to avoid underregularization (Kaipio and Somersalo, 2005). Let be corrupted with an unknown additive noise vector , i.e.,

(6)

where is the noiseless observation vector. Assuming that is a random vector, then it is convenient to define the square of the noise level as

(7)

Let the noise vector be a sample from a multivariate Gaussian distribution. Without loss of generality, assume that

with

denoting the standard deviation. Then

follows as chi-squared distribution with

degrees of freedom. Therefore, the noise level becomes

(8)

3 Regularization for Nonlinear Inverse Problems

Consider the problem of computing an approximate solution for the nonlinear least-squares problem

(9)

with being a nonlinear vector function. In order to solve (9), it is necessary to derive an iterative scheme. Consider the first-order Taylor expansion around the current guess

(10)

where

(11)

is the sensitivity matrix. The linearized version of (9) becomes

(12)

Let and , then (12) becomes

(13)

which has the same form of the linear least-squares problem (2). This problem is ill-posed and instead of (13), consider the following regularized problem

(14)

which has the same solution presented in Eq. 4. The resulting iterative procedure corresponds to the Levenberg-Marquardt (LM) algorithm for nonlinear least squares; see, e.g. (Nocedal and Wright, 2006).

Hanke (1997) used the MDP to propose a regularized LM scheme (see also (Hanke, 2010)). In his scheme, the regularization parameter should be selected such that

(15)

for some . This condition is the base for the works presented by (Iglesias and Dawson, 2013; Luo et al., 2015; Le et al., 2016) and it is the starting point for the procedure proposed in (Rafiee and Reynolds, 2017), which is discussed in the next sections.

4 Es-Mda

The ES-MDA analysis equation to update a vector of model parameter, , is

(16)

where is the size of the ensemble. is the vector containing the observation and is the vector of data predicted by the th model realization at the th data assimilation step of ES-MDA. is a random vector drawn from , where is the data-error covariance matrix. Here, is used to denote the number of data points. The covariance matrices and are computed as

(17)

and

(18)

where

(19)

and

(20)

In ES-MDA, is referred to as data-error covariance inflation factor. In the standard implementation of this method, Eq. 16 is applied a pre-defined number of times, , and the set must satisfy the condition

(21)

This condition was originally derived by Emerick and Reynolds (2013a) using linear algebra to ensure that ES-MDA is consistent with the standard ES for the linear-Gaussian case. However, the same condition can be derived directly from Bayes’ rule as discussed in (Stordal and Elsheikh, 2015; Emerick, 2016a).

4.1 Geometric Selection of ES-MDA Inflation

Here, the objective is to derive a scheme to select for to be used in the ES-MDA method. Specifically, the goal is to select a geometrically decreasing sequence such that

(22)

for . For a fixed , which can be selected with one of the procedures discussed in the next sections, the computation of is straightforward. Recall that

(23)

Hence, defining the function

(24)

it is easy to find such that using, for example, the bisection method (Press et al., 2007).

Before introducing the procedures to compute , it is convenient first to consider the ES-MDA equation to update the ensemble mean, , i.e.,

(25)

Eq. 25 uses the fact that . Using (17) and (18) in (25) and setting

(26)

or

(27)

where denotes a prior realization and is the pseudo-inverse of . Note that there will be no need to compute as the objective is only to derive a scheme to compute . Now call

(28)
(29)

and

(30)

so that Eq. 27 can be written as

(31)

Here, the dimension of the vectors and are and , respectively. Eq. 31 assumes the same form of (4), which is the solution of the regularized linear least-squares problem (3). Clearly the inflation factor plays the role of the Tikhonov regularization parameter. It is worth mentioning that the matrix defined in Eq. 30 is scaled by the square-root of the data-error covariance matrix. This has numerical benefits because the matrix may be constructed with data of different orders of magnitude. Hence, scaling avoids losing relevant information from data when truncating small singular values. Moreover, corresponds to the ensemble equivalent of the dimensionless sensitivity matrix (Zhang et al., 2002) as discussed in the Appendix A of (Emerick and Reynolds, 2012). Tavakoli and Reynolds (2010) showed that for a linear-Gaussian problem the singular values of the dimensionless sensitivity govern the reduction in uncertainty and that the smallest singular values have negligible influence on the reduction of uncertainty. In this sense, defining as in Eq. 30 is optimal.

It is important to note that the data-error covariance matrix, , plays a crucial role in the data assimilation process. The construction of this matrix, however, is not very well discussed in the petroleum literature. Here, it is assumed that the measurement errors in data are independent, such that is a diagonal matrix. This is a very common assumption and it is present in the large majority of the history-matching publications. Perhaps the main exceptions are the papers that discuss assimilation of seismic data; see, e.g., (Emerick, 2016a; Luo and Bhakta, 2017), and the work from Evensen and Eikrem (2018)

which discusses the use of rate an cumulative production data. However, even the selection of the variance (diagonal elements of

) is not very well documented. A common procedure for production data is to estimate the data variance using a smoothed version of the time series (Zhao et al., 2006; Emerick and Reynolds, 2011). However, the noise in the measured data are not the only source of errors. There are also the “model errors,” which includes errors because of coarse discretizations, simplified physics and insufficient parameterizations, just to mention a few. Model errors are difficult to be correctly accounted for and traditionally they are neglected in the majority of the history-matching applications. However, if we make the crude assumption that model errors are additive and follow the Gaussian distribution , then it is straightforward to include their effect in the data assimilation by simply replacing by . In fact, some recent works (Sun et al., 2017; Oliver and Alfonzo, 2017) show that selecting a larger variance for the data noise tends to partially compensate for model errors and improve the data assimilation results, especially in terms of uncertainty quantification. For this reason and based on previous experiences with field history-matching problems, in the last two test problems presented in this paper the selected data-error variances are overestimated compared to the variance that would be estimated based solely on the production time series.

The next sections present two procedures to compute . The first procedure was originally presented by Rafiee and Reynolds (2017) and the second is proposed in this paper.

4.1.1 Procedure Proposed by Rafiee and Reynolds

Rafiee and Reynolds (2017) used the Hanke’s regularization condition (15) as a starting point to define a criterion to select . First, they wrote Eq. 15 as

(32)

where . Then they noted that

(33)

with denoting the

th eigenvalue of

and the th singular value of . Instead of selecting that satisfies the condition (33), they use

(34)

The equality occurs for

(35)

where

(36)

with denoting the number of non-zero singular values of .

Rafiee and Reynolds (2017) proposed to use , which leads to . The resulting procedure is summarized in the Algorithm 1. Here, this method is referred to as GEO1.

  1. Specify .

  2. Compute

    where

  3. Compute by solving

    using the bisection method.

  4. Apply ES-MDA with .

Algorithm 1 ES-MDA-GEO1

4.1.2 Procedure Proposed in this Work

The procedure to define the inflation factors proposed here is inspired by the ideas in (Rafiee and Reynolds, 2017) in the sense that it is also based on a geometrically decreasing sequence. However, there are two main differences between the two procedures.

The first difference is that instead of computing independently from the value of , the new proposal computes based on a specified inflation factor for the last step of the data assimilation, . Then the value of is checked against the MDP. The idea is to ensure that the geometric scheme will not result in a sequence that decays too fast. In fact, during the computational experiments that conducted to develop this procedure, it has been noticed that the procedure proposed in (Rafiee and Reynolds, 2017) tends to generate good results as long as an appropriate value for is selected in advance. This is particularly evident for a small number of data assimilations, say , and a large initial inflation, say . This situation generates a small value for and consequently the last inflation is often very close to one. This fact is illustrated in the test cases presented in the next sections. The idea to circumvent this problem is to specify the value of the inflation at the last data assimilation and compute using

(37)

Now, define

(38)

and compute by solving using the bisection method.

The second difference between Rafiee and Reynolds (2017) and the new proposal is that Rafiee and Reynolds (2017) used the Hanke’s regularization condition (15) to derive a procedure to compute . In the new proposal the value of computed using Eq. 37 is tested against the MDP. If does not satisfy the MDP, then the number of data assimilation steps is increased until a suitable value of is obtained.

In order to check the MDP, it suffices to define a discrepancy function as

(39)

and compute such that by solving a root-finding problem.

Consider the singular value decomposition of

(40)

where and are orthogonal matrices and is a matrix with diagonal element given by

(41)

Using (40) and the properties and in (39) results in

(42)

Recall that is given by

(43)

where can be chosen arbitrarily as long as . In particular, let in the following. Using (43) in (42) results in

(44)

The value of such that can be computed using the Newton-Raphson method (Press et al., 2007), in which case it is necessary to compute the derivative

(45)

Note that for any . Therefore, is a strictly increasing function of . Moreover, note that

(46)

and

(47)

Therefore, has a unique solution as long as the noise level satisfies the following condition

(48)

The matrix forms a basis for . The left hand side of the inequality (48) can be interpreted as that one would expect that the components of orthogonal to are due to noise (Kaipio and Somersalo, 2005). However, in the specific case of ensemble data assimilation, is given by Eq. 30; i.e., the components of are computed from an ensemble of size . Therefore, . If , which is often the case for history matching with smoothers, then . Hence, it is possible to have the cases where . This effectively means that the MDP would not be violated for any . To avoid this case, a safe simplification was added to the proposed method by simply dropping the term in Eq. 44. This simplification is based on the fact that in practice we always use localization which expands the space for reconstruction of (the localized version of is full rank). The right hand side of the inequality (48) can be interpreted as that one would expect the noise level not to exceed the signal (Kaipio and Somersalo, 2005).

Note that the MDP was proposed as a procedure to regularize linear inverse problems. The condition proposed by Hanke (1997), on the other hand, was developed to generate a sequence of that stabilizes the LM algorithm in nonlinear least-squares problems. However, as discussed in the introduction section of this paper, the use of the sequence of from Eq. 15 with ES-MDA generates a method that takes too many iterations for practical applications. This fact motivated Rafiee and Reynolds (2017) to use Eq. 15 to define only . However, one might argue that there is no reason to enforce Eq. 15 to compute if a different sequence is to be used. In the new proposal, on the other hand, the idea is to use the MDP to check if the selected value of is enough to ensure a stable linear inversion in the first data assimilation step. But, similarly to Rafiee and Reynolds (2017), a geometric decreasing sequence is imposed to keep the computational cost affordable for large-scale problems.

Algorithm 2 summarizes the proposed method, which is referred to as GEO2. Besides the algorithm requires the specification of a maximum inflation, , which is used as a safeguard to ensure that the method will not take too many data assimilations. The fourth step of Algorithm 2 requires to solve the discrepancy function which is done using the Algorithm 3. The cases presented in this paper use the following values for the parameters of Algorithms 2 and 3: , , , and .

  1. Specify , and .

  2. Compute by solving

    using the bisection method.

  3. Compute

  4. Compute by solving using the Newton-Raphson method (Algorithm 3).

  5. If set and return to step 2.
    Else apply ES-MDA with .

Algorithm 2 ES-MDA-GEO2
  1. Specify , , and .

  2. Compute

  3. If set and exit.

  4. Compute

  5. For to do

    1. Compute

      and

    2. Compute

    3. If set and exit.
      Else if set and exit.

    End for.

  6. Set .

Algorithm 3 Newton-Rapshon to compute

5 Test Cases

5.1 2D Model

The first test case is a small synthetic history-matching problem. The reservoir corresponds to a 2D model discretized into a uniform grid. The only uncertain property is the permeability field. The reference (true) permeability field (Fig. 1) was generated using a spherical covariance function with isotropic correlation length of 40 gridblocks. The value of the prior mean of the natural logarithm of permeability (log-permeability) is constant and equal to 5.5 and the prior standard deviation is equal to one. No hard data were used at well locations to make the problem more challenging. The same setup was used to generate 200 prior realizations for history matching. The reservoir produces with four five-spots (Fig. 1). All wells operate at constant bottom-hole pressure of 20000 kPa for producers and 30000 kPa for water injectors. The observed data correspond to values of oil and water rates every 90 days for a period of five years. Random noise with standard deviation corresponding to 3% of the actual data were added to form the synthetic measurements. A small value for the measurement errors was intentionally selected because it tends to make the problem more prone to overcorrection. Also note that in this simple problem there are no model errors, so a diagonal matrix is used with variance equal to the square of the standard deviation of the random noise used to generated the synthetic measurements.

Figure 1: True log-permeability (2D model). Circles correspond to oil producing and triangles to water injection wells.

The Algorithm 3 was applied to this test problem resulting in , which is illustrated in the plot of the discrepancy function presented in Fig. 2. Using in the Algorithm 2 resulted in and . For comparisons, was also used in the cases with constant inflation and the method proposed in (Rafiee and Reynolds, 2017). During all data assimilations a Schur product-based localization (Houtekamer and Mitchell, 2001) was applied to the Kalman gain using a constant correlation length corresponding to the size of 50 gridblocks in the Gaspari-Cohn correlation function. The matrix inversion in Eq. 16 was done using the subspace inversion procedure (Evensen, 2004) keeping 99% of the sum of the singular values.

Figure 2: Discrepancy function (2D model).

Table 1 summarizes the inflation factors for the three cases. The case with constant inflation is labeled as -CONST, the method from (Rafiee and Reynolds, 2017) is labeled as -GEO1 while the method proposed in this paper is labeled as -GEO2. The method -GEO1 resulted in a large initial inflation () and because the number of data assimilations is fixed, the procedure resulted in an inflation of only 1.11 in the last data assimilation step. Fig. 3 shows the mean log-permeability fields obtained by each method. In all cases, the data assimilation recovered the main permeability features of the true model. However, the permeability resulting from the case -CONST shows some regions with overcorrection (very low permeability values). Table 2 presents the average values of the root mean squared error (RMSE) and the squared data-mismatch norm divided by the number of data points obtained by each method. The RMSE of the th posterior realization, , is computed as

(49)

where is the true model. The lowest RMSE was obtained by -GEO2 indicating a better recovering the true values of log-permeability. The data-mismatch norm is computed as

(50)

The case -CONST resulted in the lowest average data-mismatch norm indicating that this case obtained the best agreement with the observations. This fact is also illustrated in Fig. 4 which shows the predicted water rate for one well of the field. In addition to the historical period, the plots in Fig. 4 include a five-years forecast period. The results in this figure show that although the case -CONST obtained the best data match, it resulted in a narrow span of forecasted water rate, which do not cover the forecast from the true model. The case -GEO1 clearly do not match the observations sufficiently well and the case

-GEO2 resulted in a reasonable compromise between data-match and forecast. Ideally, a successful history-matching method should be able to provide a good sampling of the posterior probability density function of the model parameters. However, it is difficult to compare methods in terms of their sampling performance because generating a reference sampling is very computationally demanding. Even though it is not correct to claim that

-GEO2 obtained a better sampling than -CONST and -GEO1, the results show evidences that it resulted in better model estimates and better uncertainty assessment of the production forecast.

-CONST -GEO1 -GEO2
1 7 1442941.18 1087.48
2 7 138031.75 362.83
3 7 13204.12 121.05
4 7 1263.11 40.39
5 7 120.83 13.48
6 7 11.56 4.50
7 7 1.11 1.50
1 0.0957 0.3336
Table 1: Inflation factors (2D model)
(a) -CONST
(b) -GEO1
(c) -GEO2
-CONST
Figure 3: Mean log-permeability (2D model).
-CONST -GEO1 -GEO2
RMSE 1.714 2.227 1.680
Squared data-mismatch norm 3.390 63.359 14.828
Table 2: Average metrics (2D model)
(a) -CONST
(b) -GEO1
(c) -GEO2
Figure 4: Water rate at well P2. Red dots are the observations, red curve is the prediction from the true model and blue curves are the predictions from the posterior ensemble (2D model).

5.2 Unisim-I-H

The second test case is the benchmark problem UNISIM-I-H (Avansi and Schiozer, 2015), which is available in (Maschio et al., 2013). The benchmark was constructed with actual data from Namorado Field (Campos Basis, Brazil). The UNISIM-I-H model has gridblocks, with dimensions  meters and 37,000 active gridblocks. A total of 500 realizations of petrophysical properties (porosity, permeability in the three orthogonal directions and net-to-gross ratio) is provided with the dataset. There are 14 oil producing and 11 water injection wells. Fig. 5 shows the position of the wells projected in the first layer of the model. Note that the actual wells are perforated in different layers of the model, typically water injectors are perforated in bottom while oil producing wells are perforated in the top of the reservoir. The synthetic observed data provided with the dataset correspond to monthly “measurements” of water cut and bottom-hole pressure (BHP) for the oil producing wells and BHP for the water injection wells for a period of 4018 days. These observations were generated based on a fine-scale model with 3.5 million active gridblocks and are corrupted with noise. Here, the noise in each individual datum was assumed to be an independent sample from a Gaussian distribution with zero mean and standard deviation corresponding to 10% of the water cut data and 10  for pressure measurements. Besides the petrophysical properties, the UNISIM-I-H case also includes six global parameters, whose prior uncertainties were modeled as independent triangle distributions with values presented in Table 3. A detailed description of the UNISIM-I-H case can be found in (Avansi and Schiozer, 2015; Maschio et al., 2013).

Figure 5: Position of the wells (UNISIM-I-H)
Parameter Mode Min. Max.
Rock compressibility in
Depth of the oil-water contact at the East block (Fig. 5) 3174 3169 3179
Maximum water relative permeability 0.35 0.15 0.55
Corey exponent of water relative permeability 2.4 2.0 3.0
Multiplier for vertical permeability 1.5 0.0 3.0
Table 3: Prior distribution of global parameters (UNISIM-I-H)

ES-MDA was applied to update the prior ensemble considering 13 different configurations of inflation factors. Table 4 summarizes the cases considered. In this table, -CONST and -CONST stand for ES-MDA with constant inflation factors with four and eight data assimilations, respectively. -GEO(1E2), -GEO(1E3), -GEO(1E4) and -GEO(1E5) stand for ES-MDA with four data assimilations and geometric selection of the inflations, where the values of the first inflation are given between the braces, i.e., , , and , respectively. The same initial values for are also used for the cases with eight data assimilations. The method proposed in (Rafiee and Reynolds, 2017) was tested with four and eight data assimilations (-GEO1 and -GEO1). The method proposed in this paper required eight data assimilations to satisfy the MDP and it is labeled as -GEO2. Fig. 6 shows a plot of the discrepancy function indicating that the inflation value such that is 1172. Using in the Algorithm 2 resulted in , and . In all cases, Kalman gain localization was applied using the Gaspari-Cohn correlation function with correlation length of 2000 meters. Subspace inversion keeping 99% of the sum of the singular values was applied.

-CONST -GEO(1E2) -GEO(1E3) -GEO(1E4) -GEO(1E5) -GEO1
1 4 100 1000 10000 100000 4010.30
2 4 23.54 103.71 471.69 2172.79 258.07
3 4 5.54 10.76 22.25 47.21 16.61
4 4 1.30 1.12 1.05 1.03 1.07
1 0.2354 0.1037 0.0472 0.0217 0.0644
-CONST -GEO(1E2) -GEO(1E3) -GEO(1E4) -GEO(1E5) -GEO1 -GEO2
1 8 100 1000 10000 100000 4010.30 3273.79
2 8 58.64 401.08 2812.60 19929.85 1296.25 1091.58
3 8 34.39 160.87 791.07 3971.99 418.99 363.96
4 8 20.17 64.52 222.50 791.61 135.43 121.36
5 8 11.83 25.88 62.58 157.77 43.77 40.46
6 8 6.94 10.38 17.60 31.44 14.15 13.49
7 8 4.07 4.16 4.95 6.27 4.57 4.50
8 8 2.39 1.67 1.39 1.25 1.48 1.50
1 0.5864 0.4011 0.2813 0.1993 0.3232 0.3334
Table 4: Inflation factors (UNISIM-I-H)
Figure 6: Discrepancy function (UNISIM-I-H).

Even though a rigorous comparison in terms of the sampling performance is not possible because the extreme computational requirements to generate a reference distribution, it is possible to compare history-matching methods in terms of some desirable properties. One important quality for a good history-matching method is the ability to achieve a good data match at a reasonable computational cost. Here, the data match quality is evaluated in terms of the data-mismatch norm . Another important quality of a good history-matching method is the ability to preserve the geological realism. In this sense, it is desirable that the method makes the smallest changes possible in the prior realizations to match data if these realizations were generated from a rigourous geomodeling process. The model changes can be evaluated with the norm , where with denoting the th prior realization and denoting a diagonal matrix containing the inverse of the prior standard deviation of . History matching practitioners usually give more emphasis on the data match quality. However, minimizing the amount of changes in the geological model is also important. Unfortunately, these objectives are often conflicting. In this sense, it is desirable to obtain models with a good balance between the two norms.

Fig. 7 plots the values of the average of the square of the data-mismatch norm, divided by the number of data points, against the average squared-norm of the model change, divided by the number of model parameters. The model-change norm presented in Fig. 7 was computed only for the log-permeability in the -direction of the grid. Similar plots were obtained for the other petrophysical properties, but they are omitted here. The dashed curve in Fig. 7 was drawn connecting the “best” points considering the two norms to generate the analogous of a Pareto frontier. This procedure is closely related to the L-curve method; see, e.g. (Hansen and O’Leary, 1993), which is another method described in the literature to select regularization parameters for linear inverse problems. For example, note that the case -CONST resulted in the smallest data-mismatch norm, but at a cost of a large change in the model. The cases closer to the region of maximum curvature correspond to the ones with the best balance between data match and model change. The results in Fig. 7 indicate that and must be selected in conjunction. For example, for , it seems that the best choice of is on the order of . However, for , the best choice of is on the order to . In fact, the proposal from Rafiee and Reynolds (2017) obtained a particularly well-balanced result for the case with eight data assimilations (-GEO1). However, the case -GEO1 resulted in a larger data-mismatch norm, making this case distant from the maximum-curvature region. The procedure proposed in this paper (-GEO2) obtained very similar results to the case -GEO1 showing a good compromise between data match and model change.

Figure 7: Average data-mismatch and model-change norms (UNISIM-I-H).

Fig. 8 shows the ensemble mean of log-permeability for a intermediate layer for the cases -CONST, -GEO1 and -GEO2. This figure shows that -GEO1 and -GEO2 resulted smoother fields with smaller changes in the prior realization. Fig. 9 shows the normalized variance of log-permeability for the same three cases. The normalized variance is computed dividing the variance of the posterior ensemble by the variance of the prior ensemble. Clearly the methods -GEO1 and -GEO2 preserved more the variability in the posterior ensemble. Even though the correct level of variance is unknown, previous experience with ensemble data assimilation methods show a tendency of variance underestimation even when localization is used; see, e.g., (Emerick, 2016b). Therefore, the larger normalized variances of the cases -GEO1 and -GEO2 are a very positive result, especially if the posterior ensemble is to be used for uncertainty quantification of production forecast and to assimilated future data. Fig. 10 shows the predicted data for a selected well of the field. This figure illustrates the fact that even though -GEO1 and -GEO2 resulted in larger data-mismatch norm than -CONST, the three procedures results in predicted data with a reasonable agreement with the observations.

(a) Prior
(b) -CONST
(c) -GEO1
(d) -GEO2
Prior
Figure 8: Mean log-permeability (UNISIM-I-H, layer 12).
(a) -CONST
(b) -GEO1
(c) -GEO2
-CONST
Figure 9: Normalized variance of log-permeability (UNISIM-I-H, layer 12).
(a) -CONST, Water cut
(b) -CONST, BHP
(c) -GEO1, Water cut
(d) -GEO1, BHP
(e) -GEO2, Water cut
(f) -GEO2, BHP
Figure 10: Predicted production data for the well NA1A. Red dots are the observations, gray curves are the predictions from the prior and blue curves are the predictions from the posterior ensemble (UNISIM-I-H).

5.3 Field Case

The last test case corresponds to the same field problem presented in (Emerick, 2016a). This is heavy-oil reservoir in Campos Basis operating with 36 wells. The main recovery mechanism is waterflooding. The model contains 130,000 active gridblocks and the prior ensemble includes 200 realizations of porosity and horizontal permeability. These realizations are the same used in (Emerick, 2016a). Besides porosity and permeability, the history matching parameters include the transmissibility across several faults, the maximum water relative permeability curve and the Corey exponent of the water-oil relative permeability. The production history includes monthly measurements of water cut, gas-oil ratio (GOR) and BHP. Compared to the original study presented in (Emerick, 2016a), the production history has been updated with five additional years. The noise in the data was assumed Gaussian with zero mean and standard deviation of 10% of the water cut, 20% of the GOR and a constant value of 5 kgf/cm for pressure measurements. During data assimilations, Schur-product localization was applied to the Kalman gain using the Gaspari-Cohn correlation function with a constant correlation length of 2000 meters. Similarly to the other cases, subspace inversion was applied retaining 99% of the singular values.

Unlike the UNISIM-I-H case, for this field problem, the discrepancy function (42) was greater than zero any . For comparisons, the cases with constant inflations and the method proposed in (Rafiee and Reynolds, 2017) were also performed using four data assimilations. Table 5 presents the inflation factors for each case. Note that -GEO1 started with relatively large inflation (), as a result the inflation at the last data assimilation is only 1.04. For -GEO2, on the other hand, because the geometric sequence is selected such that at the last data assimilation , the initial inflation was only .

-CONST -GEO1 -GEO2
1 4 16986.84 37.33
2 4 670.47 12.79
3 4 26.46 4.38
4 4 1.04 1.50
1 0.0395 0.3425
Table 5: Inflation factors (Field Case)

Fig. 11 shows the plot of the average data-mismatch and model-change norms for the three cases. The error bars in this figure correspond to one standard deviation for each side. The results in this figure shows that -CONST resulted in the best data match and -GEO1 resulted in the lowest model changes. The proposed approach obtained a compromise between both norms. Fig. 12 shows the ensemble mean log-permeability indicating that all cases resulted in reasonable model estimates, which preserve the main characteristics of the prior realization. In terms of the posterior variance, the case -GEO1 resulted in the highest values as illustrated in Fig. 13. This results is consistent with the fact that this case resulted in the smallest changes in the prior realizations. However, the data matches after -GEO1 are clearly of inferior quality as indicated in Fig. 14, which shows the predicted data for three different wells of the field.

Figure 11: Average data-mismatch and model-change norms (Field Case). The error bars correspond to one standard deviation.
(a) Prior
(b) -CONST
(c) -GEO1
(d) -GEO2
Prior
Figure 12: Mean log-permeability (Field Case, layer 17).
(a) -CONST
(b) -GEO1
(c) -GEO2
-CONST
Figure 13: Normalized variance of log-permeability for (Field Case, layer 17).
(a) -CONST, Water cut
(b) -CONST, GOR
(c) -CONST, BHP
(d) -GEO1, Water cut
(e) -GEO1, GOR
(f) -GEO1, BHP
(g) -GEO2, Water cut
(h) -GEO2, GOR
(i) -GEO2, BHP
Figure 14: Predicted production data for three wells. Red dots are the observations, gray curves are the predictions from the prior and blue curves are the predictions from the posterior ensemble (Field Case).

6 Conclusions

This paper presented the results of an investigation on the use of ES-MDA with geometric inflation factors. The results warranted the following conclusions:

  • For the three test cases, the use of constant inflation factors resulted in the lowest values for the data-mismatch norm and the highest values for the model-change norm.

  • The procedure proposed by Rafiee and Reynolds (2017) obtained a good balance between data match and model change as longs as an appropriate number of data assimilation is select a priori. For the three test cases, this method selected a large value for the initial inflation. As a result, it is necessary to use a larger number of data assimilations to obtain acceptable results.

  • The procedure introduced in this paper avoided this problem by selecting the geometric sequence based on the inflation at the last data assimilation. Overall, the initial inflation required to satisfy the MDP is significantly smaller than the initial inflation obtained by (Rafiee and Reynolds, 2017).

  • The procedure introduced in this paper resulted in a trade-off between data match and model change for the three test problems.

  • The use of geometric sequences also resulted in posterior ensembles with larger variances, which is typically desirable for uncertainty quantification, especially considering the fact that ensemble data assimilation methods tend to underestimate the posterior variance.

  • Overall, the results are in agreement with the conjecture that starting with large inflations and gradually decreasing during the data assimilation steps improves the final model estimates.

Acknowledgement

The author would like to thank Petrobras for supporting this research and for the permission to publish this paper.

References

  • Avansi and Schiozer (2015) Avansi, G. D. and Schiozer, D. J. UNISIM-I: synthetic model for reservoir development and management applications. International Journal of Modeling and Simulation for the Petroleum Industry, 9(1):21–30, 2015. URL http://www.ijmspi.org/ojs/index.php/ijmspi/article/view/152.
  • Emerick (2016a) Emerick, A. A. Analysis of the performance of ensemble-based assimilation of production and seismic data. Journal of Petroleum Science and Engineering, 139:219–239, 2016a. doi: 10.1016/j.petrol.2016.01.029.
  • Emerick (2016b) Emerick, A. A. Estimating uncertainty bounds in field production using ensemble-based methods. Journal of Petroleum Science and Engineering, 145:648–656, 2016b. doi: 10.1016/j.petrol.2016.06.037.
  • Emerick and Reynolds (2011) Emerick, A. A. and Reynolds, A. C.

    History matching a field case using the ensemble Kalman filter with covariance localization.

    SPE Reservoir Evaluation & Engineering, 14(4):423–432, 2011. doi: 10.2118/141216-PA.
  • Emerick and Reynolds (2012) Emerick, A. A. and Reynolds, A. C. History matching time-lapse seismic data using the ensemble Kalman filter with multiple data assimilations. Computational Geosciences, 16(3):639–659, 2012. doi: 10.1007/s10596-012-9275-5.
  • Emerick and Reynolds (2013a) Emerick, A. A. and Reynolds, A. C. Ensemble smoother with multiple data assimilation. Computers & Geosciences, 55:3–15, 2013a. doi: 10.1016/j.cageo.2012.03.011.
  • Emerick and Reynolds (2013b) Emerick, A. A. and Reynolds, A. C. History matching of production and seismic data for a real field case using the ensemble smoother with multiple data assimilation. In Proceedings of the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 18–20 February, number SPE-163675-MS, 2013b. doi: 10.2118/163675-MS.
  • Evensen (2004) Evensen, G. Sampling strategies and square root analysis schemes for the EnKF. Ocean Dynamics, 54(6):539–560, 2004. doi: 10.1007/s10236-004-0099-2.
  • Evensen and Eikrem (2018) Evensen, G. and Eikrem, K. S. Conditioning reservoir models on rate data using ensemble smoothers. Computational Geosciences, First Online, 2018. doi: 10.1007/s10596-018-9750-8.
  • Hamarik et al. (2012) Hamarik, U., Palm, R., and Raus, T. A family of rules for parameter choice in Tikhonov regularization of ill-posed problems with inexact noise level. Journal of Computational and Applied Mathematics, 236:2146–2157, 2012. doi: 10.1016/j.cam.2011.09.037.
  • Hanke (1997) Hanke, M. A regularizing Levenberg-Marquardt scheme, with applications to inverse groundwater filtration problems. Inverse Problems, 13(1):79–95, 1997. doi: 10.1088/0266-5611/13/1/007.
  • Hanke (2010) Hanke, M. The regularizing Levenberg-Marquardt scheme is of optimal order. Jounal of Integral Equations Applications, 22(2):259–283, 2010. doi: 10.1216/JIE-2010-22-2-259.
  • Hansen and O’Leary (1993) Hansen, P. C. and O’Leary, D. P. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific Computing, 14(6):1487–1503, 1993. doi: 10.1137/0914086.
  • Houtekamer and Mitchell (2001) Houtekamer, P. L. and Mitchell, H. L. A sequential ensemble Kalman filter for atmospheric data assimilation. Monthly Weather Review, 129(1):123–137, 2001. doi: 10.1175/1520-0493(2001)129¡0123:ASEKFF¿2.0.CO;2.
  • Iglesias and Dawson (2013) Iglesias, M. A. and Dawson, C. The regularizing Levenberg-Marquardt scheme for history matching of petroleum reservoirs. Computational Geosciences, 17(6):1033–1053, 2013. doi: 10.1007/s10596-013-9373-z.
  • Kaipio and Somersalo (2005) Kaipio, J. P. and Somersalo, E. Statistical and Computational Inverse Problems. Springer New York, 2005. doi: 10.1007/b138659.
  • Le et al. (2016) Le, D. H., Emerick, A. A., and Reynolds, A. C. An adaptive ensemble smoother with multiple data assimilation for assisted history matching. SPE Journal, 21(6):2195–2207, 2016. doi: 10.2118/173214-PA.
  • Luo and Bhakta (2017) Luo, X. and Bhakta, T. Estimating observation error covariance matrix of seismic data from a perspective of image denoising. Computational Geosciences, 21(2):205–222, 2017. doi: 10.1007/s10596-016-9605-0.
  • Luo et al. (2015) Luo, X., Stordal, A. S., Lorentzen, R. J., and Nævdal, G. Iterative ensemble smoother as an approximate solution to a regularized minimum-average-cost problem: Theory and applications. SPE Journal, 20(5), 2015. doi: 10.2118/176023-PA.
  • Mannseth (2015) Mannseth, T. Comparison of five different ways to assimilate data for a simplistic weakly nonlinear parameter estimation problem. Computational Geosciences, 19(4):791–804, 2015. doi: 10.1007/s10596-015-9490-y.
  • Maschio et al. (2013) Maschio, C., Avansi, G. D., Santos, A. A., and Schiozer, D. J. UNISIM-I-H: case study for history matching. Dataset, 2013. URL www.unisim.cepetro.unicamp.br/benchmarks/br/unisim-i/unisim-i-h.
  • Morozov (1984) Morozov, V. Methods for Solving Incorrectly Posed Problems. Springer-Verlag New York, 1984. doi: 10.1007/978-1-4612-5280-1.
  • Nocedal and Wright (2006) Nocedal, J. and Wright, S. J. Numerical Optimization. Springer, New York, 2006.
  • Oliver and Alfonzo (2017) Oliver, D. S. and Alfonzo, M. Calibration of imperfect models to biased observations. Computational Geosciences, First Online, 2017. doi: 10.1007/s10596-017-9678-4.
  • Oliver et al. (2008) Oliver, D. S., Reynolds, A. C., and Liu, N. Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge University Press, Cambridge, UK, 2008.
  • Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, Cambridge, England, 3rd edition, 2007.
  • Rafiee and Reynolds (2017) Rafiee, J. and Reynolds, A. C. Theoretical and efficient practical procedures for the generation of inflation factors for ES-MDA. Inverse Problems, 33(11):115003, 2017. doi: 10.1088/1361-6420/aa8cb2.
  • Stordal and Elsheikh (2015) Stordal, A. S. and Elsheikh, A. H. Iterative ensemble smoothers in the annealed importance sampling framework. Advances in Water Resources, 86:231–239, 2015. doi: 10.1016/j.advwatres.2015.09.030.
  • Sun et al. (2017) Sun, W., Vink, J. C., and Gao, G. A practical method to mitigate spurious uncertainty reduction in history matching workflows with imperfect reservoir models. In Proceedings of the SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 20-22 February, number SPE-182599-MS, 2017. doi: 10.2118/182599-MS.
  • Tavakoli and Reynolds (2010) Tavakoli, R. and Reynolds, A. C. History matching with parameterization based on the SVD of a dimensionless sensitivity matrix. SPE Journal, 15(12):495–508, 2010. doi: 10.2118/118952-PA.
  • Zhang et al. (2002) Zhang, F., Reynolds, A. C., and Oliver, D. S. Evaluation of the reduction in uncertainty obtained by conditioning a 3D stochastic channel to multiwell pressure data. Mathematical Geology, 34(6):713–740, 2002. doi: 10.1023/A:1019805310025.
  • Zhao et al. (2006) Zhao, Y., Li, G., and Reynolds, A. C. Characterizing the measurement error with the EM algorithm. In Proceedings of 10th European Conference on the Mathematics of Oil Recovery, Amsterdam, 4–7 September, 2006. doi: 10.3997/2214-4609.201402491.