1. Introduction

Broadly speaking, the measured data used for description of reservoir porosity and permeability characterization consist of static and dynamic data. Static data such as well logs and core samples can resolve heterogeneity at a scale of a few inches or feet with high reliability. However, dynamic data such as fractional flow or water cut (neglecting any pre-existing mobile water in the reservoir, this could be defined as the ratio of the injection fluid to the total fluid produced at the production wells), pressure transient, and tracer test data typically scan the length scales comparable to the interwell distances. Additional dynamic data such as time-lapse seismic images [1] can provide improved spatial sampling but at a lower precision. The ensemble Kalman filter (EnKF) is now being used in a number of studies for reservoir history matching. Some of the recent studies are listed in Evensen [2]; also see Nævdal et al. [3], Wen and Chen [4], Gu and Oliver [5], and Jafarpour and McLaughlin [6]. In general, reservoir data is often scale-dependent and exhibits multiscale features, and integration of additional multiscale data could further reduce the uncertainty (see Lee et al. [7], Efendiev et al. [8, 9] and references therein). Also, it is important to resolve fine-scale heterogeneity for various purposes such as enhanced oil recovery, environmental remediation, etc. With that perspective, integration of data at coarse and fine scales is an important objective. Computationally efficient assimilation of multiscale data using EnKF to estimate fine-scale fields for subsurface characterization is the main topic of this study. The main reason we used EnKF in this paper is because it requires fewer ensemble members than the particle filters (where, rather than updating the ensemble members model state, we update the probability assigned to each ensemble member based on model data misfit), e.g., see [10] and references therein for further details.

In this paper, apart from the water-cut data, we consider two kinds of coarse-scale measured data as well. The coarse-scale data are assumed to be permeability and/or saturation at some specified level of precision. The unknown variables (permeability, at the fine scale), are estimated using a modification to the EnKF algorithm, linking the data at different scales via upscaling (from the finest to the coarsest scales). The main idea behind upscaling is to obtain an effective coarse-scale permeability which yields the same average response as that of the underlying fine-scale field, locally. First we consider coarse-scale permeability data, which could be obtained either from geological considerations or coarse-scale inversion of dynamic, fractional flow data on a coarse grid [7, 9] or also using Markov Chain Monte Carlo (MCMC) techniques [8]. This coarse-scale, static data could be viewed as prior information regarding the permeability or in other words, a constraint which is to be satisfied up to the prescribed variance while obtaining the fine-scale estimates in every data assimilation cycle using the EnKF. Upscaling methods relate the solution at the fine scale to the coarse scale; therefore, in the Kalman filtering context, it amounts to modeling a nonlinear observation operator. In this paper we study two ways to assimilate the coarse-scale data using the EnKF. The standard EnKF [2] could be used for assimilating all the available data in one step, or alternatively, the measured data could be used in batches. For example, the estimate with one data becomes a prior while assimilating the other measured data; further details are given in Section 3.

The second kind of coarse-scale observed data we consider is dynamic and is motivated based on the increasing availability of time-lapse seismic images (or 4d seismic data). Integration of inverted 4d seismic data (at fine scale) using the EnKF has been addressed in Dong et al. [11] and Skjervheim et al. [12]. In this article we consider the seismic data, not to correspond to the finest scale but to a coarse scale, since time-lapse seismic data typically have a lower spatial resolution compared to the fine-scale geologic models [13]. Since the time-lapse seismic data is collected only at specific time intervals, we used coarse-scale fluid saturation as measured data to be available at a prescribed level of precision (which accounts for the inaccuracies involved in inversion of 4d seismic data) and only for certain assimilation cycles. Therefore, unlike the coarse-scale static permeability data considered earlier, the coarse-scale saturation data is assimilated only in certain assimilation cycles (see Section 4.3 for details).

Following is the plan of this paper. For the paper to be self-contained and for notational clarity, we briefly review the governing equations and sequential data assimilation using the EnKF in Section 2. This is followed by a description of the EnKF for assimilation of coarse-scale data in Section 3. For our numerical results in Section 4, we consider a five-spot pattern, with the injection well placed in the middle of a rectangular domain and four production wells located at the vertices of the rectangle. A reference case is used to provide true data, which is randomly perturbed to obtain synthetic measurements in a twin experiment setup. After presenting the assimilation results with both coarse-scale permeability and saturation data, we conclude with some directions for future work in Section 5.

2. Preliminaries

2.1 Fine-Scale Model

In this paper we consider two-phase flow in a subsurface formation under the assumption that the displacement is dominated by viscous effects. For simplicity, we neglect the effects of gravity, compressibility, and capillary pressure, although our proposed approach is independent of the choice of physical mechanisms. Also, porosity is considered constant. The two phases are referred to as water and oil, designated by subscripts w and o, respectively. We write Darcy’s law for each phase as follows:

(1)

(2a)

(2b)

The above descriptions are henceforth referred to as the fine-scale model of the two-phase flow problem. Here κ_f is the (fine-scale) permeability of the medium, λ(S) is the total mobility, μ_j denotes phase viscosity, pr is the pressure, h is the source term, and ϕ and S denote porosity and water saturation (volume fraction), respectively.

2.2 Sequential Estimation Using EnKF

Using dynamic measured data such as water cut, we can sequentially estimate the unknown parameters (permeabi- lity, porosity, etc.) and state variables such as pressure, water saturation (two-phase flow) and production data at well locations using the EnKF [3, 5, 6, 14]. The combined state-parameter to be estimated is given by = [ln(κ_f),pr,S, W_c]^T, where ln(·) is natural logarithm of the permeability field, W_c denotes water cut, and porosity is assumed to be known.

The EnKF introduced by Evensen [15] is a sequential Monte Carlo method where an ensemble of model states evolves in state-space, with the mean as the best estimate and spread of the ensemble as the error covariance, as summarized in the following steps. Each of the ensemble members is forecasted independently,¹

(3)

where F[·] is the forecast operator [fine-scale model Eqs. (1) and (2b)], superscript (i) denotes the ith ensemble member; from this point we on drop the time subscript. The ensemble mean and covariance are defined as

(4a)

(4b)

