%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a long AMS abstract for the 2004 Seattle meeting.
\documentclass{AMSabs}
%\documentstyle[preprint,aps,epsf]{revtex4}
%\documentclass[twocolumn]{article}
\usepackage{graphicx,amsmath,AMShelv}
\newcommand{\papernumber}{J5.4}

\title{4D ENSEMBLE KALMAN FILTERING FOR ASSIMILATION \\OF ASYNCHRONOUS OBSERVATIONS}

\author{T. Sauer
 \thanks{{\it Corresponding author address:}
 George Mason University, Fairfax, VA 22030, USA} \\George Mason University, Fairfax, VA 22030\\ B.R. Hunt, J.A. Yorke, A.V. Zimin, E. Ott, E.J. Kostelich\thanks{{\it Current address:} Arizona State University, Tempe, AZ 85287, USA}, I. Szunyogh, G. Gyarmati, E. Kalnay,  D.J. Patil, \\University of Maryland, College Park, MD 20742}

\begin{document}

\maketitle \large
\begin{abstract}
A 4-dimensional ensemble Kalman filter method (4DEnKF), which
adapts ensemble Kalman filtering to the assimilation of
observations that are asynchronous with the analysis cycle, is
discussed. In the ideal case of linear dynamics between
consecutive analyses, the algorithm is equivalent to Kalman
filtering assimilation at each observation time. Tests of 4DEnKF
on the Lorenz 40 variable model are conducted.
\end{abstract}

\maketitle

\section{INTRODUCTION}
In ensemble Kalman filtering, a set of background trajectories is
integrated by the dynamical model, and used to estimate the
background covariance matrix. Numerical experiments have shown that
ensemble Kalman filters (EnKF, e.g., Evensen 1994; Evensen and van Leewen
1996, Houtekamer and Mitchell 1998, 2001; Hamill and Snyder 2000) are efficient ways to carry out data
assimilation from simple models to state-of-the-art operational
numerical prediction models.  The ensemble square-root Kalman
filter approach (Tippett et al.~2002; Bishop et al.~2001; Anderson
2001; Whitaker and Hamill 2002; Ott et al.~2002) has attracted
much recent attention.

A further advantage of the ensemble square root Kalman filter is
that it allows asynchrononous observations to be naturally
assimilated. The Four-Dimensional Ensemble
Kalman Filter (4DEnKF), first proposed in Hunt et al. (2003), 
is a
practical way of unifying the ensemble Kalman filter and the
four-dimensional variational approach. Instead of treating
observations as if they occur only at assimilation times, we can
take observations times into account in a natural way, even if
they are different from the assimilation times. The observational
increments are propagated at intermediate time steps using the
ensemble of background forecasts. This extension of the EnKF to a
4DEnKF can be considered analogous to the extension of the
three-dimensional variational technique (3D-Var) to the four
dimensional variational technique (4D-Var).  The idea  is to infer
the linearized model dynamics from the ensemble instead of the
tangent-linear map, as done in conventional 4D-Var schemes.
Furthermore, in the case of linear dynamics, our technique is
equivalent to instantaneous assimilation  of measured data.

\section{ENSEMBLE KALMAN FILTERS}
To set notation, we recall the EnKF method when the observations
are synchronous with the analysis.  Let
\begin{equation} \label{e1}
\dot{x}_m =G_m(x_1, \ldots, x_M)
\end{equation}
for  $m=1,\ldots, M$ be a continuous dynamical system representing
the background vector field, where
 $x = (x_1, \ldots, x_M)$.
The Ensemble Kalman Filter is designed to track the evolution,
under this dynamical system, of an $M$-dimensional Gaussian
distribution centered at $\overline{x}(t)$ with covariance matrix
$P(t)$.

