\documentclass[draft,jgrga]{AGUTeX}

\begin{document}


\title{Kalman Filters Implemented in FVCOM}

\authors{Changsheng Chen, Paola Malanotte-Rizzoli, Robert C. Beardsley,
Zhigang Lai, Pengfei Xue, Jun Wei and Geoffrey W. Cowles}

\begin{article}

1. The Kalman Filters

Let ${\bf x^t}$ be an array of the true values, ${\bf x^f}$ be an
array of the forecast values, ${\bf x^a}$ an array of analysis
values, and ${\bf y}$ an array of observational values. We can
define that

Analysis error:
${\bf e_a}={\bf x^a}-{\bf x^f}$         (1.1)

Forecast error: ${\bf e}_f={\bf x}^t-{\bf x^f}$         (1.2)

Observational error: ${\bf e}_o={\bf y}-{\bf x^f}$           (1.3)

The forecast error covariance ${\bf P}$ and the observational error
covariance ${\bf R}$ thereafter can be defined as

${\bf P^f}=\overline{{\bf e^fe_f^T}}$, ${\bf R}=\overline{{\bf
e}_o{\bf e}_o^T}$           (1.4)

In a forecast model, the value of ${\bf x^f}$ at time step $i$ can
be predicted by

${\bf x^f}(i)={\bf M}_{i-1\rightarrow i}[{\bf x^a}(i-1)]$ (1.5)

where ${\bf M}$ presents the linearized model operator. The forecast
error covariance is equal to

${\bf P^f}={\bf M}_{i-1\rightarrow i}{\bf P^a}(i-1){\bf
M}_{i-1\rightarrow i}^T+{\bf Q}(i-1)$       (1.6)

where ${\bf Q}(i-1)$ is the system error covariance matrix. In the
Kalman filter forecast model system, the analysis values ${\bf x^a}$
is calculated by

${\bf x^a}(i)={\bf x^f}(i)+{\bf K}(i)[{\bf y}(i)-{\bf Hx^f}(i)]$
(1.7)

where ${\bf K}$ is the Kalman gain matrix, which is equal to

${\bf K}(i)={\bf P}^f(i){\bf H}^T[{\bf HP}^f(i){\bf H}^T+{\bf
R}(i)]^{-1}$        (1.8)

and ${\bf H}$ is an observation operator that functions as an
objective map to interpolate the model data onto the observational
points. The analysis error covariance is given as

${\bf P^a}(i)=[{\bf I}-{\bf K}(i)]{\bf P}^f(i)$        (1.9)

In general, the size of the covariance matrix is huge. For example,
$P^f:[{\bf N}\times{\bf N}]$ and $M:[{\bf N}\times{\bf N}]$, and $N$
can be of $O(10^6-10^7)$, which makes it impractical to use this
method on most computers. For this reason, a family of Kalman
Filters (Reduced Rank Kalman Filter, Ensemble Kalman Filter,
Ensemble Square Root Kalman Filter and Ensemble Transform Kalman
Filter) have been developed that require less computational power
than the original filter. A brief description of these other Kalman
Filters and how they are implemented in FVCOM is given below.

2. Reduced Rank Kalman Filter (RRKF)

Let ${\bf x^f}$ be an array with a dimension of $N$ (FVCOM output
and RRKF input); ${\bf x^a}$ an array with a dimension of $N$ (that
are output from RRKF to use to refresh the initial conditions for
FVCOM input); ${\bf y}$ an array with a dimension of ${\bf N_o}$
(that is an input for RRKF); ${\bf E}_r$ the resolved empirical
orthogonal functions (EOFs) with dimensions of $N \times {\bf N_e}$
(the stationary input of RRKF); ${\bf K}_r$ the stationary Kalman
gain in the reduced space of ${\bf N_e \times N_o}$; ${\bf N_e}$ the
dimension of the EOFs subspace; ${\bf i}$ and  ${\bf i + \Delta T}$
two subsequent assimilation time steps of FVCOM; ${\bf \Delta T}$
the assimilation interval; ${\bf M_i}$ the forward (nonlinear) FVCOM
model with initial conditions specified by the analysis solution;
and ${\bf D}$ is the spatially averaged standard deviation of each
variable.

A detailed description of RRKF was given in Buehner and
Malanotte-Rizzoli [2003]. S. Lyu, who worked with P. Rizzoli as a
postdoctoral investigator at MIT, helped us implement RRKF into
FVCOM. He worked together with P. Xue, Q. Xu and Z. Lai in testing
the RRKF code in idealized cases.