where A = (b⁽¹⁾,b⁽²⁾,...,b^(N_e)), b⁽ⁱ⁾ = Ψ⁽ⁱ⁾ - Ψ, and N_e is the number of ensemble members.

In a twin experiment, the observed water cut W_c^o is related to the truth via W_c^o = H[ Ψ^t], where H[ Ψ^t] is the true water cut. For each ensemble member, we randomly perturb W_c^o to generate observational samples,

(5)

where ν⁽ⁱ⁾ simulates observational error sampling, obtained as independent and identically distributed (iid) samples [16] from a normal distribution with zero mean and variance R. We note that if only the water-cut data is being measured, the mapping from model to observational space H is trivially equal to [0 0 0 I], since Ψ = [ln(κ),pr,S,W_c]^T .

The forecasted ensemble Eq. (3) is updated by assimilating the observed data,

(6)

where K is the Kalman gain, given by

Computationally efficient implementation of the EnKF is discussed, for example, in [2] and [17] (note that the ensemble error covariance, before or after assimilation, is not explicitly computed and we instead use the ensemble members for obtaining the covariance information). We use the above set of assimilated ensemble states, {⁽ⁱ⁾}_i=1^N_e, in the fine-scale simulation model Eq. (3) for prediction until the next set of observational data is available.

3. Coarse-scale data assmilation

The EnKF presented so far used only the dynamic production data (water cut) y with error ν = y - H[ Ψ^t], ν ~ (0,R) to update the ensemble Eq. (6). In addition to y, we now consider another independently measured data with independent errors z, which is static and is on a coarser scale compared to the fine-scale variables in Eqs. (1) and (2b). We assume that the corresponding measurement error is given by ω = z - U[ Ψ^t], with zero mean and Q covariance, ω ~ (0,Q), and U is a mapping of fine-scale variables Ψ to coarse-scale data, z, i.e., U : Ψ z. Then the likelihood of z is given by

(7)

If this static data z corresponds to coarse-scale permeability data (as considered in [9] and [7]), then U = [0 0 0], where : κ_f   κ_c is a nonlinear mapping that maps the fine-scale permeability field (κ_f) to coarse-scale field (κ_c) via an upscaling procedure (e.g., Durlofsky [18] and Durlofsky [19]). (Details are provided in Section 3.3.) Alternatively, if z corresponds to coarse-scale saturation inverted from 4d seismic data (as mentioned in the Introduction), then U = [0 0 0], such that is a mapping of fine-scale saturation S_f to coarse-scale saturation S_c = S_f. (Here we consider a simple volume averaging for ; further details are provided in Section 4.3).

Now, our goal is to obtain an estimate which is based on both water-cut and available coarse-scale data. The likelihood of water-cut data y is given by

(8)

The probability distribution function (pdf) of the predicted ensemble,

(9)

where Ψ and P^f are the predicted ensemble mean and covariance, respectively Eqs. (4a) and (4b). Then, using Bayes theorem, we obtain

(10)

Based on the above equation, following the (*) term, all the available data (z and y) could be assimilated in one step (details follow in Section 3.1), whereas based on the (†) term, the measured data y and z can be assimilated in a sequential manner. First assimilate the fractional flow (y) to obtain an ensemble conditioned on y, i.e., p( Ψ|y), which could then be used to assimilate the coarse-scale data z (further explained in the following Section 3.2).

3.1 Coarse-Scale Data Assmilation: In One-Step

All the available data, y and z, could be assimilated in one assimilation step by a modification to the model-to-observation space operator, H. If coarse-scale permeability data at a single coarse-scale is to be assimilated, H = [0 0 I]. Alternatively, if coarse-scale saturation data is available, H = [0 0 I]. The fine-scale estimated ensemble is obtained as in Section 2.2, with Eq. (6) modified to account for the additional coarse-scale data,

(11)

where . From now we on refer to this one-step assimilation of coarse-scale data procedure as reg EnKF. A consequence of the particular form of H is that it introduces contributions for cross-correlations between upscaled variables and production data. From a computational point of view, it has been known that EnKF is not very efficient for assimilation of large amounts of data [20], which could arise in our case, in complex three-dimensional cases, and also if data at multiple coarse scales is to be assimilated. In such a situation, different kinds of data could be assimilated in batches [21], which is described in the Section 3.2.

3.2 Coarse-Scale Data Assmilation: In a Sequence

We obtain an intermediate ensemble by assimilating y, denoted by {⁽ⁱ⁾}_i=1^N_e,

(12)

as discussed in Section 2.2. This intermediate ensemble and likelihood in Eq. (7) can then be combined [† term in Eq. (10)] to obtain the final estimate {⁽ⁱ⁾}_i=1^N_e,

(13)

Therefore, in a least-squared sense, the final estimate maximizes the posterior pdf p( Ψ|z,y) and corresponds to the minimum of = _z + _y + _f. See Appendix A for further details (where we show that the solution ⁽ⁱ⁾ corresponds to the minimum of , for any i^th ensemble member).

If coarse-scale data is available at only one coarse scale, then the fine-scale estimated ensemble is obtained by first assimilating production data followed by assimilation of the coarse-scale data,

(14a)

(14b)

^f is approximated using the intermediate ensemble ⁽ⁱ⁾; henceforth we refer to this sequential, coarse-scale EnKF data assimilation procedure as cs-EnKF. Note that data at multiple coarse scales can be sequentially assimilated by suitable repetition of Eq. (14b), with corresponding upscaling operators. For the coarse-scale saturation data, which may be available at only certain times, for only those assimilation cycles is Eq. (14b) applicable, whereas in the case of permeability data, considering it to be prior information regarding the fine-scale permeability, it is always to be honored; hence, both of the above steps (14a) and (14b) are to be always applied. The cs-EnKF algorithm is detailed in Appendix B, and a flow chart is given in Fig. 1. Implementation of this algorithm entails upscaling of each ensemble member at every assimilation step, i.e., N_e times the upscaling operator U[·] needs to be applied. In addition, if the dimension of the coarse-scale grid is N_c = n_c × n_c, then we need to perform an Singular Value Decomposition (SVD) of a rectangular matrix of size N_c × N_e. Hence, the total computation expense involves N_e upscales and SVD of the N_c × N_e matrix. Note that a similar upscaling is involved in the case of the reg EnKF, but the size of the matrix to compute SVD is now (N_c + N_wc) × N_e, where N_wc is the dimension of the water-cut data. In addition, if there are a number of coarse scales then the size of the matrix whose SVD is to be computed will grow for the reg EnKF, since all the data is assimilated in one step. Whereas for the cs-EnKF, coarse-scale data is assimilated in a sequence, the estimate from one coarse scale being used as a prior for the next scale, the matrix size to compute SVD is always N_c × N_e. (See also remark 3 of the coarse-scale EnKF algorithm in Appendix B.)

