JANUARY 2001 123HOUTEKAMER AND MITCHELL
q 2001 American Meteorological Society
A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation
P. L . H OUTEKAMER AND HERSCHEL L. MITCHELL
Direction de la Recherche en Me´te´orologie, Meteorological Service of Canada, Dorval, Quebec, Canada
(Manuscript received 6 March 2000, in final form 12 June 2000)
ABSTRACT
An ensemble Kalman filter may be considered for the 4D assimilation of atmospheric data. In this paper, an
efficient implementation of the analysis step of the filter is proposed. It employs a Schur (elementwise) product
of the covariances of the background error calculated from the ensemble and a correlation function having local
support to filter the small (and noisy) background-error covariances associated with remote observations. To
solve the Kalman filter equations, the observations are organized into batches that are assimilated sequentially.
For each batch, a Cholesky decomposition method is used to solve the system of linear equations. The ensemble
of background fields is updated at each step of the sequential algorithm and, as more and more batches of
observations are assimilated, evolves to eventually become the ensemble of analysis fields.
A prototype sequential filter has been developed. Experiments are performed with a simulated observational
network consisting of 542 radiosonde and 615 satellite-thickness profiles. Experimental results indicate that the
quality of the analysis is almost independent of the number of batches (except when the ensemble is very small).
This supports the use of a sequential algorithm.
A parallel version of the algorithm is described and used to assimilate over 100 000 observations into a pair
of 50-member ensembles. Its operation count is proportional to the number of observations, the number of
analysis grid points, and the number of ensemble members. In view of the flexibility of the sequential filter and
its encouraging performance on a NEC SX-4 computer, an application with a primitive equations model can
now be envisioned.
1. Introduction
The standard Kalman filter explicitly propagates the
error covariances from one assimilation time to the next.
This expensive computation is approximated in the en-
semble Kalman filter (Evensen 1994) by performing an
ensemble of short-range forecasts. The forecast-error
covariances are calculated directly from the ensemble
when they are needed to assimilate data.
The ensemble Kalman filter technique has been em-
ployed to assimilate data in a number of different con-
texts. For example, Evensen and van Leeuwen (1996)
used it for the assimilation of gridded Geosat data, Ev-
ensen (1997) examined its performance in the context
of the Lorenz equations, and Houtekamer and Mitchell
(1998, hereafter HM) used it to assimilate data into a
quasigeostrophic three-level T21 model. As shown in
HM, the approximation improves as the ensemble size
increases. Furthermore, for linear dynamics, the method
converges to the standard Kalman filter as the number
of ensemble members becomes very large. In recent
Corresponding author address: Dr. P. L. Houtekamer, Direction de
la Recherche en Me´te´orologie, 2121 Route Trans-Canadienne, Dor-
val, PQ H9P 1J3, Canada.
E-mail: Peter.Houtekamer@ec.gc.ca
pilot studies, both Anderson and Anderson (1999) and
Miller et al. (1999) considered generalizations of the
ensemble Kalman filter in which the probability distri-
bution underlying the ensemble is considered to consist
of a sum of distributions rather than of a single Gaussian
distribution. Although good results were obtained, these
nonlinear filtering methods are too computationally de-
manding to be considered for operational atmospheric
data assimilation in the near future.
In an operational application, model error has to be
estimated and accounted for. A way of doing this, in
the context of the ensemble Kalman filter, has recently
been proposed by Mitchell and Houtekamer (2000) in
a follow-up study to HM. Also required for the filter to
be feasible in an operational meteorological environ-
ment is a flexible and computationally efficient analysis
algorithm. The purpose of this paper is to propose and
test such an algorithm.
As will be shown, the algorithm has a number of
desirable properties. (i) It uses ensembles to estimate
flow-dependent forecast-error covariances and does not
require a parameterized multivariate correlation model
of forecast error. (ii) It is independent of the forecast
model, so that any (ensemble of) forecast model(s), pos-
sibly including sophisticated nonlinear parameteriza-
tions, can be used to generate the ensemble of forecast
fields. (iii) It is able to utilize nonconventional obser-