The RRKF is developed based on linear theory. In a linear system, if
the observational network, observational and model error covariances
are stationary and all neutrally stable and unstable modes are
measurable, the forecast error covariance reaches an asymptotically
stationary result and then a stationary Kalman gain (${\bf K_r}$)
can be derived efficiently by the doubling algorithm [Anderson and
Moore, 1979]. This stationary ${\bf K_r}$ is calculated from the
control model run and it is applied in the data assimilation process
of the RRKF.

The RRKF in FVCOM is operated following the procedure given below:

Step 1:  Determination the number of resolved EOFs ($E_r$) from the
control model run:

${\bf D}^{-1}{\bf \hat X_f \hat X_f^T D^{-1}}={\bf E \Lambda E^T}$
(2.1)

where ${\bf \hat {X}_f}={\bf x^f}-{\bf \overline{x}^f}$, ${\bf
\overline{x}^f}$ the mean of ${\bf x^f}$, $E$ the EOF matrix; and
$\Lambda$ is the diagonal covariance eigenvalue matrix.

Step 2: Linearization of the model in the resolved EOF subspace
(build ${\bf M}_r$ from ${\bf E}_r$):

${\bf M_{ri}}=\frac{1}{\alpha}{\bf E_r^TD^{-1}}[{\bf M}({\bf
x}_0+\alpha{\bf De_i})-{\bf M}({\bf x}_0)]$       (2.2)

where ${\bf M}$ presents the nonlinear model; subscripts ``$r$'' and
``$i$'' of ${\bf M}$ denote the linearized model in the resolved
subspace and the $i$th column of ${\bf M_r}$; $\alpha$ is the
perturbation size; ${\bf e_i}$ the $i$th retained EOF; and ${\bf
x}_0$ is the specified time mean of a long model run without
assimilation.

Step 3: Projection of the error covariance into the resolved EOFs
subspace and estimation of the model and observation errors:

${\bf P_r^f}={\bf E_r^fP^fE_r}; {\bf P_r^a} = {\bf E_r^TP^\alpha
E_r}; {\bf M_r}={\bf E_r^TME_r}; {\bf Q_r}=\gamma \Lambda$     (2.3)

${\bf R}={\bf R_m}+{\bf H_uP_u^fH_u^T}$       (2.4)

where ${\bf P_r^a}$ is the analysis error covariance matrix in the
resolved subspace; ${\bf Q_r}$ the ‘pesudo’ model error covariances;
${\bf R}$ the observational error covariance; ${\bf R_m}$ the actual
measurement error; and ${\bf P_r^f}$ the forecast error covariance
matrix in the resolved subspace.

Step 4: Calculation of the stationary Kalman gain $K_r$ in the
reduced subspace by the doubling algorithm, estimation of the
difference between the observations and forecast and projection to
the full space by multiplying $E_r$

${\bf K_r}({\bf t})={\bf  P_r^f}({\bf t})({\bf H_rP_r^f}({\bf
t}){\bf H_r^t}+{\bf R})^{-1}$   (2.5)

${\bf P_r^a}({\bf t})=({\bf I}-{\bf K_r}({\bf t}){\bf H_r}){\bf
P_r^f}({\bf t})$       (2.6)

${\bf P_r^f}({\bf t}+1)={\bf M_rP_r^a}({\bf t}){\bf M_r^T}+{\bf
Q_r}$ (2.7)

where ${\bf H_r}={\bf HDE_r}$ and ${\bf P_r^f}$ is asymptotically
stationary as $t \rightarrow \infty$ if ${\bf H}$, ${\bf R}$ and
${\bf Q}$ is stationary with linear dynamics.

Step 5:  Data assimilation by using a stationary Kalman gain $K_r$

${\bf x}^a({\bf t})={\bf x}^f({\bf t})+{\bf DE_rK_r}[{\bf y}({\bf
t})-{\bf Hx}^f({\bf t})]$         (2.8)

RRKF works efficiently in a linear system but not for a nonlinear
system. We tested it for various idealized cases such as tidal waves
in the circular lakes, the flooding/drying process in the
rectangular shape estuary. It produces a fast convergence solution
for the linear tidal wave case, but never converges in the estuarine
case characterized with nonlinear dynamics.

3. Ensemble Kalman Filter (EnKF)

Evensen (1994) suggested that the error covariance relative to the
mean of ensemble model results could provide a better estimation of
the error covariance defined in the classical Kalman Filter. The
EnKF is constructed by running a forecast model driven by a set of
initial conditions and then estimate the error covariance relative
to the ensemble mean to determine the ensemble analysis values for
the next time step forecast.

