Formulating the history matching problem with consistent error statistics

It is common to formulate the history-matching problem using Bayes’ theorem. From Bayes’, the conditional probability density function (pdf) of the uncertain model parameters is proportional to the prior pdf of the model parameters, multiplied by the likelihood of the measurements. The static model parameters are random variables characterizing the reservoir model while the observations include, e.g., historical rates of oil, gas, and water produced from the wells. The reservoir prediction model is assumed perfect, and there are no errors besides those in the static parameters. However, this formulation is flawed. The historical rate data only approximately represent the real production of the reservoir and contain errors. History-matching methods usually take these errors into account in the conditioning but neglect them when forcing the simulation model by the observed rates during the historical integration. Thus, the model prediction depends on some of the same data used in the conditioning. The paper presents a formulation of Bayes’ theorem that considers the data dependency of the simulation model. In the new formulation, one must update both the poorly known model parameters and the rate-data errors. The result is an improved posterior ensemble of prediction models that better cover the observations with more substantial and realistic uncertainty. The implementation accounts correctly for correlated measurement errors and demonstrates the critical role of these correlations in reducing the update’s magnitude. The paper also shows the consistency of the subspace inversion scheme by Evensen (Ocean Dyn. 54, 539–560 2004) in the case with correlated measurement errors and demonstrates its accuracy when using a “larger” ensemble of perturbations to represent the measurement error covariance matrix.


