This preview shows page 1. Sign up to view the full content.
Unformatted text preview: An Iterative Ensemble Kalman Filter for
Multiphase Fluid Flow Data Assimilation
Yaqing Gu, SPE, BP; and Dean S. Oliver, SPE, U. of Oklahoma
Summary The dynamical equations for multiphase ﬂow in porous media are
highly nonlinear and the number of variables required to characterize the medium is usually large, often two or more variables per
simulator gridblock. Neither the extended Kalman ﬁlter nor the ensemble Kalman ﬁlter is suitable for assimilating data or for characterizing uncertainty for this type of problem. Although the ensemble Kalman ﬁlter handles the nonlinear dynamics correctly during
the forecast step, it sometimes fails badly in the analysis (or updating) of saturations.
This paper focuses on the use of an iterative ensemble Kalman
ﬁlter for data assimilation in nonlinear problems, especially of the
type related to multiphase ﬂow in porous media. Two issues are key:
(1) iteration to enforce constraints and (2) ensuring that the resulting
ensemble is representative of the conditional pdf (i.e. that the uncertainty quantiﬁcation is correct). The new algorithm is compared
to the ensemble Kalman ﬁlter on several highly nonlinear example
problems, and shown to be superior in the prediction of uncertainty. {Dobs,k−1 , dobs,k } is the collection of all data through time tk . Bayes
theorem relates the probability density for the state variables, fk ,
and model variables, mk , after assimilation of data dobs,k at time tk
to the prior probability density at tk as follows:
p( fk , mk Dobs,k ) ∝ p(dobs,k  fk , mk )p( fk , mk Dobs,k−1 )
∝ p(dobs,k  fk , mk )p( fk mk , Dobs,k−1 )
×p(mk Dobs,k−1 )
∝ p(dobs,k mk )p(mk Dobs,k−1 ) . . . . (1) The equivalence is a result of the fact that the state variables fk can
be computed from the model variables mk . Although the ﬁrst line
of Eq. 1 is fundamental for the traditional Kalman ﬁlter because it
explains how to update the model and state variables directly, the
third line simply points out the possibility of using the Kalman ﬁlter to update only the model variables (and initial conditions), then
computing the state from the model variables. It seems that the uncertainty in initial conditions is relatively small in most petroleum
reservoir engineering applications because the reservoir is typically
Introduction
in a state of static equilibrium at the time of ﬁrst production. ExFor linear problems, the Kalman ﬁlter is optimal for assimilating ceptions include uncertainty in the initial ﬂuid contacts, uncertainty
measurements to continuously update the estimate of state vari in irreducible water saturation (Tjølsen et al. 1994), and the posables. Kalman ﬁlters have occasionally been applied to the problem sibility of tilted initial ﬂuid contacts. If important, these types of
of estimating values of petroleum reservoir variables (Eisenmann uncertainty can be treated by including additional variables in the
et al. 1994; Corser et al. 2000), but they are most appropriate state vector (Evensen et al. 2007; Thulin et al. 2007).
when the problems are characterized by a small number of variables
The posterior pdf from which samples should be drawn (Eq. 1)
and when the variables to be estimated are linearly related to the is the product of two terms. The prior is represented by the enobservations. Most data assimilation problems in petroleum reser semble of model and state vectors. As in the traditional EnKF,
voir engineering are highly nonlinear and are characterized by many for purposes of updating the variables, we approximate the prior
variables, often two or more variables per simulator gridblock.
by a Gaussian whose mean and covariance are estimated from the
The problem of weather forecasting is in many respects simi ensemble. The other term is the likelihood; when both terms are
lar to the problem of predicting future petroleum reservoir perfor Gaussian, the product is Gaussian and the Kalman ﬁlter can be used
mance. The economic impact of inaccurate predictions is substan to compute the updated model and state variables. When the liketial in both cases, as is the difﬁculty of assimilating very large data lihood is not Gaussian (e.g. when the relationship between model
sets and updating very large numerical models. One method that variables and observation variables is nonlinear), the product is not
has been recently developed for assimilating data in weather fore a Gaussian, but sampling from an approximation to the posterior
casting is ensemble Kalman ﬁltering (Evensen 1994; Houtekamer can still be accomplished fairly efﬁciently.
and Mitchell 1998; Anderson and Anderson 1999; Hamill et al.
For a problem in which the relationship between the state vari2000; Houtekamer and Mitchell 2001; Evensen 2003). It has been ables, the model parameters, and the data is linear, both the model
demonstrated to be useful for weather prediction over the North At parameters and the state variables can be updated simultaneously
lantic. The method is now beginning to be applied for data assim using the Kalman ﬁlter. The result is an improved estimate of the
ilation in groundwater hydrology (Reichle et al. 2002; Chen and (nonvarying) model parameters and also an improved estimate of
Zhang 2006) and in petroleum engineering (Nævdal et al. 2002, the current value of the state variables.
2005; Gu and Oliver 2005; Liu and Oliver 2005a; Wen and Chen
For a nonlinear problem, it may be impossible to update the state
2006, 2007; Zafari and Reynolds 2007; Gao et al. 2006; Lorentzen variables to be consistent with the updated model parameters withet al. 2005; Skjervheim et al. 2007; Dong et al. 2006), but the out resolving the nonlinear forward problem to obtain state variapplications to state variables whose density functions are bimodal ables. In many applications of the ensemble Kalman ﬁlter, the obhas proved problematic (Gu and Oliver 2006).
jective is primarily to estimate the current state of the system. For
For applications to nonlinear assimilation problems, it is useful petroleum reservoir applications, however, it is generally important
to think of the ensemble Kalman ﬁlter as a least squares method to estimate not only the current state of the system (the pressures
that obtains an averaged gradient for minimization not from a vari and saturations), but also the values of permeability, porosity, and
ational approach but from an empirical correlation between model fault transmissibility.
variables (Anderson 2003; Zafari et al. 2006). In addition to proModel variables, m, are variables that are uncertain but are not
viding a mean estimate of the variables, a Monte Carlo estimate time varying. These include rock properties such as absolute perof uncertainty can be obtained directly from the variability in the meability and porosity. The estimates of these properties change
ensemble.
as data are assimilated, but the parameter itself should not be interAssume that at time tk , we have an ensemble of samples of state preted to be changing with time. If initial conditions and boundand model vectors from the posterior p( fk , mk Dobs,k ) where Dobs,k = ary conditions are uncertain, they can be included in the list of
model variables. State variables, f , are uncertain, timedependent
variables that deﬁne the state of the system. For petroleum reserCopyright © 2007 Society of Petroleum Engineers
voirs, these could include pressure in each ﬂuid phase, saturations
Original SPE manuscript received for review 11 May 2007. Revised manuscript received 13
June 2007. Paper (SPE 108438) peer approved 9 July 2007. 438 December 2007 SPE Journal of phases, or mass fraction. The state variables are frequently solutions of systems of differential (or difference) equations. If the
model is valid, and the initial conditions and model variables are
known, then it is possible to compute the state variables for any
time. Data, d, are observable quantities related to the state of the
model and indirectly to the model parameters. For petroleum reservoirs, data might include bottomhole pressure (possibly at several
locations in the wellbore), surface ﬂow rates, and amplitude of seismic reﬂection. Theoretical data can be computed from the state and
model variables. Observations always have some level of measurement error or noise associated with them.
Application of the ensemble Kalman ﬁlter to petroleum reservoir ﬂow problems is much simpliﬁed by the introduction of a joint
modelstateobservation vector consisting of model variables, state
variables, and theoretical data (Tarantola 1987; Anderson 2001).
In a typical application of the ensemble Kalman ﬁlter to a twophase
petroleum reservoir ﬂow and transport problem, we might deﬁne an
augmented state vector, Y , of the form,
T
Y = φT , ln kT , PT , Sw ,WORT T , . . . . . . . . . . . . . . . . . . . (2) which has been partitioned into model variables that do not change
with time (porosity, φ, and logpermeability, ln k), state variables
that change substantially (pressure, P, and water saturation, Sw ),
and predictions of observations or theoretical data (producing wateroil ratio, WOR). In our history matching (and some geophysical inverse theory) terminology, we denote the static model variables by
the symbol m, the dynamic state variables at time tk by fk (m), and
the predictions of observations at time tk by gk (m). In this nomenclature, the augmented state vector at time tk is
m
Yk = fk (m) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (3)
gk (m)
Because the primary focus of this paper is on the assimilation step,
we will omit the time subscript, understanding that the functions
f (·) and g(·) are generally functions of time. The relationship between the observations and the true static model variables is
dobs = g(mtrue ) + ε . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (4)
where ε is the measurement error is assumed to be Gaussian and
E[εεT ] = CD . The relationship between the observations and the
true augmented state vector can also be written as
dobs = HYtrue + ε, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (5)
where H can be represented as a matrix whose elements are all ones
or zeroes.
Linear Dynamic System. If the relationships between the model
variables, the state variables, and the theoretical data are linear, then
the augmented state vector can be written as
I
m
Y = Fm = F m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (6)
G
Gm
where F and G are linear operators. If the model variables m are
multivariate Gaussian with autocovariance CM , then the autocovariance of the augmented state vector is
I
F CM
G
⎡
CM
= ⎣FCM
GCM CY = IT FT CM F T
FCM F T
GCM F T GT
⎤
CM GT
T ⎦ . . . . . . . . . . . . . . (7)
FCM G
GCM GT After acquisition of new data dobs , the estimate of the model variables with minimum variance is
m = mpr +CM GT (GCM GT +CD )−1 (dobs − Gmpr ). . . (8)
December 2007 SPE Journal This solution can also be written in terms of the augmented state
vector
Y = Ypr +CY H T (HCY H T +CD )−1 (dobs − HYpr ) . . . . (9)
because HCY H T = GCM GT . Eq. 9 is the analysis or updating step
in the Kalman ﬁlter. It simultaneously updates the model variables,
the state variables, and the estimate of the data. Approximations to
the products HCY H T and CY H T can be efﬁciently computed from
the ensemble of state vectors (Evensen 2003).
Nonlinear Dynamic System. For nonlinear dynamic system, we
can consider a series of linear approximations to the nonlinear functions, f (m) and g(m), by linearizing them at point m . Suppose that
g(m) ≈ g(m ) + G (m − m ) and f (m) ≈ f (m ) + F (m − m ), so
the approximation to the covariance based on linearization at m is
⎡
⎤
CM
CM GT
CM F T
⎣ F CM F CM F T F CM GT ⎦ . . . . . . . . . . . . (10)
CY ≈
G CM G CM F T G CM GT
The vector of model variables that maximizes the conditional
probability density also minimizes the following objective function.
S(m) = 1
T −1
g(m) − d CD g(m) − d
2
1
T −1
+ m − mpr CM m − mpr . . . . . . (11)
2 In this notation, mpr denotes the estimate of m at the end of the
forecast step (before the assimilation or analysis step). When the
number of data is smaller than the number of model variables, the
most appropriate form of the GaussNewton method for ﬁnding the
+ 1 iterative estimate of the model vector that minimizes the objective function in Eq. 11 is
m +1 = mpr −CM GT (CD + G CM GT )−1
× g(m ) − dobs − G (m − mpr ) . . . . (12) When the number of data is larger than the number of model variables, the most appropriate form is
m +1 −1
−1
−1
= mpr − (GTCD G +CM )−1 GTCD × g(m ) − dobs − G (m − mpr ) . . . . (13)
On the other hand, if the problem is sufﬁciently nonlinear that
a reduced step length is required, the GaussNewton formula for
iteration is
m +1 = β mpr + (1 − β )m − β CM GT (CD + G CM GT )−1
× g(m ) − dobs − G (m − mpr ) . . . . (14) where β is an adjustment to the step length whose optimal value
can be determined by standard methods (Dennis and Schnabel 1983).
These formulas have been the basis for most of the GaussNewton
or LevenbergMarquardt approaches for automatic history matching
(Gavalas et al. 1976; Tan and Kalogerakis 1992; Li et al. 2003).
Note that CM in Eqs. 12 and 14 is the model covariance prior to
assimilation of the current data but after assimilation of all data before the current time. It does not change during the GaussNewton
iteration, although the linear approximation to the measurement operator, G, may change with each iteration. As a result, the computation of the product G CM GT is not as straightforward as in
the EnKF where it is only necessary to compute (H∆Y )(H∆Y )T
(Evensen 2003).
Iterative forms of the Kalman ﬁlter are not entirely new. The
iterated Kalman ﬁlter (Jazwinski 1970) is in fact similar in character to what we proposed here, Ensemble Randomized Maximum
Likelihood Filter (EnRML), although without the ensemble. The
standard Kalman ﬁlter attempts to estimate the conditional mean,
but for highly nonlinear problems, the conditional mode is often a 439 more appropriate measure of the central tendency. Bell and Cathey
(1993) discussed the equivalency of the iterated Kalman ﬁlter to the
GaussNewton method for approximating the maximum likelihood
estimate. Also, Zupanski (2005) discusses a maximum likelihood
ensemble ﬁlter (MLEF) that is viewed as a maximum likelihood
approach to the ensemble transform Kalman ﬁlter of Bishop et al.
(2001). In the MLEF method, the maximum likelihood model and
the Hessian are estimated.
Several ad hoc iterative applications of the ensemble Kalman ﬁlter have been proposed in the petroleum engineering literature. Gu
and Oliver (2006) used iteration to reduce the magnitude of the nonlinear effects on the saturation correction. Wen and Chen (2007)
proposed a conforming EnKF method in which the ﬂow simulator
was rerun from a previous step to generate plausible saturation values. Liu and Oliver (2005b) iterated on the Kalman correction to
enforce nonlinear constraints for facies observations with a truncated plurigaussian model. Recently, Zafari et al. (2006) review the
EnKF ﬁlter through the lens of optimization, to derive an iterative
EnKF procedure for nonlinear problems.
This paper includes ﬁve test problems designed to validate various aspects of the proposed iterative ﬁlter, including veriﬁcation
that it gives correct results on linear problems.
Problem 1 veriﬁes that the results with reduced and full step length
are identical to results from the analysis step of the EnKF for
a linear problem.
Problem 2 veriﬁes that the estimates of the mean and the variance
from the new iterative analysis step are much better than the
results from the analysis step in EnKF for a nonlinear problem.
Problem 3 tests the impact of the ensemble size on the estimation
of the “sensitivity matrix” from the ensemble for a nonlinear, 10parameter problem. The ensemble is smaller than
the number of model variables in some cases so the estimate
of G is underdetermined, as it will be in most real cases.
Problem 4 is a linear dynamic problem in which a passive tracer
is injected and the time of breakthrough is measured. Both
the dynamics and the observation operator are linear so the
results can be compared with EnKF, which should be correct.
Problem 5 is a nonlinear dynamic problem on which the methods can be tested, but for which the correct results are not
known. We can, however, show that the state variables (saturation) from the EnRML are physically plausible and match
the data.
Implementation of the EnRML. Let Mpr be the matrix whose
columns consist of the ensemble of model vectors after assimilation of all previous data. There are Ne of these vectors, and hence
Ne columns of Mpr . Denote the vector of means of the prior variables by mpr and the matrix of deviations from the means by ∆Mpr .
¯
The ensemble estimate of the prior model variable covariance (after
T
assimilation of all previous data) is CM = ∆Mpr ∆Mpr /(Ne − 1). One
feature that makes the implementation of the traditional ensemble
Kalman ﬁlter so efﬁcient is that it is never necessary to compute
CM , only the products HCY H T and CY H T . This computation is not
as straightforward in an iterative ﬁlter, because it is important to
maintain the distinction between the model covariance matrix estimate, which should be based on the prior models, and the sensitivity
matrix, which should be based on the current values. At the th iteration, let ∆D represent the deviation of each vector of computed
data from the mean vector of computed data and let ∆M represent
the deviation of each vector of model variables from the current
mean. The ensemble average sensitivity matrix G is the coefﬁcient
matrix relating the changes in model parameters to the changes in
computed data,
∆D = G ∆M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (15)
where ∆M is NM × Ne ; ∆D is ND × Ne ; G is ND × NM . NM is the
number of the model parameters and ND is the number of data. As
440 ∆M is not generally invertible (or even square), we use the singular value decomposition (SVD) to solve the system. Although the
dimension of the model variables is large, the size of the ensemble
is typically fairly small so the effort required to compute the SVD
is reasonable (Golub and van Loan 1989, p. 239).
Except for the iterative aspects, the EnRML procedure is very
similar to the EnKF procedure and we use the ﬁrst step of the recursive process to illustrate the EnRML procedure in the following.
1. Compute the reservoir state variables using the updated model
parameters, mj,0 , from the initial time 0 to the ﬁrst measurement time t1
ψ(mj,0 ,t : 0 → t1 ) j = 1, 2, . . . , Ne . . . . . . . . . (16) where ψ(·) denotes the reservoir simulator, the ﬁrst argument
speciﬁed in the parentheses are model parameters used to
reinitialize the ﬂow equations at initial time 0. The state
variables for the reinitialization are determined by the initial
conditions. The second subscript on m is the time index.
2. Generate an approximation of the sensitivity matrix by “solving” Eq. 15 for G using a pseudoinverse based on singular
value decomposition of ∆M .
3. Apply the GaussNewton update formula in Eq. 14 to im+1
prove estimates of model variables, m j,0 ( j = 1, 2, . . . , Ne ).
+1
4. Evaluate the data mismatch term for both m j,0 and m j,0 S(M0 ) = Ne ∑ g(mj,0 ) − dobs,j,1 T j=1
−1
×CD g(mj,0 ) − dobs,j,1 . . . . . (17)
+1
Note the computation of g(mj,0 ) involves solving the for+1
ward ﬂow equations again from initial time 0 with mj,0
+1
ψ(mj,0 ,t : 0 → t1 ) j = 1, 2, . . . , Ne . . . . . . . . . (18) 5. If S(M0+1 ) < S(M0 ), overwrite m with m
β ; otherwise, keep m and decrease β . +1 and increase 6. Check if the convergence criteria are satisﬁed. If not, go to
Step 1 and iterate the procedure; otherwise, exit the assimilation step and begin the next forecast.
The following criteria are used to determine if the solution has
converged:
• MAX1≤i≤NM ;1≤ j≤Ne  mi,+1 − mi, j < ε1 or
j
• S(M +1 ) − S(M • S(M +1 ) ≤ N
D ) < ε2 S(M ) or or • Iteration exceeds the preset maximum number of iterations,
IMAX .
In this paper, ε1 = 10−5 , ε2 = 10−4 , IMAX = 20 for the linear problem, and IMAX = 6 for the nonlinear problem.
A comprehensive description of the noniterative EnKF procedure, including reﬁnements that can be used to decrease the effects
of small sample size can be found in Evensen (2003).
Test Problems Because the key feature of the EnRML approach is the iteration in
the analysis step, the ﬁrst several examples do not involve integration of dynamical systems. Subsequent tests apply the method to
linear and nonlinear dynamical ﬂow and transport problems.
Problem 1: SingleVariable Linear, Analysis Only. The ensemblebased implementation of the Randomized Maximum Likelihood method
of sampling (Kitanidis 1995; Oliver et al. 1996) should perform
quite well on a linear problem with large ensembles, because both
the variance and the crosscovariance can be estimated accurately. December 2007 SPE Journal TABLE 1—SUMMARY STATISTICS PROBLEM 1 (LINEAR)
Estimate EnKF EnRML (0.5) EnRML (1.0) mean
var
obj
iterations 0.000
0.498
99.97
1 0.000
0.498
99.97
12.2 0.000
0.498
99.97
2 2.2 True
0
0.5
100
NA 2 EnKF EnRML (0.5) EnRML (1.0) 2.04
0.033
1 2.80
0.069
22.7 2.80
0.070
9.8 8 10 6 8 10 1.6 True mean
var
iterations 6 1.8 TABLE 2—SUMMARY STATISTICS FOR PROBLEM 2 (NONLINEAR)
Estimate 4 2.84
0.067
NA (a) EnRML 2 obj = (g(mi,j ) − dobs
Cε
j=1 100 ∑ )2 are computed for each ensemble. The values mean, var, and obj
reported in Table 1 are the averages over 10,000 ensembles. Not
surprisingly, results from the EnKF method are quite good. The
iterated ﬁlter also gives the correct results, at a somewhat higher
price. In particular, note that the iterated ﬁlter (EnRML) with restricted step length (β = 0.5) took an average of 12 iterations to
converge, but obtained results that were identical to those of the
EnKF. This veriﬁes that the iterated update, even with a restricted
step, gives correct results for a linear problem.
Problem 2: SingleVariable Nonlinear, Analysis Only. In this
example, a single measurement is made of a random variable drawn
from a standard normal distribution. The observed value is g(−3),
but the measurement is assumed to be contaminated with Gaussian
noise with mean 0 and variance 0.01. An ensemble of 100 states is
generated to test the methods with relatively small sampling error.
We tested the analysis step from both the EnKF and the EnRML.
Two values of β were tested in the iterative update, simply to verify
consistency. For this nonlinear test problem,
g(m) = m + (m/3)2 .
Prior to acquisition of data, the best estimate is mpr = 0 and the variance in the estimate is σ2 = 1. In this case, the results for EnRML
m
with the full step length (β = 1) and for EnRML with restricted step
length (β = 0.5) are nearly identical, except that the restricted step
length took twice as many iterations to converge (Table 2). Results from the EnRML are much better than results from the EnKF
for this problem. In particular, note that the variance estimate from
EnRML is very close to the correct result, while the estimate from
EnKF is too small by a factor of 2.
Problem 3: Ten Variable Nonlinear, Analysis Only. One important aspect of the iterative updating method is the need to compute
(and update) the sensitivity matrix for the data G when the number
of model variables is larger than the number of ensemble members. December 2007 SPE Journal 4 2.8
2.6
2.4
2.2 (b) EnKF
Mean of Ensemble Means In this simple linear example, a single measurement is made of a
random variable drawn from a standard normal distribution. The
goal is to estimate the value of the variable and the uncertainty
of the estimate. The observed value is 0, but it is assumed to be
contaminated with Gaussian measurement noise with mean 0 and
variance 1. To reduce the effect of sampling error, an ensemble of
100 states is generated to test the performance of both the EnKF
and the EnRML methods. For this linear test problem, the variable itself is observed so g(m) = m, but the observation is noisy, so
dobs = g(mtrue ) + ε.
Prior to acquisition of data, the best estimate of m is 0, and the
uncertainty is completely characterized by the variance, σ2 = 1.
m
After assimilation of the data, the best estimate is still 0, but the
posterior variance is 1/2. For each ensemble, after the analysis step,
the mean, the variance, and the data objective function, 3
2.8
2.6 EnRML 40 2.4
2.2 EnKF 40 2
1.8 Rejectio 1.6
0 2 4
6
8
Model Variable Index 10 (c) Comparison Fig. 1—Mean of ensemble means from 10,000 trials with ensemble size varying from 5 (diamonds), 10 (stars), 20 (squares), 40
(triangles), 80 (dashed). The test problem has a larger number of variables than the ensemble
size, to test the impact of the estimation of G from a small ensemble. The variables form a 1D gaussian random ﬁeld on a uniform
lattice. The prior expectation for the variables is mpr = {0, . . . , 0}
and the prior covariance is (CM )i, j = exp(−3i − j/4). (Note that
the prior variance is uniform, but that the covariance is not.)
A single nonlinear measurement is made of the quantity
g(m) = m + 0.2m2
¯
¯
where m = (m1 + · · · + m10 )/10. We assume that the measurement
¯
is d = g({2, . . . , 2}) = 2.8 and that the standard deviation of the
measurement error is 0.01. (This implies that the mean of the model
variables must be approximately 2.0.)
Means and variances for each of the 10 variables were estimated
using standard implementations of the analysis step in EnKF and
the EnRML. Results from the EnKF and EnRML were compared
with results from an acceptance/rejection algorithm (Ripley 1987)
in which 10 million random samples were proposed from a multigaussian approximation to the posterior pdf. 22,000 of the proposals were accepted.
Ensembles for EnKF and EnRML varied in size from 5 to 40
and the experiment was repeated 10,000 times for each ensemble
size. Fig. 1 shows the means of the ensemble estimates for each
of the variables from the EnKF and EnRML. Note that the results
for the mean do not seem to depend signiﬁcantly on the size of the
ensemble, but results from the EnKF are much different from those
of EnRML, which are nearly identical to results from the rejection
method (Fig. 1c). The estimate of variance from the ensemble after 441 Mean Variance 0.85
0.8
RML 5 0.75
0.7 RML 10 0.65 RML 20 0.6 RML 80 0.55 Rejection
4
2
6
8
Model Variable Index 0 10 i Normalized Objective (a) Mean variance as a function of ensemble size. txk = txk−1 + 60
55
50
45
40
35
30
25
20 30 40 50 60
Ensemble Size 70 80 Fig. 2—Results from 10,000 trials with ensemble size varying
from 5 (diamonds), 10 (stars), 20 (squares), 40 (triangles), 80
(dashed). the EnRML iterative update does depend on the size of the ensemble, but for an ensemble size of 80, the results are nearly identical to
the results from the rejection algorithm (Fig. 2a). The normalized
objective function for an ensemble of realizations mi with perturbed
data dobs,i is the mean of the following over the 10,000 trials:
1 Ne
2 −1
∑ g(mi ) − dobs,i CD + mi − mpr,i
2Ne i=1 T −1
×CM mi − mpr,i . . . . . . . . . . . . . . . (19) The magnitude of the mean normalized objective function initially
decreases rapidly as the size of the ensemble increases, but increasing the size above 40 results in little further improvement for this
problem (Fig. 2b).
Problem 4: Linear Dynamics and Observation Operator. In this
problem, dynamics are included so that the result of iterating from
time 0 can be compared with the standard EnKF without iteration,
and with another iterative ﬁlter that does not require a return to the
initial time. Assume that a passive tracer is injected at one end of a
core and that the times for transport of the tracer to different locations of the core is measured. The ﬂow is single phase; the length
of the core is L; the crosssectional area is A; the permeability is
k(x); the porosity is φ(x); and the viscosity of the ﬂuid is µ. For
simplicity, assume also a consistent set of units, so Darcy’s Law is
q = uA = kA∆p(x)/(µx), where q is the ﬂow rate, u is the superﬁcial
velocity, and p(x) is the pressure drop from the inlet to location
x. The velocity of the tracer front is v = u/φ.
For ﬁxed ﬂow rate q, the arrival time of the tracer at location x is
x A
=
v
q Z x
0 φ(x )dx , . . . . . . . . . . . . . . . . . . . . . . . . . (20) or for uniformly discretized grids,
tx = 442 φi
k−1 i txk = (b) Mean normalized objective function as a function of ensemble size. tx = xk
A x
∑
q i=1+ix and to propagate the system from time 0 (location 0) is 10 S(m) = where ix is the discretized grid index corresponding to the location
x. We can see from Eq. 21 that the travel time of the tracer to the grid
ix depends on the summation of the porosity of the core from the
inlet up to the grid ix . It does not depend on the porosity elsewhere,
or on the permeability or the ﬂuid viscosity.
For this selected problem, the core is discretized into 20 uniform
grids. The arrival times of the tracer at the downstream ends of
grids 4, 7, 12, and 20 are used to update the porosity distribution,
and the tracer concentrations of the core. To propagate the system
from previous assimilation time (i.e. previous measured location) is A x ix
∑ φi , . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (21)
q i=1 A x xk
∑ φi ,
q i=1 where ixk−1 and ixk are the integer grid indices corresponding to the
location xk−1 and xk , respectively.
Results from three methods were compared for this problem:
EnKF, EnRML, and conforming EnKF (Wen and Chen 2007). Conforming EnKF is an iterative form of the ensemble Kalman ﬁlter
in which only the model variables are updated using the Kalman
gain. The reservoir simulator is used to update the state variables,
by restarting the reservoir simulator from the previous assimilation
time with the updated model variables. By using the simulator
to update the state variables, the values are physically plausible.
The EnRML method is somewhat similar, except that the iterative
scheme is developed from the GaussNewton formulation of the
Randomized Maximum Likelihood method, which has been extensively tested by comparison with Markov chain Monte Carlo results
(Liu and Oliver 2003). EnRML does demand greater computation as it requires the state variables to be recomputed from initial
conditions. Differences between the methods were accentuated by
assuming that the ﬁrst measurement of arrival time was highly inaccurate (σε = 100.0), while subsequent measurements were more
accurate (σε = 0.25). Histograms of the computed data after assimilation of the ﬁrst and third observation are shown in Fig. 3. There
is very little difference between the results from the two methods
after the ﬁrst data assimilation (Fig. 3a), but the differences are pronounced (Fig. 3b) after the assimilation of the third data. Residuals
from EnRML at the third assimilation time are consistent with the
assumption of small noise in the observations.
Fig. 4 compares the mean of the ensemble porosity estimates
with 30 members from the three methods at measurement time 4.
The vertical green line indicates the measurement location. The
results from the EnKF and EnRML overlay each other, indicating
that for a problem with linear dynamics and a linear observation
operator, EnRML gives results that are identical to the results from
EnKF.
Problem 5: Nonlinear Dynamics and Observation Operator.
This example deals with data assimilation for 1D, twophase, immiscible ﬂow without capillary pressure. The problem is chosen because the saturation shock results in a bimodal probability density
for saturation that is difﬁcult to handle with the standard Kalman ﬁlters. The test case has 32 grid cells in a 1D grid. Water is injected at
a constant rate in grid 1 and ﬂuid is produced at a constant pressure
from grid 32. Water saturations are measured at the observation
well in grid 16.
An exponential covariance model is used to generate 65 initial
realizations of porosity and log permeability. The practical correlation range (the distance at which the covariance drops to 5% of
the variance value (Journel and Huijbregts 1978) of the covariance
model is approximately 15 grids. The mean and standard deviation of porosity ﬁelds are 0.2 and 0.04, respectively. The mean and
standard deviation of ln k ﬁelds are 5.5 and 0.7. The correlation
coefﬁcient between porosity and ln k is 0.6.
December 2007 SPE Journal TRUE
ConEnKF (a=0.5)
EnRMLF (B=0.5)
EnKF 0.32
0.30 Synthetic data: 24.0553 Mean: 26.1051
STD : 3.5160 Freqency 8 6 4 Porosity Estimates 0.28 EnRMLF(B=0.5) 0.26
0.24
0.22
0.20
0.18
0.16 2 0 2 4 6 8 16 18 20 22 24 26 28 30 Freqency 32 Mean: 25.7721
STD : 3.6493 6 14 16 18 20 Fig. 4—Ensemble mean porosity with 30 ensemble members
from the EnKF, conforming EnKF and EnRML both with half
step length at the end of their iterations at measurement time 4. 4 2 0
16 18 20 22 24 26 28 30 32 34 Computed Tracer Travel Time (s) (a) At measurement time 1
Synthetic data: 70.8425
EnRMLF(B=0.5)
Mean: 70.8453
STD : 0.1693 8 Freqency 12 34 ConEnKF(a=0.5) 8 6 4 2 0
68 70 72 74 76 78 ConEnKF(a=0.5) 8 Freqency 10 Grid 0 Mean: 72.7645
STD : 1.6791 6 4 2 0
68 70 72 74 76 78 In a previous investigation, it was found that the ensemble Kalman
ﬁlter worked fairly well at updating the water saturation ﬁeld when
corrections to the saturation were small, but that the updated values
became nonphysical when the measurements started at a late time
(Gu and Oliver 2006). For this test, water saturation observations at
days 40, 90, 140, 190, 210, 230, 250, and 270 are used to reﬁne the
reservoir models. None of the models have water breakthrough at
the measurement location before 190 days, so there are no changes
to the model or state variables until 190 days. At that time, the
correction is fairly large. The standard deviation of the error contained in the data is 0.5%, a small value that was chosen to make
the problem more nonlinear.
Fig. 5 shows the ensemble of water saturation proﬁles after the
assimilation of data at day 190. The vertical line indicates the observation location (grid 16). The red curve is the true saturation proﬁle
at this time. The multiple black lines are saturation proﬁles from
the ensemble. The EnKF gives saturation proﬁles that are completely implausible; in some cases water saturations exceed 1 − Sor
and in others the saturation proﬁles oscillate between high and low
values (Fig. 5a). Results from the conforming EnKF are plausible
but, even after six iterations, inconsistent with the observed saturation (Fig. 5b). Results from the EnRML are both plausible and
consistent with the observed saturation at the measurement location
(Fig. 5c).
In each of these methods, both the model variables (porosity and
logpermeability) and the state variables (saturation and pressure)
are updated to be consistent with observations. The saturation proﬁles in Fig. 5 show, however, that the errors in saturation in the
neighborhood of the front can be quite large for both EnKF and
conforming EnKF. The errors in the porosity are generally smaller
and more consistent between the methods. We compare the estimates of porosity quantitatively using the rootmeansquare error
(RMSE) in the porosity variables. Computed Tracer Travel Time (s) (b) At measurement time 3 Fig. 3—Histograms of the computed tracer travel times with
30 ensemble members from the Conforming EnKF and EnRML
both with half step length at the end of their iterations at measurement times 1 and 3. December 2007 SPE Journal RMSEk = 1
Ne Ne true
∑ (xk, j − xk )2 . . . . . . . . . . . . . . . . . (22) j=1 Fig. 6 plots the individual gridblock RMSE in the porosity values at
190 days (Fig. 6a) and the maximum value of the gridblock RMSE
at each of the measurement times (Fig. 6b). Although there is no
particular reason that the RMSE in model parameters should always
decrease as data are assimilated, we can make a couple of general
observations. Despite the large differences in saturation estimates
between EnKF and EnRML, the errors in the porosity estimates are
very similar. None of the methods make any changes in the model
until 190 days, at which time all of the methods show a substantial 443 EnRMLF
EnKF
Conforming EnKF 0.11
0.10 Gridblock RMSE for Porosity
Estimates at Day 190 0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0 5 10 15 20 25 30 35 Grid (a) After EnKF correction.
(a) RMSE in porosity at 190 days.
EnRMLF
EnKF
Conforming EnKF Maximum Gridblock RMSE
for Porosity Estimates 0.120 0.115 0.110 0.105 0.100 0.095 0.090
0 50 100 150 200 250 300 Time (Days) (b) After 6 conforming EnKF iterations.
(b) Maximum RMSE in porosity at all assimilation times. Fig. 6—Comparison of porosity values for EnKF, conforming
EnKF, and EnRML. (c) After 6 EnRML iterations. Fig. 5—Assimilation of Sw observation at 190 days. reduction in the RMS error. In this example, the maximum RMSE
in porosity after 200 days is smallest for the EnKF.
Discussion and Conclusions Although the EnKF has generally proven to be effective for history
matching and assimilation of data in multiphase ﬂuid ﬂow problems, there are times when the resulting state estimates or realizations are inconsistent with the data, or are inconsistent with the dynamical equations and physical limits. The solutions that we have 444 proposed resolve the issue of inconsistency in the state variables by
generating the updated (history matched) state variables from the
dynamical equations instead of from the ﬁlter — only the static variables are updated from the ﬁlter. Because of the nonlinearity in the
relationship between model and state variables, and between model
variables (e.g. porosity and permeability) and observable variables
(e.g. water saturation, pressure, rates), it is necessary to iterate to
enforce the measurement constraints. We demonstrated through a
number of examples that the iterative ﬁlter (EnRML) gave the correct results for small linear static problems, and agreed with the
EnKF for a linear dynamical problem involving the ﬂow of tracer in
a porous medium.
We also compared the results on two static nonlinear problems,
and showed that the estimates of the mean and the variance from EnRML were very close to the correct values, while results from EnKF
were not. For the dynamic problem of twophase ﬂow in a porous
medium, with observations of saturation, we showed that EnRML,
EnKF, and conforming EnKF gave similar results for porosity (and
permeability) but that the estimate of the state variable (saturation)
was much better from EnRML.
The new EnRML method is more expensive than EnKF because
of the requirement to rerun the simulation models from the beginning time to update the estimates of saturation and pressure. In
practice (Gu 2006), we have found that it is best to use the EnKF
correction when the corrections to the state variables are small, and
to use the EnRML correction when the changes are large. In this December 2007 SPE Journal way, estimates and realizations of the state variables and the model
variables are consistent with all of the data and the physical constraints.
References Anderson, J. L. 2001. An Ensemble Adjustment Kalman Filter for
Data Assimilation. Monthly Weather Review 129 (12): 2884–
2903. DOI
Anderson, J. L. 2003. A Local Least Squares Framework for Ensemble Filtering. Monthly Weather Review 131 (4): 634–642.
Anderson, J. L. and Anderson, S. L. 1999. A Monte Carlo Implementation of the Nonlinear Filtering Problem to Produce Ensemble Assimilations and Forecasts. Monthly Weather Review
127 (12): 2741–2758.
Bell, B. M. and Cathey, F. W. 1993. The iterated Kalman Filter
Update as a GaussNewton Method. IEEE Transactions on Automatic Control 38 (2): 294–297.
Bishop, C. H., Etherton, B. J., and Majumdar, S. J. 2001. Adaptive
Sampling With the Ensemble Transform Kalman Filter. Part I:
Theoretical Aspects. Mon. Wea. Rev. 129: 420–436.
Chen, Y. and Zhang, D. 2006. Data Assimilation for Transient
Flow in Geologic Formations via Ensemble Kalman Filter. Advances in Water Resources 29 (8): 1107–1122.
Corser, G. P., Harmse, J. E., Corser, B. A., Weiss, M. W., and
Whitﬂow, G. L. 2000. Field Test Results for a RealTime Intelligent Drilling Monitor. Paper SPE 59227 presented at the
IADC/SPE Drilling Conference, New Orleans, 23–25 February.
DOI: 10.2118/59227MS.
Dennis, J. E. and Schnabel, R. B. 1983. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. PrenticeHall, Englewood Cliffs.
Dong, Y., Gu, Y., and Oliver, D. S. 2006. Sequential Assimilation
of 4D Seismic Data for Reservoir Description Using the Ensemble Kalman Filter. Journal of Petroleum Science and Engineering 53 (1–2): 83–99. Gu, Y. and Oliver, D. S. 2005. History Matching of the PUNQS3 Reservoir Model Using the Ensemble Kalman Filter. SPEJ
10 (2): 51–65. SPE89942PA. DOI: 10.2118/89942PA.
Gu, Y. and Oliver, D. S. 2006. The Ensemble Kalman Filter for
Continuous Updating of Reservoir Simulation Models. Journal
of Energy Resources Technology 128 (1): 79–87.
Hamill, T. M., Snyder, C., Baumhefner, D. P., Toth, Z., and Mullen,
S. L. 2000. Ensemble Forecasting in the Short to Medium
Range: Report From a Workshop. Bull. Amer. Meteor. Soc. 81:
2653–2664.
Houtekamer, P. L. and Mitchell, H. L. 1998. Data Assimilation
Using an Ensemble Kalman Filter Technique. Monthly Weather
Review 126 (3): 796–811.
Houtekamer, P. L. and Mitchell, H. L. 2001. A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation. Monthly
Weather Review 129 (1): 123–137.
Jazwinski, A. H. 1970. Stochastic Processes and Filtering Theory.
Academic Press, New York.
Journel, A. and Huijbregts, C. J. 1978. Mining Geostatistics. Academic Press, New York. 600 p.
Kitanidis, P. K. 1995. QuasiLinear Geostatistical Theory for Inversing. Water Resour. Res. 31 (10): 2411–2419.
Li, R., Reynolds, A. C., and Oliver, D. S. 2003. History Matching
of ThreePhase Flow Production Data. SPEJ 8 (4): 328–340.
SPE87336PA. DOI: 10.2118/87336PA.
Liu, N. and Oliver, D. S. 2003. Evaluation of Monte Carlo Methods
for Assessing Uncertainty. SPEJ 8 (2): 188–195. SPE84936PA. DOI: 10.2118/84936PA.
Liu, N. and Oliver, D. S. 2005a. Critical Evaluation of the Ensemble Kalman Filter on History Matching of Geologic Facies.
SPEREE 8 (6): 470–477. SPE 92867PA. DOI: 10.2118/92867PA.
Liu, N. and Oliver, D. S. 2005b. Ensemble Kalman ﬁlter for
Automatic History Matching of Geologic Facies. Journal of
Petroleum Science and Engineering 47 (3–4): 147–161. Eisenmann, P., Gounot, M.T., Juchereau, B., and Whittaker, S. J.,
1994. Improved Rxo Measurements Through SemiActive Focusing. Paper SPE 28437 presented at the SPE Annual Technical Conference and Exhibition, New Orleans, 25–28 September.
DOI: 10.2118/28437MS. Lorentzen, R. J., Nævdal, G., Vàlles, B., Berg, A. M., and AA.
Grimstad 2005. Analysis of the Ensemble Kalman Filter for Estimation of Permeability and Porosity in Reservoir Models. Paper
SPE 96375 presented at the SPE Annual Technical Conference
and Exhibition, Dallas, 9–12 October. DOI: 10.2118/96375MS. Evensen, G. 1994. Sequential Data Assimilation With a Nonlinear
Quasigeostrophic Model Using Monte Carlo Methods to Forecast Error Statistics. Journal of Geophysical Research 99 (C5):
10143–10162. Nævdal, G., Johnsen, L. M., Aanonsen, S. I., and Vefring, E. H.
2005. Reservoir Monitoring and Continuous Model Updating
Using Ensemble Kalman Filter. SPEJ 10 (1): 66–74. SPE84372PA. DOI: 10.2118/84372PA. Evensen, G. 2003. The Ensemble Kalman Filter: Theoretical Formulation and Practical Implementation. Ocean Dynamics 53:
343–367. Nævdal, G., Mannseth, T., and Vefring, E. H. 2002. NearWell
Reservoir Monitoring Through Ensemble Kalman Filter. Paper
SPE 75235 presented at the SPE/DOE Improved Oil Recovery
Symposium, Tulsa, 13–17 April. DOI: 10.2118/75235MS. Evensen, G., Hove, J., Meisingset, H. C., Reiso, E., Seim, K. S.,
and Espelid, O. 2007. Using the EnKF for Assisted History
Matching of a North Sea Reservoir Model. Paper SPE 106184
presented at the SPE Reservoir Simulation Symposium, Houston, 26–28 February. DOI: 10.2118/106184MS.
Gao, G., Zafari, M., and Reynolds, A. C. 2006. Quantifying Uncertainty for the PUNQS3 Problem in a Bayesian Setting with
RML and EnKF. SPEJ 11 (4): 506–515. SPE93324PA. DOI:
10.2118/93324PA.
Gavalas, G. R., Shah, P. C., and Seinfeld, J. H. 1976. Reservoir
History Matching by Bayesian Estimation. SPEJ 16 (6): 337–
350. SPE5740PA. DOI: 10.2118/5740PA.
Golub, G. H. and van Loan, C. F. 1989. Matrix Computations. The
Johns Hopkins University Press, Baltimore, second edition.
Gu, Y. 2006. History Matching Production Data Using the Ensemble Kalman Filter. PhD dissertation, University of Oklahoma.
December 2007 SPE Journal Oliver, D. S., He, N., and Reynolds, A. C. 1996. Conditioning
Permeability Fields to Pressure Data. In European Conference
for the Mathematics of Oil Recovery, V 1–11.
Reichle, R. H., McLaughlin, D. B., and Entekhabi, D. 2002. Hydrologic Data Assimilation With the Ensemble Kalman Filter.
Monthly Weather Review 130 (1): 103–114.
Ripley, B. D. 1987. Stochastic Simulation. John Wiley & Sons,
New York.
Skjervheim, J.A., Evensen, G., Aanonsen, S. I., Ruud, B. O., and
Johansen, T. A. 2007. Incorporating 4D Seismic Data in Reservoir Simulation Models Using Ensemble Kalman Filter. SPEJ
12 (3): 282–292. SPE95789PA. DOI: 10.2118/95789PA.
Tan, T. B. and Kalogerakis, N. 1992. A ThreeDimensional ThreePhase Automatic History Matching Model: Reliability of Parameter Estimates. Journal of Canadian Petroleum Technology
31 (3): 34–41. 445 Tarantola, A. 1987. Inverse Problem Theory: Methods for Data
Fitting and Model Parameter Estimation. Elsevier, Amsterdam,
The Netherlands.
Thulin, K., Li, G., Aanonsen, S. I., and Reynolds, A. C. 2007. Estimation of Initial Fluid Contacts by Assimilation of Production
Data With EnKF. Paper to be presented at the SPE Annual Technical Conference and Exhibition, Anaheim, California, 11–14
November.
Tjølsen, C. B., Damsleth, E., and Bu, T. 1994. The Effect
of Stochastic Relative Permeabilities in Reservoir Simulation.
Journal of Petroleum Science and Engineering 10: 273–290.
Wen, X.H. and Chen, W. H. 2006. Realtime Reservoir Model
Updating Using Ensemble Kalman Filter. SPEJ 11 (4): 431–
442. SPE92991PA. DOI: 10.2118/92991PA.
Wen, X.H. and Chen, W. H. 2007. Some Practical Issues on Realtime Reservoir Model Updating Using Ensemble Kalman Filter.
SPEJ 12 (2): 156–166. SPE111571PA. DOI: 10.2118/111571PA.
Zafari, M., Li, G., and Reynolds, A. C. 2006. Iterative Forms of the
Ensemble Kalman Filter. In Proceedings of the 10th European
Conference on the Mathematics of Oil Recovery — Amsterdam
A030.
Zafari, M. and Reynolds, A. C. 2007. Assessing the Uncertainty
in Reservoir Description and Performance Predictions With the
Ensemble Kalman Filter. SPEJ 12 (3): 378–387. SPE95750PA. DOI: 10.2118/95750PA.
Zupanski, M. 2005. Maximum Likelihood Ensemble Filter: Theoretical Aspects. Monthly Weather Review 133 (6): 1710–1726.
Dean Oliver is the Mewbourne Chair Professor in Petroleum and
Geological Engineering at the U. of Oklahoma, where he has
been a faculty member since 2002. From 1997 to 2002, he was
on the faculty at the University of Tulsa. He has seventeen years
of experience with Chevron prior to joining academia. Oliver
is the Executive Editor of SPEJ. He holds a BS degree in physics
from Harvey Mudd College and a PhD degree in geophysics
from the U. of Washington.
Yaqing Gu is a reservoir engineer and member of Reservoir Management Team of EPTG,
BP America Inc. Email: [email protected] Her working interests include reservoir simulation and characterization as well
as uncertainty analysis. She holds a Ph. D. degree in Petroleum
Engineering from the University of Oklahoma and B. S. degree
in Applied Geophysics from the Guilin Institute of Technology. 446 December 2007 SPE Journal ...
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
 tan

Click to edit the document details