124 VOLUME 129MONTHLY WEATHER REVIEW
vations such as satellite radiances. Such observations
have become an essential component of the available
observation set (Andersson et al. 1994; Derber and Wu
1998; McNally et al. 2000) and their importance is ex-
pected to grow in the future. (iv) It functions reasonably
even if the number of ensemble members is modest.
This is important in an operational context, where one
might be limited to O(100) ensemble members. (v) It
produces spatially smooth analysis increments. Discon-
tinuous analysis increments can be expected to be a
significant source of noise and imbalance (Derber et al.
1991; Cohn et al. 1998). (vi) The proposed algorithm
can efficiently assimilate a large number of observa-
tions. In an operational environment, the number of
available observations may exceed 10
5
(6-h period)
21
.
(vii) The proposed algorithm is suitable for parallel
computation. [Some of these properties (in particular,
ii, iii, and implicitly v and vi) overlap with the design
objectives of the Physical-space Statistical Analysis
System (PSAS) described by Cohn et al. (1998) and are
discussed in more detail in section 2 of that paper.]
Recently Keppenne (2000) implemented an ensemble
Kalman filter for a two-layer spectral T100 shallow-
water model to assimilate 2775 gridded observations
per analysis time. To render the algorithm computa-
tionally feasible, it was run on a parallel computer. That
algorithm would become expensive if it were used to
assimilate a significantly larger number of observations
(say, 10
5
) using a fairly small number of processors
(say, 10), as is the operational context at our center.
In the present paper, we address the main weaknesses
of the HM algorithm: its inability to deal efficiently with
a large number of observations and the possibility of
imbalance in the analyses resulting from the imposition
of a cutoff radius. The proposed algorithm, a sequential
ensemble Kalman filter, is described in the next section
and its performance validated in section 3. A parallel
version of the filter is presented in section 4 and its
computational performance in a multiprocessor envi-
ronment is discussed in section 5. Section 6 consists of
a summary and concluding discussion.
2. The algorithm
The algorithm is based on the one described by Ev-
ensen (1994; see also Burgers et al. 1998) as modified
by HM. Its main elements are first presented in terms
of the properties enumerated in section 1. Only the first
six properties are discussed in this section; the seventh,
suitability for parallel computation, will be addressed
in section 4. The presentation begins with the features
inherited from the earlier algorithm.
a. Inherited features
The algorithm uses a pair of ensembles, as in HM,
to deal with a problem of inbreeding. [See van Leeuwen
(1999) and Houtekamer and Mitchell (1999) for further
discussion of this problem.] This strategy allows rep-
resentative ensembles to be maintained even when the
ensemble size is rather small. Let C
f
and C
a
denote
the forecast (i.e., background or first guess) and analysis
vectors, respectively. Here, these are defined on a grid.
The basic analysis equation is [cf. HM Eq. (10)]
51K
j9
(d
i,j
2 H ), i 5 1,...,N, (1)
af f
CC C
i, j i, j i,j
where all quantities apply at the analysis time and we
assume two N-member ensembles, denoted j 5 1 and
j 5 2. Here, K
j9
represents a Kalman gain, d
i,j
a set of
perturbed observations, and H the forward interpolation
operator from the first guess to the observations. Here,
j9 represents the ensemble that is complementary to en-
semble j, that is, j952 for j 5 1 and j951 for j 5
2. The Kalman gain used for the assimilation of ensem-
ble j is thus computed from the complementary ensem-
ble j9.
The background-error covariances are calculated di-
rectly from the ensembles, as in HM. In particular, the
two terms H
T
and H
T
, which occur in the ex-
ff
P H P
jj
pression for the Kalman gain, are defined as [cf. HM
Eqs. (13)–(15)]
N
1
T
f ff f f
P H [ (C 2 C )(H C 2 H C ) (2)
O
ji,ji,jj
N 2 1
i51
and
N
1
f ffff
H P H [ (H C 2 H C )(H C 2 H C ),
O
j i,j j i,j j
N 2 1
i51
(3)
where
NN
11
f ff f
C 5 C and H C 5 H C .
OO
ji,jj i,j
i51 i51
Equation (2) is for the covariances between C and H C,
while (3) is for the covariance of H C with itself.
While H is linear for the observations used in the
current study, it need not be restricted in this way. Since
H is applied to each background field individually (rath-
er than to the covariance matrix , which summarizes
f
P
j
the ensemble statistics), it is possible to use nonlinear
operators. For example, H could be a radiative transfer
model that yields radiances or a convection scheme that
yields precipitation rates, if radiance or precipitation rate
observations were available for assimilation. However,
it should be noted that the Kalman filter equations have
been derived for linear interpolation operators. If the
nonlinearity is large, the improvement due to the cor-
responding observations may be small.
With an ensemble Kalman filter, the first-guess fields,
, are obtained directly by integrating the complete
f
C
i,j
(nonlinear) forecast model(s). Letting M
i,j
denote the
model used to integrate ensemble member i, j, the prop-
agation of the ensemble from time t to time t 1 1 can
be written as