Unformatted text preview: SPE 125101
Data Assimilation Using The Constrained Ensemble Kalman Filter
Hemant A. Phale, Dean S. Oliver, SPE, University of Oklahoma
Copyright 2009, Society of Petroleum Engineers
This paper was prepared for presentation at the 2009 SPE Annual Technical Conference and Exhibition held in New Orleans, LA, USA, 4–7 October 2009.
This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s).
Contents of the paper have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not
necessarily reﬂect any position of the Society of Petroleum Engineers, its oﬃcers, or members. Electronic reproduction, distribution, or storage of any part of this
paper without the written consent of the Society of Petroleum Engineers is prohibited. Permission to reproduce in print is restricted to an abstract of not more
than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of SPE copyright. Abstract
When the ensemble Kalman ﬁlter (EnKF) is used for history matching, the resulting updates to reservoir properties
sometimes exceed physical plausible bounds, especially when the problem is highly nonlinear. Problems of this type
are often encountered during history matching compositional models using EnKF. In this paper, we illustrate the
problem using an example in which the updated molar density of CO2 in some regions is observed to take negative
values while molar densities of the remaining components are increased. Standard truncation schemes avoid negative
values of molar densities, but do not address the problem of increased molar densities of other components. The
results can include a spurious increase in reservoir pressure with a subsequent inability to maintain injection. In
this paper, we present a method for constrained EnKF which takes into account the physical constraints on the
plausible values of state variables during the data assimilation. The proposed method can be implemented in two
diﬀerent approaches, both of which convert inequality constraints to a small number of equality constraints. The
ﬁrst approach uses Lagrange multipliers to apply the active constraints. In the second approach, the constraints
are used as virtual observations for calibrating the model parameters within plausible ranges.
The constrained EnKF method is tested on a 2D compositional model and on a highly heterogeneous 3-phase
ﬂow reservoir model. The eﬀect of the constraints on mass conservation are illustrated using a 1D Buckley-Leverett
ﬂow example. Results show that the constrained EnKF is able to enforce the nonnegativity constraints on molar
densities and bound constraints on saturations (all phase saturations must be between 0 and 1), and achieve a
better estimation of reservoir properties than is obtained using only truncation with EnKF. Introduction
The reliability of reservoir models generally increases as more data are included in their construction. Geostatistical methods have traditionally been used to incorporate sparsely located static data (e.g. geological, geophysical,
and well log or core data) into the construction of reservoir geological models (Deutsch and Journel, 1998). These
methods, however, are typically limited by the requirement that the observations are linearly related to the model
variables. On the other hand, dynamic production data (e.g. oil production rate, water production rate, gas oil ratio, and cumulative ﬂuid production rates) are available in abundant quantity but are nonlinearly related to model
variables. Methods for incorporating both static as well as dynamic data in geological reservoir model construction are often termed history matching. Assisted or automatic history matching is formulated as a minimization
problem in which parameters of the model are adjusted to minimize the mismatch between actual observations
and observations predicted from the model. In most cases a regularization term is included in the minimization
function to penalize deviations from plausible values (Tarantola, 1987; Oliver et al., 2001; Li et al., 2003). If the
data are assimilated sequentially in time, the process is called data assimilation.
In recent years, the ensemble Kalman ﬁlter (EnKF) technique has established itself as a viable data assimilation
method for models with large numbers of variables (Evensen, 1994). The EnKF relies on the Monte Carlo approach
to forecast the error statistics and to compute an approximate Kalman gain matrix for updating model and state
variables. In a recent review paper, Aanonsen et al. (2009) described a number of applications where the EnKF has
been used for problems speciﬁc to the petroleum industry. The ﬁrst application of the EnKF within the petroleum
industry was presented by Lorentzen et al. (2001) and the ﬁrst application to reservoir monitoring was presented by 2 SPE 125101 Nævdal et al. (2002). Since its introduction to petroleum engineering ﬁeld, EnKF technique has gained increasing
attention for history matching and continuous reservoir model updating (Nævdal et al., 2005; Haugen et al., 2006;
Bianco et al., 2007; Evensen et al., 2007; Skjervheim et al., 2007; Zhang and Oliver, 2009).
Despite the general success, one problem that can occur when using the EnKF to assimilate reservoir data
is that the updated saturations may take values which are negative or are greater than unity. Gu and Oliver
(2006) showed examples in which this problem occurred and proposed three diﬀerent methods for preventing nonphysical water saturations. Using normal score transformed values of the saturations as state variables ensures
that saturation values can never exceed plausible limits but it did not prevent some nonphysical oscillations in
the values. A normal score transformation of one set of state variables (termed Gaussian anamorphosis) was
reported to give good results in an oceanographic application (Simon and Bertino, 2009). A much diﬀerent type
of reparameterization was applied by Chen et al. (2009) in an attempt to limit problems with the updates of
saturation. Although the method was successful in improving the quality of the updates, it required substantial
additional computation. The most generally successful method are iterative forms of the ensemble Kalman ﬁlter
(Gu and Oliver, 2007; Li and Reynolds, 2009) in which the Kalman gain is used only to update model variables.
Because state variables are updated using the simulator they naturally satisfy the inequality constraints. The cost
of iteration can be substantial, however, so the goal should be to avoid recourse to this solution.
One approach that has sometimes proved successful with the Kalman ﬁlter is explicit incorporation of the
constraints into the analysis step. Equality constraints have been addressed by a number of authors in the context
of the Kalman ﬁlter (Qureshi et al., 1989; De Geeter et al., 1997; Simon and Chia, 2002; Ko and Bitmead, 2007) but
while these approaches could generally be easily extended to the ensemble Kalman ﬁlter, the problem of inequality
constraints is generally more important for data assimilation and certainly more challenging.
The problem of applying inequality constraints to estimates of state variables has also been addressed (Massicotte et al., 1995; Rao et al., 2003; Gupta and Hauser, 2007; Rotea and Lana, 2008). Ungarala et al. (2007)
presented a formulation for incorporating inequality constraints on state variables in an extended Kalman ﬁlter
technique. In their approach, active constraints were identiﬁed by an exhaustive search of all combinations of state
variables. Their approach was successfully applied to a non-linear batch reactor with small number of state variables and limited number of constraints, but direct application to reservoir problems with few hundred thousands
of state variables and constraints is not feasible. Thacker (2007) discussed two diﬀerent approaches for incorporating inequality constraints into a Kalman ﬁlter framework. In his investigation, it is determined that when the
number of state variables is relatively large and correlated, it may not always be important to determine the active
set of constraints precisely. As a result, the extension of the approach to the ensemble Kalman ﬁlter is relatively
In this paper, we describe a new method based on the work by Thacker (2007) for updating reservoir models
using the ensemble Kalman ﬁlter while avoiding the problem of non-physical values for saturations and densities.
This paper includes three test problems designed to validate the performance of the proposed method against the
standard implementation of the ensemble Kalman ﬁlter. The ﬁrst test problem shows the beneﬁts of the constrained
solution in a compositional model for which truncation of negative densities results in signiﬁcant problems with
the mass of ﬂuid in the model. The second test problem is a three-phase black oil example for which the state
variables have both lower and upper bounds. The ﬁnal example is a simple 1D Buckley-Leverett ﬂow example that
investigates the eﬀects of constraints on mass conservation during the update step. Methodology
The ensemble Kalman ﬁlter The EnKF methodology is a two step process. The ﬁrst step is a forecast step in
which the system is evolved forward in time. The second step is a data assimilation step in which the model and
state variables that describe the system are updated to reduce the observation data mismatch.
At a data assimilation time tk , each state vector of the ensemble is updated using the ensemble approximation
to the Kalman gain in the following manner:
yk,j = yk,j + Kk (dobs,k,j − Hk yk,j ), for j = 1, 2, . . . , Ne (1) where, the superscripts f and u stand for forecast and update, respectively; the subscript k stands for the time
index; and H is the measurement operator that maps the the state vector yf to the simulated data. The forecast
state is obtained by running the simulator forward in time up to the time at which observations are available. In
Eq. 1, Kk is the ensemble approximation to the Kalman gain and is computed as:
Kk = Cf HT (Hk Cf HT + CD,k )−1
Y,k k (2) The collection of all ensemble state vectors is denoted by Y and is a matrix of dimensions Ny × Ne . In Eq. 2, SPE 125101 3 Cf denotes the prior covariance matrix for the state variables which can be estimated from the ensemble using
the following expression:
(Yf − Yk )(Yk − Yk )T
Ne − 1 k
where, Yk is the mean of state variables calculated across the ensemble members and is a vector of dimensions Ny .
The state covariance matrix can be very large, but it is never necessary to compute the entire matrix. Instead,
we compute products of the form (Yk − Yk )T HT which are typically much smaller in dimension than the state
covariance matrix. The Constrained ensemble Kalman ﬁlter The kinds of constraints that are of particular interest in this paper
are linear inequality constraints of the form
ymin ≤ yi ≤ ymax
which is of the type satisﬁed by phase saturations and molar fractions, or of the form
yi ≥ ymin
which is of the type satisﬁed by pressures and densities. It is not easy to determine which of the inequality
constraints should be converted to equality constraints and which might be ignored when the problem has many
variables. A relatively easy approach would be to ﬁrst assimilate the data without any constraints (data assimilation
using the standard EnKF) and then check for constraint violations. In this paper, we apply an iterative procedure
for identifying the constraints that are active in a particular data assimilation step. Once the active set is identiﬁed,
equality constraints are applied to the problem in the update step. If constraints are still violated, the procedure is
repeated. When it is determined that a state variable is on the boundary of the feasible region, then the inequality
constraint can be replaced by equality constraints. The system of active constraints can be written in the following
Py = c
where, the vector c contains the constraint values.
Thacker (2007) discussed two diﬀerent approaches for incorporating equality constraints into a Kalman ﬁlter
framework. The ﬁrst approach uses the undetermined Lagrange multipliers and it is shown that, the resulting
formulation has a form very close to the Kalman ﬁlter updating equation. The second approach incorporates
the constraints along with data in an expanded data vector and subsequently, the entire system of equations is
expanded to reﬂect the constraints. Following is a brief description of the two approaches mentioned above.
Lagrange multipliers This approach uses the undetermined Lagrange multipliers (Thacker and Long, 1988) for
incorporating equality constraints into the data assimilation formalism. While the standard Kalman ﬁlter computes
variables that minimize
J(y) = 1
(y − yf )T (Cf )−1 (y − yf ) + (Hy − d)T C−1 (Hy − d)
2 (5) the constrained Kalman ﬁlter will compute variables that minimize
L(x, y) = J(y) − (Py − c)T x. (6) The vector x is a vector of the Lagrange multipliers that enforce the constraints. The values of the Lagrange
multipliers must be determined as a part of the solution.
The constrained minimum solution can be written as a correction to the unconstrained updated minimum
solution yu (Thacker, 2007):
= yj + Cu PT (Pj Cu PT )−1 (cj − Pj yj )
and yj are the constrained and unconstrained updated solutions for the jth ensemble member,
respectively; CY is the posterior covariance matrix; Pj is the coeﬃcient matrix for the constraint equations for the
jth ensemble member; and cj is the vector of constraint values for jth ensemble member.
The constrained solution (Eq. 7) is clearly similar in form to the unconstrained solution from the EnKF (Eq. 1).
The largest diﬀerence is that, while in the standard EnKF the same Kalman gain is used to update each ensemble
member, here each ensemble member may have a diﬀerent set of active constraints, so terms such as Pj , Cu PT ,
Pj Cu PT , cj , and Pj yj must be computed for each ensemble member at each data assimilation time step in order
to compute the constrained solution. 4 SPE 125101 Augmented data This approach also starts by computing the unconstrained solution. Once the unconstrained
solution is obtained, the next step includes checking for any constraint violations that need to be enforced. The
constraint data cj are then appended to the data vector dobs,j . The data-error covariance matrix CD should be
enlarged with zero values on the diagonal corresponding to the constraint errors, and H should be enlarged to
reﬂect the constraints. With these modiﬁcations, an augmented Kalman gain that includes columns corresponding
to data and active constraints can be constructed for computing the constrained updates to the original forecast
yj . The modiﬁed form of Eq. 1 for implementing the constrained EnKF method can be given as:
= yj + Cf HT (HCf HT + CD )−1 (dobs,j − Hyj )
Y (8) ˜
where, HT is the modiﬁed form of the data extraction matrix which extracts the state variables that are subject
to constraint enforcement along with the simulated data; CD is the modiﬁed form of the data-error covariance
˜ obs,j is the modiﬁed vector of the observations appended with the
matrix which also reﬂects the constraints; and d
As we can notice from Eq. 7 and Eq. 8, the primary diﬀerence between the two approaches is in the sequence of
incorporating the equality constraints. In the Lagrangian approach, the constraints are applied after assimilating
data while in the augmented data approach, the constraints are assimilated simultaneously with data. For linear
observation operators and Gaussian probability densities, the two approaches are equivalent. In the examples that
follow, we use the augmented data approach.
Algorithm description for the constrained EnKF Following is a brief generalized algorithm for implementing
the augmented data approach of the constrained EnKF method. This algorithm describes various steps which need
to be carried out in order to enforce the constraints on any type of state variable.
1. At each data assimilation time step, compute the unconstrained solution for each ensemble member using the
standard EnKF methodology as given by Eq. 1.
2. Identify the state variable with the largest degree of constraint violation.
3. Append the constraint to the vector of actual observations and modify the observation operator matrix to
include the constraint relation. Also expand the data-error covariance matrix, CD for the constraint.
4. With the newly expanded matrices H and CD , compute the modiﬁed Kalman gain for the particular ensemble
member under consideration.
5. Update the state vector for that ensemble member using Eq. 8.
6. If the state variables of the updated realization still violate the inequality constraints by a signiﬁcant amount,
repeat steps 2 through 5. During each additional iteration, expand all the matrices and vectors in order to
reﬂect the new constraints while keeping the old constraints from the previous iterations.
7. Repeat steps 2 through 6 for all ensemble members.
8. Repeat steps 1 through 7 for all the data assimilation time steps in the history matching period.
Application to a Compositional Model In this section, we describe the application of constrained EnKF to
a data assimilation problem of CO2 injection into a compositional reservoir model with 12 components, 11 in the
hydrocarbon phase. The number of components is set to be large in order to test the robustness of the method for
diﬃcult problems. Results obtained from the implementation of the standard EnKF will be discussed ﬁrst. These
results illustrate problems that can occur with the application of the standard EnKF when the state variables must
satisfy inequality constraints. Results from the implementation of the constrained EnKF method are presented
Reservoir Model Description The synthetic 2-dimensional (x-z cross-section) compositional reservoir model
has immobile water saturation with no initial gas phase. The 11 components in the oil phase range from C1 to C7+
along with CO2 and N2 . There is one producer and one injector in the ﬁeld. The producer is located on the east
boundary of the ﬁeld and is operated with a bottom hole pressure (BHP) constraint of 1600 psia. The CO2 injector
is located on the west boundary of the ﬁeld with a BHP constraint of 2200 psia. The reservoir model is divided
into a uniform grid of dimensions 50 × 1 × 30. Grid block dimensions are 50 ft × 30 ft × 20 ft. The true porosity
and log permeability are realizations of a correlated anisotropic Gaussian random ﬁeld. The mean porosity is 0.32 SPE 125101 5 with standard deviation of 0.06. The mean log permeability is 4.0 with standard deviation of 1.4. Figs. 1(a) and
1(b) show the porosity and log permeability ﬁelds for the reference case, respectively. Completion locations are
also shown in Fig. 1(b). The compositional suite from a commercial reservoir simulator ECLIPSE (Schlumberger,
2007) is used for the forward modeling of ﬂow and transport.
0 Porosity 5 Z-direction Producer 0.48
-10 0 10 20 30 40 50 60 X-direction
(a) Reference: Porosity (b) Reference: Log Perm Figure 1: The porosity and log permeability ﬁelds for the compositional model (reference). Ensemble Kalman ﬁlter Forty initial realizations of porosity and log permeability are generated using the same
geostatistical parameters as the reference model. None of the realizations are conditioned to static well data.
Measurements used in the data assimilation include oil production rate, producing gas-oil ratio, CO2 injection rate,
gas-phase CO2 mole fraction at the producer, and liquid-phase CO2 mole fraction at the producer. An additional
observation of CO2 molar density at grid block (16,22) is assimilated at the ﬁrst data assimilation time.
The state vector for model updating includes the molar densities of 12 components from the reservoir ﬂuid
system, pressure, porosity and log permeability. Thus, there are a total of 15 variables in each grid block. The
total simulated production history time is 1106 days and production data are assimilated approximately every 100
days. The noise associated with the production data is sampled from a normal distribution with a mean of 0 and
standard deviation which varies by type of observation. Following is a list of the standard deviation values used
for diﬀerent types of production data:
• oil production rate: 10.0 STB/day
• CO2 injection rate: 20.0 lb moles/day
• producing gas-oil ratio: 0.50 MSCF/STB
• producing gas-phase CO2 mole fraction: 0.0005
• producing liquid-phase CO2 mole fraction: 0.00005
Figure 2 shows the distribution of CO2 throughout the reservoir at day 102 before data assimilation. Because
of the relatively large degree of heterogeneity in permeability, the molar densities are highly variable.
Results from standard EnKF At the ﬁrst data assimilation time (day 102), in addition to the observations of
production data, we assume that the CO2 molar density is measured at one location in the reservoir at that time.
Fig. 3 shows the CO2 molar densities for the 22nd layer of the reservoir model before and after data assimilation.
The updated CO2 molar density from the analysis step in EnKF (Fig. 3(b)) is negative in some locations for most
of the realizations.
The standard solution in the EnKF is to truncate state variables that violate constraints to the value at the
boundary of the feasible region. Molar density cannot be negative so all CO2 densities in the model that are
less than 0 would be set equal to 0 at the end of the analysis step. When that approach is taken in this case,
the consequences on reservoir behavior are striking and are most easily understood through investigation of the
behavior of one particular realization. We will examine the eﬀect of the use of the standard EnKF update with
truncation on the reservoir behavior of realization #31.
At the end of the ﬁrst data assimilation step, the values of all state variables are updated and returned to
the reservoir simulator. Fig. 4 shows the updated pressure and molar densities of two of the components. The 6 SPE 125101 (a) Realization 6 (b) Realization 18 (c) Realization 32 Figure 2: Forecast CO2 molar density distributions at t=102 days for various ensemble members. 3
Molar density Molar density 2.0
10 20 30 40 50 2
1 Grid block index 10 20 30 40 50 Grid block index (a) Forecast CO2 molar density (b) Updated CO2 molar density (without truncation) Figure 3: The prior and posterior ensemble of molar density of CO2 for the 22nd layer.
pressure after updating is quite reasonable with high pressure on the left and low pressure on the right near the
producing completions. The molar density of NC4 and C7 are somewhat odd, however, as they show anomalously
high densities to the right of the measurement location. These densities are higher than at any other location in
the reservoir. (a) Pressure: Day 102 (at restart) (b) NC4 molar density: Posterior (c) C7 molar density: Posterior Figure 4: Realization #31 - Pressure and molar density distributions in the reservoir after the ﬁrst data assimilation
The model is restarted using all the static and dynamic variables updated at the ﬁrst data assimilation time
step (day 102). At day 102.5, a region of high pressure can be seen between the injector and the producer.
The magnitude of the maximum pressure in this region is about 2700 psia which is substantially higher than the
maximum allowable bottom hole pressure for the injector (2200 psia). The anomalously high pressure due to the
excess mass is signiﬁcant enough that it results in an inability to inject ﬂuid into the reservoir (Fig. 5(b)). The SPE 125101 7 CO2 injection rate lb moles day regions with higher molar densities of components NC4 and C7 in Figs. 4(b) and 4(c) overlap quite well with the
region of pressure build-up shown in Fig. 5(a). Similar observation of increased molar density is made for all of the
other components from the oil phase with the exception of CO2 . 5000
0 200 400 600 800 1000 1200 Time Days
(a) Day 102.5251 (immediately after restart) (b) CO2 Injection Rate (c) Legend Figure 5: Pressure in Realization #31 at the end of the ﬁrst forecast time step after data assimilation (left) and
the ensemble of injection rates after data assimilation.
All of the components in the oil phase for realization 31, except CO2 , exhibited an anomalous increase in the
molar density after the update step. The increase in mass resulted in a subsequent increase in the pressure which
then resulted in other problems. Truncation of the negative CO2 molar densities contributed to the mass excess
problem by reducing only the amount of negative mass change without aﬀecting the positive changes. In the next
section, we show that the constrained EnKF solves the problem of state variables that violate constraints, without
introducing the problem of excess mass.
Results from constrained EnKF The constrained EnKF uses a diﬀerent Kalman gain matrix and a diﬀerent
observation vector for each ensemble member. Here we describe in detail the application of the constrained EnKF
to the problem of data assimilation for the same realization studied in the previous section. For the ﬁrst data
assimilation time step, the unconstrained solution for realization 31 showed a value of -0.5632 for the updated CO2
molar density at the grid block (27,23) as the largest violation of the positivity constraint on CO2 molar density
value. Fig. 6 shows the results for the forecast, unconstrained and the constrained solutions for the components
CO2 , NC4 , and C7 obtained from iteration 1 for the 23rd layer in the model. The green vertical line indicates
the location of the grid block (27,23) at which the constraint is enforced. At subsequent iterations, the largest
constraint violation occurs in other layers, but there is still a small eﬀect from the application of constraints on
CO2 molar density in other layers on the density in layer 23 as can be seen in Fig. 6.
In the regions where the CO2 molar density takes negative values, the unconstrained EnKF analysis step
gives anomalously high molar densities for all other components in the hydrocarbon phase. Fig. 6 shows that the
consequence of applying the constraints on the molar density of CO2 is to reduce the increase in molar densities of
components NC4 and C7 . The constrained EnKF technique helps to ensure that the magnitudes of the increases
and decreases in molar densities are applied in a consistent manner. This consistency is missing in case of the
simple ad hoc truncation scheme used during the implementation of the standard EnKF. The results showed that
the constrained EnKF method is able to solve the problem of increased molar density (and thereby, additional
mass) of components in a particular region of the reservoir ﬁeld which was responsible for the high pressure build
up in that region ultimately shutting the injector oﬀ.
Fig. 7 shows the CO2 molar density values with maximum degree of constraint violation within the entire
reservoir model for realization 31 at diﬀerent iterations of the constrained EnKF method for some of the data assimilation time steps. With the application of the constrained EnKF method in an iterative manner, the magnitude
of the CO2 molar density with the maximum degree of constraint violation moves in the direction of the enforced
constraint value. SPE 125101 2.5
Molar density 0.4
Molar density Molar density 8 0.3
0.0 0 10 20 30 40 50 0.3
0.0 0 10 Grid block index 20 30 40 50 0 10 20 30 40 Grid block index (b) NC4 (a) CO2 50 Grid block index (c) C7 Figure 6: Layer 23: The forecast, unconstrained (standard EnKF), and constrained (CEnKF) solutions obtained at
the 1st data assimilation time step using iteration 1 for the 31st ensemble member. The black curve stands for the
forecast solution, red curve denotes the unconstrained EnKF update, and the blue curve represents the constrained
EnKF update. 0.00
CO2 molar density CO2 molar density 0.1
0.6 0 2 4 6 iteration index
(a) assimilation time step 1 8 10 0.05
0.20 0 2 4 6 8 10 iteration index
(b) assimilation time step 4 Figure 7: Realization #31: The minimum value of CO2 molar density at diﬀerent iterations of the constrained
EnKF method for some of the data assimilation time steps.
Matching Production Data The ﬁnal updated ensembles of porosity and permeability ﬁelds obtained from
the standard EnKF and the constrained EnKF method are run from the initial time (day 0) up to the end of the
prediction period (day 1735). Fig. 9 and Fig. 10 show the comparison of the history match of diﬀerent types of
production data from the two methods. As the standard EnKF method suﬀered from the issue where most of the
ensemble members had problem with continuous ﬂuid injection during the data assimilation process, it had adverse
eﬀect on updating the static model parameters. This is reﬂected in the history matching results obtained from the
standard EnKF which are shown in Fig. 9. The results for CO2 injection rate obtained from the standard EnKF
method are biased and do not match the historical injection rate data. Similarly, the results for the oil production
rate obtained from the standard EnKF (Fig. 9(a)) show bias in the early production period where most of the
ensemble members produce signiﬁcantly more amount of oil than the observations. The constrained EnKF method
gave signiﬁcantly better history matching results compared to the standard EnKF.
The results for the history match of the gas-phase CO2 mole fraction (WYMF) data are shown in Fig. 10(c)
and Fig. 10(d) for the two methods. The results indicate that the constrained EnKF resulted in a better match
to the CO2 breakthrough time, compared to the standard EnKF. The root mean squared error (RMSE) values for
the production data from the two methods are given in Table 1. For all types of data, results from the constrained
EnKF are better than the results from the standard EnKF.
Application to a Black Oil Model In this section, the results from the implementation of the ensemble Kalman
ﬁlter in its standard form and the implementation of the constrained EnKF method (augmented data approach)
to a data assimilation problem of black oil model are discussed. The reservoir model description is given below.
Reservoir Model Description The synthetic 3-phase, 2-dimensional (x-y areal) black oil reservoir model has
connate water saturation of 0.15 with no initial gas saturation. The initial reservoir pressure is 3600 psia, which is 5000
0 200 400 600 800 1000 1200 CO2 injection rate lb moles day 9 CO2 injection rate lb moles day SPE 125101 0
20 000 0 200 400 Time Days 600 800 1000 1200 Time Days (a) Standard EnKF (b) Constrained EnKF Figure 8: The CO2 injection rate proﬁles for all the ensemble members during data assimilation obtained from the
standard EnKF and the constrained EnKF 1000
Oil Prod Rate STB day Oil Prod Rate STB day 1000
0 0 500 1000 800
0 1500 500 Time Days 0
500 1000 1500 (b) Constrained EnKF: OPR 1500 CO2 Injection rate lb moles day CO2 Injection rate lb moles day (a) Standard EnKF: OPR 0 1000
Time Days Time Days 0
0 500 1000 1500 Time Days (c) Standard EnKF: WCMIR (inj) (d) Constrained EnKF: WCMIR (inj) Figure 9: Comparison of the oil production rate (OPR) and CO2 injector rate (WCMIR) data obtained from
rerunning the ﬁnal updated ensembles obtained from the standard EnKF and the constrained EnKF. The black
lines represent ensemble members, the red dots represent the observations from the reference model, and the vertical
line indicates the end of the history matching period (day 1106). Table 1: Compositional model - RMSE values for diﬀerent types of observations
WXMF Standard EnKF
0.00185 Constrained EnKF
0.00106 10 SPE 125101 Gas Oil Ratio MSCF STB Gas Oil Ratio MSCF STB 140
0 500 1000 0 1500 500 1.0
0.2 500 1000 1500 (b) Constrained EnKF: GOR
Gas phase CO2 mole fraction Gas phase CO2 mole fraction (a) Standard EnKF: GOR 0.0
Time Days Time Days 1500 Time Days (c) Standard EnKF: WYMF 1.0
0.0 0 500 1000 1500 Time Days (d) Constrained EnKF: WYMF Figure 10: Comparison of the gas oil ratio (GOR) and gas-phase CO2 mole fraction (WYMF) data obtained from
rerunning the ﬁnal updated ensembles obtained from the standard EnKF and the constrained EnKF. The black
lines represent ensemble members, the red dots represent the observations from the reference model, and the vertical
line indicates the end of the history matching period (day 1106).
also the bubble point pressure. The reservoir is divided into a uniform cartesian grid of dimensions 50 × 50 × 1.
The grid block dimensions are 30 ft × 30 ft × 20 ft. There are 4 producers and 1 injector in the ﬁeld arranged in a
ﬁve-spot well pattern. The porosity ﬁeld has a mean of 0.25 and a standard deviation of 0.03. The log permeability
ﬁeld is generated with a mean and a standard deviation of 4.8 and 1.4, respectively. An anisotropic Gaussian
variogram model with practical ranges of 30 and 5 gridblocks in the two principal directions was used to generate
the ﬁelds. Fig. 11(a) and 11(b) show the porosity and log permeability ﬁelds for the reference case, respectively.
All the producers are operated with a maximum oil production rate of 1500 STB/day as the primary constraint
and a minimum bottom-hole pressure of 1000 psia as the secondary constraint. The injector, located in the center
of the ﬁeld, is operated with a maximum water injection rate of 8000 STB/day as the primary constraint and a
maximum bottom hole pressure of 6000 psia as the secondary constraint.
Data Assimilation The porosity and log permeability ﬁelds for the reference model and the initial ensemble are
generated using the sequential Gaussian simulation algorithm with a correlation coeﬃcient of 0.75. An ensemble
of 100 realizations is conditioned to values of permeability and porosity at well locations from the reference model.
The state vector includes static model parameters such as porosity and log permeability along with dynamic state
variables such as pressure, water saturation, and gas saturation. The simulated production data, including oil
production rate, producing gas-oil ratio, water production rate, and water injection rate, are used in the data
assimilation. The total simulated production history is 280 days. The production data are assimilated at day 10,
20, 40, and thereafter at every 40 days. The measurement noise associated with the production data is sampled
from a normal distribution with a mean of 0 and a standard deviation which varies by type of observation. Following
is a list of the standard deviation values used for diﬀerent types of production data:
• oil production rate: 4.0 STB/day
• producing gas-oil ratio: 0.003 MSCF/STB
• water production rate: 0.0001 STB/day when the true rate is less than 1.0 STB/day; otherwise 4.0 STB/day
• water injection rate: 5.0 STB/day SPE 125101 11 0 y-direction 10 0.34
0.18 PROD3 20 INJ 30
40 PROD1 0
PROD4 10 y-direction PROD4 20 INJ 30
40 PROD2 50 PROD3 PROD1 8.00
-1.00 PROD2 50 0 10 20 30 40 50 0 10 20 30 40 50 x-direction x-direction (a) Reference: Porosity (b) Reference: Log Perm Figure 11: The porosity and log permeability ﬁelds for the black oil model (reference).
In the implementation of the constrained EnKF, the updated water saturation is constrained to the connate
water saturation value of 0.15 whenever the updated value is smaller than that. The gas saturation is constrained
to a value of zero upon taking a negative value after data assimilation. The water and gas saturations from the
updated state vector (unconstrained solution) are analyzed simultaneously for any constraint violation and the
constraints are enforced together.
Results from Data Assimilation Fig. 12 shows the prior and posterior water saturations at the 3rd data
assimilation time step for realization #68 using the standard EnKF. Fig. 12(b) indicates that the posterior water
saturation from the data assimilation step is less than the connate water saturation (0.15087) in some part of the
reservoir. As the water saturation cannot go lower than the connate water saturation, the truncation approach
used in the implementation of the standard EnKF method would set all such water saturations in the model to the
connate water saturation. Truncation or constraints had also been used at previous assimilations steps. Although
the standard truncation scheme handles the issue of constraint violation, it does so by addressing the state variables
(water saturation) only in the aﬀected region and does not guarantee the consistency between the state variables
in the remaining part of the model.
50 40 50 PROD1 PROD2 40 PROD1 PROD2
0.2 0.5 0.3 0.5 0.6
0.4 0.6 30 0.2 20 30 0.3
0.1 20 0.6
PROD4 0 PROD3 10 20 30 (a) Prior 40 PROD4 50 0 PROD3 10 20 30 40 50 (b) Posterior Figure 12: The prior and posterior water saturations at the 3rd data assimilation time step for realization #68
using the standard EnKF
The application of the constrained EnKF method ensures that the inequality constraints on the water and gas
saturations are honored at each data assimilation time step. Fig. 13 shows the prior and posterior water saturations
for realization #68 using the constrained EnKF. The posterior water saturations from the standard EnKF and the
constrained EnKF are shown in Fig. 12(b) and Fig. 13(b), respectively, and they both match the production data
from the reference model. The posterior water saturations from the constrained EnKF, however, also honor the 12 SPE 125101 inequality constraints.
50 40 50 PROD1 PROD2 40 PROD1 PROD2 0.3
30 20 30 0.4 0.7
20 0.2 0.6 0.3 0.5
0.4 10 10
PROD4 0 PROD3 10 20 30 (a) Prior 40 PROD4 50 0 PROD3 10 20 30 40 50 (b) Posterior Figure 13: The prior and posterior water saturations at the 3rd data assimilation time step for realization #68
using the constrained EnKF Matching Production Data The production data are continually assimilated using the standard EnKF and the
constrained EnKF. Fig. 14 shows the production data match for the oil production rate for producers PROD1,
PROD2, and PROD3 from the application of the two methods. The constrained EnKF performed better in terms
of matching the production data compared to the standard EnKF. Even though the ensemble obtained from the
standard EnKF appears to cover the observations from the reference model, the spread of the ensemble members
around the measurements is greater in case of the standard EnKF than the constrained EnKF. The ensemble
obtained from the constrained EnKF appears to have enough variability in the prediction phase which is important
for assimilating additional production data in the future.
The results of the production data match for the producing gas-oil ratio for some of the producers from the
two methods are shown in Fig. 15. The constrained EnKF method generally performed better than the standard
EnKF in terms of matching the observations of gas-oil ratio for producers PROD2 and PROD3. Similar results
were observed for water production and injection rates. Mass conservation
The most signiﬁcant problem that was observed with the application of the standard EnKF to the compositional
data assimilation problem was that excess mass was added during the update step at the ﬁrst data assimilation
time. In general, mass should not be expected to be conserved in data assimilation or history matching, as the
mass of a component is typically an unknown quantity that must be estimated. On the other hand, when a
substance is injected into the reservoir at a known rate, it would be desirable for that quantity to be conserved in
the updating process. In this section, we will investigate experimentally the consequences of the analysis step on
mass conservation in a simple model.
Consider one-dimensional, constant-rate injection of water into a reservoir of length 50 that is initially at
irreducible water saturation (0.2) and assume that porosity is a correlated random ﬁeld with mean 0.25 and
standard deviation 0.07. The covariance for porosity is nearly Gaussian with a practical range of 25. At some
particular time (560 in dimensionless units), we measure the saturation at the center of the reservoir. Fig. 16(a)
shows the distribution of forecast water saturation from a large number of reservoir realizations at the measurement
time. If the true water saturation at x = 25 is 0.2, then the result of a standard ensemble Kalman ﬁlter analysis
step is as shown in Fig. 16(b). Approximately 10% of the updated realizations contain negative saturations after
updating. Fig. 16(c) shows the same saturation proﬁles after updating using the constrained EnKF. Constraining
the saturation values to lie in the feasible region is not suﬃcient to ensure that the proﬁles are monotonic and the
advantage over truncation is not obvious because in the two-phase incompressible model, the problems with volume
conservation are not so signiﬁcant as in the compositional model. Truncation of negative water saturations results
in a corresponding reduction in the updated oil saturation, which does not occur in the compositional model.
While the total volume will not be eﬀected by truncation in this case, the volume of the individual phases will
be. Fig. 17 compares the volume of moveable water at the end of the forecast with the volume of moveable water
after EnKF with truncation and with the constrained EnKF. The small variability in the volume at the end of the 1000
100 200 300 400 500 1000 500 0
0 100 Time Days 600
300 400 500 Oil production rate STB day Oil production rate STB day 800 200 400 500 100 1000 500 100 Time Days 200 300 400 500 1400
0 100 Time Days (d) Constrained EnKF: PROD1 300 400 500 (c) Standard EnKF: PROD3 1500 0
0 200 Time Days (b) Standard EnKF: PROD2 1000 100 300 1400
0 Time Days (a) Standard EnKF: PROD1 0
0 200 Oil production rate STB day 0
0 1500 Oil production rate STB day 13 Oil production rate STB day Oil production rate STB day SPE 125101 200 300 400 500 Time Days (e) Constrained EnKF: PROD2 (f) Constrained EnKF: PROD3 2.4
0 100 200 300 400 500 Gas oil ratio MSCF STB 4
Gas oil ratio MSCF STB Gas oil ratio MSCF STB Figure 14: Comparison of the oil production rate (OPR) for producers PROD1, PROD2, and PROD3 obtained
from rerunning the ﬁnal updated ensemble from the standard EnKF and the constrained EnKF 3
0 100 Time Days 400 500 2.4 1.8
0 100 400 Time Days (d) Constrained EnKF: PROD1 500 3
0 100 200 300 200 300 400 500 (c) Standard EnKF: PROD3
Gas oil ratio MSCF STB 2.0 200 1.0 Time Days 4 2.2 100 1.5 (b) Standard EnKF: PROD2
Gas oil ratio MSCF STB Gas oil ratio MSCF STB 300 2.0 Time Days (a) Standard EnKF: PROD1 1.0
0 200 2.5 400 Time Days (e) Constrained EnKF: PROD2 500 2.5
0 100 200 300 400 500 Time Days (f) Constrained EnKF: PROD3 Figure 15: Comparison of the gas oil ratio (GOR) for producers PROD1, PROD2, and PROD3 obtained from
rerunning the ﬁnal updated ensemble from the standard EnKF and the constrained EnKF 14 SPE 125101 0.8 0.8 0.8 0.7 0.7
0.6 0.6 0.6
0.4 0.2 0.3 0.3
10 20 30 40 20 30 40 50 50 10 (a) Forecast (b) Standard EnKF 20 30 40 50 (c) Constrained EnKF Figure 16: Water saturation proﬁles
forecast is a result of discretization of the reservoir and the use of grid block center saturation values. The standard
EnKF with truncation shows a bias towards higher water volume because the negative values are truncated. The
constrained EnKF does not remove the bias, but does reduce the problem substantially. 12.5 12.0 11.5 11.0 Forecast EnKF constrained EnKF Figure 17: Distributions of volume of moveable water in the reservoir at the end of several operations: the forecast
step, the standard EnKF update with truncation of saturations, and the constrained EnKF. If mass is explicitly included in the state vector, and all forecast realizations have the same mass, then the mass
in all analyzed realizations will also be identical. The mass of a component in each gridblock must of course be
non-negative, but mass will not necessarily be non-negative after an EnKF update. The consequence of truncating
negative values of mass is that mass will not be conserved in the analysis step. Application of constraints to the
mass variable would necessarily conserve the mass of each component if all realizations initially contained the same
In general, however, the variables that are included in the state vectors are saturations, porosities, molar
densities, or mass fractions. When saturations or molar densities are used in the state vector, mass is not necessarily
conserved even when truncation is not used. When truncation is used, the results are typically much worse. Fig. 18
shows results from the same Buckley-Leverett displacement experiment whose results were shown in Fig. 16. In
this case, however, the molar densities and porosities are included in the state variable instead of saturations
and porosities. One consequence of this choice of parameterization is that unlike saturation, there is no simple
upper limit on the molar density. Note that truncation of the low values of molar density of water (solid curve
in Fig. 18(a)), does not eﬀect the molar density of oil (solid curve in Fig. 18(b)). Also, truncation of the molar
density of water has no eﬀect on the porosity ﬁeld (solid curve in Fig. 18(c)). As a result, there is 9% more total
moles of both components in the realization after truncation than in the realization after the analysis step. If the
constrained EnKF is used, the molar density of water is quite similar to the molar density after truncation (dashed
curve in Fig. 18(a)) while the molar density of oil and the porosity ﬁeld are both reduced substantially (dashed
curves in Fig. 18(b) and (c)). The net eﬀect is that the total number of moles is only 2% less than the forecast, so
the mass is much more nearly conserved, even though the mass (or number of moles) was not used directly in the
state vector. SPE 125101 15 45
10 0.40 8 0.35 6 0.30 4 0.25 40
20 10 20 30 40 50
10 (a) Molar density of water 20 30 (b) Molar density of oil 40 50 10 20 30 40 50 (c) Porosity Figure 18: Comparison of molar densities and porosity for Realization 1, with simple truncation (solid) and with
constrained EnKF (dashed). Conclusions
In this paper we describe a new method for sequential data assimilation that ensures that the updated state
variables satisfy inequality constraints. The method is a constrained version of the ensemble Kalman ﬁlter that is
eﬃcient in that it does not require additional model evaluations for constraint application. The constrained EnKF
method is demonstrated on a compositional model for which state variables have lower limits, and on a threephase black oil model for which saturations are constrained with both upper and lower limits. When the standard
EnKF approach is used for history matching with a compositional reservoir model having the molar densities of the
components as a part of the state vector, it resulted in problems with the updated molar densities. The use of simple
truncation in this case solved the issue of non-physical (negative) molar density values for the reservoir simulator,
but it did not address the issue of the increased molar densities of some of the other components. The ultimate
result was a substantial pressure build-up in the reservoir and subsequent loss of injection. The application of the
constrained EnKF with a constraint on the CO2 molar density reduced the molar densities of other components in
a systematic manner and solved the issues of excess mass and pressure. We also showed that in cases where mass
conservation of a component is important, the constrained EnKF gives superior results to the standard truncation
The proposed constrained EnKF method can be easily implemented with some minor modiﬁcations to the
pre-existing codes for the standard EnKF. The additional computational cost and time for using the constrained
EnKF is signiﬁcantly lower compared to the time required for running the reservoir simulator. Acknowledgment
Support for Hemant Phale was provided by the OU Consortium on Ensemble Methods (OUCEM). The authors
would like to acknowledge the donation of multiple licenses of ECLIPSE by Schlumberger and the computational
resources provided by the OU Supercomputing Center for Education and Research (OSCER). References
Aanonsen, S. I., G. Nævdal, D. S. Oliver, A. C. Reynolds, and B. Vall`s, Ensemble Kalman ﬁlter in reservoir
engineering — a review, SPE Journal, approved for publication, 2009.
Bianco, A., A. Cominelli, L. Dovera, G. Nævdal, and B. Vall`s, History matching and production forecast
uncertainty by means of the ensemble Kalman ﬁlter: A real ﬁeld application (SPE-107161), in EUROPEC/EAGE
Conference and Exhibition, 2007.
Chen, Y., D. S. Oliver, and D. Zhang, Data assimilation for nonlinear problems by ensemble Kalman ﬁlter with
reparameterization, Journal of Petroleum Science and Engineering, 66, 1–14, 2009.
De Geeter, J., H. Van Brussel, and J. De Schutter, A smoothly constrained Kalman ﬁlter, IEEE Transactions
on Pattern Analysis and Machine Intelligence, 19, 1171–1177, 1997.
Deutsch, C. V. and A. G. Journel, GSLIB: Geostatistical Software Library and User’s Guide, second edn., Oxford
University Press, New York, 1998.
Evensen, G., Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods
to forecast error statistics, Journal of Geophysical Research, 99(C5), 10,143–10,162, 1994. 16 SPE 125101 Evensen, G., J. Hove, H. C. Meisingset, E. Reiso, K. S. Seim, and Ø. Espelid, Using the EnKF for assisted history
matching of a North Sea reservoir model (SPE 106184), in Proceedings of the 2007 SPE Reservoir Simulation
Gu, Y. and D. S. Oliver, The ensemble Kalman ﬁlter for continuous updating of reservoir simulation models,
Journal of Energy Resources Technology, 128(1), 79–87, 2006.
Gu, Y. and D. S. Oliver, An iterative ensemble Kalman ﬁlter for multiphase ﬂuid ﬂow data assimilation, SPE
Journal, 12(4), 438–446, 2007.
Gupta, N. and R. Hauser, Kalman ﬁltering with equality and inequality state constraints, Tech. Rep. Report
No. 07/18, Numerical Analysis Group, Oxford Univ. Computing Laboratory, Univ. Oxford, Oxford, UK, 2007.
Haugen, V., L.-J. Natvik, G. Evensen, A. Berg, K. Flornes, and G. Nævdal, History matching using the ensemble
Kalman ﬁlter on a North Sea ﬁeld case (SPE-102430), in SPE Annual Technical Conference and Exhibition, 2006.
Ko, S. and R. R. Bitmead, State estimation for linear systems with state equality constraints, Automatica, 43(8),
Li, G. and A. C. Reynolds, Iterative ensemble Kalman ﬁlters for data assimilation, SPE Journal, approved for
Li, R., A. C. Reynolds, and D. S. Oliver, History matching of three-phase ﬂow production data, SPE Journal,
8(4), 328–340, 2003.
Lorentzen, R. J., K. K. Fjelde, J. Frøyen, A. C. V. M. Lage, G. Nævdal, and E. H. Vefring, Underbalanced
drilling: Real time data interpretation and decision support, in SPE/IADC Drilling Conference, 2001.
Massicotte, D., R. Z. Morawski, and A. Barwicz, Incorporation of a positivity constraint into a Kalman ﬁlterbased algorithm for correction of spectrometric data, IEEE Transactions on Instrumentation and Measurement,
44(1), 2–7, 1995.
Nævdal, G., L. M. Johnsen, S. I. Aanonsen, and E. H. Vefring, Reservoir monitoring and continuous model
updating using ensemble Kalman ﬁlter, SPE Journal, 10(1), 66–74, 2005.
Nævdal, G., T. Mannseth, and E. H. Vefring, Near-well reservoir monitoring through ensemble Kalman ﬁlter:
SPE 75235, in Proceedings of SPE/DOE Improved Oil Recovery Symposium, 2002.
Oliver, D. S., A. C. Reynolds, Z. Bi, and Y. Abacioglu, Integration of production data into reservoir models,
Petroleum Geoscience, 7(SUPP), 65–73, 2001.
Qureshi, A. G., M. P. Res, and B. C. Burnaby, Constrained Kalman ﬁltering for image restoration, in Acoustics,
Speech, and Signal Processing, 1989. ICASSP-89., 1989 International Conference on, pp. 1405–1408, 1989.
Rao, C. V., J. B. Rawlings, and D. Q. Mayne, Constrained state estimation for nonlinear discrete-time systems:
Stability and moving horizon approximations, IEEE Transactions on Automatic Control, 48(2), 246–258, 2003.
Rotea, M. and C. Lana, State estimation with probability constraints, International Journal of Control, 81(6),
Schlumberger, Eclipse Technical and Reference Manual, 2007.
Simon, D. and T. L. Chia, Kalman ﬁltering with state equality constraints, IEEE Transactions on Aerospace
and Electronic Systems, 38(1), 128–136, 2002.
Simon, E. and L. Bertino, Application of the Gaussian anamorphosis to assimilation in a 3-D coupled physicalecosystem model of the North Atlantic with the EnKF: a twin experiment, Ocean Science Discussions, 6(1),
Skjervheim, J.-A., G. Evensen, S. I. Aanonsen, B. O. Ruud, and T. A. Johansen, Incorporating 4D seismic data
in reservoir simulation models using ensemble Kalman ﬁlter, SPE Journal, 12(3), 282–292, 2007.
Tarantola, A., Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation, Elsevier,
Amsterdam, The Netherlands, 1987. SPE 125101 17 Thacker, W. C., Data assimilation with inequality constraints, Ocean Modelling, 16(3-4), 264–276, 2007.
Thacker, W. C. and R. B. Long, Fitting dynamics to data, Journal of Geophysical Research, 93(1), 1227–1240,
Ungarala, S., E. Dolence, and K. Li, Constrained extended Kalman ﬁlter for nonlinear state estimation, in
Proccedings of the 8th International IFAC Symposium on Dynamics and Control of Process Systems, 2007.
Zhang, Y. and D. S. Oliver, History matching using a hierarchical stochastic model with the ensemble Kalman
ﬁlter: A ﬁeld case study, SPE-118879, in Proceedings of the 2009 SPE Reservoir Simulation Symposium, The
Woodlands, February 2–4, 2009. ...
View Full Document
This note was uploaded on 02/06/2012 for the course ECON 101 taught by Professor Tan during the Spring '11 term at The Petroleum Institute.
- Spring '11