Data-Space Inversion with Ensemble Smoother

TL;DR

This paper presents an improved data-space inversion method using an iterative ensemble smoother to enhance reservoir production forecasting, offering faster and more robust predictions without updating the geological models.

Contribution

It introduces a new DSI implementation with iterative ensemble smoothing, demonstrating improved speed and robustness over previous PCA-based methods in reservoir history matching.

Findings

Faster and more robust DSI implementation.

Comparable forecast accuracy to traditional history matching.

Effective in complex reservoir models with long production histories.

Abstract

Reservoir engineers use large-scale numerical models to predict the production performance in oil and gas fields. However, these models are constructed based on scarce and often inaccurate data, making their predictions highly uncertain. On the other hand, measurements of pressure and flow rates are constantly collected during the operation of the field. The assimilation of these data into the reservoir models (history matching) helps to mitigate uncertainty and improve their predictive capacity. History matching is a nonlinear inverse problem, which is typically handled using optimization and Monte Carlo methods. In practice, however, generating a set of properly history-matched models that preserve the geological realism is very challenging, especially in cases with complicated prior description, such as models with fractures and complex facies distributions. Recently, a new data-space inversion (DSI) approach was introduced in the literature as an alternative to the model-space inversion used in history matching. The essential idea is to update directly the predictions from a prior ensemble of models to account for the observed production history without updating the corresponding models. The present paper introduces a DSI implementation based on the use of an iterative ensemble smoother and demonstrates with examples that the new implementation is computationally faster and more robust than the earlier method based on principal component analysis. The new DSI is also applied to estimate the production forecast in a real field with long production history and a large number of wells. For this field problem, the new DSI obtained forecasts comparable with a more traditional ensemble-based history matching.

Tables

Table 1. Table 1: Geostatistical parameters used to create the reference model and the prior ensemble. Test case 1

Parameter	Porosity	Log-permeability
Mean	0.22	7.2 ln-mD
Standard deviation	0.05	0.60 ln-mD
Variogram type	Spherical	Spherical
Variogram maximum range	1124.0 m	1124.0 m
Variogram minimum range	281.0 m	281.0 m
Azimuth	45∘	45∘

Table 2. Table 2: CPU time for inversion (70 data points). Test case 1

Method	CPU time (seconds)
DSI	1506
DSI-ESMDA	2

Table 3. Table 3: Prior distribution of global parameters. UNISIM-I-H

Parameter	Mode	Min.	Max.
Rock compressibility in $\left({\text{cm}}^2/\text{kgf}\right)$	$5.3\times {10}^{-5}$	$1\times {10}^{-5}$	$9.6\times {10}^{-5}$
Depth of the oil-water contact at the East block (Fig. 7)	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 4. Table 4: Normalized data-mismatch objective function. UNISIM-I-H

Case	Mean	Standard deviation
Prior	456.7	333.462
DSI	26.6	7.379
DSI-ESMDA	2.7	0.003
DSI-ESMDA with localization	1.9	0.008
ES-MDA	3.0	0.109

Table 5. Table 5: Normalized data-mismatch objective function. Field case

Case	Mean	Standard deviation
Prior	32.19	11.309
DSI-ESMDA (with localization)	0.55	0.004
ES-MDA	4.11	0.581

Equations

$$\mathbf{d}=\left[\begin{array}[]{c}\mathbf{d}_{\textrm{h}}\\ \mathbf{d}_{\textrm{f}}\end{array}\right].$$

$${\mathbf{d}}_{\text{obs}}={\mathbf{d}}_{\text{true}}+{\mathbf{e}}_{\text{d}},$$

$${\mathbf{d}}_{\text{true}}={\mathbf{d}}_{\text{h}}\left({\mathbf{m}}_{\text{true}}\right)+{\mathbf{e}}_{\text{m}},$$

$$\mathcal{L}(\mathbf{m}|{\mathbf{d}}_{\text{obs}})=\text{const}\times \exp \left(-\frac{1}{2}{\left({\mathbf{d}}_{\text{obs}}-{\mathbf{d}}_{\text{h}}(\mathbf{m})\right)}^{\top }{\mathbf{C}}_{\mathbf{e}}^{-1}\left({\mathbf{d}}_{\text{obs}}-{\mathbf{d}}_{\text{h}}(\mathbf{m})\right)\right),$$

$${\mathbf{C}}_{\mathbf{e}}=\text{cov}\left[{\mathbf{e}}_{\text{d}}\right]+\text{cov}\left[{\mathbf{e}}_{\text{m}}\right]$$

$$p(\mathbf{m}|{\mathbf{d}}_{\text{obs}})$$

$$\mathcal{O}(\mathbf{m})=\frac{1}{2}{\left({\mathbf{d}}_{\text{obs}}-{\mathbf{d}}_{\text{h}}(\mathbf{m})\right)}^{\top }{\mathbf{C}}_{\mathbf{e}}^{-1}\left({\mathbf{d}}_{\text{obs}}-{\mathbf{d}}_{\text{h}}(\mathbf{m})\right)+\frac{1}{2}{\left(\mathbf{m}-{\mathbf{m}}_{\text{pr}}\right)}^{\top }{\mathbf{C}}_{\mathbf{m}}^{-1}\left(\mathbf{m}-{\mathbf{m}}_{\text{pr}}\right).$$

$${\mathcal{O}}_{\text{rml}}(\mathbf{m})=\frac{1}{2}{\left({\mathbf{d}}_{\text{obs}}^{\ast }-{\mathbf{d}}_{\text{h}}(\mathbf{m})\right)}^{\top }{\mathbf{C}}_{\mathbf{e}}^{-1}\left({\mathbf{d}}_{\text{obs}}^{\ast }-{\mathbf{d}}_{\text{h}}(\mathbf{m})\right)+\frac{1}{2}{\left(\mathbf{m}-{\mathbf{m}}^{\ast }\right)}^{\top }{\mathbf{C}}_{\mathbf{m}}^{-1}\left(\mathbf{m}-{\mathbf{m}}^{\ast }\right),$$

$${\mathbf{d}}_{\text{pca}}=\overline{\mathbf{d}}+{\mathbf{C}}_{\mathbf{d}}^{\text{1/2}}\mathbf{x},$$

$$\overline{\mathbf{d}}=\frac{1}{N_e}\sum _{j=1}^{N_e}{\mathbf{d}}_j$$

$${\mathbf{C}}_{\mathbf{d}}$$