For nonlinear upscaling operators, such as , the final estimates from the reg EnKF would be different from those obtained using the cs-EnKF. When both the coarse-scale and water-cut data are assimilated together as in the above reg EnKF, it would imply fitting a multivariate normal likelihood to the different measured data together, whereas when the different kinds of data are assimilated one after another, as in the cs-EnKF, we fit each data separately, with a different pdf. (For further details on this topic, please see Dance [21] and references therein; also see Section 5.)

3.3 Upscaling Methods

In brief, the main idea behind upscaling of absolute fine-scale permeability is to obtain effective coarse-scale permeability for each coarse-grid block.Upscaling techniques in conjunction with the upscaling of absolute permeability have been used in groundwater applications (see, e.g., [19, 22, 23]). The link between the coarse- and the fine-scale permeability fields is usually nontrivial, because one needs to take into account the effects of all the scales present at the fine level. In the past simple arithmetic, harmonics, or power averages have been used to link properties at various scales. These averages can be reasonable for low heterogeneities or for volumetric properties such as porosity. For permeabilities, simple averaging can lead to inaccurate and misleading results. In this paper we use the flow-based upscaling methods.

Consider the fine-scale permeability that is defined on a domain with underlying fine grid as shown in Fig. 2. On the same plot a coarse-scale partitioning of the domain is also illustrated. To calculate the coarse-scale permeability field at this coarse level, we need to determine it for each coarse block, Ω_c. The coarse-scale permeability is computed so that it delivers the same average response as that of the underlying fine-scale problem, locally. The calculation of the coarse-scale permeability based on local solutions is schematically shown in Fig. 2. For each coarse domain Ω_c we solve the local problems

(15)

with some coarse-scale boundary conditions. An example of such boundary conditions is given by ϕ_j = 1 and ϕ_j = 0 on the opposite sides along the direction e_j and no flow boundary conditions on all other sides, alternatively, ϕ_j = x_j on ∂Ω_c. For these boundary conditions the coarse-scale permeability (κ_c) is given by

(16)

where ϕ_j is the solution of Eq. (15) with prescribed boundary conditions. Various boundary conditions such as periodic, Dirichlet, etc. can have some influence on the accuracy of the calculations. These issues have been discussed, e.g., in [24]. In particular, for determining the coarse-scale permeability field one can choose local domains that are larger than target coarse block, Ω_c, in Eq. (15). Furthermore Eq. (16) is used in the domain Ω_c, where ϕ_j are computed in the larger domains with correct scaling (see [24]). This way one reduces the effects of the artificial boundary conditions imposed on Ω_c (for details see [24]).

We denote by the local operator that maps the local fine-scale permeability field κ_f onto κ_c, defined on the coarse grid as in the above Eq. (16). For our computations we assume

(17)

where ϵ are some random fluctuations that represent inaccuracies in the coarse-scale permeability. One source of these fluctuations is the errors associated with solving inverse problems on the coarse grid. The other source of the inaccuracies include the fact that the inversion on the coarse grid does not take into account the adequate form of the coarse-scale models. Indeed, the inversion on the coarse grid for flow problems often involves the same flow equations as the underlying fine ones, for example, the same relative permeabilities are used for the coarse-scale problems as those for the fine-scale problems or the effects of macrodispersion are neglected. It is known that the flow equations at the coarse level may have a different form than the underlying fine-scale equations [19, 25-27]. In general, this form depends on the detailed nature of the heterogeneities, which are very difficult to obtain in solving inverse problems. In our paper we use Gaussian errors in Eq. (17) and consider the impact of coarse-scale data precision (i.e., nature of ϵ) by varying the variance of ϵ (see Section 4.2 for more details).

4. Numerical results

For our numerical tests, we use a 50 × 50 fine grid (dimensionless domain size 50 × 50) and two kinds of coarse-scale data in a twin experiment setup. First we consider coarse-scale permeability, which in reality, could be obtained by coarse-scale inversion of fractional flow data on a coarse grid [9, 28]. In this study we upscaled the reference fine-scale permeability (described below) to a 5 × 5 grid to obtain a coarse-scale permeability using flow based upscaling (Section 3.3). This coarse-scale field could be thought of as static data, which is to be honored as constraint (up to the prescribed measurement data variance) in Eq. (7); hence, we need to always assimilate it in every assimilation cycle. In reality we never know the reference field; therefore, this experimental setting is unrealistic. However, it allows us to compare and contrast a variety of test cases.

For the second set of results, a coarse-scale saturation is used which in practice could be obtained from inversion of 4d seismic measurements (see Section 1). Here, the coarse-grid saturation was obtained by volume averaging of true fine-scale saturation at some specific observation times (further details are given in Section 4.3). Therefore, unlike coarse-scale permeability, static data constraint, which is to be always satisfied, the coarse-scale saturation data is assumed to be available at only a few observation times. Following the flowchart in Fig. 1, we always have coarse-scale data if it is coarse permeability, and only at those few observation times for coarse-scale saturation data.

An initial ensemble with different permeability realizations was generated using the sequential Gaussian simulation (Deutsch and Journel [29]). We specified a Gaussian variogram model with a correlation length of 20 gridblocks in the x direction and 5 gridblocks in the y direction. One of the realizations is used as the “true” field (shown in Fig. 3) and was removed from the ensemble. Porosity (ϕ) is assumed to be equal to 0.15 for all grid blocks. For simplicity, relative permeabilities k_rj are assumed to be linear functions of water saturation (S): k_rw(S) = S,k_ro(S) = 1 - S. One injection well at the center of the field (injection rate: 71.4 m³/day) and four producing wells at the four corners (all with equal rate of 17.85 m³/day) were considered. The fine-scale model Eqs. (1)-(2b) are solved with no flow boundary conditions, zero initial water saturation, and by discretizing the transport equation using the first-order upwind finite volume method. In the top panel of Fig. 4 we provide the predicted fractional flow for 256 initial ensemble members along with the true fractional flow (obtained from the true permeability field).

