\svnidlong
{$HeadURL: $}
{$LastChangedDate: $}
{$LastChangedRevision: $}
{$LastChangedBy: $}

\odachapter{Data assimilation methods available in \oda}

\begin{tabular}{p{4cm}l}
\textbf{Origin:} & CTA memo200802\\
\textbf{Last update:}    & \svnfilemonth-\svnfileyear\\
\end{tabular}

\section{Introduction to Data assimilation}
The terminology of 'data assimilation' originates from the field of meteorology
where in 1950's and 1960's new methods were developed to improve weather
forecasts. Due to rapid development of computers in the 1950's the applications
of large numerical models to weather forecasting became possible. In this early
phase the initial state of the model was estimated directly from the
measurement. The numerical model was then used to produce the forecast. It was
soon recognized that the forecasts could be improved if the initial state was
not only based on the measurements but also on the forecast produced by the
previous model run. The first data assimilation method was direct insertion
which was based on replacement of a model variable by its measured value.
Although it was simple to implement, it suffered from the lack of smoothness in
the assimilation results. Much of the research at the time was concentrated on
finding a proper way to introduce the measurements without introducing
oscillations.

In recent years more and more complex methods are being used. These methods are
either based on the minimization of a criterion or on statistical principles.
This difference divides the data assimilation methods into two classes: the
variational and the sequential methods.

The following data assimilation methods are implemented in \oda :
\begin{itemize}
\item Reduced-rank square-root filter (RRSQRT)
\item Ensemble KF (EnKF)
\item Ensemble Square-Root filter (EnSRF)
\item Complementary Orthogonal subspace Filter For Efficient Ensembles (COFFEE)
\item Steady State KF
\item Particle Filter
\item 3DVar
\end{itemize}

\subsection{Variational methods}
Variational methods aim at adjusting a model solution to all the observations
available over the assimilation period (\cite{Talagrand1997}). In variational
methods, one first defines a scalar function which, for any model solution over
the assimilation interval, measures the ``distance'' or ``misfit'' between that
solution and the available observations. The so-called \emph{cost function}
will typically be a sum of squared differences between the observations and the
corresponding model values. One will then look for the model solution that
minimizes the cost function.

A common use of variational method in meteorology is to obtain an optimal
initial condition for a model forecast. A suitable cost function for this is
the following:
\begin{equation}\label{eq.costfunction1}
J(x_o) = \sum_{k=1}^N (y^o(k)-Hx(k))R^{-1}(y^o(k)-Hx(k))'
\end{equation}
where $x_o$ is the initial value of the variable to be determined, $x(k)$ the
variable at time $t_k$, $H$ the observation operator, $y^o(k)$ the observation
at time $t_k$ and $R$ the representation covariance. The optimal initial
condition, $x_o$, is the one that minimizes $J$.

If a suitable initial state $x_o$ has been obtained, the analysed states are
formed with a model forecast:
\begin{eqnarray}
x^a_o = x_o \\
x^a(k) = M[x^a(k-1)]
\end{eqnarray}
for $k=1, \cdots, N$. The final analysed state $x^a_N$ is optimised given data
from spatial different locations and from different times in the interval
$(t_o,t_N)$. This approach is then referred to as \emph{4D-VAR}.
%Data from the period before $t_o$ is used implicitly if it was used to form the background state $x^b$.

The minimization of the cost function is often based on quasi-Newton methods.
These methods require computation of the gradient of the cost function. In most
situations, it is impossible to establish explicit analytical expressions for
the gradient. It is possible to numerically and approximately determine the
gradient through explicit finite perturbations of the initial state. But this
would be much too costly for practical implementation since it requires to
compute the cost function, i.e. to effectively integrate the model over the
assimilation period, as many times as there are independent components in the
initial states. Therefore to compute the gradient efficiently an adjoint model
should be used.

\subsection{Sequential Methods}
While variational method is based on minimization of the cost function within a
time interval, sequential method assimilates the data each time the observation
becomes available. In sequential method the adjusted model solution is
expressed as a linear combination of the forecast state and the data elements
following the equation:
\begin{equation}
x^a(k) = x^f(k) + K(k)(y^o(k)-Hx^f(k))
\end{equation}
Here $x^f(k)$ represents the forecast state, while $x^a(k)$ the analysis state,
i.e. the adjusted state. The gain matrix $K(k)$ describes how elements of the
state should be adjusted given the difference between the measurement and the
forecast. Different methods are available to determine $K$.