$$\Delta \mathbf{D}=\frac{1}{\sqrt{N_e-1}}\left[{\mathbf{d}}_1-\overline{\mathbf{d}},\dots ,{\mathbf{d}}_{N_e}-\overline{\mathbf{d}}\right].$$

$$\Delta \mathbf{D}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\top },$$

$${\mathbf{C}}_{\mathbf{d}}^{\text{1/2}}=\mathbf{U}\mathbf{\Sigma }.$$

$$\frac{{\sum }_{i=1}^{N_r}{\sigma }_i}{{\sum }_{i=1}^{\min \{N_d,N_e-1\}}{\sigma }_i}\ge \xi$$

$${\hat{d}}_i={\text{cdf}}_1^{-1}\left({\text{cdf}}_2\left(d_{\text{pca},i}\right)\right),$$

$${\mathcal{O}}_{\text{rml}}(\mathbf{x})=\frac{1}{2}{\left({\mathbf{d}}_{\text{obs}}^{\ast }-{\hat{\mathbf{d}}}_{\text{h}}(\mathbf{x})\right)}^{\top }{\mathbf{C}}_{\mathbf{e}}^{-1}\left({\mathbf{d}}_{\text{obs}}^{\ast }-{\hat{\mathbf{d}}}_{\text{h}}(\mathbf{x})\right)+\frac{1}{2}{\left(\mathbf{x}-{\mathbf{x}}^{\ast }\right)}^{\top }\left(\mathbf{x}-{\mathbf{x}}^{\ast }\right),$$

$${\mathbf{d}}_j^{k+1}={\mathbf{d}}_j^k+\mathbf{R}\circ {\mathbf{K}}^k\left({\mathbf{d}}_{\text{obs}}+\sqrt{{\alpha }_k}{\mathbf{e}}_j^k-{\mathbf{d}}_{\text{h}}^k\right),$$

$${\mathbf{K}}^k=\Delta {\mathbf{D}}^k{\left(\Delta {\mathbf{D}}_{\text{h}}^k\right)}^{\top }{\left(\Delta {\mathbf{D}}_{\text{h}}^k{\left(\Delta {\mathbf{D}}_{\text{h}}^k\right)}^{\top }+{\alpha }_k{\mathbf{C}}_{\mathbf{e}}\right)}^{-1}.$$

$$\Delta {\mathbf{D}}_{\text{h}}=\frac{1}{\sqrt{N_e-1}}\left[{\mathbf{d}}_{\text{h},1}-\overline{{\mathbf{d}}_{\text{h}}},\dots ,{\mathbf{d}}_{\text{h},N_e}-\overline{{\mathbf{d}}_{\text{h}}}\right].$$