To compare and contrast our results using coarse-scale data and different versions of EnKF, we use the following mean L₂-norm error. Since we know the true (fine- and coarse-scale) field for our synthetic problem, denoting the true permeability field by κ^true, the error for any ensemble member is given by

Consider the L₂ norm of the error for each member, ||e⁽ⁱ⁾||₂ = Σ_j[e⁽ⁱ⁾_j]₂, by which we define the mean L₂ error as

(18)

so that e gives us an indication of the distance of the entire ensemble from the true solution κ^true. Since after assimilating any observation we updated all the ensemble members, we can monitor the variation of e over the time of assimilation; the success of assimilation can therefore be related to the decrease in e.

4.1 EnKF with Water Cut Data Only

We start with a presentation of results obtained with assimilation of water-cut data only. Next we discuss results with coarse-scale data.

The water-cut data from the reference field is assumed to be available every 200 days, with mean zero and standard deviation of 0.01 (therefore R^1/2 = 0.01I₄, where I₄ is the unit matrix of size 4 × 4, since there are four producing wells). The observed data is assumed to be available up to 2400 days; hence, we performed assimilation between 200 and 2400 days. A prediction beyond the interval of data assimilation, up to 4000 days, is also provided.

The choice of ensemble size (N_e) is very important for successful data assimilation using EnKF. This is because a finite size ensemble prediction is used to estimate the prior error covariance P^f Eq. (4b). For small sample sizes, sampling errors in the covariance estimates result in insufficient variance for P^f, so that observations which lie outside the small ensemble spread are completely ignored [17, 30]. (We are trying to sample a covariance matrix for unknown variables: ln(κ),pr,S,W_c, i.e., an unknown of size 3 × 2500 plus four fractional flow data in this case, using sample sizes that are far lesser, resulting in severely reduced rank covariance matrices.) Different approaches such as covariance inflation and localization have been proposed to alleviate this problem of ensemble inbreeding, which is discussed elsewhere (see [31-34] and references therein for further details). An ensemble with sufficiently large number of members needs to be selected so that the assimilation system would not severely suffer from the above-described problem.² Here we present our data assimilation results with N_e = 256, and in Sections 4.2 and 4.3 we briefly discuss some important characteristics of the error covariance matrix such as variance and eigenspectrum in the context of coarse-scale data assimilation. The issue of coarse-scale data assimilation with smaller ensemble sizes would be tackled in the future.

In the bottom panel of Fig. 4 we plot the ensemble and true water-cut data after assimilation of only water-cut data with the EnKF. Comparing with the initial forecast (top panel), we observe that the assimilated ensemble better envelopes the true data. Also, the ensemble mean saturation field after 500, 1000, 2000, 3000, and 4000 days of simulation better compares with the true saturation than with no assimilation in Fig. 7. The final permeability field after assimilation for the ensemble mean and a few members is compared with the true field in Fig. 9. Note that the central, southeast-northwest channel is prominent, but the features at the southwest and northeast corners are not well captured, which is reflected in the plot of mean saturation (Fig. 7), where many fine-scale features present in the true saturation field are not present in the ensemble mean. Therefore, assimilation of only water-cut data helps in identifying some of the important features.

4.2 EnKF with Water-Cut and Coarse-Scale Permeability Data

In addition to water-cut production data, the coarse-scale permeability data, as described in Section 3.3 was used as additional measured data. Flow-based upscaling of the reference permeability field was used as a proxy for the inverted coarse field. Following our previous notation, this coarse-scale permeability data is denoted by z Eq. (7). The mapping between state variables (at fine-scale) and observations (at coarse-scale) as given by U = [0 0 0], , denotes flow-based upscaling. For the reg EnKF, H = [0 0 I] in Eq. (11) of Section 3.1.

Exactly as in the previous section, we prescribed the same frequency (of availability) and precision R for the water-cut data. Since we use coarse-scale permeability as additional data, it is to be assimilated whenever we assimilate water-cut data. For 5 × 5 coarse-scale data with mean zero and variance, Q = qI₂₅ (we present results with q = 4,2,1, and 0.1), so that we can consider the impact of coarse-scale data precision. In the left panel of Fig. 5 we plot the variation of mean L₂ error e Eq. (18) with observation time at the coarse scale for different values of q and using reg EnKF as well as cs-EnKF. In the right panel of the same figure we show the correlation between coarse-scale ensemble mean and true fields for q = 1. The values of correlation coefficients for different values of q are provided in Table 1. Note that as the precision of coarse-scale data is increased, i.e., for smaller value of variance, we observe a larger decrease in coarse-scale mean L₂ error and higher correlation with true coarse-scale field. This would be expected because smaller variance Q implies more strict coarse-scale data constraint in Eq. (7) and hence, the coarse-scale data is more accurately assimilated as it is made more precise. The water-cut data prediction using the final permeability field after assimilation for different coarse-scale data precisions is plotted in Fig. 6. (The nature of results with q = 4 is similar to those with q = 2,1,0.5; thus, we drop it.) Notice the improved fit of ensemble prediction to the true data for more precise coarse-scale data and also when compared to the assimilation of only water cut in Fig. 4, which is a consequence of the additional coarse-scale data being available. However, the water-cut prediction with the cs-EnKF compares better with the truth than that with the reg EnKFfor higher values of q = 2 and 1; with q = 0.5, the reg EnKF prediction is highly improved. In Fig. 7 we compare the ensemble mean with the true saturation. Once again, as the coarse-scale data constraint is more precisely imposed, the ensemble mean saturation captures most of the features in the true field. Also, notice that the reg EnKF saturation prediction improves more markedly as q is lowered when compared to the cs-EnKF, which could explain the better water-cut fit in Fig. 6 for q = 0.5 with the reg EnKF.

