Noname manuscript No.
(will be inserted by the editor)
The Ensemble Kalman Filter: Theoretical Formulation and Practical
Implementation
Geir Evensen
Nansen Environmental and Remote Sensing Center, Bergen Norway
The date of receipt and acceptance will be inserted by the editor
Abstract The purpose of this paper is to provide a com-
prehensive presentation and interpretation of the Ensemble
Kalman Filter (EnKF) and its numerical implementation. The
EnKF has a large user group and numerous publications have
discussed applications and theoretical aspects of it. This pa-
per reviews the important results from these studies and also
presents new ideas and alternative interpretations which fur-
ther explain the success of the EnKF. In addition to providing
the theoretical framework needed for using the EnKF, there is
also a focus on the algorithmic formulation and optimal nu-
merical implementation. A program listing is given for some
of the key subroutines. The paper also touches upon specific
issues such as the use of nonlinear measurements, in situ pro-
files of temperature and salinity, and data which are avail-
able with high frequency in time. An ensemble based optimal
interpolation (EnOI) scheme is presented as a cost effective
approach which may serve as an alternative to the EnKF in
some applications.
Key words Data assimilation – Ensemble Kalman Filter
1 Introduction
The Ensemble Kalman Filter has been examined and applied
in a number of studies since it was first introduced by Evensen
(1994b). It has gained popularity because of it’s simple con-
ceptual formulation and relative ease of implementation, e.g.,
it requires no derivation of a tangent linear operator or ad-
joint equations and no integrations backward in time. Fur-
ther, the computational requirements are affordable and com-
parable with other popular sophisticated assimilation meth-
ods such as the representer method by Bennett (1992); Ben-
nett et al. (1993); Bennett and Chua (1994); Bennett et al.
Send offprint requests to: Geir Evensen, Nansen Environmental and
Remote Sensing Center, Edvard Griegsvei 3A, 5059 Solheimsviken,
Norway, e-mail: Geir.Evensen@nersc.no
(1996) and the 4DVAR method which has been much stud-
ied by the meteorological community (see e.g. Talagrand and
Courtier, 1987; Courtier and Talagrand, 1987; Courtier et al.,
1994; Courtier, 1997).
This paper gives a comprehensive presentation of the EnKF,
and it may serve as an EnKF reference document. For a user
of the EnKF it provides citations to hopefully all previous
publications where the EnKF has been examined or used. It
also provides a detailed presentation of the method both in
terms of theoretical aspects and the practical implementation.
For experienced EnKF users it will provide a better under-
standing of the EnKF through the presentation of a new and
alternative interpretation and implementation of the analysis
scheme.
In the next section, an overview is givenof previous works
involving the EnKF. Further, in Section 3, an overview of
the theoretical formulation of the EnKF will be given. There-
after the focus will be on implementation issues starting with
the generation of the initial ensemble in Section 4.1 and the
stochastic integration of the ensemble members in Section 4.2.
The major discussion in this paper relates to the EnKF anal-
ysis scheme which is given in Section 4.3. Section 5 dis-
cusses particular aspects of the numerical implementation.
Appendix A presents an approach for examining the con-
sistency of the EnKF based on comparisons of innovations
and predicted error statistics. In Appendix B an optimal in-
terpolation algorithm is presented. It uses a stationary en-
semble but is otherwise similar to the EnKF, and it can thus
be denoted Ensemble Optimal Interpolation (EnOI). In Ap-
pendix C we have given an algorithm which is currently used
for assimilation of observations of subsurface quantities. In
Appendix D the Ensemble Kalman Smoother (EnKS) is pre-
sented in terms of the terminology developed in this paper.
It is illustrated how the smoother solution can be very ef-
ficiently computed as a reanalysis following the use of the
EnKF. Finally in Appendix E we have reviewed and detailed
the presentation of the algorithm used for the generation of
pseudo random fields.
G. Evensen: The ensemble Kalman Filter...
2 Chronology of ensemble assimilation developments
This section attempts to provide a complete overview of the
developments and applications related to the EnKF. In addi-
tion it also points to other recently proposed ensemble based
methods and some smoother applications.
2.1 Applications of the EnKF
Applications involving the EnKF are numerous and includes
the initial work by Evensen (1994b) and an additional exam-
ple in Evensen (1994a) which showed that the EnKF resolved
the closure problems reported from applications of the Ex-
tended Kalman Filter (EKF).
An application with assimilation of altimeter data for the
Agulhas region was discussed in Evensen and van Leeuwen
(1996) and later in an intercomparison with the Ensemble
Smoother (ES) by van Leeuwen and Evensen (1996).
An example with the Lorenz attractor was given by Evensen
(1997) where it was shown that the EnKF could track the
phase transitions and find a consistent solution with realistic
error estimates even for such a chaotic and nonlinear model.
Burgers et al. (1998) reviewed and clarified some points
related to the perturbation of measurements in the analysis
scheme, and also gave a nice interpretation supporting the
use of the ensemble mean as the best estimate.
Houtekamer and Mitchell (1998) introduced a variant of
the EnKF where two ensembles of model states are integrated
forward in time, and statistics from one ensemble is used to
update the other. The use of two ensembles was motivated
by claiming that this would reduce possible inbreeding in the
analysis. This has, however, lead to some dispute discussed
in the comment by van Leeuwen (1999a) and the reply by
Houtekamer and Mitchell (1999).
Miller et al. (1999) included the EnKF in a comparison
with nonlinear filters and the Extended Kalman Filter, and
concluded that it performed well, but could be beaten by a
nonlinear and more expensive filter in difficult cases where
the ensemble mean is not a good estimator.
Madsen and Ca
˜
nizares (1999) compared the EnKF and
the reduced rank square root implementation of the Extended
Kalman filter with a 2–D storm surge model. This is a weakly
nonlinear problem and good agreement was found between
the EnKF and the extended Kalman filter implementation.
Echevin et al. (2000) studied the EnKF with a coastal ver-
sion of the Princeton Ocean Model and focussed in particular
on the horizontal and vertical structure of multivariate covari-
ance functions from sea surface height. It was concluded that
the EnKF could capture anisotropic covariance functions re-
sulting from the impact of coastlines and coastal dynamics,
and had a particular advantage over simpler methodologies
in such areas.
Evensen and van Leeuwen (2000) rederived the EnKF as
a suboptimal solver for the general Bayesian problem of find-
ing the posterior distribution given densities for the model
prediction and the observations. From this formulation the
general filter could be derived and the EnKF could be shown
to be a suboptimal solver of the general filter where the prior
densities are assumed to be Gaussian distributed.
Hamill and Snyder (2000) constructed a hybrid assimila-
tion scheme by combining 3DVAR and the EnKF. The esti-
mate is computed using the 3DVAR algorithm but the back-
ground covariance is a weighted average of the time evolving
EnKF error covariance and the constant 3DVAR error covari-
ance. A conclusion was that with increasing ensemble size
the best results were found with larger weight on the EnKF
error covariance.
Hamill et al. (2000) report from working groups in a work-
shop on ensemble methods.
Keppenne (2000) implemented the EnKF with a two layer
shallow water model and examined the method in twin ex-
periments assimilating synthetic altimetry data. He focused
on the numerical implementation on parallel computers with
distributed memory and found the approach efficient for such
systems. He also examined the impact of ensemble size and
concluded that realistic solutions could be found using a mod-
est ensemble size.
Mitchell and Houtekamer (2000) introduced an adaptive
formulation of the EnKF where the model error parameter-
ization was updated by incorporating information from the
innovations during the integration.
Park and Kaneko (2000) presented an experiment where
the EnKF was used to assimilate acoustic tomography data
into a barotropic ocean model.
Gronnevik and Evensen (2001) examined the EnKF for
use in fish stock assessment, and also intercompared it with
the Ensemble Smoother (ES) and the Ensemble Kalman Smoother
(EnKS).
Heemink et al. (2001) have been examining different ap-
proaches which combine ideas from RRSQRT filtering and
the EnKF to derive computationally more efficient methods.
Houtekamer and Mitchell (2001) have continued the ex-
amination of the two-ensemble approach and introduced a
technique for computing the global EnKF analysis in the case
with many observations, and also a method for filtering of
eventual long range spurious correlations caused by a lim-
ited ensemble size. As will be seen below the current paper
presents a much more efficient way to compute the global
analysis and also argues against filtering of covariances.
Pham (2001) reexamined the EnKF in an application with
the Lorenz attractor and intercompared results with those ob-
tained from versions of the Singular Evolutive Extended Kalman
(SEEK) filter and a particle filter. Ensembles with very few
members were used and this favoured methods like the SEEK
where the “ensemble” of EOFs is selected to best possible
represent the model attractor.
Verlaan and Heemink (2001) applied the RRSQRT and
EnKF filters in test examples with the purpose of classifying
and defining a measure of the degree of nonlinearity of the
model dynamics. Such an estimate may have an impact on
the choice of assimilation method.
Hansen and Smith (2001) proposed a method for produc-
ing analysis ensembles based on integrated use of the 4DVAR
Submitted to Ocean Dynamics 2 December 18, 2002
G. Evensen: The ensemble Kalman Filter...
method and the EnKF. A probabilistic approach was used and
lead to high numerical cost, but an improved estimate could
be found compared to 4DVAR and the EnKF used separately.
Hamill et al. (2001) examined the impact of ensemble size
on noise in distant covariances. They evaluated the impact
of using an “inflation factor” as introduced by Anderson and
Anderson (1999), and also the use of a Schur product of the
covariance with a correlation function to localize the back-
ground covariances as previously discussed by Houtekamer
and Mitchell (2001). The inflation factor is used to replace
the forecast ensemble according to
ψ
j
= ρ(ψ
j
− ψ) + ψ, (1)
with ρ slightly greater than one (typically 1.01). The purpose
is to account for a slight under representation of variance due
to the use of a small ensemble.
Bishop et al. (2001) used an implementation of the EnKF
in an observation system simulation experiment. Ensemble
predicted error statistics were used to determine the optimal
configuration of future targeted observations. The application
typically looked at a case where additional targeted measure-
ments could be deployed over the next few days and the de-
ployment could be optimized to minimize the forecast errors
in a selected region. The methodology was named Ensem-
ble Transform Kalman Filter and it was further examined by
Majumdar et al. (2001).
Reichle et al. (2002) give a nice discussion of the EnKF
in relation to the optimal representer solution. They find good
convergence of the EnKF toward the representer solution with
the difference being caused by the Gaussian assumptions used
in the EnKF at analysis steps. These are avoided in the rep-
resenter method which solves for the maximum likelihood
smoother estimate.
Evensen (2002) provided an intercomparison and review
of sequential assimilation methods including some simple ex-
amples.
Bertino et al. (2002) applied the EnKF and the Reduced
Rank Square Root (RRSQRT) filter with a model for the Odra
estuary. The two methods were compared and used to assim-
ilate real observations to assess the potential for operational
forecasting in the lagoon. This is a relatively linear model and
the EnKF and the RRSQRT filter provided similar results.
Eknes and Evensen (2002) examined the EnKF with a 1–
D three component marine ecosystem model with focus on
sensitivity to the characteristics of the assimilated measure-
ments and the ensemble size. It was found that the EnKF
could handle strong nonlinearities and instabilities which oc-
cur during the spring bloom.
Allen et al. (2002) takes the Eknes and Evensen (2002)
work one step further by applying the method with a 1–D
version of ERSEM for a site in the Mediterranean Sea. They
showed that even with such a complex model it is possible to
find an improved estimate by assimilating in situ data into the
model.
Haugen and Evensen (2002) applied the EnKF to assimi-
late sea level anomalies and sea surface temperature data into
a version of the Miami Isopycnic Coordinate Ocean Model
(MICOM) by Bleck et al. (1992) for the Indian Ocean. The
paper provided an analysis of regionally dependent covari-
ance functions in the tropics and subtropics and also the mul-
tivariate impact of assimilating satellite observations.
Brusdal et al. (2002) discussed a similar application as
Haugen et al. (2002), but focussed on the North Atlantic.
In addition, this paper included an extensive intercompari-
son of the theoretical background of the EnKF, EnKS and
the SEEK filter, and also compared results from these meth-
ods using the same model and measurements. This paper, to-
gether with Haugen et al. (2002), are the first applications of
the EnKF with a full-blown OGCM in a realistic application,
and have proved the feasibility of the assimilation system for
real oceanographic problems.
Natvik and Evensen (2002a,b) presented the first realis-
tic 3–D application of the EnKF with a marine ecosystem
model. These papers proved the feasibility of assimilating
SeaWiFS ocean colour data to control the evolution of a ma-
rine ecosystem model. In addition several diagnostic methods
were introduced which can be used to examine the statistical
and other properties of the ensemble.
Mitchell et al. (2002) examined the EnKF with a global
atmospheric general circulation model with simulated data
resembling realistic operational observations. They assimi-
lated 80 000 observations a day. The system was examined
with respect to required ensemble size, and the effect of lo-
calization (local analysis at a grid point using only nearby
measurements). It was found that severe localization could
lead to imbalance, but with large enough ratio of influence
for the measurements, this was not a problem and no digital
filtering was required. In the experiments they also included
model errors and demonstrated the importance of this to avoid
filter divergence. This work is a significant step forward and it
shows promising results with respect to using the EnKF with
atmospheric forecast models.
Keppenneand Rienecker (2002) implemented a massively
parallel version of the EnKF with the Poseidon isopycnic co-
ordinate ocean model for the tropical Pacific. They demon-
strated the assimilation of in situ observations and focussed
on the parallelization of the model and analysis scheme for
computers with distributed memory. They also showed that
regionalization of background covariances has negligible im-
pact on the quality of the analysis.
2.2 Other ensemble based filters
The EnKF can also be related to some other sequential filters
such as the Singular Evolutive Extended Kalman (SEEK) fil-
ter by Pham et al. (1998); Brasseur et al. (1999); Carmillet
et al. (2001) (see also Brusdal et al., 2002, for an intercompar-
ison of the SEEK and the EnKF); the Reduced Rank Square
Root (RRSQRT) filter by Verlaan and Heemink (2001); and
the Error Subspace Statistical Estimation (ESSE) filter by
Lermusiaux and Robinson (1999a,b); Lermusiaux (2001) which
can be interpreted as an EnKF where the analysis is computed
in the space spanned by the EOFs of the ensemble.
December 18, 2002 3 Submitted to Ocean Dynamics
G. Evensen: The ensemble Kalman Filter...
Anderson (2001) proposed a method denoted the “En-
semble Adjustment Kalman Filter” where the analysis is com-
puted without adding perturbations to the observations. This
still gives the correct mean of the analyzed ensemble but in
the EnKF it would give too low variance as explained by
Burgers et al. (1998). This is accounted for by deriving a lin-
ear operator which replaces the traditional gain matrix and re-
sults in an updated ensemble which is consistent with theory.
A drawback may be the required inversion of the measure-
ment error covariance when this is nondiagonal. This method
becomes a variantof the square root algorithm used by Bishop
et al. (2001). It is demonstrated that for small ensembles (10–
20 members) the EAKF performs better than the EnKF.
Whitaker and Hamill (2002) proposed another version of
the EnKF where the perturbation of observations are avoided.
The scheme provides a better estimate of the analysis vari-
ance by avoiding the sampling errors of the observation per-
turbations. The scheme was tested for small ensemble sizes
(10–20 members) where it had a clear benefit on the results.
The scheme is based on a redefinition of the Kalman gain
derived from the equation
P
a
e
= (I − KH)P
f
(I − H
T
K
T
) + KRK
T
= (I − KH)P
f
.
(2)
where the term KRK
T
= 0 without perturbations of mea-
surements. A solution of this equation is
K = P
f
H
T
p
HP
f
H
T
+ R
−1
T
×
h
p
HP
f
H
T
+ R +
√
R
i
−1
.
(3)
This is essentially a Monte Carlo implementation of the square
root filter and was named (EnSRF).
2.3 Ensemble smoothers
Some publications have focussed on the extensionof the EnKF
to a smoother. The first formulation was given by van Leeuwen
and Evensen (1996) who introduced the Ensemble Smoother
(ES). This method has later been examined in Evensen (1997)
with the Lorenz attractor; applied with a QG model to find a
steady mean flow by van Leeuwen (1999b) and for the time
dependent problem in van Leeuwen (2001); and for fish stock
assessment by Gronnevik and Evensen (2001). Evensen and
van Leeuwen (2000) re-examined the smoother formulation
and derived a new algorithm with better properties named the
Ensemble Kalman Smoother (EnKS). This method has also
been examined in Gronnevik and Evensen (2001) and Brus-
dal et al. (2002).
2.4 Nonlinear filters and smoothers
Another extension of the EnKF relates to the derivation of an
efficient method for solving the nonlinear filtering problem,
i.e., taking non-Gaussian contributions in the predicted error
statistics into account when computing the analysis. These
are discarded in the EnKF (see Evensen and van Leeuwen,
2000), and a fully nonlinear filter is expected to improve the
results when used with nonlinear dynamical models with multi-
modal behaviour where the predicted error statistics are far
from Gaussian. Implementations of nonlinear filters have been
proposed by Miller et al. (1999), Anderson and Anderson
(1999), Pham (2001) and Miller and Ehret (2002), although
they are still not practical for high dimensional systems due to
large numerical cost. A promising exception is the approach
taken by van Leeuwen (2002) based on importance sampling.
3 Sequential data assimilation
This section gives a brief introduction to sequential data as-
similation methodologies such as the Kalman Filter (KF) and
the Extended Kalman Filter (EKF) and outline the general
theory of the EnKF.
3.1 A variance minimizing analysis
The Kalman Filter is a sequential filter method, which means
that the model is integrated forward in time and whenever
measurements are available these are used to reinitialize the
model before the integration continues. We neglect the time
index and denote a model forecast and analysis as ψ
f
and ψ
a
respectively and the measurements are contained in d. Fur-
ther, the respective covariances for model forecast, analysis
and measurements are denoted P
f
, P
a
and R, and the anal-
ysis equation becomes
ψ
a
= ψ
f
+ P
f
H
T
(HP
f
H
T
+ R)
−1
(d − Hψ
f
), (4)
with the analysis error covariances given as
P
a
= P
f
− P
f
H
T
(HP
f
H
T
+ R)
−1
HP
f
. (5)
Here H is the measurement operator relating the true model
state ψ
t
to the observations d allowing for measurement er-
rors , i.e.
d = Hψ
t
+ . (6)
The reinitialization, ψ
a
, is determined as a weighted lin-
ear combination of the model prediction, ψ
f
, and the mea-
surements, d. The weights are the inverses of the error co-
variances for the model prediction and the measurements, and
the optimal linear-combination becomes the Best Linear Un-
biased Estimator (BLUE).
The error covariances for the measurements, R, need to
be prescribed based on our best knowledge about their accu-
racy and the methodologies used to collect them. The error
covariances for the model prediction is computed by solving
an equation for the time evolution of the error covariance ma-
trix of the model state.
A derivation of these equations can be found in several
publications (see e.g. Burgers et al., 1998).
Submitted to Ocean Dynamics 4 December 18, 2002
G. Evensen: The ensemble Kalman Filter...
3.2 The Kalman Filter
Given a linear dynamical model written on discrete form as
ψ
k+1
= F ψ
k
, (7)
the error covariance equation becomes
P
k+1
= F P
k
F
T
+ Q, (8)
where the matrix Q is the error covariance matrix for the
model errors. The model is assumed to contain errors, e.g.
due to neglected physics and numerical approximations. The
equations (7) and (8) are integrated to produce the forecasts
ψ
f
and P
f
, used in the analysis equations (4) and (5).
3.3 The Extended Kalman Filter
With a nonlinear model
ψ
k+1
= f(ψ
k
), (9)
the error covariance equation would still be (8) but with F
being the tangent linear operator (Jacobian) of f(ψ). Thus,
a linearized and approximate equation is used for the predic-
tion of error statistics in the Extended Kalman Filter (EKF).
A comprehensive discussion of the properties of the EKF
can be found in the literature, but for a convenient summary
which intercompares the EKF with Ensemble Kalman Filter
(EnKF), to be discussed next, see Evensen (2002).
3.4 The Ensemble Kalman Filter
The ensemble Kalman filter as proposed by Evensen (1994b)
and later clearified by Burgers et al. (1998) is now introduced.
We will adapt a three stage presentation starting with the rep-
resentation of error statistics using an ensemble of model
states, then an alternative to the traditional error covariance
equation is proposed for the prediction of error statistics, and
finally a consistent analysis scheme is presented.
3.4.1 Representation of error statistics The error covariance
matrices for the forecasted and the analyzed estimate, P
f
and
P
a
, are in the Kalman filter defined in terms of the true state
as
P
f
= (ψ
f
− ψ
t
)(ψ
f
− ψ
t
)
T
, (10)
P
a
= (ψ
a
− ψ
t
)(ψ
a
− ψ
t
)
T
, (11)
where the overline denotes an expectation value, ψ is the
model state vector at a particular time and the superscripts
f, a, and t represent forecast, analyzed, and true state, respec-
tively. However, since the true state is not known, it is more
convenient to consider ensemble covariance matrices around
the ensemble mean,
ψ,
P
f
' P
f
e
= (ψ
f
− ψ
f
)(ψ
f
− ψ
f
)
T
, (12)
P
a
' P
a
e
= (ψ
a
− ψ
a
)(ψ
a
− ψ
a
)
T
, (13)
where now the overline denote an average over the ensemble.
Thus, we can use an interpretation where the ensemble mean
is the best estimate and the spreading of the ensemble around
the mean is a natural definition of the error in the ensemble
mean.
Since the error covariances as defined in (12) and (13) are
defined as ensemble averages, there will clearly exist infini-
tively many ensembles with an error covariance equal to P
f
e
and P
a
e
. Thus, instead of storing a full covariance matrix, we
can represent the same error statistics using an appropriate
ensemble of model states. Given an error covariance matrix,
an ensemble of finite size will always provide an approxima-
tion to the error covariance matrix. However, when the size
of the ensemble N increases the errors in the representation
will decrease proportional to 1/
√
N.
Suppose now that we have N model states in the ensem-
ble, each of dimension n. Each of these model states can be
represented as a single point in an n-dimensional state space.
All the ensemble members together will constitute a cloud
of points in the state space. Such a cloud of points in the
state space can be approximately described using a proba-
bility density function
φ(ψ) =
dN
N
, (14)
where dN is the number of points in a small unit volume and
N is the total number of points. With knowledge about either
φ or the ensemble representing φ we can calculate whichever
statistical moments (such as mean, covariances etc.) we want
whenever they are needed.
The conclusion so far is that the information contained
by a full probability density function can be represented by
an ensemble of model states.
3.4.2 Prediction of error statistics The EnKF was designed
to resolve two major problems related to the use of the EKF
with nonlinear dynamics in large state spaces. The EKF ap-
plies a closure scheme where third- and higher order mo-
ments in the error covariance equation are discarded. This
linearization has been shown to be invalid in a number of ap-
plications, e.g., Evensen (1992) and Miller et al. (1994). In
fact, the equation is no longer the fundamental equation for
the error evolution when the dynamical model is nonlinear.
In Evensen (1994b) it was shown that a Monte Carlo method
can be used to solve an equation for the time evolution of
the probability density of the model state, as an alternative to
using the approximate error covariance equation in the EKF.
For a nonlinear model where we appreciate that the model
is not perfect and contains model errors, we can write it as a
stochastic differential equation (on continuous form) as
dψ = f (ψ)dt + g(ψ, dq). (15)
This equation states that an increment in time will yield an
increment in ψ, which in addition, is influenced by a random
contribution from the stochastic forcing term, g(ψ, dq), rep-
resenting the model errors. The dq describe a vector Brown-
ian motion process with covariance Qdt. Because the model
December 18, 2002 5 Submitted to Ocean Dynamics