Let $k$ denote the $k$th ensemble member and $N_E$ the total number
of the ensemble members selected in the forecast model run. The
forecast value at time step $i$ for the $k$th ensemble model run can
be estimated by

${\bf x_k^f}(i)={\bf M}_{i-\Delta t \rightarrow i}[{\bf
x_k^a}(i-1)], k=1,2,...,N_e$       (3.1)

and the analysis values at time step $i$ for the $k$th ensemble
model run are calculated by

${\bf x_k^a}(i)={\bf x_k^f}(i)+{\bf K}(i)[{\bf y}_k(i)-{\bf
Hx_k^f}(i)], k=1,2,...N_e$       (3.2)

Define that

${\bf X_f}=\{[{\bf x_k^f}-\overline{{\bf x}}^{\bf f}]/\sqrt{{\bf
N_e}-1}\}, k=1,2,...N_e$      (3.3)

then the forecast error covariance can be estimated by

${\bf P}^f(i)\approx {\bf X}_f(i){\bf X}_f^T(i): [N\times
N_e][N_e\times N]$         (3.4)

The Kalman gain is equal to

${\bf K}(i)={\bf X}_f{\bf X}_f^T{\bf H}^T(H{\bf X}_f{\bf X}_f^T{\bf
H}^T+{\bf R})^{-1}$         (3.5)

To conduct the EnKF, we need to create an ensemble of the
observational data constructed with the perturbation relative to the
real value, i. e.,

${\bf y_k}={\bf y}+\delta _{\bf k}, \delta _k=N(0,\sqrt{R})$ (3.6)

where

${\bf R}=\overline{\delta \delta^{\bf T}}$              (3.7)

In the situation with a sufficiently large number of ensembles, the
ensemble analysis error covariance matrix can be updated with a
relationship as

${\bf P_e^a}({\bf t})=({\bf I}-{\bf K_r}({\bf t}){\bf H}){\bf
P_e^f}({\bf t})$        (3.8)

With assumption of Gaussian distribution, this will give us optimal
analysis values in the maximum likelihood sense and in the minimum
variance sense. In the situation with a small number of ensembles,
the perturbed observations required by EnKF may cause a rank
deficiency problem for the estimation of ${\bf P^f}$ and an
underestimate of ${\bf P^a}$, which leads to filter divergence
[Whitaker and Hamill, 2002]. The filter divergence problem can be
controlled by using the covariance localization [Houterkamer and
Mitchell, 2001] and covariance inflation [Wang and Bishop, 200].
EnKF is suitable for both linear and nonlinear systems.  For a
linear system, RRKF works well and also fast, but it sometimes fails
to resolve linear waves in the idealized, linear coastal ocean
system. In such a case, EnKF works well.

4. Ensemble Square-Root Kalman Filter (EnSKF)

The EnKF is a filter that requires perturbed sets of observational
values. In this system, the perturbed observational values are
usually constructed by a control observation in addition to random
noise sampled from the assumed observational error distribution.
There are derivatives of EnKF that are conducted with deterministic
observational ensembles. These filters are known as the EnSKF
[Whitaker and Hamill, 2002], Ensemble Transform Kalman Filter
(EnTKF) [Bishop et al., 2001] and the Ensemble Adjustment Kalman
Filter (EnAKF) [Anderson, 2001]. Actually, these filters are a
family of the deterministic square root filters [Tippett et al.,
2003]. A brief description of one type of EnSKF is given below.

Assuming that the forecast and observational error covariance is
Gaussian distributed, the ensemble Kalman Filter provides optimal
analysis values which satisfy the classical Kalman Filter covariance
form as

${\bf P_e^a}=({\bf I}-{\bf K_rH}){\bf P_e^f}$      (4.1)

and the Kalman gain is

${\bf K}_e={\bf P}_e^f{\bf H}^T[{\bf HP}_e^f{\bf H}^T+{\bf
R}_e]^{-1}.$                                           (4.2)

Eq. (4.1) can be rewritten into

${\bf P_e^a}={\bf X_aX_a^T}$         (4.3)

where

${\bf X_a}={\bf X_f\hat{T}}$         (4.4)

and ${\bf \hat{T}}$ is the square root matrix defined as

${\bf \hat{T}}=\{ {\bf I}-{\bf X_fH^T}[{\bf HX_fX_f^TH}^T+{\bf
R}]^{-1}{\bf HX_f}\}$     (4.5)

where ${\bf I}$ is $N_e \times N_e$ unit matrix. In this case, an
ensemble of the analysis deviation ${\bf X}_a$ can be estimated
deterministically from an ensemble of the forecast deviation, and
also the analysis ensemble has a desired error covariance.