In the left panel of Fig. 8 we plot the fine-scale mean L₂ error for the permeability field with different values of q, and in the right panel of the same figure we plotted the correlation between the ensemble mean and the true fine-scale permeability for q = 1. The correlation coefficients are given in Table 1. (For assimilation of only water-cut data, we obtained a correlation coefficient equal to 0.3074.) Though the mean L₂ error is lower with the cs-EnKF, a slightly higher correlation is obtained with the reg EnKF. We observe that higher precision, i.e., lower q, does not necessarily imply highest correlation, whereas we obtained a lower mean L₂ error. The final permeability field after assimilation for the ensemble mean and a few members is shown in Fig. 9. We note that the low-permeability region at the northeast and high permeability at the southwest corners are well captured. Also, the ensemble mean permeability with the cs-EnKF has some features not very well present, in that with the reg EnKF, particularly with q = 2,1, e.g., the low-permeability region to the left of the central southwest-northeast channel. As we noted in Section 3.1, the estimates obtained using the reg EnKF are expected to be different from those obtained with the cs-EnKF. We observed that for higher values of the coarse-scale data variance (q = 2 and 1), the cs-EnKF yielded better water cut, ensemble mean saturations, and permeability estimates when compared to the reg EnKF, whereas with q = 0.5, the results are similar (but not exactly same) with both versions of the EnKF. These results indicate that an optimal value for the coarse-scale data variance is important, particularly for the different versions of the EnKF (reg EnKF or cs-EnKF), and it could be obtained by a prior calculation of the uncertainty in the coarse-scale data which can be addressed in a future study. Also, more complex, studies in coarse-scale data at more than one scale and for three-dimensional models are needed to understand the merits and demerits of each version of the EnKF.

Regardless of the version of EnKF being used, after every assimilation cycle we can obtain an estimate of the analysis error covariance matrix for the ensemble fine-scale permeability. Using the fine-scale assimilated permeability fields κ_f⁽ⁱ⁾ we define the following ensemble mean permeability and error covariance,

(19a)

(19b)

Minimization of variance of the covariance matrix P_κf,κf is desired in various Kalman filtering applications [2] as it provides a measure of uncertainity. In the top row (left panel) of Fig. 10 we plotted the normalized variance of P_κf,κf during assimilation with various precisions of coarse-scale data and for both reg EnKF and cs-EnKF, and also for assimilation of only fractional flow data (we normalized using the variance from the initial ensemble). We observed that a higher reduction in variance is obtained as coarse-scale data precision is increased, and the cs-EnKF obtains more reduction than the reg EnKF. This trend is seen even for an ensemble of much larger size, for example, with 1000 members as shown in the bottom row (left panel) of Fig. 10.

Another important property of interest, particularly for ensemble data assimilation is the eigenspectrum of P_κf,κf, which is plotted in the right panels of Fig. 10; the top row represents 256 ensemble members, whereas the bottom row represents 1000 members. This eigenspectrum is from the last assimilation cycle at the end of 2400 days; for other assimilation cycles, the spectrum was similar. A smooth distribution of eigenvalues is desired, and indeed, we observe that, but the tails of the spectrum are steep for only water-cut data assimilation and cs-EnKF. As the coarse-scale data precision is improved, the trailing eigenspectrum gets smoothed out, which may indicate that the eigenvectors associated with the small eigenvalues are able to resolve some small-scale correlation features. Also, the larger variance associated with the leading eigenvectors as obtained for only water-cut data assimilation seems not to be the case with coarse-scale data assimilation (some of the desirable properties regarding eigenspectrum of covariance matrices are given, e.g., by Hamill et al. [32]). We note that for an ensemble size of N_e = 1000, the spectra are almost identical for both versions of the EnKF. In Fig. 11 we compare side-by-side the spectra with N_e = 256 for the reg EnKF and cs-EnKF, with the 1000 ensemble member eigenspectrum. (Since the spectra are similar for the 1000 ensemble member case, we plot and compared only with the reg EnKF.) As noted in Fig. 10, with higher precision of coarse-scale data, the tail of the spectrum gets smoothened, but also the magnitude of trailing eigenvalues is decreasing, which suggests that for smaller ensemble sizes there could be issues with loss of rank of the error covariance matrix. This problem seems to be slightly more aggrevated for the cs-EnKF. When compared to the 1000 member case, the leading eigenvalues seem to be slightly larger, as observed in [32]. In any case, in this study we did not apply any covariance inflation or localization (and always used an ensemble of size 256, which seems to preserve the rank of ensemble up to about 245 eigenvalues; see Fig. 11). These topics in the context of coarse-scale data assimilation will be covered in a future study.

4.3 EnKF with Water-Cut and Coarse-Scale Saturation Data

As mentioned in the Introduction, by coarse-scale inversion of 4d-seismic data we could obtain dynamic data such as coarse-scale pressure and saturation. In this section we attempt to assimilate such a coarse-scale saturation in a twin experiment along with the fractional flow data using the EnKF algorithms discussed in Sections 3.2 and 3.1. To this end, the saturation obtained by using the reference permeability is saved at three different times: 200,1200, and 2400 days, which respectively correspond to the start, middle, and end of the time window of data assimilation. This saved fine-scale saturation field is then upscaled (see Section 3.3) by volume averaging to a 5 × 5 coarse-scale grid and used as observed coarse-scale saturation data. If we denote the volume averaging by operator , acting on fine-scale saturation S_f, to give coarse-scale saturation S_c = S_f, then the mapping between state variables at fine scale and measured data at coarse-scale is given by U = [0 0 0]. Therefore for the reg EnKF (Section 3.1), the measured data is related to the fine-scale variables via H = [0 0 I] in Eq. (11). For the cs-EnKFwe use above U operator to compute the misfit: z-U[ ] in Eq. (14b) (i.e., steps 2.1 and 2.4 of the cs-EnKF algorithm in Appendix B). Unlike the coarse-scale permeability data which was taken into account at every assimilation step, by construction, the coarse-scale saturation data is available only at a few assimilation steps, in this particular case, assimilation after 200,1200, and 2400 days.