In the implementation of (Ott et al.~2003; Tippett et al.~2003),
$k+1$ trajectories of (\ref{e1}) are followed starting from
initial conditions $x^{a(1)}, \ldots, x^{a(k+1)}$ over a time
interval $[t_a,t_b]$. Since the system is typically
high-dimensional, assume that $k+1\leq M$. The $k+1$ initial
conditions are chosen so that their sample mean and sample
covariance are $ \overline{x}(t_a)$ and $P(t_a)$, respectively.
After running the system over the time interval, we denote the
trajectory points at the end of the interval by $x^{b(1)}, \ldots,
x^{b(k+1)}$, and compute a new sample mean $\overline{x}^b$ and
sample covariance $P^b$ from these $k+1$ vectors. Define the mean
vector $$\overline{x}^b= \frac{1}{k+1}\sum_{i=1}^{k+1}x^{b(i)}$$
and $$\delta x^{b(i)} = x^{b(i)} - \overline{x}^b.$$ Set  the
matrix
$$X^b = \frac{1}{\sqrt{k}}[\delta x^{b(1)} | \cdots | \delta
x^{b(k+1)}]$$ to get the $M\times M$ covariance matrix
\begin{eqnarray}
 P^b&=& X^b(X^b)^T.
\end{eqnarray}
Since the sum of the columns of $X^b$ is zero, the maximum
possible rank of $P^b$ is $k$.

