About#

A simple 2D reservoir simulator and a straight-forward implementation of the basic Ensemble Smoother with Multiple Data Assimilation (ESMDA) algorithm.

ESMDA#

In the following an introduction to the Ensemble Smoother with Multiple Data Assimilation (ESMDA) algorithm following [EmRe13]:

In history-matching problems, it is common to consider solely the parameter-estimation problem and thereby neglecting model uncertainties. Thus, unlike with the ensemble Kalman filter (EnKF), the parameters and states are always consistent (Thulin et al., 2007). This fact helps to explain the better data matches obtained by ESMDA compared to EnKF. The analyzed vector of model parameters \(z^a\) is given in that case by

\[z_j^a = z_j^f + C_\text{MD}^f \left(C_\text{DD}^f + \alpha C_\text{D} \right)^{-1}\left(d_{\text{uc},j} - d_j^f \right) \qquad \text{(1)}\]

for ensembles \(j=1, 2, \dots, N_e\). Here,

  • \(^a\): analysis;

  • \(^f\): forecast;

  • \(z^f\): prior vector of model parameters (\(N_m\));

  • \(C_\text{MD}^f\): cross-covariance matrix between \(z^f\) and \(d^f\) (\(N_m \times N_d\));

  • \(C_\text{DD}^f\): auto-covariance matrix of predicted data (\(N_d \times N_d\));

  • \(C_\text{D}\): covariance matrix of observed data measurement errors (\(N_d \times N_d\));

  • \(\alpha\): ESMDA coefficient;

  • \(d_\text{uc}\) : vector of perturbed data, obtained from the vector of observed data, \(d_\text{obs}\) (\(N_d\));

  • \(d^f\): vector of predicted data (\(N_d\)).

The prior vector of model parameters, \(z^f_j\), can in reality be \(j\) possible models \(z^f\) given from an analyst (e.g., the geologist). In theoretical tests, these are usually created by perturbing the prior \(z^f\) by, e.g., adding random Gaussian noise.

The ESMDA algorithm follows [EmRe13]:

  1. Choose the number of data assimilations, \(N_a\), and the coefficients \(\alpha_i\) for \(i = 1, \dots, N_a\).

  2. For \(i = 1\) to \(N_a\):

    1. Run the ensemble from time zero.

    2. For each ensemble member, perturb the observation vector using \(d_\text{uc} = d_\text{obs} + \sqrt{\alpha_i} C_\text{D}^{1/2} z_d\), where \(z_d \sim \mathcal{N}(0,I_{N_d})\).

    3. Update the ensemble using Eq. (1) with \(\alpha_i\).

The difficulty is the inversion of the large (\(N_d \times N_d\)) matrix \(C=C_\text{DD}^f + \alpha C_\text{D}\), which is often poorly conditioned and poorly scaled. How to compute this inverse is one of the main differences between different ESMDA implementations.

Also note that in the ESMDA algorithm, every time we repeat the data assimilation, we re-sample the vector of perturbed observations, i.e., we recompute \(d_\text{uc} \sim \mathcal{N}(d_\text{obs}, \alpha_i C_\text{D})\). This procedure tends to reduce sampling problems caused by matching outliers that may be generated when perturbing the observations.

One potential difficultly with the proposed MDA procedure is that \(N_a\) and the coefficients \(\alpha_i\)’s need to be selected prior to the data assimilation. The simplest choice for \(\alpha\) is \(\alpha_i = N_a\) for all \(i\). However, intuitively we expect that choosing \(\alpha_i\) in a decreasing order can improve the performance of the method. In this case, we start assimilating data with a large value of \(\alpha\), which reduces the magnitude of the initial updates; then, we gradually decrease \(\alpha\).

Reservoir Model#

The implemented small 2D Reservoir Simulator was created by following the course material of AESM304A - Flow and Simulation of Subsurface processes at Delft University of Technology (TUD); this particular part was taught by Dr. D.V. Voskov, https://orcid.org/0000-0002-5399-1755.