To be consistent with our previous results, the frequency (of availability) and precision R for the water-cut data has been kept the same. For the coarse-scale saturation data we prescribed zero mean and variance, Q = q_sI₂₅, with q_s = 0.1,0.01, such that the precision is varied from low to high. Since the saturation ranges between 0 and 1, and the fractional flow data is usually more accurately measured than 4d-seismic data, we picked q_s to be always larger than the variance in fractional flow data. In the left panel of Fig. 12 we plotted the variation of mean L₂ error for the coarse-scale saturation (while assimilating) vs observation time. (For our test case, we had assumed zero initial water saturation; therefore, the water saturation increases in time and hence, the inherent, increasing trend in this figure.) Note that whenever the coarse-scale saturation is assimilated (200,1200, and 2400 days) the error decreases for both values of q_s considered. The water-cut data prediction using the assimilated ensemble members is given in Fig. 13; the fit of ensemble water-cut with the truth gets better as the coarse-scale saturation is prescibed higher precision.

We discuss the fine-scale results starting with fine-scale saturation and then the fine-scale permeability. The ensemble mean saturation is compared to the true field at certain times in Fig. 14. By assimilating the coarse-scale saturation data we are able to capture many of the subtle features that are present in the true saturation field, such as the fingers that develop off the center toward the northeast corner and sharp contrast between different levels of saturation, throughout the entire time interval (up to 4000 days) considered. Regarding the fine-scale permeability, a comparison of the mean L₂ error is shown in the right panel of Fig. 12. Note the marked reduction in error when coarse-scale saturation data is assimilated. Also, the correlation of the ensemble mean permeability with the truth is improved as q_s is decreased, as shown in Table 2. The fine-scale permeability fields for a few ensemble members and the mean are shown in Fig. 15. Based on the above results, we observed that unlike in the case of coarse-scale permeability data assimilation, coarse-scale saturation data assimilation does not yield significant reduction in error (comparing the left panel in Fig. 8 and the right panel in Fig. 12); neither are the correlations of fine-scale ensemble mean with the truth (second column of Tables 1 and 2). The results (ensemble mean saturation and permeability) are improved when compared to assimilation of only water-cut data, particularly with q_s = 0.01, but again, not as much improved as with coarse-scale permeability data assimilation. This could be anticipated, since the fine-scale permeability is more correlated to coarse-scale permeability than to the coarse-scale saturation. Regarding the two versions of the EnKF considered for assimilating the coarse-scale saturation data, results are very similar to each other, even those for the variance of error covariance matrix P_κf,κf and its eigenspectrum at the end of 2400 days (Fig. 16). Note the steepness of the trailing eigenvalues with q_s = 0.1 (for both reg EnKF and cs-EnKF), which is similar to the assimilation of only fractional flow data; however, for q_s = 0.01, this undesirable effect has been smoothned out. This is similar to the result obtained with coarse-scale permeabilty data assimilation in Fig. 10. Our observation that both versions of EnKF performed similarly could be due to the linearity of the fine-scale to coarse-scale saturation mapping ; however, identical results would not be possible due to the different analysis equations and sampling of errors used in Eqs. (14b) and (11).

5. Conclusions

The EnKF is increasingly being used for subsurface characterization in various geological and groundwater applications to identify fine-scale state and parameters. Recently, dynamic data other than production data has been considered in the EnKF context [11, 12]; nevertheless, the observed data to be assimilated was assumed to be at the finest scale. For a number of reasons, it is widely recognized that usage of additional multiscale data could further reduce the uncertainty at the fine scale. Also, it is often important to preserve large-scale features of the permeability field. These are coarse-scale features that can typically represent connectivity of the media. For example, facies consisting of high permeable regions are important for flow and transport. On the other hand, these features alone are not sufficient to match water-cut data. By performing inversion on the coarse grid, we can capture these large-scale features with more accuracy and certainty, and in turn, the coarse-scale inverted permeability could be used as a prior constraint in EnKF data assimilation of the fine-scale fields.

Here we proposed assimilation of coarse-scale data along with water-cut production data. We showed that the modifications to the EnKF for multiscale data assimilation are completely recursive and easily implementable. The relation between fine and coarse scales was modeled via flow-based upscaling, which could be thought of as a nonlinear observation operator linking the coarse-scale data to the unknown fine-scale variables. In addition, the proposed methodology could be used in any other sequential data assimilation method, as well to assimilate data at multiple coarse scales. Two versions of EnKF were suggested: (i) reg EnKF, where all the data (coarse-scale and water-cut) were assimilated together, and (ii) cs-EnKF, where the data were assimilated sequentially in batches. Ensemble members obtained after assimilating water-cut data are used as a prior to assimilate coarse-scale data. Though in our current paper we used only one coarse scale, the proposed method can be easily implemented to integrate as many coarse scales as required by the available data. Also the methodology is independent of the upscaling operator.

The assimilation setup was tested and compared for a two-dimensional synthetic 50 × 50 heterogenous true field. Two kinds of coarse-scale data were considered. In the first implementation, coarse-scale permeability data was considered and in the second, coarse-scale saturation. In our twin experiment setup, both of these data were derived from the reference field and, in both cases a 5 × 5 coarse grid was used. The coarse-scale data variance was varied from low to high in order to study its impact on fine-scale assimilated fields and water-cut predictions. In all cases we observed that the assimilated, ensemble mean coarse-scale (permeability/saturation) field for all variances was highly correlated to the true coarse-scale field. In addition, lower variance in the coarse-scale data yielded higher correlation. The water-cut data was better honored, for higher precision of coarse data. When assimilating coarse-scale permeability we observed that the cs-EnKF gave better fit with the true saturation, water-cut, and fine-scale permeability than the reg EnKF. It also yielded less error in an averaged L₂ norm error taken with regard to the reference field, whereas both versions of EnKF performed similarly when assimilating coarse-scale saturation. As shown in Appendix A, for a linear observation operator, assimilation of coarse-scale data in batches (i.e., as in cs-EnKF) or in shot (as in reg EnKF) would yield the same estimate. Then the cause for the difference in the performance of the two versions of EnKF for assimilation of coarse-scale permeability could be either due to ensemble size or linearity/nonlinearity of the upscaling (observation) operator. This issue has been outlined by Dance [21] and references therein. As far as the number of ensembles is concerned, we observed that even with a larger ensemble size of 1000 members, we noticed different performance of the two versions of EnKF (Fig. 10, reduction in normalized variance of P_κf,κf while having similar ranks in terms of the eigenspectrum). Certainly further study with different upscaling operators and coarse-scale data at multiple levels would be needed to study this aspect of the two versions of EnKF considered here. Over all, inclusion of coarse-scale data replicated many of the subtle features present in the fine-scale permeability and saturation fields which were not present after assimilating only water-cut data.