In the Ensemble Kalman Filter, data assimilation is done using
observations assumed to have been taken at time $t_b$. The
observations are used to replace the dynamics-generated pair
$\overline{x}^b, P^b$ at time $t_b$ with a revised pair
$\overline{x}^a, P^a$ that are used as $\overline{x}(t_a')$ and
$P(t_a')$ on the next time interval $[t_a',t_b']$ where $t_a'
\equiv t_b$.

In the typical case the rank of $P^b$ is $k$. Then the column
space $S$ of $P^b$ is $k$-dimensional, and equals the row space,
since $P^b$ is a symmetric matrix. The orthonormal eigenvectors
$u^{(1)}, \ldots, u^{(k)}$ of $P^b$ that correspond to nonzero
eigenvalues span this space. Since the variation of the ensemble
members occurs in the directions spanning the vector space $S$, we
look there for corrections to $\overline{x}^b$ in the data
analysis step. Set $Q=[u^{(1)} | \cdots | u^{(k)}]$ to be the
$M\times k$ matrix whose columns form a basis of $S$. To represent
$P^b$ in this basis, define the $k\times k$ matrix $\hat{P}^b = Q^TP^bQ$.

The data analysis step for EnSQKF uses observations $(y_1, \ldots,
y_l)$  measured at assimilation time $t_b$ that we assume are linearly related to the dynamical state $x$ by
$y=Hx$, where $H$ is the observation operator. This assumption simplifies the following discussion, but the method can be easily extended to the case of a nonlinear observation operator.
Denote by $R$ the
error covariance matrix of the observations. Define $\hat{H} = HQ$
to restrict the action of $H$ to the subspace $S$. The formula for
the solution to recursive weighted least squares with current
solution $\overline{x}^b$ and error covariance matrix $\hat{P}^b$
is
\begin{eqnarray} \label{e4}
\hat{P}^a&=& \hat{P}^b(I+\hat{H}^TR^{-1}\hat{H}\hat{P}^b)^{-1}\nonumber\\
\Delta \hat{x} &=& \hat{P}^a\hat{H}^TR^{-1}(y-H\overline{x}^b)\nonumber\\
\overline{x}^a&=&\overline{x}^b+Q\Delta \hat{x}
\end{eqnarray}
The corrected most likely solution is $\overline{x}^a$, with error
covariance matrix $\hat{P}^a$.

To finish the step and prepare for a new step on the next time
interval, we must produce a new ensemble of $k+1$ initial
conditions $x^{a(1)},\ldots, x^{a(k+1)}$ that have the analysis
mean $\overline{x}^a$ and analysis covariance matrix $\hat{P}^a$.
This can be done in many ways. One approach (Ott et al.~2002) is
to define the positive square root matrix
\begin{equation} \label{e6}
Y = \{I+(\hat{X}^b)^T(\hat{P}^b)^{-1}
(\hat{P}^a-\hat{P}^b)(\hat{P}^b)^{-1}\hat{X}^b\}^{1/2},
\end{equation}
where $\hat{X}^b = Q^TX^b$. Define the matrix ${X}^a = {X}^bY$ and
$${X}^a = \frac{1}{\sqrt{k}}[\delta x^{a(1)} | \cdots | \delta x^{a(k+1)}].$$
Next define the vectors
$$x^{a(i)} =\delta x^{a(i)}+ \overline{x}^a .$$
It can be checked that
\begin{eqnarray}\label{e7}
\overline{x}^a&=& \frac{1}{k+1}\sum_{i=1}^{k+1}x^{a(i)}\nonumber\\
{P}^a&=& {X}^a({X}^a)^T
\end{eqnarray}
satisfy (\ref{e4}).

 \section{ASSIMILATION OF ASYNCHRONOUS DATA}
 The above description assumes that the
data to be assimilated was observed at the assimilation time
$t_b$. The 4DEnKF method adapts EnKF to handle asynchronous
observations, those that have occurred at non-assimilation times.
The key idea is to mathematically treat the observation as a
slightly modified observation of the current state at the
assimilation time. The method of (Hunt et al. 2003) consists of
using the dynamics contained in the ensemble members to carry this
out. In this way we avoid the need to linearize the original
equations of motion, as is necessary in standard implementations
of 4D-Var.

Notice that Eqs.~(\ref{e4},\ref{e6},\ref{e7}) result in analysis
vectors $x^{a(1)}, \ldots, x^{a(k+1)}$ that lie in the space
spanned by the background ensemble $x^{b(1)}, \ldots, x^{b(k+1)}$.
Consider model states of the form
\begin{equation}\label{e2}
x_b = \sum_{i=1}^{k+1} w_ix^{b(i)}.
\end{equation}
The goal of the analysis is to find the appropriate set of weights
$w_1^{a(j)}, \ldots, w_{k+1}^{a(j)}$ for each analysis vector
$x^{a(j)}$.

Now let $y=h(x)$ be a particular observation made at time $t_c
\neq t_b$. We associate to the state $x_b$ in (\ref{e2}) at time
$t_b$ a corresponding state
\begin{equation} \label{e3}
x_c = \sum_{i=1}^{k+1} w_ix^{c(i)},
\end{equation}
where $x^{c(i)}$ is the state of the $i$th ensemble solution at
time $t_c$. We assign the observation $h(x_c)$ at time $t_c$ to
the state $x_b$ given by (\ref{e2}). Eqn.~(\ref{e3}) was utilized
by (Bishop et al.~2001) and (Majumdar et al.~2002)  to predict the
forecast effects of changes in the analysis error.  Here, we use
this property to propagate the dynamical information within the
analysis time window.

It remains to express the asynchronous observations $h(x_c)$ as
functions of $x_b$, the state at the analysis time. This
functional relationship is needed to apply the standard recursive
least squares equation as in (\ref{e4}). Let
$$E_b=[ x^{b(1)} | \cdots | x^{b(k+1)}]$$
and
$$E_c=[ x^{c(1)} | \cdots | x^{c(k+1)}]$$
be the matrices whose columns are the ensemble members at the
times $t_b$ and $t_c$, respectively. Then (\ref{e2}) and
(\ref{e3}) say that $E_bw=x_b$ and $E_cw=x_c$, respectively, where
$w=[w_1, \ldots, w_{k+1}]^T$. The orthogonal projection to the
column span of $E_b$ is given by the matrix
$E_b(E_b^TE_b)^{-1}E_b^T$, meaning that the coefficients $w$ in
(\ref{e2}) can be defined by $w=(E_b^TE_b)^{-1}E_b^Tx_b$. The
linear combination (\ref{e3}) is
$x_c=E_cw=E_c(E_b^TE_b)^{-1}E_b^Tx_b$. Therefore the observation
$h(x_c)$, expressed as a function of the background state $x_b$ at
the time of assimilation, is
\begin{equation}\label{e8}
h(E_cw)=h(E_c(E_b^{T}E_b)^{-1}E_b^Tx_b).
\end{equation}
The latter expression can be substituted directly into the
ensemble filter equations (\ref{e4}). For example, a set of
observations denoted by the matrix $H$ and time-stamped at $t_c$
can be represented at time $t_b$ by the matrix
$HE_c(E_b^{T}E_b)^{-1}E_b^T$. Therefore the innovation $y-Hx_c$
learned from the observations is treated instead as
$y-HE_c(E_b^{T}E_b)^{-1}E_b^Tx_b$ in the assimilation step.

This technique is equivalent to the computation of the forcing of
the observational increments at the correct time in 4D-Var;
however, it propagates the increments forward or backward in time
without the need for the linear tangent model or its adjoint. We
point out that a more efficient implementation than (\ref{e8}) is
possible using the expression $h(E_cw)$, but using analysis
equations that solve directly for the weight vector $w$, similar
to the Ensemble Transform Kalman Filter (Bishop et al.~2001).

Multiple observations are handled in the same manner. Assume the
observation matrix is $H=(h_1^T|\cdots|h_l^T)^T$, where the
observation row vectors $h_1,\ldots, h_l$ correspond to times
$t_{c_1}, \ldots, t_{c_l}$, respectively. Then the observation
matrix $H$ in (\ref{e4}) is replaced with  the matrix
\begin{equation} \label{e9}
\left(
\begin{array}{c}
h_1E_{c_1}\\
\vdots\\
h_lE_{c_l}
\end{array}\right)
(E_b^TE_b)^{-1}E_b^T.
\end{equation}
In addition, it should be noted that the $t_{c_i}$ can be smaller
or larger than $t_b$, allowing for observations to be used at
their correct observational time even after the nominal analysis
time.
 In the case of linear system dynamics, the 4DEnKF technique is
 equivalent to assimilating data at the time it is observed.

%\begin{figure}[ht]
%\resizebox{2.5in}{!}{\includegraphics{f1}} \caption[]{A solution
%of the Lorenz40 model (solid line) with noisy observations
%(squares) and analysis state from 4DEnKF (circles). Only the %first
%component $x_1$ is shown. Observations are generated at each
%$\Delta t = 0.05$ time units. Assimilation is done every four
%$\Delta t$ steps using (\ref{e4}) with $H$ replaced by
%(\ref{e9}).} \label{fig0}
%\end{figure}

\section{COMPUTER EXPERIMENTS}  The Lorenz model is
a manageable spatio-temporal dynamical model that is useful for
illustrating our results. Consider the vector field defined by
\begin{equation}
\dot{x}_m = (x_{m+1}-x_{m-2})x_{m-1}-x_m+F
\end{equation}
for $m=1,\ldots,M$ and with periodic boundary conditions $x_1 =
x_{M+1}$. Setting the forcing parameter $F=8$ and the system
dimension $M=40$, the attractor of the so-called Lorenz40 system
has information dimension approximately 27.1 (Lorenz, 1998).

We start with a long integration of the model, creating a long
background trajectory $x^*$, to be considered as the ''true''
trajectory. The average root mean square deviation from the mean
is approximately $3.61$ for the true trajectory. We produce
artificial noisy observations at each time interval $\Delta t$ by
adding uncorrelated Gaussian noise with variance $1$ to the true
state at each spatial location. 
%Fig.~\ref{fig0} shows a typical
%trajectory projected to one space dimension.

\begin{figure}[ht]
\resizebox{2.5in}{!}{\includegraphics{f2}}
\caption[]{Root mean square error of proposed 4DEnKF method (circles) compared to standard EnSQKF (diamonds) and EnSQKF with time interpolation (triangles). Variance inflation is set at 6\%. Symbols showing RMSE = 1 actually represent values $\geq 1$. RMSE is averaged over several runs of 40,000 steps.}
\label{fig1}
\end{figure}

Figure \ref{fig1} shows that if we use 4DEnKF, assimilations can be skipped with little loss of accuracy in tracking the system state. The system is advanced in steps of size $\Delta t = .05$, but instead of assimilating the observations at each step, assimilation of past data is done only every every $s$ steps. The resulting root mean square error (RMSE) is plotted as circles in Figure \ref{fig1} as a function of $s$. For $s\leq 6$, it appears that little accuracy is lost. 
This shows the
ability of 4DEnKF to take asynchronous observations into account
without carrying out analysis at each observation step. The fact
that the circles in Fig.~\ref{fig2} stay constant as $s$ increases
verifies this capability for the Lorenz40 example. As mentioned
above, it can be shown analytically that the 4DEnKF method is
equivalent to assimilating at each observation time in the case of
linear background dynamics. This experiment shows that the
property can hold as well for nonlinear dynamics, at least for
small values of $s$.

The RMSE of two other methods are shown in Fig.~\ref{fig1} for comparison.
The diamonds plotted in Fig.~\ref{fig1} are the RMSE found by using EnSQKF, allowing $s$ steps of length $\Delta t$ to elapse between assimilations.  Only those observations occurring at the assimilation time were used for assimilation. 
The triangles refer to time-interpolation of the data since the last assimilation. In this alternative, linear interpolation of individual observations as a function of the ensemble background state evolved by the model is used to create an improved observation $y_\Delta(t_b)$ at the assimilation time. In other words, the innovation at time $t_c$ is added instead at assimilation time $t_b$. For the Lorenz example, where the observations are noisy states, this amounts to replacing the observation at time $t_c$ with $y_\Delta(t_b) \equiv y(t_c)+\overline{x}_b-\overline{x}_c$. Assimilation is done by EnSQKF. The idea behind this technique is widely used in operational 3D-Var systems to assimilate asynchronous observations (e.g., Huang et al.~2002; Benjamin et al.~2003). Our implementation provides somewhat optimistic results for this technique, since our background error covariance matrix is not static (independent of time) and homogeneous (independent of location) as it is assumed in a 3D-Var.  As Figure \ref{fig1} shows, for the latter two methods, the accuracy of the assimilated system state becomes considerably worse compared to 4DEnKF as the steps per assimilation $s$ increases.

\begin{figure}[ht]
\resizebox{2.5in}{!}{\includegraphics{f4}} \caption[]{Root mean
square error of proposed 4DEnKF method (circles) compared to
standard EnSQKF (asterisks), when both pre- and post-assimilation data is used. Six percent variance inflation is
used. Symbols showing RMSE = 1 actually represent values $\geq 1$.
RMSE is averaged over several runs of 40,000 steps.} \label{fig2}
\end{figure}

In another computer simulation we tested the ability of 4DEnKF to
assimilate a combination of asynchronous observations from before
and after the analysis time.
 The Lorenz40 system is advanced in steps of size $\Delta t = .05$,
but instead of assimilating the observations at each step,
analysis is done only every $s$ steps, using observations
from (the integer part of) $s/2$ time steps before and after the
analysis time. If $t_b$ is the analysis time, observations from
times $t_b-(\Delta t)[s/2], \ldots,t_b, \ldots, t_b+(\Delta
t)[s/2]$ were used, where the square bracket notation is used to mean ''integer part''. As observations in this computer simulation,
we used all 40 spatial state variables.

The root mean square error (RMSE) of the analysis trajectory with
respect to the original true trajectory is plotted (open circles) as a function
of steps per assimilation $s$ in Fig.~\ref{fig2}. 
Also, we compare 4DEnKF to the strategy of treating the
observations at times $t_b-(\Delta t)[s/2], \ldots,t_b, \ldots,
t_b+(\Delta t)[s/2]$ as if they occurred at the assimilation time
$t_b$, that is, ignoring the precise timing information as frequently done in operational forecasts (Kalnay, 2003). The RMSE
under this alternative strategy, using EnSQKF without the 4D
improvements we have outlined, is plotted in Fig.~\ref{fig2} as
asterisks. Clearly, a penalty is paid for not taking the time
stamp of the observations into account, as 4DEnKF does.

Variance inflation was used in the experiment described above,
meaning that the analysis covariance matrix was artificially
inflated by adding $\epsilon I$ to $\hat{P}^a$ for small
$\epsilon$. In Figures \ref{fig1} and \ref{fig2}, $\epsilon=0.06$ was used for all methods. Variance inflation helps to compensate for
underestimation of the uncertainty in the background state due to
nonlinearity, limited ensemble size, and model error.

The results in Fig.~\ref{fig1} show that for pre-assimilation data, 4DEnKF is superior to straightforward EnKF as well as an alternative form where observational increments were computed with the background at the observing time, a method also used in operational centers. The additional use of post-assimilation observations as in Fig.~\ref{fig2} decreases the RMSE and extends the range of viable $s$, steps per assimilation, to an even greater degree.

We have also achieved similar results by applying the 4DEnKF
methodology to the Local Ensemble Kalman Filter (LEKF), as
developed in (Ott et al.~2003). The local approach is based on the
hypothesis that assimilation can be done on moderate-size spatial
domains and reassembled. The 4D treatment of the asynchronous
local observations can be exploited in the same way as shown in
this article.

The computational savings possible with the 4DEnKF technique arise
from the ability to improve the use of asynchronous observations
without more frequent assimilations. The extra computational cost
of 4DEnKF is dominated by inverting the $(k+1)\times (k+1)$ matrix
$E_b^TE_b$ in (\ref{e8}), which is comparatively small if the
ensemble size $k+1$ is small compared to the number of state
variables $M$. Moreover, applying this technique in conjunction
with local domains as in LEKF allows $k$ to be greatly  reduced in
comparison with $M$. We will report on a combination of the two
ideas in a future publication.

\section{ACKNOWLEDGEMENTS} This work was supported by a James S.
McDonnell 21st Century Research
Award, by the Office of Naval Research (Physics), by the Army
Research Office, by the NPOESS Integrated Program Office, and by
the National Science Foundation (Grants DMS 0104087, DMS 0208092,
and PHYS 0098632).


%\begin{bibliography}

%\end{bibliography}
\section{REFERENCES}
\vspace{.1in} \noindent Anderson, J. L., 2001: An ensemble
adjustment filter for data assimilation. \emph{Mon. Wea. Rev.}
{\bf 129}, 2884-2903.

\vspace{.1in} \noindent Anderson, J. L., and S. L. Anderson, 1999:
A Monte Carlo implementation of the nonlinear filtering problem to
produce ensemble assimilations and forecasts. \emph{Mon. Wea.
Rev.} {\bf 127}, 2741-2758.

\vspace{.1in} \noindent Benjamin, S. G. et al., 2003: An hourly
assimilation/forecast cycle: The RUC. \emph{Mon. Wea. Rev.}, (in
press).

\vspace{.1in} \noindent Bishop, C. H., B. J. Etherton, and S.
Majumdar, 2001: Adaptive sampling with the Ensemble Transform
Kalman Filter. Part I: Theoretical aspects. \emph{Mon. Wea. Rev.}
{\bf 129}, 420-436.

%\vspace{.1in} \noindent Daley, R., 1991: \emph{Atmospheric data
%analysis.} Cambridge University Press, New York.

\vspace{.1in} \noindent Evensen, G., 1994: Sequential data
assimilation with a nonlinear quasi- \linebreak geostrophic model
using Monte Carlo methods to forecast error statistics. \emph{J.
Geophys. Res.} {\bf 99}(C5), 10 143-10 162.

\vspace{.1in} \noindent Evensen, G., and P. J. van Leeuwen, 1996:
Assimilation of Geosat altimeter data for the Agulhas current
using the ensemble Kalman Filter with a quasi-geostrophic model.
\emph{Mon. Wea. Rev.} {\bf 124}, 85-96.

\vspace{.1in} \noindent Hamill, T. M., and C. Snyder, 2000: A
hybrid ensemble Kalman Filter-3D variational analysis scheme.
\emph{Mon. Wea. Rev.} {\bf 128}, 2905-2919.

\vspace{.1in} \noindent Houtekamer, P. L., and H. L. Mitchell,
1998: Data assimilation using an ensemble Kalman Filter technique.
\emph{Mon. Wea. Rev.} {\bf 126}, 796-811.

\vspace{.1in} \noindent Houtekamer, P. L., and H. L. Mitchell,
2001: A sequential ensemble Kalman Filter for atmospheric data
assimilation. \emph{Mon. Wea. Rev.} {\bf 129}, 123-137.

\vspace{.1in} \noindent Huang, X.-Y., K. S. Morgensen, and X.
Yang, 2002: First-guess at the appropriate time: the HIRLAM
implementation and experiments. Proceedings, HIRLAM Workshop on
Variational Data Assimilation and Remote Sensing, Helsinki,
Finland, 28-43.

\vspace{.1in} \noindent Hunt, B.R.,  E. Kalnay, E.J. Kostelich, E.
Ott, D.J. Patil, T. Sauer, I. Szunyogh, J.A. Yorke, A.V. Zimin,
2003: Four-Dimensional Ensemble Kalman Filtering. Preprint.

\vspace{.1in} \noindent Kalnay, E., 2003: \emph{Atmospheric
modeling, data assimilation, and predictability.} Cambridge
University Press, Cambridge, 341 pp.

%\vspace{.1in} \noindent Lorenc, A., 1986: Analysis methods for
%numerical weather prediction. \emph{Quart. J. Roy. Meteor. Soc.}
%{\bf 112}, 1177-1194.

\vspace{.1in} \noindent Lorenz, E.N. 1998: In Proc. Seminar on
Predictability, Vol. 1 (ECMWF, Reading, Berkshire, UK 1996); E. N.
Lorenz and K.A. Emanuel, 1998: \emph{J. Atmos. Sci.} {\bf 55}, 399.

\vspace{.1in} \noindent Majumdar, S., C.H. Bishop, J. Etherton,
and Z. Toth, 2002: Adaptive sampling with the ensemble transform
Kalman filter. Part II: Field program implementation. \emph{Mon.
Wea. Rev.} {\bf 130}, 1144-1165.

\vspace{.1in} \noindent Ott, E., B. R. Hunt, I. Szunyogh, A.V.
Zimin, E.J. Kostelich, M. Corazza, E. Kalnay, D.J. Patil, J.A.
Yorke, 2002: A local ensemble Kalman filter for atmospheric data
assimilation. Submitted.

\vspace{.1in} \noindent Ott, E., B. R. Hunt, I. Szunyogh, A.V.
Zimin, E.J. Kostelich, M. Corazza, E. Kalnay, D.J. Patil, J.A.
Yorke, 2003: Estimating the state of large spatio-temporal chaotic
systems. Submitted to Phys. Rev.

\vspace{.1in} \noindent Tippett, M.K., J.L. Anderson, C.H. Bishop,
T.M. Hamill, J.S. Whitaker, 2003: Ensemble square-root filters.
Preprint.

\vspace{.1in} \noindent Whitaker, J. S., and T. H. Hamill, 2002:
Ensemble Data Assimilation without perturbed observations.
\emph{Mon. Wea. Rev.} {\bf 130}, 1913-1924.



\end{document}

The Kalman filter is a key feature of data assimilation in
spatio-temporal systems, in particular for numerical weather
forecasting. Its adaptive nature is the cornerstone of efforts to
optimize state estimation on the basis of a physical model and
measured observations of related variables. The Kalman filter is
an optimal solution in the case of linear systems dynamics.

The difficulty of data assimilation in numerical weather
forecasting is tied to the large number of variables in the model
and the nonlinearity of the dynamics being modelled. The method of
ensemble Kalman filters (EnKF) has developed as a means of
attacking both problems. An ensemble of background vector fields
is integrated by the model and used to estimate the current
covariance matrix, a key part of the Kalman computation. Research
has shown that ensemble filters enhance the assimilation of
observations where system nonlinearities are critical. More
recently, the ensemble square root Kalman filter (EnSQKF) has
proved to be especially useful (Tippett 2003; Ott et al. 2003).

In this article we point out a further benefit of the ensemble
square root Kalman filter: it can be easily adapted to naturally
integrate the time course of system variables and observations, in
what we call a Four-Dimensional Ensemble Kalman Filter (4DEnKF).

\bibitem{prl} E. Ott, B. R. Hunt, I. Szunyogh, A.V. Zimin, E.J. Kostelich, M. Corazza, E. Kalnay, D.J. Patil, J.A. Yorke.
Preprint.

\bibitem{Ott}  E. Ott, B. R. Hunt, I. Szunyogh, A.V. Zimin, E.J. Kostelich, M. Corazza, E. Kalnay, D.J. Patil, J.A. Yorke. A local ensemble Kalman filter for atmospheric data assimilation. Preprint.
\bibitem{Tippett} M.K. Tippett, J.L. Anderson, C.H. Bishop, T.M. Hamill, J.S. Whitaker. Ensemble square-root filters. Preprint.

\bibitem{Lorenz} E.N. Lorenz in Proc. Seminar on Predictability, Vol. 1 (ECMWF, Reading, Berkshire, UK 1996); E. N. Lorenz and K.A. Emanuel, J. Atmos. Sci. {\bf 55}, 399 (1998).