Introduction
Reservoir simulators are typically forced by measured oil, gas, and water rates during the historical simulation to generate the predicted production [1]. This forcing can be applied directly using individual measured rates (e.g., the oil rate) or maybe more commonly by the produced reservoir volume as computed from the fluid rate data (i.e., the total volume of fluids produced by a well). Thus, one attempts to prescribe the observed production and injection from each well to the model, and if it is physically possible, the model will produce the prescribed rates. Sometimes the model Geir Evensen geir.evensen@norceresearch.no 1 NORCE (Norwegian Research Center), Bergen, Norway 2 Nansen Environmental and Remote Sensing Center, Bergen, Norway is unable to satisfy the imposed production constraints. In such cases, the purpose of the conditioning is to update the reservoir description so that the model can produce the observed rates.
The approach described above is standard in most history-matching studies based on ensemble methods, see e.g., [4,6,21,25]. Additionally, the inclusion of seismic information complements the information contained in the historical rates [5,19]. Lately, the conditioning on seismicinformation has become more efficient with adaptive localization schemes, such as the one proposed by Luo and Bhakta [20]. Petroleum companies now use ensemble conditioning methods operationally to manage their reservoirmodel workflows [13,23]. This paper aims to enhance the understanding of the formulation and solution of ensemble history-matching methods.
We will now formulate the history matching problem while allowing the historical rates to contain errors. Start by storing all the observed production and injection rates of oil, gas, and water, for the whole history of the model, in the control vector, u ∈ m u , and define the x ∈ n , to include all the uncertain parameters that characterize the reservoir. We can then write a model for the simulated rates of produced oil, water, and gas, i.e., the predicted measurements, y ∈ m , as y = g x, u . ( Thus, the model depends on the uncertain parameters in x ∈ n and the uncertain rates in the vector of controls, u. Here m u is the dimension of the control variable, u, m is the dimension of the possible subset of rate data that we condition on, and n is the dimension of the vector of uncertain model parameters. Now, characterize the random variables x and u by a prior joint probability density function (pdf), f (x, u). The model in Eq. 1 is then a stochastic equation starting from an uncertain set of parameters, x, and forced by the uncertain controls, u. A joint pdf for the model prediction and the poorly known inputs is f (y, x, u) = f x, u f y − g x, u = f x, u δ y − g x, u , with δ being the Dirac-delta function. One traditionally neglects the stochastic noise in the measured rates, which results in an underestimated prediction uncertainty for y. The impact of this noise can be significant, particularly in the case with time-correlated errors in the forcing u, which yields significantly more substantial error-variance growth. In the current study, we neglect other model errors, such as those that may arise from numerical discretization. However, it is possible to include these additional errors using the algorithm from Evensen [10]. Define the relation between the measured rate data used in the conditioning, d ∈ m , and the actual production, y, by where e ∈ m denote the errors in the measured rates, d.
We differentiate between the data conditioned on, d, and the controls used to force the model, u, (although d may typically include a subset of u). Now, write the following joint pdf: When deriving the ensemble methods, it is not required to include the normalizing denominator, and one can instead write (5) as f y, x, u | d ∝ f d | y, x, u f y, x, u .
The prior-term defined in Eq. 2 is straightforward to evaluate. The measurements, d, only depends on y, and the likelihood becomes f d | x, y, u = f d | y , as usual. Still, it is necessary to separate the control variables used to force the model and the data on which we condition the posterior. Then, using (2), it is possible to write (6) as The marginal is obtained by integration over y and becomes and, together with Eq. 1 for the model prediction, constitutes the reformulated history-matching problem. The formulation in Eq. 8 relates to the one discussed by Evensen [10], who formulated the history-matching problem when including model errors. Here the "model errors" correspond to the errors in the rate data. As in Evensen [10], we end up having to estimate both the parameters in x and the rate data in u.
Note that when using ensemble methods, the parameterspace dimension is not of crucial importance. The ensemble methods compute the sensitivities between measurements and parameters, and parameters are updated only if the sensitivity is unequal to zero. Unimportant parameters will retain their prior variance. Thus, the addition of uncertain controls to the estimation does not increase the complexity of the problem.

Solution method with normal priors
Define the new random variable z T = x T , u T ∈ n+m u , (9) and rewrite model as y = g x, u = g(z). (10) It is then possible to write the problem stated by Eq. 8 as Now, assume normal priors where and Eq. 11 becomes Thus, Eq. 15 has the same form as in the traditional historymatching problem. The difference is that we solve for both the model parameters and the rate errors in z rather than only the model parameters, x.

Cost funtion and ensemble formulation
Maximizing f (z | d) is equivalent to minimizing the cost function We can approximately sample the posterior density by minimizing an ensemble of cost functions [17,22], Thus, we have an ensemble of models forced by uncertain controls, representing the prior information, and we have perturbed measurements, representing the information from the likelihood.

Gradient and Hessian
The gradient of the cost function in Eq. 17 becomes and the approximate Hessian is The above equations are similar to the ones obtained when deriving the Iterative Ensemble Smoother with model errors included [10, compare (20) and (21) with Eqs. 71 and 72 in].

Introduce ensemble averaged model sensitivity
Next, define a joint or averaged model sensitivity by the linear regression which represents the best least-squares-fit to the ensemble of model sensitivities One can then approximate the gradient and the Hessian from Eqs. 20 and 21 with Algorithm 1 Subspace EnRML algorithm. For the case with estimation and use of correlated rate errors one need to include the lines three and sixteen, and also the E i in the model in lines 9 and 12, all marked in red. : 12: 15: By replacing the model sensitivities, G j , with the average least-squares fit, G, we introduce an approximation in the gradient in Eq. 20. We, therefore, also change the location of the minimum of the cost function. With the model adjoint's availability, it would be possible to compute the correct model sensitivity and thereby eliminate this approximation. The approximation introduced by using the averaged model sensitivity comes on top of the approximate use of the ensemble of cost functions in Eq. 17 to sample the posterior.

Subspace EnRML solution
For computing the Ensemble Smoother (ES) or Iterative Ensemble Smoother (IES) solutions, see [3,9,10,12,24], where the two last papers present new subspace-EnRML algorithm used in this work (see Algorithm 1). The method computes for each ensemble member z j , the following iteration for iterate number i, which uses the definitions of the gradient and Hessian from Eqs. 24 and 25. With the step length, γ = 1.0, and only one iteration, the solution obtained is the ES solution.
Otherwise, one needs to specify a value of γ ∈ (0, 1], typically γ = 0.5 is a reasonable choice, and iterate until convergence. The subspace Algorithm 1 takes as input the ensemble of model parameters, X = (x 0 , . . . , x N ) ∈ n×N , and the ensemble of perturbed measurements, D ∈ m×N . If we include rate errors in the estimation, we must also input the ensemble of rate perturbations E 0 ∈ m u ×N (red line number 3).
From the input ensemble one can compute the ensemble of model predictions, g(X) ∈ m×N . Furthermore, the algorithm uses the scaled matrices of model predictions anomalies, Y = g(X) , and scaled measurement perturbations, E = D . The projection , as defined in line number 5 of Algorithm 1, subtracts the ensemble mean and then scales the resulting ensemble anomalies by √ N − 1. In line number 11 of Algorithm 1, S i is most easily computed by writing Ω T S T i = Y T and then computing the LU factorization of Ω T i followed by m back substitutions to a cost of O(mN 2 ) floating point operations. In line number 12, one needs to compute the matrix multiplication S i W i to a cost of mN 2 operations.
It is possible to approximately compute the inversion in line number 13 to a cost of O(mN 2 ) operations by using the ensemble subspace algorithm discussed in [7,8,12] (and see also the discussion in the following section).
Finally, in line number 15, we obtain the updated model ensemble in nN 2 floating-point operations. In the case of uncertain rate errors, we must also update E i in the red line number 16 to a cost of m u N 2 , and use E i to force the model. Thus, the most expensive computation in Algorithm 1 is likely the evaluation of the model-ensemble prediction Y i = g(X i ) in each iteration. The cost of integrating the model is proportional to the number of time steps, n t . With each time step requiring n c n computations, the total cost of the ensemble integration is Nn t n c n. For nonlinear models, the number of timesteps is likely much larger than N, and the number of computations per time step can be several times n. All other computations in this algorithm are linear in the number of measurements, m, and the number of state variables, n. For the detailed derivation of the subspace EnRML algorithm, see [12,24].
Thus, in the case where we include rate errors with uncertainty, we include the lines in red in Algorithm 1. The only required modification is that we need an unscaled version of the matrix of measurement perturbations stored in E 0 . Then the algorithm will use E 0 to compute updated estimates of the rate errors, E i , in each iteration, and use these updated rates to force the model.

Ensemble subspace inversion
In Algorithm 1, we have represented the measurement errorcovariance matrix by an ensemble of measurement perturbations, E. The inversion is in this case computed using ensemble subspace scheme as was proposed by Evensen [7] and further discussed in Evensen [8] and recently in Evensen [12]. This scheme projects the measurement error perturbations onto the ensemble subspace and computes the pseudo inverse of the following factorization where we have defined the singular-value decomposition and the I N ∈ N×N is just the identity matrix. The eigenvalue decomposition in Eq. 30 is of the matrix product in Eq. 29. Note that this eigenvalue decomposition is most efficiently computed by a singular value decomposition of the product + U T E. The left singular vectors will then equal the eigenvectors in Z and the squares of the singular values will equal the eigenvalues in . Thus, the inversion becomes The subspace inversion algorithm's properties will be discussed in more detail in Section 3.

EnKF update example
Before testing the proposed algorithm on a reservoir model, it is instructive to examine a single update step using a toy model. The purpose of this example is to illustrate the properties of the update scheme, in particular, when one uses the measurement perturbations to represent the measurement error covariance matrix. This example verifies the robustness of the projection of the measurement error covariance matrix onto the ensemble of predicted measurements. The update is identical to the EnRML algorithm's solution in the linear case, and as such, the learnings are representative for the iterative smoother update as well.

Example description
The test example uses a one-dimensional periodic domain with 1024 grid points and x = 1. In this domain, we simulate a smooth pseudo-random function with mean μ = 4, variance σ 2 = 1, and decorrelation length r d = 40, representing the unknown truth, The constant μ = 4, is just added for plotting purposes.
The first guess solution is generated by simulating another realization x ← N (0, 1, 40) and adding it to the truth, i.e., The factor √ 2 ensures that the variance of x fg is equal to one.
The initial ensemble is created by adding random realizations x i ← N (0, 1, 40) to the first guess x fg , The measurements are distributed uniformly over the domain and sampled from a perturbed true solution according to where either uncorrelated, e i ← N (0, 0.25, 0), or Gaussian, e i ← N (0, 0.25, r d ), errors are used. Here, M is the linear measurement operator that extracts the measurements from the functions x true + e i .
The following experiments use the same referencetruth, measurements, initial ensemble, and random seed. When it is essential to eliminate sampling errors, we use extended ensemble sizes. Thus, although there are seed dependencies in the obtained solutions, the different methods should produce the same answer. We can, in most of the experiments, attribute differences in the solutions to the methods used. Thus, our approach is different from running multiple data assimilation experiments with varying seeds while attempting to estimate the correct answer.

Solution methods
In the following, we study two prominent cases, one with a diagonal measurement error-covariance matrix and another with correlated errors, e i , simulated from Eq. 37 and with r d = 40. For the two cases of uncorrelated and correlated measurement errors, the EnKF computes the analysis using both an exactly specified measurement error-covariance matrix C dd and while representing the measurement error covariance by the perturbations in E.
The case with a full-rank measurement error-covariance matrix solves where the ensemble perturbation matrix is A = X . The matrix C = SS T + C dd is formed and then inverted by computing an eigenvalue decomposition C = Z Z T . The inverse is just C = Z + Z T where the use of a pseudo inverse is needed in case the matrix C is poorly conditioned.
When using an ensemble representation for the measurement error-covariance matrix C dd = EE T , one can solve for the update from In the examples below, the line labels used in the figures indicate the inversion scheme used to compute (SS T + C dd ) −1 . The line label Cdd denotes using the standard EnKF analysis equation with a full rank measurement errorcovariance matrix C dd , as explained above. The curves with line label EE correspond to the EnKF update when the samples in E replace the "exact" analytic measurement error covariance matrix C dd , and we use the ensemble subspace scheme.
Using E to represent the measurement error-covariance matrix introduces additional sampling errors. However, we will see below how it is possible to reduce these sampling errors to a negligible level with a simple algorithm modification. I.e., one uses a larger number of realizations in E to represent C dd better. The code used is the test case from https://github.com/geirev/EnKF Analysis.git.

Example 1 (large ensemble size)
The first example uses 50 measurements and a large ensemble size of 2000 to reduce sampling errors. Figure 1 shows the results for the two cases with either diagonal or correlated measurement errors.
The upper-left plot shows the EnKF estimates for the case with uncorrelated measurement errors. The two schemes, represented by the lines labeled Cdd and EE, give similar results in this case. The upper-right plot shows the prior and posterior error variances for the two updates, and again they are nearly identical. The lower-left plot presents the EnKF estimates for the case with correlated measurement errors. In this case, we also see that the results using the exact and approximate schemes (Cdd and EE) are nearly identical.
An apparent difference between the two cases is that, with uncorrelated errors, the measurements are scattered randomly about the correct solution. In contrast, with correlated measurement errors, successive measurements will have similar error values, and they follow a smooth curve.
The nonzero measurement correlations' role is to reduce the strength of the update, and the result is an update with a more substantial variance. By taking the measurement error correlations into account, we inform EnKF that neighboring measurements make the same error, and we reduce their accumulated impact. Table 1 shows that for this case with a large ensemble, the approximation we introduce by using EE T in the subspace inversion scheme, instead of inverting the full C dd , is negligible. Thus, the learning from these two cases is that correlated measurement errors reduce the measurements' impact on the update. Furthermore, the two schemes used to compute the inversion are mutually consistent and give approximately the same result with only minor sampling errors. At least, this result is accurate when the scales of the measurement perturbations are the same as in the predicted ensemble anomalies.
The lower plots in Fig. 1 include two inconsistent updates conditioned on the same measurements with correlated errors. In the first inconsistent update (named ICA), C dd = I, thus neglecting the measurement-error correlations in C dd , and uncorrelated measurement error perturbations are sampled and used to form D. This update is inconsistent, but it is the procedure that most practical EnKF applications use. I.e., one neglects the measurement-error correlations, both because they are challenging to prescribe but also because it allows for efficient computation of the update using (41) discussed below. We notice that the ICA estimate slightly overfits the measurements. The reason is that by neglecting the measurement error correlations in C dd , EnKF gives "full weight" to all the measurements. Additionally, the ensemble variance becomes way too low and similar to the case with uncorrelated measurement errors. Thus, neglecting measurement error correlations leads to a too close fit to the observations and a too low variance estimate.
In the second inconsistent update (named ICB), we still set C dd = I, but we retain the original correlated perturbations when forming D. The surprising result is that we obtain a too large variance in this case. One can explain this larger variance: A set of measurements of a smooth "true" field, with correlated errors, will have similar values. By creating the measurement perturbation with correlated errors, one generates perturbed measurements of similar value for each realization. Then, using a diagonal C dd , each realization is pulled too strongly towards the common value of the perturbed measurements. Thus, the updated ensemble members will more closely resemble the variance of the measurement perturbations.
In Table 1, we see that cases ICA and ICB lead to significant errors, both for the estimated mean and the variance. Here N is the ensemble size, m is the number of measurements, n e dimension inflation of the E-matrix, and r d is the decorrelation length used for the measurement perturbations. The RMSE is the root mean square difference between the solutions (ensemble means and variances) computed using an exact inversion of C dd and an approximation EE T These two inconsistent updates stress the importance of specifying the correct measurement error statistics in the ensemble update scheme.

Example 2 (ensemble size of 100)
We now repeat the previous experiment from Example 1 using a more common ensemble size of 100. The purpose is to illustrate the impact of sampling errors when using the measurement error perturbations in E to represent C dd . In Fig. 2, we observe that with 100 realizations, the additional sampling errors introduced by scheme EE lead to a slight deviation between the two estimates. More problematic is the underestimation of the ensemble variance. In a sequential data-assimilation context, this underestimation would have to be compensated for, e.g., by using inflation, to avoid possible filter divergence. In the next example, we will learn how to reduce these sampling errors to a negligible level.

Example 3 (augmenting realizations to E)
The benefit of using (39) over (38) is the reduced computational cost, but also the fact that it is easier to sample perturbations with accurate statistics than constructing a full-rank measurement error covariance matrix. An approach for reducing the sampling errors in scheme EE is to augment columns of new realizations of measurement perturbations to E. This modification only slightly increases the computational cost of the algorithm when computing + U T E in Eq. 29 and is a simple modification of the code. In Fig. 3, we see the results using 100 realizations and correlated measurement errors, and when using 1000 samples in E.
From Table 1

Example 4 (different scales in E)
Example 4 repeats the experiment from Example 3 using measurements with two different error decorrelation lengths r d = 20 and r d = 80. These length scales are in contrast to the length scale of r d = 40 used in the previous experiments. From the results shown in Fig. 4, we confirm that the subspace inversion algorithm works well also when the measurement error correlations have a different length scale than the model ensemble (i.e., r d = 40). The result from the case with smooth measurement-error correlations, r d = 80, is slightly more accurate than the r d = 20 result. If the measurement perturbations in E are smoother than the measured ensemble anomalies in S, then E should be well represented when projected on the space spanned by S.
However, in the case with r d = 20, the measurement perturbations include small scales that are not represented by the predicted measurements' ensemble. Thus, when conditioning on measurements with small-scale noise in the measurements that cannot be well represented by the ensemble anomalies in S, the subspace projection introduces an approximation. The truncation of small scales in the  Table 1 support this conclusion, where it is clear that the residual of the estimated mean increases for the case r d = 20 and decreases for the case r d = 80 compared to the r d = 40 case.

Example 5 (large number of measurements)
The final example, shown in Fig. 5, increases the number of measurements from 50 to 200, i.e., twice the ensemble size. In this case, we apply a truncation at 99% of the variance when computing the inversion, which retained 29 singular values when computing the singular value decomposition of S. Again the results obtained are very similar using the two algorithms. It is also interesting to see how the measurements' impact reduces at the grid indices 400-500 in this plot. Note also that there is no indication of the so-called "ensemble collapse," and the analysis ensemble retains a significant variance. The posterior variance using 200 measurements is similar to the one obtained using only 50 observations. This result indicates that including additional measurements does not introduce much new information in this example.

A note on the EnKF analysis equation
Most operational ensemble-based schemes apply an assumption of uncorrelated measurement errors and uses a diagonal C dd = I, see, e.g., the reviews on data assimilation in the geosciences [2], weather prediction [15], and petroleum applications [1]. This assumption is emploied for simplicity for two reasons. First, the measurement error covariances are often not well known, and additionally, the update scheme (38) simplifies considerably. With C dd = I, Eq. 38 becomes which makes it possible to use an efficient algorithm from Hunt et al. [16] where, by using a Woodbury identity, the EnKF update becomes dd being a symmetrical square root of a full rank C dd . E.g., write the eigenvalue decomposition and define the symmetrical square root and its' inverse Now, by scaling the predicted measurement anomalies and the innovations according to: and with some algebra, Eq. 38 becomes and using the Woodbury identity, There are, however, significant numerical costs associated with the establishment of C 1 2 dd and the associated rescalings in Eqs. 45 and 46 which are both O(m 2 N) operations.
Additionally, C 1 2 dd needs to be of full rank or a formulation based on pseudo inverses must be employed. Thus, the discussion in this section justifies the use of the ensemble subspace projection scheme in Eq. 33 for computing consistent updates at the cost of O(mN 2 ), while taking measurement error-correlations into account.

Reservoir experiments
We will now discuss the application of the subspace EnRML method in Algorithm 1 for history matching a reservoir cells on a 40-times-64 grid with 14 layers. The uncertain parameters include the three-dimensional porosity field and six fault multipliers F2-F7. There are five producing wells, OP1-OP5 and three injectors I1-I3. The model is the same as was used in Evensen et al. [12], and we will not describe the model at length here but rather focus on the impact of introducing correlated measurement errors for a reservoir case. The model was previously used also by Leeuwenburgh et al. [18],Hanea et al. [14], and Wang and Oliver [26], and we refer to these papers for a detailed presentation.

Overview of experiments
In the cases below, we condition on the monthly production rates. Thus, we have up to 36 triplets of the oil production rate (OPR), the gas production rate (GPR), and the water production rate (WPR) for each well (depending on the date the well is started). The error standard deviations for the rate measurements are five percent of the measured value of the OPR and 10 percent for the GPR and WPR. We didn't condition the model on the injection rates, and we didn't have access to pressure data. The measurements are real historical production rates from the reservoir.
We will present results from the experiments listed in Table 2. The tree first "PRED" experiments are ensemble predictions without conditioning the model on observations. We randomly picked one model realization from the ensemble and used this model's porosity and fault multipliers for all the ensemble members. The difference between the three ensemble-prediction experiments is the stochastic controls used to force the ensemble integration (white, red, or "bias"). In the experiment PREDW, the model forcing is the historical rate data contaminated by white noise, i.e., the measurement errors are uncorrelated in time. In PREDR, we set the time-decorrelation length of the measurement errors to 15 months, resulting in measurement errors that are smooth (red) in time. Finally, in PREDB, the measurement error perturbations are constant in time. Thus, the measurement errors are infinitely correlated, and for each realization, the errors act as a bias. These experiments' purpose is to assess the random forcing's impact on the prediction uncertainty while eliminating any contributions from model-parameter uncertainties. Case IES0 is the default or standard case using a diagonal measurement error covariance matrix C dd = I. Thus, there are no measurement error correlations and no stochastic noise added to the historical rate data used to force the reservoir model. Accordingly, there is no updating of the measurement perturbations. This case represents the traditional approach for solving the historymatching problem using ensemble methods. We computed the solution using Algorithm 1, with the red lines excluded.
The second case, IESnd, includes the effect of a nondiagonal C dd represented by the product of correlated measurement perturbations EE T . Like in IES0, there is no stochastic forcing of the model. This case allows us to examine the impact of accounting for the measurement error correlations in the analysis update.
The third case, IESR, is the full new Algorithm 1, including the red lines, using correlated measurement perturbations in E to perturb the measurements conditioned on, as well as taking measurement error correlations into account in the inversion. The algorithm also updates the ensemble of perturbations in E i and uses these to force the model realizations. The initial forcing ensemble is the red one used in PREDR.
The last case, IESB, is similar to IESR, except that it uses the "bias" initial forcing ensemble used in PREDB.
The ensemble size is 100 for all experiments. Like in Section 3, the experiments also use the same initial ensemble and random seed to make it easier to compare cases. Figure 6 presents the resulting oil production for the OP2 well when forcing the model ensemble with a white, red, and a "biased" noise model. The uncertainty in the produced rates will depend on the spread in the historical rates used to force the model realizations. Thus, the prediction uncertainty is more apparent in the first half of the simulation than in the second half since the forcing-uncertainty scales with the data value. For PREDW with white-noise forcing, we notice the instant response in the model ensemble's predicted rates. Also, we see that the uncertainty in the accumulated production remains low throughout the integration. This "random-walk" process, where the rate uncertainties nearly cancel out over time, yields a nearlinear variance increase in time (there is a nonlinearity contribution introduced by the model dynamics). In PREDR and PREDB, we notice an increasing spread in the predicted ensembles' accumulated production. This uncertainty increase is most significant in PREDB. We can explain it by the constant-in-time errors added to the historical rates, which will continuously drive the different realizations towards different productions. In fact, for the PREDB case, the uncertainty in the total accumulated production for OP2 is almost of the same magnitude as the uncertainty of The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments the prior ensemble in the IES0 case (which has no stochastic forcing and only reflects the reservoir uncertainty). The inclusion of stochastic historical production rates is essential since it increases the predicted rates' variance, thereby enhancing the impact of observations on the update. Figure 7 presents the prior and posterior ensembles of scalar fault multipliers for four of the six faults in the model, for the cases IES0, IESnd, IESR, and IESB. The typical behavior is that IES0 results in the most substantial update, as we would expect from the discussion of the analysis scheme in Section 3. The introduction of correlated measurement errors in IESnd leads to a slight increase in the spread of the fault multipliers associated with a slightly weaker update. The additional stochastic forcing used in IESR and IESB leads to further weakening of the parameter updates. The reason is that the additional stochastic forcing reduces the correlations between the model parameters and the predicted production. Contrary, the stochastic forcing The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments leads to an increased variance in the model-predicted rates, leading to an increased impact from the measurements. Thus, the total effect of including correlated rate errors is quite complicated to assess. The estimated porosity fields (not shown here) have similar behavior between the four examples, as seen for the fault multipliers.

Discussion of results
For the experiments IES0, IESnd, IESR, and IESB, the Figs. 8, 9, 10 and 11 present the production rates and the accumulated production for oil, gas, and water for OP2, and correspondingly the Figs. 12, 13, 14 and 15 give the results from well OP4. Additionally, Fig. 16 shows the prior and posterior estimates of the uncertain forcing rates for IESR. The producers, OP1, OP2, and OP3, had somewhat similar behavior while the match obtained for OP5 was nearly perfect in all the cases but with more substantial uncertainty in IESR and IESB. Thus, the figures include results from the wells OP2 and OP4, which appropriately illustrate the results. For OP2, the posterior ensemble has very low uncertainty in IES0. The ensemble spread increases slightly in IESnd, but it is still too small compared to the measurement uncertainty. The small ensemble spread is particularly the case early in the simulation. However, in IESR, the ensemble spread is significantly increased, and the posterior ensemble represents better the historical oil production.
For OP4, the prior ensemble covers the measurements well, but the posterior from IES0 overestimates the oil production, and the posterior ensemble does not cover the historical rates. Thus, the posterior uncertainty is too low. The solution from IESnd slightly increases the ensemble spread. On the other hand, in IESR, the posterior ensemble has a more considerable variance, and the ensemble mean shifts towards lower production, in better agreement with the measurements.
For the case, IESR, the posterior ensemble of uncertain historical rates in Fig. 16, has a lower variance than the prior Case IESB: Results from producer OP2. The upper plots show oil production rate OPR (left), and accumulated oil production OPT (right). Similarly, the middle plots are for gas production GPR and GPT, and the bottom plots show water production WPR and WPT. The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments one for all wells and variables. Also, we observe a shift in the historical rates, e.g., towards lower gas production for both wells OP2 and OP4, and lower water production for OP2. The magnitude of the update applied to the rates is modest. The update of the model parameters balances the rates' update, based on their relative impact on the solution and their prior variances.
A general observation is that case IES0 results in a significant improvement when compared to the prior ensemble.
Thus, the traditional history-matching approach, using the subspace EnRML, can generate a posterior ensemble solution with predictive skills. However, we notice that the posterior IES0 ensemble has too low uncertainty for most of the wells, with an ensemble spread that is much less than the variance in the historical measurements we condition on. This underestimation of the posterior variance results from using a diagonal measurement error-covariance matrix and ignoring the measurement errors' stochastic contribution in The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments the forward simulations. In several instances, the posterior ensemble is also unable to cover the observed history.
On the other hand, case IESR provides a slightly weaker update for the model parameters but compensates by updating the ensemble of perturbed historical rates. The result is a posterior prediction ensemble with significantly higher uncertainty, which in most instances, covers the historical data. The collective impact of updated rates and model parameters gives a better agreement with the observed production. The result is a posterior ensemble consistent with the historical rates and more substantial and realistic uncertainty than in case IES0.
The final experiment, IESB, assumes "biased" rate errors for each realization. The numerical implementation of the analysis scheme, using the formulation from Section 2.5, handles this case with a rank-one E matrix. IESB has qualitatively similar behavior to IESR, but the more considerable spread of the prior IESB ensemble leads to higher show oil production rate OPR (left), and accumulated oil production OPT (right). Similarly, the middle plots are for gas production GPR and GPT, and the bottom plots show water production WPR and WPT. The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments variance in the posterior IESB ensemble. We find that this ultimate increase of time correlations leads to a further weakening of the solution update. On the other hand, it seems that the overall match to the historical rate data may be better in IESB than IESR. Thus, maybe the assumption of infinite decorrelation time in the historical rate errors is accurate, e.g., like, if one extracts all a well's rate data using the same allocation table.
The uncertainty in the model parameters leads to biases in the individual realizations, and the posterior ensemble variance will grow with time, as we see from experiment IES0. However, the uncertainty in the controls in experiment IESR leads to more considerable uncertainty in the ensemble prediction early in the production period. However, after some time, the accumulated effect of the time-correlated rate-control errors starts canceling out. Thus, in the total The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments production, we don't see a massive increase in the production uncertainty. This situation changes in experiment IESB where the rate-control errors act as a constant bias in time. An effect of the increased production uncertainty early in the simulation is that the first measurements will have a more substantial impact on the update.
In our experiments, we do not impose constraints on the total production of oil, gas, and water, as these quantities are uncertain. Thus, we allow the ensemble realizations to produce different reservoir fluid quantities within their uncertainty range to span their uncertainty.

Summary and conclusions
The paper revisits the Bayesian formulation of the historymatching problem. It shows that we need to include stochastic errors in the rates used to force the ensemble  Fig. 15 Case IESB: Results from producer OP4. The upper plots show oil production rate OPR (left), and accumulated oil production OPT (right). Similarly, the middle plots are for gas production GPR and GPT, and the bottom plots show water production WPR and WPT. The error bars on the rate data indicate the plus-minus two standard deviations of the rate errors assumed in the assimilation experiments simulations for a statistically consistent formulation. We must further augment these rates to the model state vector where they, together with the uncertain model parameters, are being updated.
The current manuscript complements the discussion in Evensen and Eikrem [11] on the impact of correlated rate errors and redundancy in rate measurements. They claimed that the negligence of correlated measurement errors would lead to an underestimate of the posterior variance. In extreme cases, with vast data sets, this underestimation may appear as a so-called "ensemble collapse" where all realizations become nearly identical, and the variance approaches zero. In Evensen and Eikrem [11] this problem was eluded by conditioning only on the total accumulated production for each fluid in each well. The current paper has used an alternative and statistically sound approach retaining all the data. However, it prescribes the correct measurement error statistics and implements a consistent The paper presents the new formulation's properties on a simple reservoir case that illustrates the benefits of using a complete measurement error-statistics treatment. The algorithm used to solve the conditioning problem is the subspace EnRML method described by Evensen et al. [12], combined with the procedure for correcting stochastic model errors in iterative smoothers presented by Evensen [10].
A section in the paper is dedicated to further examining and illustrating the consistency of the ensemble-subspace inversion algorithm [7,8,12] for cases with correlated measurement errors. This paper uses a larger ensemble of measurement perturbations to represent the measurement error-covariance matrix with reduced sampling errors. The article also examines the impact of the measurement errors' different correlation length scales compared to the ensemble members' smoothness. The examples show that in the case with shorter length scales in the measurement perturbations than in the ensemble realizations, the shorter scales are truncated away when projecting the measurement perturbations onto the predicted measurements' ensemble. Thus, this is the only case where the new scheme introduces a significant approximation.
A practical implication of this work is that it resolves the reappearing question concerning the sampling frequency of rate data used in the conditioning. Using the current method, the sampling frequency of production rates does not matter as long as we specify realistic measurementerror correlations. The possibility of taking measurement error correlations into account will also be of vital importance when conditioning on sizeable seismic data sets. Seismic data often consists of gridded maps, and through the gridding procedure, one introduces additional error correlations. The proper formulation of measurement error statistics will also resolve the issues noted by Lorentzen et al. [19], related to the balanced impact of seismic data versus production data.
Algorithm 1, with the subspace inversion, is a computationally efficient approach for consistently solving the history matching problem. The algorithm's computational cost is linear in both the state dimension and the number of measurements, and the primary computation remains the integration of the model ensemble. Thus, Algorithm 1 is suitable for the conditioning of big models on big data.
The new algorithm's practical use will require reservoir engineers to model the rate data's error statistics properly. The common practice has been to set the error variance of the rate data to a fixed but sometimes arbitrary or subjective value and neglecting the presence of any measurement error correlations. In some situations, one has also inflated the measurement error variance to reduce the measurements' disproportionate impact when ignoring the measurement error correlations. Alternatively, it is possible to subsample the rate data in time to lessen their impact on the ensemble update. The recommendation is to specify as realistic as possible error variances for all the rate data, and there is no need for subsampling the data. The time correlations should reflect the procedure used to obtain the time series of rates. Unless one has access to a statistically sound method for estimating the time correlations of the measurement errors, the recommendation is to assume a time correlation with several months' decorrelation time. If we use allocation tables to obtain the production rates, then more than a year's decorrelation time is likely. Thus, for the reservoir engineers, the new algorithm doesn't pose a significant complication of use. In fact, by specifying one number, i.e., the time decorrelation of the errors in the production rates, one avoids previous issues with subsampling of data or inflation of measurement error variances to remedy formulation inconsistencies.
Finally, while the HM problem fits the model to past observations, the optimization problem maximizes future production or net-present value. The ensemble methods adapt well into a robust ensemble-based optimization framework [14,18] where one can optimize future production and net-present value, and one can assess new optimal well-placements while taking the geological uncertainty into account.

Appendix: Implementation in ERT
The reservoir example adopts the Ensemble Reservoir Tool (ERT) available from Github https://github.com/ equinor/ert.git for the management and conditioning of the model ensemble. The forward integration of the reservoir models is computed using the ECLIPSE reservoir simulator provided by Schlumberger. To be able to test the new formulation, some modifications were implemented in ERT and in the simulation job that runs ECLIPSE. The current version of ERT lacks functionality for simulating or specifying correlated measurement errors, and as most other assimilation implementations it works with a diagonal C dd . On the other hand, the actual implementation of the update scheme uses the same methods as discussed in the previous section, and has already full functionality for handling a nondiagonal measurement error-covariance matrix, C dd , or directly using the measurement perturbations, E, see Evensen et al. [12]. The following sections briefly discuss how ERT was modified to accommodate for correlated measurement errors and stochastic forcing of the reservoir simulations. For now, the implementation uses a casespecific implementation with communication through files, but the plan is to upgrade ERT to generally support correlated measurement errors.

A.1 Convergence of the subspace EnRML
The step-length scheme used in Algorithm 2 (i.e., the value of γ ), as implemented in ERT, is the following: One defines the maximum step length t 1 = 0.5, the minimum step length, t 2 = 0.2, and a step-length decline from one iteration to the next t 3 = 2.5, and we compute the value of γ i , for iteration i, from Here, γ i follows a geometrical series starting from the value t 1 in the first iteration and then reducing geometrically with the number of iterations towards t 2 . The formula in Eq. 49 allows us to define and test different steplength schemes easily. In real applications, it is desirable to start with a conservative step-length that works well "in most applications." If it turns out that the selected scheme becomes unstable, then the experiment should be restarted with a new and more conservative step-length. There is no exact theory on how much we would need to reduce the step size, but a good starting point can be to reduce t 1 = 0.5 to t 1 = 0.4.
The termination of the iterations can be based on the magnitude of the gradient, or one can use the relative change of the cost functions from one iteration to the next. For real reservoir applications, the time required per iteration can vary from hours to days, and the affordable number of iteration steps will always be limited. Thus, a practical procedure is to manually stop the iterations when one sees that the cost functions for the realizations are "almost" identical from one iteration to the next.
The rapid convergence of the assimilation algorithm for case IESR is illustrated in Fig. 17. Ten iterations were run, but after about five to six iterations there is only a marginal change in the parameters, and for practical applications five or six iterations should suffice for models of similar nonlinearity to the one used here.

A.2 Simulation of measurement perturbations
First, there is need of an application for simulating an ensemble of measurement perturbations for the time series of oil-production-rate (OPR), gas-production-rate (GPR), and water-production-rate (WPR), for each of the production wells. In the examples shown below, possible errors in the injection rates are ignored. The measurement perturbations can be white in time, red in time, or just a constant bias. The simulated ensemble of rate perturbations is then stored in a file. The current implementation computes perturbations for all the rates given in the ECLIPSE schedule file, no matter if they are used or not in the conditioning. The code used for simulating the measurement perturbations is available from the Github repository https://github. com/geirev/EnKF sampling/tree/ERTOBS/ERTsamp.

A.3 Importing E in the EnRML algorithm
In the current implementation, line 5 of Algorithm 2 illustrates how the subspace EnRML subroutine in ERT reads all the simulated rate perturbations from a file and store them in E 0 . The new measurement perturbations replace the original measurement perturbations supplied by ERT in D (line 6) and are input to the inversion in the iteration of W i (line 14). The product I e E 0 in lines 6 and 7, extracts the new simulated perturbations corresponding to the measurements contained in D, which can typically be a subset of all the rates represented in E 0 . Thus, with this simple modification, it is possible to compute the EnRML update while allowing for correlated measurement errors.

A.4 Updating measurement perturbations E i
As soon as one has computed the transition matrix, T i , for iteration i, in line 15 of Algorithm 2, it is possible to augment the measurement perturbations to the model state vector and update them according to The beautiful property of the subspace EnRML algorithm is its independence of the model state vector, which only enters as input to the forward model. Thus, augmenting the measurement errors to the state vector does not impact the computational procedure used to obtain the transition matrix. However, when the measurement errors are input to the forward model, they will change the model predicted measurements, and thus lead to a different transition matrix.

A.5 Updating historical rates in the schedule file
In iteration i, the forward model is forced by the updated perturbations, g(X i , E i ). Practically, one needs to update the Eclipse schedule file with new historical OPR, GPR, and WPR rates, by adding the updated perturbations to the corresponding "WCONHIST" rates. ERT handles this schedule modification by running a simple program, which is available from https://github.com/geirev/Schedule parser. git. This program reads the updated perturbations in E i from file and adds them to the historical rates defined in the original schedule file. ERT runs this program immediately before starting the Eclipse simulation for each realization.

IESR
Acknowledgements This work received support from the Research Council of Norway and the companies AkerBP, Wintershall-DEA, Vår Energy, Petrobras, Equinor, Lundin, and Neptune Energy, through the Petromaks-2 DIGIRES project (280473) (http://digires.no). We acknowledge the access to Eclipse licenses granted by Schlumberger.
Funding Open Access funding provided by NORCE Norwegian Research Centre AS.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.