As our results indicate that the inclusion of coarse-scale data enhances identification of the multi-scale reservoir characteristics, it is important to study methods to obtain coarse-scale data as well as its precision. In a realistic scenario, coarse-scale inversion [7, 28], and in the future perhaps with more computing resources, MCMC methods [8], could be used to obtain such data. The coarse-scale saturation obtained using inversion has been shown to yield improved estimates in a three-dimensional reservoir case by Devegowda et al. [35]. Our current and future work is directed toward obtaining coarse-scale data with higher precision and its assimilation using ensembles of smaller size for complex three-dimensional cases.

Acknowledgments

The work of Akhil Datta-Gupta and Yalchin Efendiev is partially supported by DOE (DE-FG03-00ER15034), NSF CMG 0724704, and KAUST award number KUS-C1-016-04.

References

1. Lumley, D. E., Time-lapse seismic reservoir monitoring, Geophysics, 66:50-53, 2001.

2. Evensen, G., Data Assimilation: The Ensemble Kalman Filter, Springer, 2006.

3. Nævdal, G., Johnson, L., Aanonsen, S., and Vefring, E., Reservoir monitoring and continuous model updating using ensemble Kalman filter, SPE J., 10:66-74, 2005.

4. Wen, X.-H. and Chen, W. H., Real-time reservoir model updating using ensemble Kalman filter, SPE Reservoir Simulation Symposium, SPE 92991, 2005.

5. Gu, Y. Q. and Oliver, D. S., The ensemble kalman filter for continuous updating of reservoir simulation models, J. Energy Resour. Technol.-Trans. ASME, 128:79-87, 2006.

6. Jafarpour, B. and McLaughlin, D. B., History matching with an ensemble Kalman filter and discrete cosine parametrization, Comput. Geosci., 12:227-244, 2008.

7. Lee, S. H., Malallah, A., Datta-Gupta, A., and Higdon, D., Multiscale data integration using Markov random fields, SPE Reservoir Eval. Eng., 5(1):68-78, 2002.

8. Efendiev, Y., Datta-Gupta, A., Ginting, V., Ma, X., and Mallick, B., An efficient two-stage Markov chain Monte Carlo method for dynamic data integration, Water Resour. Res., 41:W12423, 2005.

9. Efendiev, Y., Datta-Gupta, A., Osako, I., and Mallick, B., Multiscale data integration using coarse-scale models, Adv. Water Resour., 28:303-314, 2005.

10. Liu, Y. and Gupta, H. V., Uncertainty in hydrologic modeling: Toward an integrated data assimilation framework, Water Resour. Res., 43:W07401, 2007.

11. Dong, Y., Gu, Y., and Oliver, D. S., Sequential assimilation of 4D seismic data for reservoir description using the ensemble Kalman filter, J. Pet. Sci. Eng., 53:83-99, 2006.

12. Skjervheim, J.-A., Evensen, G., Aanonsen, S., Ruud, B., and Johansen, T., Incorporating 4D seismic data in reservoir simulation models using an ensemble Kalman filter, SPE J., 12:282-292, 2007.

13. Gosselin, O., Aanonsen, S. I., Aavatsmark, I., Cominelli, A., Gonard, R., Kolasinski, M., Ferdinandi, F., Kovacic, L., and Neylon, K., History matching using time-lapse seismic (HUTS), SPE Annual Technical Conference and Exhibition, Denver, SPE 84464-MS, Oct 2003.

14. Gu, Y. and Oliver, D. S., History Matching of the PUNQ-S3 reservoir model using the ensemble Kalman filter, SPE J., 10:217-224, 2005.

15. Evensen, G., Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, J. Geophys. Rese., 99(10):14310, 162, 1994.

16. Burgers, G., van Leeuwen, P. J., and Evensen, G., Analysis scheme in the ensemble Kalman filter, Mon. Weather Rev., 126(6):1719-1724, 1998.

17. Chen, Y. and Zhang, D., Data assimilation for transient flow in geologic formations via ensemble Kalman filter, Adv. Water Resour., 29:1107-1122, 2006.

18. Durlofsky, L. J., Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media, Water Resour. Res., 27:699-708, 1991.

19. Durlofsky, L. J., Coarse scale models of two phase flow in heterogeneous reservoirs: Volume averaged equations and their relationship to the existing upscaling techniques, Comput. Geosci., 2:73-92, 1998.

20. Lawniczak, W., Hanea, R., Hemmink, A., and McLaughlin, D., Multiscale ensemble filtering for reservoir engineering applications, Comput. Geosci., 13:245-254, 2009.

21. Dance, S. L., Issues in high resolution limited area data assimilation for qualitative precipitation forecasting, Physica D, 196:1-27, 2004.

22. Durlofsky, L. J., Behrens, R. A., Jones, R. C., and Bernath, A., Scale up of heterogeneous three-dimensional reservoir descriptions, SPE Paper 30709, 1996.

23. Durlofsky, L. J., Jones, R. C., and Milliken, W. J., A nonuniform coarsening approach for the scale up of displacement processes in heterogeneous media, Adv. Water Resour., 20:335-347, 1997.

24. Wu, X. H., Efendiev, Y. R., and Hou, T. Y., Analysis of upscaling absolute permeability, Discrete Contin. Dyn. Syst., Ser. B, 2:185-204, 2002.

25. Efendiev, Y. R. and Durlofsky, L. J., Numerical modeling of subgrid heterogeneity in two phase flow simulations, Water Resour. Res., 38(8):W1128, 2002.

26. Efendiev, Y. and Durlofsky, L., Generalized convection-diffusion model for subgrid transport in porous media, SIAM Multiscale Model. Simul., 1:504-526, 2003.

27. Efendiev, Y. and Durlofsky, L., Accurate subgrid models for two-phase flow in heterogeneous reservoirs, SPE J., 9:219-226, 2004.