One of the popular sequential methods is the \emph{optimal interpolation} (\cite{Daley1991}. This method uses a gain matrix based on an empirical covariance
function or matrix. The basic assumption is that both the forecast and the
observation error are normally distributed with mean zero and known
covariances. The idea of optimal interpolation is to set the analysed state to
the conditional mean of the true state given the observations,
$x^a(k)=E[x^t(k)|y^o(k)]$. Application of Bayes theorem to Gaussian
distribution shows that this could be achieved with a linear gain:
\begin{eqnarray}
x^a(k)=x^f(k)+K(k)(y^o(k)-H'x^f(k)) \\
K(k)=P^f(k) H [H' P^f(k) H + R(k)]^{-1} \label{eq.K_OI}
\end{eqnarray}
The gain matrix $K$ in equation (\ref{eq.K_OI}) is known as the
\emph{conditional mean gain} or the \emph{minimal variance gain}. A problem is
how to choose suitable covariance matrices $P$ and $R$. In this method the user
needs to specify the error covariance at each assimilation time.

Another development in this class is the Kalman filtering. The Kalman filter
can be seen as an extension of the optimal interpolation scheme, accounting for
the evolution of errors from previous times. The target of the Kalman filter is
to obtain a distribution for the true state in terms of a mean $\hat{x}$ and
covariance $P$, given the model and the measurements. Like in the optimal
interpolation method, the Kalman filter also assumes normally distributed
forecast and observation errors. The Kalman filter is originally derived for
linear systems, which in state-space form can be written as:
\begin{eqnarray}
x(k+1)= M(k) x(k) + \eta(k) \\
y(k) = H(k) x(k) + \nu(k)
\end{eqnarray}
where $x$ is the system state, $A$ the linear model operator, $\eta \sim
N(0,Q)$ the system noise, $y$ the predicted observation, $H$ the observation
operator, and $\nu \sim N(0,R)$ the observation error. The Kalman filter
algorithm consists of two steps:
\begin{enumerate}
\item Forecast step:
   \begin{eqnarray}
     x^f(k+1) = M(k) x^a(k) \\
     P^f(k+1) = M(k) P^f(k) M'(k) + Q(k)
   \end{eqnarray}
 \item Analysis step:
   \begin{eqnarray}
     x^a(k)=x^f(k) + K(k) (y^o(k) - H(k) x^f(k)) \\
     P^a(k)=(I-K(k) H(k)) P^f(k) \\
     K(k) = P^f(k) H(k) (H'(k) P^f(k) H(k) + R(k))^{-1}
   \end{eqnarray}
\end{enumerate}

\section{Data assimilation  with the RRSQRT method}

The Kalman filter gives optimal estimates of $x$ and $P$ for linear models. The
main problem of applying the Kalman filter directly to environmental models is
the computation of the covariance matrix $P$. Since such models usually have a
big number of states (e.g. $O(10^4)$) the covariance will also become very big,
which causes very expensive computational costs or even the impossibility to
compute. Another problem is that the real life model is usually nonlinear.
Therefore methods are proposed and developed to modify the Kalman filter to
solve these difficulties. Two most popular algorithms are the
\emph{reduced-rank square-root} filter (\cite{VerlaanHeemink1997}) and
the \emph{ensemble Kalman filter} (\cite{Evensen1994}).

The reduced-rank square-root (RRSQRT) filter algorithm is based on a
factorization of the covariance matrix $P$ of the state estimate according to
$P=LL'$, where $L$ is a matrix with the $q$ leading eigenvectors $l_i$ (scaled
by the square root of the eigenvalues), $i=1,...,q$, of $P$ as columns. The
algorithm can be summarized as follows.

\begin{enumerate}
\item Initialization
  \begin{equation}
    [L^a(0)]=[l_1^a(0),\cdots,l_q^a(0)]
  \end{equation}
  where $l_i^a$ denotes the $q$ leading eigenvectors of the initial covariance
  matrix $P_o$.

\item Forecast step
  \begin{eqnarray}
    x^f(k)=M[x^a(k-1)] \\
    l^f_i(k)=\frac{1}{\epsilon} \{ M[x^a(k-1)+ \epsilon  l^a_i(k-1)]-M[x^a(k-1)] \} \\
    \tilde{L}^f(k)=[l_1^f(0),\cdots,l_q^f(0),Q(k-1)^{1/2}] \\
    L^f(k)=\Pi^f(k)\tilde{L}^f(k)
  \end{eqnarray}
where $\epsilon$ represents a perturbation, often chosen close to 1, $\Pi^f(k)$
is a projection onto the $q$ leading eigenvectors of the matrix
$\tilde{L}^f(k)\tilde{L}^f(k)'$. Note that here $M$, the model operator, need
not be linear.

\item Analysis step
  \begin{eqnarray}
    P^f(k)=L^f(k)L^f(k)' \\
    K(k) = P^f(k)H(k)'[H(k)P^f(k)H(k)'+R(k)]^{-1} \\
    x^a_k = x^fa_k + K_k [y^o_k - H_k x^f_k] \\
    \tilde{L}^a(k)=\{[I - K_k H_k] L_f, K_k R_k^{1/2}\} \\
    L^a(k)=\Pi^a(k)\tilde{L}^a(k)
  \end{eqnarray}
where $\Pi^a(k)$ is a projection onto the $q$ leading eigenvectors of the
matrix $\tilde{L}^a(k)\tilde{L}^a(k)'$. This reduction step is again introduced
to reduce the number of columns in $L^a_k$ to $q$ in $\tilde{L}^a(k)$.
\end{enumerate}

\subsection{Using the RRSQRT method}


%HIER MOET NOG WAT!!!

\subsection{The configuration of the RRSQRT method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of modes
\item {\tt mode}: Model configuration
\end{itemize}

\subsection{XML-example}

\begin{verbatim}
<?xml version="1.0" encoding="UTF-8"?>
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <!-- Observations used for assimilation -->
  <CTA_SOBS id="obs_assim" database="obs_lorenz.sql"/>

  <!-- Used model class -->
  <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />

  <!-- Filter configuration -->
  <method name="lbfgs">
  <parameter_calibration>
    <!-- Simulatie timespan en stapgrootte via CTA_Time -->
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <output>
        <filename>results_oscill.m</filename> 
       <CTA_TIME id="times_mean" start="0" stop="50.0" step="1"/>
    </output>
    <iteration maxit="10" maxln="20" tol_step="0.0001" tol_grad="0.0002">
    </iteration>
    <method nstore="3" c1="1E-4" c2="0.5" delta="1E-8">
    <!-- nstore specifies the storage size, i.e. max number of vector s and y to store -->
    <!-- c1 and c2 are the parameters used in the 1st and 2nd Wolfe conditions-->
    <!-- delta specifies the perturbation size in computing gradient-->
    </method>
       <model>
         <xi:include href="models/lorenz_sp.xml"/>
       </model>
  </parameter_calibration>
  </method>
</costa>
\end{verbatim}


\section{Data assimilation  with the Ensemble method}

While the RRSQRT represents the covariance matrix $P$ based on the first $q$
leading eigenvectors, the ensemble Kalman filter (EnKF) is based on a
representation of the probability density of the state estimate by a finite
number $N$ of randomly generated system states. The EnKF algorithm can be
summarized as follows.

\begin{enumerate}
\item Initialization

An ensemble of $N$ states $\xi^a_i(0)$ are generated to represent the
uncertainty in $x_o$.

\item Forecast step
  \begin{eqnarray}
    \xi_i^f(k)=M[\xi_i^a(k-1)] + \eta_i(k-1) \label{eq.EnKF_forecast} \\
    x^f(k)=\frac{1}{N} \sum_{i=1}^N \xi_i^f(k) \\
    E^f(k)=[\xi_1^f(k)-x^f(k),\cdots,\xi_N^f(k)-x^f(k)]
  \end{eqnarray}

\item Analysis step
  \begin{eqnarray}
    P^f(k)=\frac{1}{N-1}E^f(k)E^f(k)' \\
    K(k) = P^f(k)H(k)'[H(k)P^f(k)H(k)'+R(k)]^{-1} \\
    \xi^a_i(k) = \xi^f_i(k) + K(k) [y^o(k) - H(k) \xi^f_i(k) + \nu_i(k)]
  \end{eqnarray}
\end{enumerate}
where the ensemble of state vectors are generated with the realizations
$\eta_i(k)$ and $\nu_i(k)$ of the model noise and observation noise processes
$\eta(k)$ and $\nu(k)$, respectively.

For most practical problems the forecast equation (\ref{eq.EnKF_forecast}) is
computationally dominant. As a result the computational effort required for the
EnKF is approximately $N$ model simulations. The standard deviation of the
errors in the state estimate are of a statistical nature and converge very
slowly with the sample size ($\approx N$). Here it should be noted that for
many atmospheric data assimilation problems the analysis step is also a very
time consuming part of the algorithm.


%\subsubsection{Using the Ensemble method}

\subsection{The configuration of the Ensemble method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of ensembles
\item {\tt model}: The configuration of a model
\end{itemize}

\subsection{XML-example}
\begin{verbatim}
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <!-- De invoer m.b.t. de observaties -->
  <CTA_SOBS id="obs_assim" database="obs_advec1d.sql"/>

  <!-- Used model class -->
  <xi:include href="models/advec1d_cls.xml"/>

  <!-- De invoer m.b.t. het filter.  -->
  <!-- Dit filter maakt zelf de model-instantiaties aan -->
  <method name="ensemble"> 
  <ensemble_filter>
    <CTA_TIME id="simulation_span" start="0" stop="5" step="0.0005"/>
    <output>
       <filename>results.m</filename>
       <CTA_TIME id="times_mean" start="0" stop="50" step="0.05"/>
    </output>
    <modes max_id="50">
       <model>
       </model>
    </modes>
  </ensemble_filter>
  </method> 
</costa>


\end{verbatim}


\section{Data assimilation  with the Ensemble Square-Root method (ENSRF)}

%   \subsubsection{Kalman filtering with Ensemble Square-Root filter}
There are two fundamental problems associated with the use of EnKF. First is
that the ensemble size is limited by the computational cost of applying the
forecast model to each ensemble member. The second one is that small ensembles
have few degrees of freedom available to represent errors and suffer from
sampling errors that will further degrade the forecast error covariance
representation. Sampling errors lead to loss of accuracy and underestimation of
error covariances. This problem can progressively worsen, resulting in filter
divergence.

In ensemble square-root filters (ENSRF), the analysis step is done
deterministically without generating any observation noise realization
(\cite{Tippettetal2003}; Evensen, 2004). Since no random sample is generated,
this extra source of sampling error is eliminated. Therefore, these methods are
expected to perform better than the ones with perturbed observations for a
certain type of applications.

Dropping the time index $k$ for simplicity, the covariance analysis update of
the ENSRF is obtained by rewriting the covariance as
\begin{eqnarray}
  P^a = E^a E^{a'}=[I-P^f H'(H P^f H' + R)^{-1} H]P^f \\
  = E^f[I-E^{f'}H'(HE^f E^{f'}H'+R)^{-1}HE^f] E^{f'} \\
  = E^f [I- V D^{-1} V'] E^{f'}
\end{eqnarray}
where $V=(HL^f)'$ and $D=V' V + R$. Then from this equation it is clear that
the analysis ensemble can be calculated from
\begin{equation}
  \label{eq.ENSRF_EX}
  E^a = E^f X
\end{equation}
where $X X'=(I-V D^{-1} V')$. Therefore one can say that the updated ensemble
$E^a$ is a linear combination of the columns of $E^f$ and is obtained by
inverting the matrix $D$ and computing a matrix square root $X$ of the matrix
$[I- V D^{-1} V']$. Note that $X$ can also be replaced by $XU$, where $U$ is an
arbitrary orthogonal matrix, so that $(XU)(XU)'=XX'$.

As we have seen that in ENSRF the analysis step consists of determining the
\emph{transformation matrix}, $X$. A number of methods are available to compute
$X$. The method implemented in this study is derived as follows:
\begin{equation}
  \label{eq.xx1}
  X X'=(I-V D^{-1} V') = (I-V (V' V + R)^{-1} V')
\end{equation}
If we write $R=SS'$ then we can rewrite Equation \ref{eq.xx1} as
\begin{equation}
  \label{eq.xx2}
  X X'= \Psi' (\Psi \Psi' + I)^{-1} \Psi)
\end{equation}
where $\Psi=S^{-1}V=S^{-1}HE^f$. When computing the singular value
decomposition of $\Psi$, i.e. $\Psi=\Gamma \Sigma \Lambda'$, equation
(\ref{eq.xx2}) can be written as
\begin{eqnarray}
  X X'= (\Gamma \Sigma \Lambda')' ((\Gamma \Sigma \Lambda') (\Gamma \Sigma \Lambda')' + I)^{-1} (\Gamma \Sigma \Lambda')) \\
  = \Lambda (I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma) \Lambda' \\
  = (\Lambda \sqrt{(I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma)})(\Lambda \sqrt{(I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma)})'
\end{eqnarray}
Thus, a solution of the transformation matrix $X$ is given by
\begin{equation}
  X = \Lambda \sqrt{(I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma)}
\end{equation}
%       Moreover, the transformation matrix, $T$, for updating the ensembles can be shown to be
%       \begin{equation}
%           
%       \end{equation}

\subsection{Using the ENSRF method}

The user specifies the number of ensembles in the input file. Like in COFFEE
filter, the ensembles in ENSRF filter are also stored as instances of COSTA
models. An array of COSTA state vectors is also used for representing matrix
$E$ in equation (\ref{eq.ENSRF_EX}). There is always a decision to make whether
to store $E$ as a COSTA matrix or as an array of COSTA state vectors. We
inclined sometimes to store it as a COSTA matrix, since the ENSRF algorithm
consists of mostly matrix operations. However, since the model state is stored
as COSTA state vector it is easier to store $E$ as an array of COSTA state
vector. Moreover, there is also a function available in COSTA for performing
linear algebra operation between COSTA matrix and COSTA state, i.e.
cta\_state\_gemm. This function is used for example in computing the correction
for all ensemble members, which is done by multiplying the transformation
matrix and the forecast ensembles.

The implementation of ENSRF filter requires a singular value decomposition
(SVD). At the moment there is no COSTA function, which can perform SVD directly
to a COSTA matrix nor to an array of COSTA state vectors. Therefore in this
study this is done by assigning the values of the COSTA variables to a Fortran
array and use the LAPACK function DGESVD to perform the SVD.

This method also requires the square root of observation error covariance
$R^{1/2}$. However, as also mentioned in the previous subsection, the function
cta\_sobs\_getcovmat gives the covariance matrix $R$. Moreover, there is no
COSTA function yet available for computing the square root of a square matrix
like for example Cholesky decomposition. However since in the example models
the observation noise is always independent, we can easily compute $R^{1/2}$ by
taking square root of its diagonal elements.

\subsection{The configuration of the ENSRF method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of ensembles
\item {\tt model}: Model configuration
\end{itemize}

\subsection{XML-example}
 \begin{verbatim}
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <CTA_SOBS id="obs_assim" database="obs_oscill.sql"/>
  <!-- Used model class -->
  <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />
  <method name="ensrf"> 
<ensrf_filter>
  <!-- Simulatie timespan en stapgrootte via CTA_Time -->
  <CTA_TIME id="simulation_span" start="0" stop="5" step="0.0005"/>
    <output>
       <filename>results.m</filename>
       <CTA_TIME id="times_mean" start="0" stop="50" step="0.05"/>
    </output>

     <modes max_id="50">
        <model>
          <xi:include href="models/oscill_sp.xml"/>
        </model>        
     </modes>
</ensrf_filter>
  
  </method> 
</costa>
 \end{verbatim}




\section{Data assimilation  with the COFFEE method}

%      \subsubsection{Kalman filtering with COFFEE filter}
COFFEE (Complementary Orthogonal subspace Filter For Efficient Ensembles) is a
hybrid filter, which combines the RRSQRT filter and the EnKF
(\cite{Heemink2001}). One problem of the RRSQRT algorithm is that repeated
projection on the leading eigenvalues leads to a systematic bias in forecast
errors. Because of the truncation, the covariance matrix is always
underestimated, which may result in a filter divergence problem. The truncated
part of the covariance matrix does not contribute to the improvements of the
state estimate. The COFFEE filter attempts to solve this problem by
representing the truncated part of the covariance matrix as random ensembles
and to add them to the EnKF part. The RRSQRT part acts as a variance reductor
for the ensemble filter, thus reducing the statistical errors of the Monte
Carlo approach. Moreover, by embedding the reduced-rank filter in an EnKF the
covariance is not underestimated, eliminating the filter divergence problems of
the reduced-rank approach (also for very small numbers of $q$).

The COFFEE algorithm can be summarized as follows:
\begin{enumerate}
\item Initialization
  \begin{equation}
    \label{eq.COFFEE_initialmatrix}
          [L^a(0) E^a(0)]=[l_1^a(0),\cdots,l_q^a(0),\xi_1^a(0),\cdots,\xi_N^a(0)]
  \end{equation}
where $l_i^a$ denotes the $q$ leading eigenvectors of the initial covariance
matrix $P_o$. The random ensemble members are generated only to represent the
truncated part of the covariance matrix $P_o-L^a(0)L^a(0)'$.
\item Forecast step

The $L^a$ is updated using the RRSQRT update as follows
\begin{eqnarray}
  x^f(k)=M[x^a(k-1)] \\
  l^f_i(k)=\frac{1}{\epsilon} \{ M[x^a(k-1)+ \epsilon  l^a_i(k-1)]-M[x^a(k-1)] \} \label{eq.COFFEE_lf} \\
  \tilde{L}^f(k)=[l_1^f(0),\cdots,l_q^f(0),Q(k-1)^{1/2}] \label{eq.COFFEE_Lf}\\
  L^f(k)=\Pi^f(k)\tilde{L}^f(k)
\end{eqnarray}
where $\epsilon$ represents a perturbation, often chosen close to 1, $\Pi^f(k)$
is a projection onto the $q$ leading eigenvectors of the matrix
$\tilde{L}^f(k)\tilde{L}^f(k)'$.

The ensemble $\xi_i^a$ is updated using similar equations with the EnKF update
\begin{eqnarray}
  \xi_i^f(k)=M[\xi_i^a(k-1)] \\
  x^f(k)=\frac{1}{N} \sum_{i=1}^N \xi_i^f(k) \\
  E^f(k)=[\xi_1^f(k)-x^f(k)+\eta_1(k),\cdots,\xi_N^f(k)-x^f(k)+\eta_N(k)]
\end{eqnarray}
where $\eta_i$ are the ensembles representing the truncated part of the
covariance matrix with $E[\eta_i(k)\eta_i(k)']=[I-\Pi^f(k)]\tilde{L}^f(k)
\tilde{L}^f(k)'[I-\Pi^f(k)]'$.

\item Analysis step

In the analysis step the gain matrix $K$ is computed using
\begin{eqnarray}
  P^f(k)=L^f(k)L^f(k)' + \frac{1}{N-1}E(k)E(k)' \\
  K(k) = P^f(k)H(k)'[H(k)P^f(k)H(k)'+R(k)]^{-1} \\
  x^a(k) = x^f(k) + K(k) [y^o(k) - H(k) x^f(k)] \\
  \tilde{L}^a(k)=\{[I - K(k) H(k)] L_f(k), K(k) R(k)^{1/2}\} \label{eq.COFFEE_La} \\
  L^a(k)=\Pi^a(k)\tilde{L}^a(k) \\
  \xi^a_i(k) = \xi^f_i(k) + K(k) [y^o(k) - H(k) \xi^f_i(k) + \nu_i(k)]
\end{eqnarray}
\end{enumerate}

\subsection{Using the COFFEE method}
The COFFEE filter requires the user to specify in the input file the number of
modes for the RRSQRT part and the number of ensembles for the EnKF part.
Moreover, the perturbation $\delta$ (i.e. $\epsilon$ in equation
(\ref{eq.COFFEE_lf})) needs also to be specified.

In this implementation, the modes and ensembles are represented as COSTA model
instances. Here each model instance represent one mode or ensemble. This is
chosen for supporting parallel computation implementation in the future.

In the implementation of COFFEE filter, the initial matrix $[L^a(0) E^a(0)]$ in
equation \ref{eq.COFFEE_initialmatrix} is always assumed to be zero. This
however should be extended to accommodate the option where the user can specify
this from the input file. Moreover, the size of $L$ is limited to a maximum of
200. This is due to the fact that its dimension varies according to the number
of noise parameters in the model as well as to the number of observations.
Since we know the number of observations only at the analysis-step and that the
number of observations may vary, it is not possible to determine the dimension
of $L$ a priori. Furthermore, in the present implementation the computation
cost increases with number of observations. More study is required to develop
the implementation for solving this problem. It should be noted also that the
present implementation does not support covariance localization using Schur
product either.

The COFFEE filter also requires the square root of the model error covariance
matrix $Q^{1/2}$ in equation (\ref{eq.COFFEE_Lf}) as well as the observation
error covariance matrix $R^{1/2}$ in equation (\ref{eq.COFFEE_La}). The
covariance matrices $Q$ and $R$ can be obtained by using the functions
cta\_model\_getnoisecovar and cta\_sobs\_getcovmat respectively. However since
we mostly work with the square root of the covariance matrices, it is better if
these functions are modified to give the square root matrices. Note that
cta\_model\_getnoisecovar gives the noise covariance matrix in terms of an
array of COSTA state vectors, while cta\_sobs\_getcovmat gives the observation
noise covariance matrix in term of COSTA matrix. This may cause confusion for
new programmers working with COSTA. Moreover, it may be better to use a similar
name for the two functions. This makes the functions easier to remember.

Note that for general nonlinear models, the additional columns of matrix $L^f$
with matrix $Q^{1/2}$ in equation (\ref{eq.COFFEE_Lf}) can be computed using
finite difference as follows
\begin{eqnarray}
  \hat{x}(k) = M[x(k-1),u(k-1),0] \\
  x^n_i(k) = M[x(k-1),u(k-1),\epsilon Q^{1/2}(:,i)] \\
  l^n_i(k) =  \frac{1}{\epsilon} (x_i(k) - \hat{x}(k))\\
  \tilde{L}^f(k)=[l_1^f(k),\cdots,l_q^f(k),l^n_1(k),...,l^n_w(k)]
\end{eqnarray}
where $M$ is the general nonlinear model operator, which is a function of the
previous time-step state vector $x$, the input forcing $u$ and the noise
$\eta$. The state vector $\hat{x}$ is obtained from the deterministic model,
i.e. when no noise is present, while $x^n_i$ is from the stochastic model, with
the realization $\eta=\epsilon Q^{1/2}(:,i)$, where $Q^{1/2}(:,i)$ refers to
the $i^{th}$ column vector of matrix $Q^{1/2}$. The implementation requires a
function which can set the model noise realization. However, since there is no
such function available yet, it is not possible to implement this in COSTA at
the moment.

\subsection{The configuration of the COFFEE method}
\begin{itemize}

\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of modes
\item {\tt delta}: perturbation, see equation \ref{eq.COFFEE_lf}
\item {\tt model}: Model configuration
\item {\tt ensembles}: the number of ensembles
\end{itemize}

\subsection{XML-example}

\begin{verbatim}

<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <CTA_SOBS id="obs_assim" database="obs_oscill.sql"/>

  <!-- Used model class -->
  <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />


  <method name="coffee"> 
<coffee_filter>
  <CTA_TIME id="simulation_span" start="0" stop=".5" step="0.0005"/>
    <output>
       <filename>results.m</filename>
       <CTA_TIME id="times_mean" start="0" stop="50" step="0.05"/>
    </output>

     <modes max_id="10">
        <delta id="1E-4"> </delta>
        <model>
          <xi:include href="models/oscill_sp.xml"/>
        </model>        
     </modes>

     <ensembles max_id="5"> </ensembles>
</coffee_filter>
  
  </method> 
</costa>
\end{verbatim}