Assuming that the observation errors are uncorrelated, the
observations can be assimilated serially. In this case, ${\bf
HX_fX_f^TH^T}+{\bf R}$ in (4.5) is simplified to be a scalar in the
case of a single observation. Then the square root matrix can be
easily calculated as

${\bf \hat{T}}=[1-\beta {\bf X_f^TH^THX_f}][1-\beta {\bf
X_f^TH^THX_f}]^T $      (4.6)

where $\beta = [a + \sqrt{Ra}]^{-1}$. EnSKF does not require a
perturbed observation set, so that the ensemble mean of analysis
values ${\bf \overline{x}^a}$ can be updated directly by the mean of
forecast ensemble produced from a traditional Kalman Filter equation
as

${\bf \overline{x}}^a(i)={\bf \overline{x}}^f(i)+{\bf
X_fX_f^TH^T}[{\bf HX_fX_f^TH^T+R}]^{-1}({\bf y}-{\bf H
\overline{x}^f})$     (4.7)

and each analysis member ${\bf x}_k^a$ can be calculated by

${\bf \overline{x}_k^a}(i)={\bf \overline{x}}^a(i)+{\bf
X}_k^a\sqrt{N_e-1}$ $k=1,2,...N_e$          (4.8)

5. Ensemble Transform Kalman Filter (EnTKF)

The EnTKF is very similar to EnSKF with a linear transformation of
the forecast perturbation into the analysis perturbation by ${\bf
\hat{T}}$

${\bf X_a}={\bf X_f\hat{T}}$        (5.1)

Bishop et al. [2001] proposed the form of the transformation matrix
$T$ as:

${\bf \hat{T}}=C({\bf \tilde{A}}+{\bf I})^{-1/2}$       (5.2)

where columns of the matrix $C$ contain the eigenvectors of ${\bf
X_f^TH^TR^{-1}HX}_f$ and ${\bf \tilde{A}}$ is the nonzero diagonal
matrix that satisfies a relationship with $C$ as

${\bf X_f^TH^TR^{-1}HX}_f={\bf C\tilde{A}C^T}$          (5.3)

In this case,

${\bf x}^a(i)={\bf \overline{x}}^f(i)+{\bf X_fC\tilde{A}}^{1/2}({\bf
\tilde{A}}+{\bf I})^{-1}{\bf E}^T\{{\bf R}^{-1/2}{\bf y}-{\bf H
\overline{x}^f})$       (5.4)

where

${\bf E}={\bf HX_fC \tilde{A}}^{-1/2}$        (5.5)

In FVCOM, EnTKF is used for adaptive observation optimization in
which ${\bf X}_a$ is used to evaluate the analysis ensemble error
covariance under different observational strategies. This approach
allows us to use the model to determine the optimal observational
network for a selected region.

References

Anderson, B.  and J. B. Moore (1979), Optimal Filtering,
Prentice-Hall, 357pp.

Anderson, J. L. (2001), An ensemble adjustment Kalman filter for
data assimilation. Mon Weather Rev 129: 2884–2903.

Bishop, C. H., B. J. Etherton and S. J. Majumdar (2001), Adaptive
sampling with the Ensemble Transform Kalman Filter. Part I:
theoretical aspects, Mon.Wea,Rev., 129, 420-436.

Buehner, M. and P. Malanotte-Rizzoli (2003), Reduced-rank Kalman
filters applied to an idealized model of the wind-driven ocean
circulation, J. Geophys. Res., 108, 3192, doi:10.1029/2001JC000873

Evensen, G. (1994), Sequential data assimilation with a nonlinear
quasi-geostrophic model using Monte Carlo methods to forecast error
statistics. J. Geophys. Res., 99, 10,143-10,162.

Houtekamer, P. and H. L. Mitchell (1998), Data assimilation using an
ensemble Kalman filter technique. Mon. Wea. Rev., 126, 796-811.

Tippett, M. K., J. L. Anderson, C. H. Bishop, T. M. Hamill and J. S.
Whitaker (2003). Ensemble square root filters.  Mon. Wea. Rev., 131,
1485-1490.

Wang, X., and C. H. Bishop (2003), A comparison of breeding and
ensemble transform Kalman filter ensemble forecast schemes. J.
Atmos. Sci., 60, 1140-1158.

Whitaker, J. S. and T. M. Hamill (2002), Data assimilation without
perturbed observations. Mon. Wea. Rev., 130, 1913-1924.

Figure Caption

Figure S1: Schematic of the implementation of RRKF into FVCOM.

Figure S2: Schematic of the implementation of EnKF (EnSKF/EnTKF)
into FVCOM.

\end{article}

\end{document}