$$r\left(\frac{h}{L}\right)=\left\{\begin{array}[]{ll}-\frac{1}{4}\left(\frac{h}{L}\right)^{5}+\frac{1}{2}\left(\frac{h}{L}\right)^{4}+\frac{5}{8}\left(\frac{h}{L}\right)^{3}-\frac{5}{3}\left(\frac{h}{L}\right)^{2}+1,&\text{if}~{}0\leq\frac{h}{L}\leq 1\\ \frac{1}{12}\left(\frac{h}{L}\right)^{5}-\frac{1}{2}\left(\frac{h}{L}\right)^{4}+\frac{5}{8}\left(\frac{h}{L}\right)^{3}+\frac{5}{3}\left(\frac{h}{L}\right)^{2}-5\left(\frac{h}{L}\right)+4-\frac{2}{3}\left(\frac{h}{L}\right)^{-1},&\text{if}~{}1\leq\frac{h}{L}\leq 2\\ 0&\text{if}~{}\frac{h}{L}>2\end{array}\right.,$$

$$\frac{h}{L}=\sqrt{{\left(\frac{\Delta x^{\prime }}{L_{x^{\prime }}}\right)}^2+{\left(\frac{\Delta y^{\prime }}{L_{x^{\prime }}}\right)}^2+{\left(\frac{\Delta t}{T}\right)}^2},$$

$$\left[\begin{array}[]{c}\Delta x^{\prime}\\ \Delta y^{\prime}\end{array}\right]=\left[\begin{array}[]{cc}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\end{array}\right]\left[\begin{array}[]{c}\Delta x\\ \Delta y\end{array}\right].$$

$${\mathbf{w}}_j^k\equiv {\left(\Delta {\mathbf{D}}_{\text{h}}^k{\left(\Delta {\mathbf{D}}_{\text{h}}^k\right)}^{\top }+{\alpha }_k{\mathbf{C}}_{\mathbf{e}}\right)}^{-1}\left({\mathbf{d}}_{\text{obs}}+\sqrt{{\alpha }_k}{\mathbf{e}}_j^k-{\mathbf{d}}_{\text{h}}^k\right).$$

$${\mathbf{d}}_j^{k+1}$$

$${\mathbf{d}}_j^{k+1}=\sum _{i=1}^{N_e}{\gamma }_{i,j}^k{\mathbf{d}}_i^k,$$

$${\mathbf{d}}_j^{k+1}$$

$$d_{n,j}^{k+1}$$

$$\widehat{\Delta \mathbf{D}}\equiv {\mathbf{C}}_{\mathbf{e}}^{-\text{1/2}}\Delta \mathbf{D}=\hat{\mathbf{U}}\hat{\mathbf{\Sigma }}{\hat{\mathbf{V}}}^{\top },$$

$${\mathbf{C}}_{\mathbf{d}}^{\text{1/2}}={\mathbf{C}}_{\mathbf{e}}^{\text{1/2}}\hat{\mathbf{U}}\hat{\mathbf{\Sigma }}.$$

Full text

Data-Space Inversion with Ensemble Smoother

Mateus M. Lima1, Alexandre A. Emerick2 and Carlos E. P. Ortiz3

11footnotetext: Petrobras and UENF (mateusmartins@yahoo.com.br)22footnotetext: Petrobras (aemerick@gmail.com)33footnotetext: UENF (capico.LENEP@gmail.com)

Abstract

Reservoir engineers use large-scale numerical models to predict the production performance in oil and gas fields. However, these models are constructed based on scarce and often inaccurate data, making their predictions highly uncertain. On the other hand, measurements of pressure and flow rates are constantly collected during the operation of the field. The assimilation of these data into the reservoir models (history matching) helps to mitigate uncertainty and improve their predictive capacity. History matching is a nonlinear inverse problem, which is typically handled using optimization and Monte Carlo methods. In practice, however, generating a set of properly history-matched models that preserve the geological realism is very challenging, especially in cases with complicated prior description, such as models with fractures and complex facies distributions. Recently, a new data-space inversion (DSI) approach was introduced in the literature as an alternative to the model-space inversion used in history matching. The essential idea is to update directly the predictions from a prior ensemble of models to account for the observed production history without updating the corresponding models. The present paper introduces a DSI implementation based on the use of an iterative ensemble smoother and demonstrates with examples that the new implementation is computationally faster and more robust than the earlier method based on principal component analysis. The new DSI is also applied to estimate the production forecast in a real field with long production history and a large number of wells. For this field problem, the new DSI obtained forecasts comparable with a more traditional ensemble-based history matching.

Keywords: Data-space inversion; uncertainty quantification; ensemble smoother; history matching.

1 Introduction

Reservoir characterization from static and dynamic data allows to create numerical models that can be used to simulate the performance of petroleum reservoirs under different operating conditions. These models are essential for efficient exploitation and management of oil and gas fields. The incorporation of static data is typically done using geostatistics while the incorporation of dynamic data is done using history-matching methods.

History matching is usually a very difficult task, which involves the integration of interdisciplinary teams and an intensive use of computational resources. A complete study may require months of work and the results are not always satisfactory. The task is even harder if one needs to provide uncertainty estimates, in which case several alternative history-matched models must be generated. In the last few decades, the advances in assisted (or semi-automatic) history-matching techniques were notorious. Yet, history matching remains one of the most time-consuming steps of a field study because the size and complexity of the models have also increased significantly in the same period. Oliver and Chen [36] present a review of the main history-matching methods proposed in the literature. Among these methods, the ones based on the ensemble Kalman filter (EnKF) [17, 19] have become quite popular, especially because of their ease of implementation and integration with commercial reservoir simulators and the ability to generate multiple models with large number of uncertainty parameters at an affordable computational cost. Despite the relative success in a number of recent field cases reported in the literature; see, for example, [16, 7, 9, 33, 1, 20, 28], generating a set of models properly conditioned to all historical data and still preserving the geological realism is very challenging, especially in cases with complicated prior description, such as models with fractures and complex facies distributions.

A new approach known as data-space inversion (DSI) [48] has drawn attention in the literature as an alternative to the model-space inversion approach used in history matching. The basic idea behind DSI is to update directly the predictions from a prior ensemble of models to account for the observed production history without updating the corresponding models. The upside of this approach is to be able to provide an ensemble of forecasts without going through the time-consuming history-matching step. Because there are no model inversions, there are no concerns about losing the geological realism. The downside is that DSI provides only the forecast estimates for a fixed production strategy. In practice, however, reservoir engineers may also be interested in having the corresponding models to study different drainage strategies. For this reason, data- and model-space inversions should be considered complementary rather than alternative (or even competing) approaches.

The DSI method introduced by Sun and Durlofsky [48] uses principal component analysis (PCA) to reparameterize the predicted data from the prior ensemble into a lower dimensional space and the randomized maximum likelihood (RML) method [37] to generate samples of the posterior distribution of predicted data given the observations. The authors used a data transformation before PCA to improve linearity. Sun et al. [49] extended the original DSI method introducing a more general data transformation procedure. They tested the method in a model of a complex fractured reservoir and obtained reasonable uncertainty estimates of the production forecast. Jiang [25] modified the DSI approach to allow changes in the well controls during the forecast period so that the method could be used for life-cycle optimization. They noted that the method required a larger prior ensemble to better represent a wider range of possibilities. They showed that the proposed method combined to a direct search optimization algorithm [2] was able to improve the expected net-present value of a reservoir model.

Similar ideas of DSI have appeared before in the literature. For example, Krishnamurti et al. [26] and Pagowski et al. [38] applied linear regression to combine the ozone forecast of ensembles of models. In the atmospheric literature these methods are referred to as aggregation methods or ensemble forecast [31]. Emerick [8] also used a DSI-type of approach to invert 4D seismic impedance data directly to pressure and water saturation. Scheidt et al. [43] proposed a method named prediction-focused analysis (PFA) based on projecting the prior predictions into a low-dimensional space and using kernel smoothing to estimate the joint distribution of historical and forecasted data. Using the joint distribution, they could predict the uncertainty estimates in a tracer transport problem. However, PFA method seems applicable only to problems with few data points because it requires the projection to very few dimensions (two or three dimension in the examples presented in the paper). Satija and Caers [41] modified the PFA method using canonical functional component analysis to improve the linearity in the projected data. Satija and Caers [42] used the same approach in a reservoir problem and concluded that the method provided uncertainty estimates of production forecast in reasonable agreement with rejection sampling. More recently, He et al. [22] used similar ideas from DSI to estimate the uncertainty reduction in a study to compute the value of information of data-acquisition plans. Jeong et al. [24] applied machine learning techniques (neural networks and support vector regression) to DSI. They concluded that the method can be a more efficient alternative to the computationally demanding history-matching methods, however the method fails to provide satisfactory forecast if the predictions from the prior ensemble (training set) are too far from the expected true response.

In the present paper, we introduce a new DSI implementation based on the use of an iterative ensemble smoother and demonstrate with examples that the new DSI is computationally faster and more robust than the procedure proposed in [48, 49]. Moreover, we apply the new DSI to a real field case with long production history and large number of wells and show that the method provides forecasts comparable with a more traditional ensemble-based history-matching process. The rest of the paper is organized as follows: Section 2 reviews the DSI method proposed in [48, 49] and the new DSI method. Section 3 presents three reservoir test problems. The first problem is a small synthetic case used to demonstrate that the proposed method provides results similar to the original DSI with a lower computational cost. The second problem is a benchmark history-matching case [3, 32] constructed with data from real reservoir in Campos Basis. This problem is used to compare the methods is a more realistic situation with a large number of data points. The last problem corresponds to a real brown-field case where the proposed method is compared against an ensemble-based history matching. Section 5 summarizes the conclusions of the paper.

2 Methodology

2.1 Preliminaries

Let $\mathbf{m}\in\mathds{R}^{N_{m}}$ denote the vector of uncertain parameters of a reservoir model with a historical production period $t_{\textrm{h}}$. Our goal is to predict the production performance for a period $t_{\textrm{f}}$ after $t_{\textrm{h}}$. Let $\mathbf{d}\in\mathds{R}^{N_{d}}$ denote the vector of predicted production data, which is a nonlinear function of $\mathbf{m}$, that is, $\mathbf{d}=\mathbf{d}(\mathbf{m})$. In the applications of interest of this paper, $\mathbf{d}$ is the result of a reservoir simulation. This vector contains predicted data from both, history, $\mathbf{d}_{\textrm{h}}$, and forecast periods, $\mathbf{d}_{\textrm{f}}$, that is,

[TABLE]

Let $\mathbf{d}_{\textrm{obs}}\in\mathds{R}^{N_{d,\textrm{h}}}$ denote the vector of field observations, which is corrupted with an additive random noise, $\mathbf{e}_{\textrm{d}}$, due to measurement errors such that

[TABLE]

where $\mathbf{d}_{\textrm{true}}$ is the true (noiseless) data. Moreover, assume that model errors are also additive and

[TABLE]

where $\mathbf{m}_{\textrm{true}}$ is the vector containing the “true” values for the model parameters and $\mathbf{e}_{\textrm{m}}$ is a vector of model errors. Under these conditions, it is straightforward to show [51] that the likelihood of $\mathbf{m}$ is

[TABLE]

where

[TABLE]

is the total data-error covariance matrix. Note that we use the term “model error” to refer to any imperfection in the model used to represent the real reservoir. Examples of sources of model errors include numerical and discretization errors, simplifications of the physics and insufficient parameterization. The assumption in Eq. 3 is that all sources of model errors can be aggregated into a random vector with $\textrm{E}\left[\mathbf{e}_{\textrm{m}}\right]=\mathbf{0}$ and known covariance. Even though these are strong assumptions which are unlike to hold in reality, it is important to note that including $\textrm{cov}\left[\mathbf{e}_{\textrm{m}}\right]$ in the matrix $\mathbf{C}_{\mathbf{e}}$ reduces the weights attributed to data helping to partially compensate for deficiencies in the models; see, for example, [50, 35].

If the prior model follows a multivariate Gaussian distribution, $\mathbf{m}\sim\mathcal{N}(\mathbf{m}_{\textrm{pr}},\mathbf{C}_{\mathbf{m}})$, then the posterior probability density function (PDF) of $\mathbf{m}$ given $\mathbf{d}_{\textrm{obs}}$ has the form

[TABLE]

where

[TABLE]

The model $\mathbf{m}$ that minimizes $\mathcal{O}(\mathbf{m})$ corresponds to the maximum a posteriori [51]. In practice, however, we are interested in sampling the posterior PDF to quantify uncertainty. In this case, one alternative is the RML method [37], which provides an approximate sampling of $p(\mathbf{m}|\mathbf{d}_{\textrm{obs}})$. Each RML sample is obtained by minimizing a modified version of Eq. 7 given by

[TABLE]

where $\mathbf{d}_{\textrm{obs}}^{\ast}\sim\mathcal{N}(\mathbf{d}_{\textrm{obs}},\mathbf{C}_{\mathbf{e}})$ and $\mathbf{m}^{\ast}\sim\mathcal{N}(\mathbf{m}_{\textrm{pr}},\mathbf{C}_{\mathbf{m}})$.

2.2 Data-Space Inversion

In this section, we review the data-space inversion (DSI) procedure as proposed in [48] and later improved in [49]. The main idea behind the method is to use PCA to write the vector of predicted data as

[TABLE]

where $\overline{\mathbf{d}}$ and $\mathbf{C}_{\mathbf{d}}$ are the mean and covariance of $\mathbf{d}$, respectively. Both are computed using a prior ensemble, that is,

[TABLE]

and

[TABLE]

where

[TABLE]

The square root of $\mathbf{C}_{\mathbf{d}}$ in Eq. 9 is computed using the singular value decomposition (SVD) of $\Delta\mathbf{D}$

[TABLE]

where $\mathbf{U}$ is a $N_{d}\times N_{d}$ orthogonal matrix containing the left singular values of $\Delta\mathbf{D}$, which are equivalent to the left eigenvectors of $\mathbf{C}_{\mathbf{d}}$; $\bm{\Sigma}$ is a $N_{d}\times N_{e}$ matrix containing as non-zero elements the singular values of $\Delta\mathbf{D}$, or, equivalently, the square root of the eigenvalues of $\mathbf{C}_{\mathbf{d}}$. The matrix $\mathbf{V}$ contains the right singular values of $\Delta\mathbf{D}$. The square-root of $\mathbf{C}_{\mathbf{d}}$ becomes

[TABLE]

The vector $\mathbf{x}$ in Eq. 9 is a sample from a standard normal distribution, that is, $\mathbf{x}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. In practice, we truncate small singular values using an energy criterium, which means that we consider only the $N_{r}\leq\max\{N_{d},N_{e}-1\}$ largest singular values such that

[TABLE]

where $\xi\leq 1$ is the energy threshold, typically selected between 0.9 and 0.99 and $\sigma_{i}$ is the $i$th singular value of $\Delta\mathbf{D}$. Truncating small singular values has two positive side effects: it reduces the dimension of $\mathbf{x}$ and introduces regularization in the inversion. Both are important because the DSI method uses RML for sampling, which requires solving several minimization problems.

Sun and Durlofsky [48] noted that the direct application of Eq. 9 may result in nonphysical values for the predicted data, for example negative production or pressure. According to the authors, this problem occurs mainly before water breakthrough time. Therefore, they proposed to apply a data transformation to the prior realizations of $\mathbf{d}$ before PCA. The transformation is based on shifting and compressing/stretching the time series. They claim and illustrate in a example that the transformed vectors, $\widehat{\mathbf{d}}$, have a more Gaussian prior distribution. However, this procedure is difficult to apply in cases with frequent changes in well controls. In [49], the authors proposed to use an inverse Gaussian anamorphosis procedure using the empirical cumulative density function (CDF) computed using the $N_{e}$ prior realizations of $\mathbf{d}$. Figure 1 illustrates the process, where each component of the transformed vector $\widehat{\mathbf{d}}$ is computed as

[TABLE]

where $\textrm{cdf}_{1}(\cdot)$ and $\textrm{cdf}_{2}(\cdot)$ are the CDF’s of $\mathbf{d}$ and $\mathbf{d}_{\textrm{pca}}$, respectively.

The final step of DSI is to use RML to generate a posterior ensemble of predicted data. The RML objective function can be written in terms of the vector of PCA coefficients as

[TABLE]

where $\mathbf{x}^{\ast}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. The minimization of Eq. 17 can be done with any optimization method. Sun and Durlofsky [48] used the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method [34], in which case it is necessary to compute the gradient of $\mathcal{O}_{\textrm{rml}}(\mathbf{x})$ with respect to $\mathbf{x}$. However, $\widehat{\mathbf{d}}_{\textrm{h}}(\mathbf{x})$ is a nonlinear function of $\mathbf{x}$ with no general analytical form because of the transformation of Eq. 16. One alternative, is to compute numerical gradients, which may be computationally expensive if we have a large number of data points. Jiang [25] proposed to ignore the transformation of Eq. 16 and use the analytical gradients computed with $\mathbf{d}_{\textrm{pca}}$ instead of $\widehat{\mathbf{d}}$. Our limited set of tests indicated that this procedure in fact improved the computational performance of our DSI implementation, which uses the limited-memory BFGS [34] implementation available in the C# library Accord.NET [45].

2.3 Data-Space Inversion with Ensemble Smoother

The ensemble smoother (ES) was introduced by van Leeuwen and Evensen [52] as an alternative to the sequential data assimilation scheme of EnKF. The first application of ES for history matching was presented by Skjervheim et al. [44], which concluded that the method is faster than EnKF with similar results. Despite the good results presented in [44], some authors [5, 6, 15, 14] observed that ES tends to result in unreasonable data matches when applied to more complex history-matching problems. The reason for the poor performance is because ES is similar to applying a single Gauss-Newton iteration to minimize a RML-type of objective function with a sensitivity matrix estimated based on the prior ensemble [40]. This fact lead to the development of several iterative forms of ES; see, for example, [5, 6, 15, 46, 47, 29]. Among the iterative forms of ES, the ensemble smoother with multiple data assimilation (ES-MDA) [15] is a popular choice. The popularity of ES-MDA can be attributed mainly to its good performance in history-matching problems [16, 9, 33, 20] and its simplicity of implementation. In fact, ES-MDA is essentially equivalent to repeat ES a few times with the data-error covariance matrix, $\mathbf{C}_{\mathbf{e}}$, multiplied by coefficients $\alpha_{k}$’s to avoid overweighing the measurements. The choice of the coefficients and the number of repetitions of ES, $N_{a}$, must obey the condition $\sum_{k=1}^{N_{a}}\alpha_{k}^{\scriptsize-1}=1$ to ensure that ES-MDA samples the correct posterior PDF in the linear-Gaussian case [15]. For nonlinear problems, the choice of the $\alpha_{k}$’s and $N_{a}$ has a major impact in the performance of the method. There are recent works proposing methods to select $\alpha_{k}$’s and $N_{a}$ [27, 9, 39, 30, 10]. However, here we use simplest choice, which consists of selecting $N_{a}$ in advance and setting $\alpha_{k}=N_{a}$, for $k=1,\ldots,N_{a}$.

The application of ES-MDA to DSI, that is, to generate samples of the predicted data vector, $\mathbf{d}=[\mathbf{d}_{\textrm{h}}^{\scriptsize\top},\mathbf{d}_{\textrm{f}}^{\scriptsize\top}]^{\scriptsize\top}$, given the vector of observations, $\mathbf{d}_{\textrm{obs}}$, is straightforward. We can write the resulting DSI-ESMDA update equation as

[TABLE]

for $k=1,\ldots,N_{a}$ and $j=1,\ldots,N_{e}$, where $\mathbf{K}$ is a modified version of the Kalman gain given by

[TABLE]

In the above equations, $\Delta\mathbf{D}$ was defined before (Eq. 12) and $\Delta\mathbf{D}_{\textrm{h}}$ includes only the predicted data corresponding to the historical period, that is,

[TABLE]

The vector $\mathbf{e}$ is a sample from $\mathcal{N}(\mathbf{0},\mathbf{C}_{\mathbf{e}})$ and $\mathbf{R}$ is the localization matrix with “$\circ$” denoting the Schur (element-wise) product. The matrix inversion required in Eq. 19 is done using the subspace inversion method [18] as described in the Appendix section of [13]. This procedure involves the truncated SVD of a rescaled matrix $\mathbf{C}_{\mathbf{e}}^{\scriptsize-\textrm{1/2}}\Delta\mathbf{D}_{\textrm{h}}$ to avoid loss of relevant information when removing small singular values. The number of singular values retained is computed using the same energy criterium of Eq. 15.

The reason for introducing the Schur product between the localization matrix and the Kalman gain in Eq. 18 is twofold: it regularizes the estimates of the Kalman gain removing long-distance spurious correlations and increases the degrees of freedom to assimilate data [23]. In fact, localization proved to improve substantially the results of ensemble methods applied to history matching; see, for example, [4, 12, 11].

The regularization introduced by the Schur product is obtained by constructing the localization matrix using a correlation function with compact support. A common choice is the fifth-order compact correlation function from [21], in which case each entry of the matrix $\mathbf{R}$ is computed using

[TABLE]

where $h$ is a “distance” and $L$ is a parameter known as “critical length,” which corresponds to the distance where the correlation function decays to approximately 0.21. Note that the matrix $\mathbf{R}$ has the same shape of the Kalman gain and that each row of the Kalman gain corresponds to the correction term applied to each variable updated with Eq. 18, that is, entries of the vector $\mathbf{d}$. Each column of the Kalman gain corresponds to an observation used in the conditioning process. Therefore, $h$ in Eq. 21 corresponds to the distance between entries of $\mathbf{d}$ and $\mathbf{d}_{\textrm{obs}}$. In the applications of interest of this paper, both $\mathbf{d}$ and $\mathbf{d}_{\textrm{obs}}$ contain data from wells in the field. Therefore, the ratio $h/L$ is computed based on the spatial distance between wells. Besides the spatial distance, we also introduced the time difference between data points in the calculation of $h/L$ to account for the fact that early data have lower correlations with the predictions than the late data. The final expression for computing $h/L$ is given by

[TABLE]

where

[TABLE]

In the above equations, $\Delta x$ and $\Delta y$ correspond to the spatial distance between the wells in the orthogonal directions $x$ and $y$, while $\Delta x^{\prime}$ and $\Delta y^{\prime}$ are the corresponding distances rotated by an angle $\theta$. This allows to consider anisotropy in the localization function as discussed in [12]. $\Delta t$ is the time difference between data points. $L_{x^{\prime}}$, $L_{y^{\prime}}$ and $T$ are the critical lengths in the “directions” $x^{\prime}$, $y^{\prime}$ and $t$, respectively.

Besides removing long-distance spurious correlations, localization also increases the degrees of freedom to assimilate data. To show that, first consider the case without localization. Using Eq. 19 in (18), we can define the vector $\mathbf{w}_{j}^{k}$ as

[TABLE]

Hence

[TABLE]

where $\beta_{i,j}^{k}\equiv\left(\mathbf{d}_{\textrm{h},i}^{k}-\overline{\mathbf{d}_{\textrm{h}}^{k}}\right)^{\scriptsize\top}\mathbf{w}_{j}^{k}$ is a scalar. The vector $\overline{\mathbf{d}^{k}}$ is a linear combination of the vectors $\mathbf{d}_{i}^{k}$’s. Hence, it is possible to write $\mathbf{d}_{j}^{k+1}$ as

[TABLE]

where $\gamma_{i,j}^{k}$’s are scalars. Eq. 26 means that $\mathbf{d}_{j}^{k+1}$ is a linear combination of the vectors $\mathbf{d}_{i}^{k}$’s. Applying this result recursively to all MDA iterations, we conclude that the final predictions from DSI-ESMDA without localization are simply linear combinations of the prior ones. This effectively means that we have at most $N_{e}$ coefficients (degrees of freedom) available to assimilate data. Note that the number of degrees of freedom of the standard DSI is also limited by the number of PCA coefficients which is $N_{r}\leq\min\{N_{d},N_{e}-1\}$.

Considering now the case with localization, we can write

[TABLE]

where $\delta\mathbf{d}_{j}^{k}\equiv\mathbf{d}_{\textrm{obs}}+\sqrt{\alpha_{k}}\mathbf{e}_{j}^{k}-\mathbf{d}_{\textrm{h}}^{k}$ and $\delta d_{i,j}^{k}$ is its $i$th entry. $\mathbf{r}_{i}$ and $\bm{\kappa}^{k}_{i}$ correspond to the $i$th columns of the localization and Kalman gain matrices, respectively. Writing the update equation for the $n$th entry of the vector $\mathbf{d}_{j}^{k}$, we have

[TABLE]

For each $n$, we have a different coefficient $\eta_{n,i,j}$. Therefore, Eq. 28 means that each component of $\mathbf{d}_{j}^{k+1}$ may be computed with a different linear combination of the $N_{d,\textrm{h}}$ columns of $\mathbf{K}^{k}$. Thus, localization expands the degrees of freedom to assimilate data.

The application of DSI-ESMDA is similar to the standard DSI, both methods require to run reservoir simulations only for the prior ensemble of models to generate the prior ensemble of predicted data. After that, DSI-ESMDA applies Eq. 18 $N_{a}$ times to generate the posterior ensemble. Note that the $N_{a}$ iterations are necessary because the prior ensemble is not Gaussian (if the prior is Gaussian, MDA is equivalent to a single ES update). Moreover, note that using Gaussian anamorphosis to transform variables only ensures that the marginal distributions of individual $d_{i}$’s are Gaussian, the joint distribution which is updated with DSI-ESMDA may not be Gaussian. In our tests, we noticed that only a few MDA iterations are required. In all cases presented in the next section, we use $N_{a}=4$. DSI-ESMDA can also be applied to the same parameterization of DSI, that is, update the PCA coefficients with the data transformation of Eq. 16. However, our initial tests showed that this procedure did improve the results, actually the results were slightly worse. Moreover, using the PCA parameterization of DSI would prevent to use localization. The only correction we applied is to truncate in zero if the final predicted data is negative, which in our tests occurred with water production data before breakthrough. Note that we apply this truncation only to the final estimates, not during the DSI-ESMDA iterations.

3 Test Cases

3.1 Test Case 1

The first test problem is a synthetic reservoir case created with rock and fluid properties typically found in the Campos Basin, Brazil. The reservoir model contains 50 $\times$ 70 $\times$ 10 gridblocks with uniform size of 50 m $\times$ 50 m $\times$ 5 m. A reference (true) case was generated using sequential Gaussian simulation for modeling porosity and sequential Gaussian co-simulation for permeability using the porosity as secondary variable and a correlation coefficient of 0.95. Table 1 summarizes the geostatistical parameters used to generate the reference model. The model has two vertical oil producing and two vertical water injection wells placed on the borders of the reservoir as illustrated in Fig. 2. The producers operated under as specified bottom-hole pressure (BHP) of 25,000 kPa while the injectors operate with a BHP of 35,000 kPa.

Using the same geostatistical parameters of the reference model, we generated an ensemble of 500 realizations for testing DSI and DSI-ESMDA. The synthetic measurements correspond to oil and water rate at both production wells and water injection rate at both injectors. These observations were corrupted by adding random Gaussian noise with zero mean and standard deviation corresponding to 10% of the data predicted by the reference model. For DSI, we kept 99% of the singular value energy (Eq. 15) which corresponded to 98 singular values. For DSI-ESMDA, we also kept 99% of the singular value energy in the subspace inversion. We used $N_{a}=4$ MDA iterations. Neither spatial nor temporal localizations were applied in this case. For comparisons, we also applied a standard history matching (that is, model-space inversion) using ES-MDA to update the same prior ensemble of 500 realizations with $N_{a}=4$ without localization. It is worth noting that the results of the history matching with ES-MDA do not correspond to the reference solution, as the method is not guaranteed to sample the posterior PDF correctly for nonlinear problems. The objective is to compare the DSI results against a history matching procedure that is used in practice. There are rigorous sampling methods, such as Markov chain Monte Carlo and rejection sampling [37] that could be used to generate the reference solution. Unfortunately, these methods are computationally prohibitive, even for a simple problem such as the one describe in this section. One alternative would be to simplify the problem by reducing the size of the model and the number of data points. However, this has already been done by [48], where it is shown that DSI obtained reasonable results compared to rejection sampling for a problem with few data points.

Figure 3 shows the predicted water production rate for both wells of the model obtained with DSI, DSI-ESMDA. For comparisons, we included in each plot the water rate predicted by the history-matched models with ES-MDA. The results in this figure indicate that both DSI methods obtained similar predictions, which are in reasonable agreement with the standard ES-MDA. Figure 4 shows the field cumulative production of oil and water predicted by the prior ensemble and by each method. The cumulative production from the reference case is also included in this figure for comparisons. This figure shows that the uncertainty ranges predicted by the methods are similar. Figure 5 shows the predicted water production rate for the first well considering three different sizes of the historical period. The main difference is observed for the case with no water breakthrough in the production history (cases with 20 data points in Fig. 3). For this case, DSI-ESMDA predicts a larger uncertainty range in the water production rate. In order to test the robustness of the methods, we selected another reference model with predictions outside the P10–P90 range from the prior ensemble and the results are presented in Fig. 6. Despite of the poor coverage of the prior ensemble, both methods obtained forecasts in reasonable agreement with the reference. However, we note that DSI-ESMDA obtained a sightly better data match for the first well (Figs. 6a and b) and a forecast range centered around the reference.

Both methods require to run reservoir simulations of the prior ensemble, which is clearly the dominating computational time. DSI also requires to solve the RML optimizations with the BFGS method, while DSI-ESMDA requires to apply the analysis (Eq. 18) $N_{a}=4$ times for each ensemble member. We measured the CPU running time of the inversion part of each method and the results are reported in Table 2. All computations were performed in the same computer (Intel Core i5-4690 CPU 3.5 GHz and 16 GB RAM). The results in this table show that DSI-ESMDA is notably faster than DSI, requiring only two seconds against approximately 25 minutes for DSI.

3.2 Test Case 2

The second test case is a more realistic history-matching problem known as UNISIN-I-H [3]. This problem is based on actual data from Namorado Field (Campos Basis, Brazil). The UNISIM-I-H model has 81 $\times$ 58 $\times$ 20 gridblocks, but only 37,000 are active. All gridblocks have a uniform size of 100 m $\times$ 100 m $\times$ 8 m. The original dataset is available for download at [32]. The dataset consists of 500 realizations of petrophysical properties (porosity, permeability in the three orthogonal directions and net-to-gross ratio). Besides the 500 petrophysical realizations, the UNISIM-I-H case also includes six global parameters, whose prior uncertainties were modeled as independent triangle distributions with values shown in Table 3. In this field there are 25 long horizontal wells (14 producers and 11 water injectors). Figure 7 shows the position of the wells projected in the first layer of the model. The oil producing wells are perforated close to the top and the water injection wells are perforated near to the bottom of the reservoir. The observed data were generated by adding random Gaussian data-error to the data predicted by a reference fine-scale model (UNISIM-I-R). The UNISIM-I-R was constructed with a higher level of geological details in a grid with 3.5 million active gridblocks [3]. The observed data correspond to monthly “measurements” of oil and water rate and the noise level was assumed as 10% of the data predicted by the model UNISIM-I-R. All wells are controlled by specified BHP during the historical and forecast periods.

For this test problem, we applied DSI, DSI-ESMDA and DSI-ESMDA with localization. We also applied a model-based inversion (history matching) with ES-MDA for comparisons. The dimension of the matrix $\Delta\mathbf{D}$ for the UNISIM-I-H case is $12825\times 500$ and the SVD retaining 99% of the singular values energy resulted in a vector of PCA coefficients with 395 elements for DSI. For DSI-ESMDA, we used $N_{a}=4$ data assimilations and kept 99% of the singular values in the subspace inversion. DSI-ESMDA with localization considered a critical length computed with $L_{x^{\prime}}=L_{y^{\prime}}=2000$ m and $T=6000$ days. These values were selected based on our previous experience testing ES-MDA for history matching this benchmark problem. We did not try different configurations for the critical length. The history matching with ES-MDA also used $N_{a}=4$ and a spatial localization with critical length of 2000 m. Note that in the history matching case, the Kalman gain is applied to update model parameters (porosity, permeability, etc.). Therefore, it does not seem appropriate to introduce the “time distance” to compute the localization matrix in this case.

In ours tests, the DSI implementation failed to achieve reasonable data matches for the UNISIM-I-H case. This result came as a surprise because the same problem was not observed with DSI-ESMDA using the same prior ensemble. We conducted a series of experiments trying to solve this problem, but the reason for the poor performance of the optimizations is still unclear. Among our experiments, we tried to use both numerical and the approximate analytical gradients during the minimizations, however, the large majority of the optimization did not converge properly. We also introduced a rescaling procedure before the truncated SVD. Basically, instead of computing the SVD of the matrix $\Delta\mathbf{D}$ (as in Eq. 13), we applied SVD to a scaled version

[TABLE]

and computed the square root of $\mathbf{C}_{\mathbf{d}}$ for PCA as

[TABLE]

The rationality for this procedure is to avoid losing relevant information during the truncation of small singular values because $\Delta\mathbf{D}$ may be poorly scaled. Unfortunately, this procedure did not result in significant improvements. We also tried the optimizations keeping all 499 singular values to preserve the maximum number of degrees freedom provided by the prior ensemble to match data. In this case, we observed an even worse performance. Conversely, we also tried to reduce the number of singular values by keeping only 95% and 90% of the energy, which resulted in 237 and 154 singular values, respectively. The idea was to check if using a small number of PCA coefficients would help to regularize the optimizations. Again, no noticeable improvements were observed.

In order to evaluate the quality of the final data matches, we computed the data-mismatch objective function normalized by the number of observed data points ($N_{d,\textrm{h}}=9900$) for 500 posterior realizations obtained with DSI, DSI-ESMDA, DSI-ESMDA with localization and the history matching with ES-MDA. The normalized data mismatch objective function was computed as

[TABLE]

where $\sigma_{e,i}$ is the data-error standard deviation of the $i$th datum. The last equality in (31) holds only for independent data errors. In the Appendix section of [35] it is shown that the expectation of $\mathcal{O}_{N,d}$ for a set of RML samples in a linear problem should be one half. Even though this value is valid only for linear problems, it serves as a reference. For example, if we have a predicted curve where the difference with all observed data points is exactly one standard deviation of the data error, then the normalized objective function of this curve is exactly 0.5. Analogously, we have $\mathcal{O}_{N,d}=2$ and $\mathcal{O}_{N,d}=4.5$ for two and three standard deviations, respectively. This effectively means that it would be difficult to justify in practice the acceptance of a data match with $\mathcal{O}_{N,d}$ significantly larger than five. Table 4 shows the mean and standard deviation of $\mathcal{O}_{N,d}$ for each case. Clearly, DSI resulted too large values for $\mathcal{O}_{N,d}$, indicating the lack of convergence of the minimizations. The objective functions from DSI-ESMDA and DSI-ESMDA with localization have mean roughly ten times smaller. The history matching case with ES-MDA also resulted in a mean objective function significantly smaller than DSI. In terms of computational cost, the DSI method required approximately 28 hours to complete, while the DSI-ESMDA (with and without localization) require around 45 minutes each.

The prior ensemble for the UNISIM-I-H case was provided by the authors of the benchmark and it is noteworthy that even though this is a synthetic problem, the predictions from the prior ensemble do not span the data from the reference case for several wells. This is clearly not an ideal situation for applying DSI or even any history-matching method. Ideally, the prior ensemble should provide a reasonable estimate of the prior uncertainty, at the very least, it should to be able to span the observations. In practice, we should revise the prior ensemble before using for data assimilation. However, here for the purposes of the investigation, we decided to test the methods with this deficient prior ensemble. Figures 8 and 9 show the results obtained by DSI, DSI-ESMDA and DSI-ESMDA with localization for three wells with good and poor coverage of the predictions from the prior ensemble, respectively. We selected these wells because they are representative of the results observed in this problem. Figure 8 indicates that DSI failed to match data even for the wells with good prior coverage. As a result, the prediction from the reference model lies outside the predicted P10–P90 range. For the wells with poor coverage (Fig. 9), DSI seems to fail to reduce the uncertainty range properly, for example, the posterior predictions for wells PROD021 and PROD025A have almost the same uncertainty range of the prior ones. DSI-ESMDA obtained reasonable data matches for all wells, however, the predictions are clearly too narrow as they do not span the reference. DSI-ESMDA with localization improved the results significantly, although it was not able to span the reference predictions, especially for the wells with poor prior coverage. Figure 10 compares the results of DSI-ESMDA with localization with the history matching with ES-MDA. It is interesting to note that the history matched models also suffer from the lack of representativeness of the prior ensemble; see, for example, the plots (a) and (d) in Fig. 10. Figure 11 shows boxplots of field cumulative oil and water production predicted by the prior ensemble and by each method. The cumulative production from the reference case is also included in this figure for comparisons. Overall, the predictions are biased compared to the reference production, which is probably explained by the problems with the prior ensemble.

4 Field Case

The field case corresponds to a large offshore turbidite reservoir in Campos Basis. Currently there are 24 wells producing oil in the field from a total of 43 vertical wells that have been drilled during 18 years of operation. The main recovery mechanism is pressure maintenance by a large aquifer. A total of 500 prior realizations of the reservoir were built by the asset team of the field using the best practices available to date in the Company. These realization include facies, porosity, permeability and net-to-gross ratio for each one of 11 producing zones. The engineers of the field also selected the oil-water relative permeability curves as parameters with relevant uncertainty. The relative permeability curves in six different regions of the model were parameterized using Corey functions with the exponents and the maximum relative permeability of water selected as uncertainty parameters with Gaussian prior distributions. The model has approximately 500,000 active gridblocks (total of 114 $\times$ 238 $\times$ 60 with an average size of 75 m $\times$ 75 m $\times$ 2 m). This field has a large number of sealing faults as illustrated in Fig. 12. All wells are controlled by specified total liquid production rate during the historical period and the forecast. The observed data corresponds to monthly measurements of oil and water rate at the well with a total of 26703 data points.

Unlike the UNISIM-I-H case, the predictions from the prior ensemble cover reasonably well the observed data for the large majority of the wells in the field, which creates a more favorable situation for the application of the methods. Once again, the optimizations in our DSI implementation failed to converge, providing unreasonable results for this field problem. DSI-ESMDA without localization showed severe reduction in the posterior variance, showing indications of ensemble colapse. Both results were not included in this paper. The only acceptable results obtained in this field are the ones from DSI-ESMDA with localization. In this case, we use the same configuration of the case UNISIM-I-H to define the critical length for localization, that is, $L_{x^{\prime}}=L_{y^{\prime}}=2000$ m and $T=6000$ days. For comparisons, we history matched the same ensemble of models using ES-MDA with $N_{a}=4$ and localization with a critical length of 2000 m. Table 5 presents the mean and standard deviation of the normalized data mismatch objective function. The results in this table indicate that DSI-ESMDA with localization obtained an excellent data match, with mean objective function of 0.55. The history matching with ES-MDA also improved significantly the data match, but with higher values of objective function. Figure 13 shows the water production rate for eight representative wells of the field. The values in this figure were normalized to preserve the confidentiality of the information. Overall, the forecasts after DSI-ESMDA with localization are comparable with the forecasts from the history-matched models, despite some clear differences observed in some wells. Figure 14 shows boxplots of field cumulative oil and water production. The main difference between DSI-ESMDA and ES-MDA is for the cumulative water production, where DSI-ESMDA predicts a higher water production in the field.

5 Conclusions

This paper introduced a new DSI implementation based on the data assimilation method ES-MDA. The new DSI-ESMDA was compared with the original DSI proposed in [48, 49], which is based on PCA to reparameterize the predicted data from a prior ensemble combined to RML for sampling. The new implementation preserves the main advantage of the original DSI, namely, it is able to provide an ensemble of production forecasts requiring reservoir simulations only for a prior ensemble of models. We compared the DSI-ESMDA with the original DSI in two synthetic reservoir problems. We also applied DSI-ESMDA to a real field case with long production history and large number of wells. Based on the results for these test problems, we can summarize the following conclusions:

• •

The proposed DSI-ESMDA is computationally faster than the original DSI. Even though the time to execute the reservoir simulations for the prior models tends to be dominant in both methods, the difference in the computational time for the inversion can be relevant for large problems. For example, for the UNISIM-I-H case, the original DSI required approximately 28 hours in a stand-alone computer while the DSI-ESMDA required only 45 minutes.

• •

The performance of both DSI implementations is highly dependent of the ability of the prior ensemble to provide reasonable estimates of the prior uncertainty. The same is also true for the more traditional model inversion with ES-MDA.

• •

The results indicate that the proposed DSI-ESMDA is more robust than the original DSI. In our tests, the optimizations required by DSI have not converged properly for the UNISIM-I-H case. As a result, DSI failed to obtained acceptable data matches and reasonable production forecasts. DSI-ESMDA, on the other hand, was able to match the observed data for all wells.

• •

The use of spatial and temporal localization improved significantly the results of the DSI-ESMDA method when the number of data point is large.

• •

The proposed DSI-ESMDA method with localization obtained forecasts of production comparable with the forecast provided by the posterior models generated with history matching using ES-MDA for a real field problem.

• •

The proposed method is very simple to implement. It does not require data transformations and it is straightforward integrate with different types of data and models.

Acknowledgement

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