28. Yoon, S., Malallah, A. H., Datta-Gupta, A., Vasco, D. W., and Behrens, R. A., A multiscale approach to production-data integration using streamline models, SPE J., 6:182-192, 2001.

29. Deutsch, C. V. and Journel, A., GSLIB Geostatistical Software Library and User‘s Guide, Oxford Univ. Press, 1992.

30. Hendricks Franssen, H. J. and Kinzelbach, W., Real-time groundwater flow modeling with the Ensemble Kalman Filter: Joint estimation of states and parameters and the filter inbreeding problem, Water Resour. Res., 44:W09408, 2008.

31. Anderson, J. L., Exploring the need for localization in ensemble data assimilation using a hierarchical ensemble filter, Physica D, 230(1-2):99-111, 2007.

32. Hamill. T., Whitaker, J., and Synder, C., Distance-dependent filtering of background error covariance estimates in ensemble Kalman filter, Mont. Weather Rev., 129:2776-2790, 2001.

33. Arroyo-Negrete, E., Devegowda, D., Datta-Gupta, A., and Choe, J., Streamline-assisted ensemble Kalman filter for rapid and continuous reservoir model updating, SPE Reservoir Eval. Eng., 11(6):1046-1060, 2008.

34. Sun, A. Y., Morrisa, A., and Mohantya S., Comparison of deterministic ensemble Kalman filters for assimilating hydrogeological data, Adv. Water Resour., 32(2):280-292, 2009.

35. Devegowda, D., Akella, S., Datta-Gupta, A., and Efendiev, Y., Interpretation of partitioning interwell tracer tests using EnKF with coarse-scale constraints, SPE Reservoir Simulation Symposium, SPE 119125-MS, The Woodlands, TX, 2-4 Feb 2009.

---

APPENDIX A. Two-Step Coarse-Scale Constrained Kalman Filter Estimate

From Section 3,

and

For notational simplicity we denote μ _Ψ as μ and denote P^f by B.

Step 1 (minimize _f + _y):

First we minimize the sum, ₁ = _f + _y. The gradient³ of the above quadratic cost functional with respect to (w.r.t.) Ψ is given by

Then the minimizer of ₁ satisfies (we assume H to be linear)

By rearranging the above equation we get

(20)

Note that the Hessian of ₁ w.r.t. Ψ is given by B^-1 + H^T R^-1H, and for linear quadratic cost functionals, the Hessian inverse is equal to the error covariance matrix. Therefore, the error covariance matrix for is given by

(21)

Step 2 (minimize _g + _z):

We use , in

Therefore, the minimum of _g + _z satisfies

Using Eqs. (21) and (20) we can rewrite the above as

It is trivial to show that also satisfies

Therefore, the two-step method to obtain the final estimate gives the same results as a one-shot approach of minimizing _f + _y + _z.

APPENDIX B. The Coarse-Scale EnKF Algorithm

---

Algorithm 1. Coarse-scale EnKF algorithm

---

Run the simulation model up to a particular observation time for the entire ensemble to get predicted samples: {Ψ⁽ⁱ⁾}_i=1^N_e, A = (Ψ⁽¹⁾, Ψ⁽²⁾,...,Ψ^(N_e)).

• Step 1: Using measured water-cut data y with variance R, get the updated ensemble: { Ψ⁽ⁱ⁾}_i=1^N_e.

  Step 1.1—Find ensemble mean [Eq. (4a)] Ψ.

  Step 1.2—Subtract deviation from the mean A' = (b⁽¹⁾, b⁽²⁾,...,b^(N_e)), b⁽ⁱ⁾ = Ψ⁽ⁱ⁾ - Ψ.

  Step 1.3—Apply H to each column of A' to get S = H A', i.e., simply pick the water-cut deviations in A'.

  Step 1.4—For i = 1,2,...,N_e,

  

  is predicted water cut for each ensemble member.

      End for

  Step 1.5—Compute SVD [S + R^1/2] = X_LΣX_R.

  Get retaining the first few singular values which explain most variability in Σ, corresponding to the left singular vectors: _L.

  Step 1.6—Update ensemble: Eq. (6),

  

• Step 2: Using coarse-scale data z with variance Q, get the updated ensemble: {⁽ⁱ⁾}_i=1^N_e.

  Step 2.1—Compute coarse-scale ensemble prediction:

  

  Step 2.2—Coarse-scale mean:

  

  Step 2.3—Coarse-scale deviations:

  

  Step 2.4—Repeat step 1.4 using coarse-scale measurement.

  For i = 1,2,...,N_e,

  

  End for

  Step 2.5—Compute SVD [S' + Q^1/2] = X_LΣX_R. Get and _L as in step 1.5.

  Step 2.6—Compute fine-scale mean:

  

  Step 2.7—Compute fine-scale deviations:

  

  Step 2.8—Update ensemble:

  

---

Remark 1:
Note that steps 2.6 and 2.7 in the above algorithm approximate the intermediate fine-scale error covariance,

Remark 2:
Steps 2.1-2.3 accomplish⁴

Note that the above algorithm is independent of the choice of upscaling procedure, and we can use the same algorithm for different kinds of coarse-scale observed data (if available).

Remark 3:
Note that the above coarse-scale constrained EnKF algorithm can be readily extended to incorporate data at multiple coarse scales with the appropriate upscaling procedure in U. To elaborate, if we had other independent data at a scale different from z, we could use the estimates ({⁽ⁱ⁾}_i=1^N_e) obtained using z. As an intermediate solution, repeating step 2 to assimilate the data at another scale.

---

¹ In this work we focus primarily on assimilation of coarse-scale data using the EnKF, its feasibility, and impact on fine-scale estimates with different kinds of coarse-scale data (Section 3); hence, we neglect modeling errors, which will be addressed in the future.

² Our choice of N_e = 256 was based on observing the eigenspectrum and variance, discussed in Sections 4.2 and 4.3, by comparing the results for the 256 ensemble with a 1000-sized ensemble; the 256-sized ensemble did not suffer from the insufficient variance problem discussed above.

³ We note in passing that B and R are covariance matrices and are positive definite by construction and hence, for our derivation purposes, are formally invertible.

⁴ As noted in [2], this approach of accounting for the nonlinear observations operator U works well as long as U is weakly nonlinear and a monotonic function of model variables.