%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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).