# Ocean Data Assimilation using Sequential Methods based on the Kalman Filter

Pierre Brasseur

► To cite this version:

Pierre Brasseur. Ocean Data Assimilation using Sequential Methods based on the Kalman Filter. E. Chassignet; J. Verron. Ocean Weather Forecasting, pp.271-316, 2006, 10.1007/1-4020-4028-8_10. hal-00232778

HAL Id: hal-00232778

https://hal.science/hal-00232778v1

Submitted on 19 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# OCEAN DATA ASSIMILATION USING SEQUENTIAL METHODS BASED ON THE KALMAN FILTER

## From theory to practical implementations

Pierre Brasseur

CNRS/LEGI, Grenoble, France

**Abstract:** The main purpose of this chapter is to review the fundamentals of the Kalman Filter for ocean data assimilation and to expose the basic ingredients of practical assimilation algorithms developed for applied ocean research and operational forecasting, focusing mainly on high-resolution applications. Important implementation issues such as the reduction in dimensionality of the estimation problem, the simplification of the schemes based on static error covariances, the formulation of low-rank filters, the problem of consistency verification, and the concepts of adaptivity and incremental analysis updating will be addressed using scientific and operational examples. Finally, the discussion will conclude with a number of key questions related to the assimilation challenges of the next decade.

**Keywords:** Data assimilation, mesoscale ocean circulation, reduced-order Kalman filters, statistical estimation, ocean forecasting.

## 1. Introduction

Operational ocean prediction systems are being developed with a variety of objectives in mind, such as ocean current hindcasting and short-range forecasting, estimation of the thermodynamic state of the ocean for seasonal and climate predictions, and production of retrospective analyses of the changing ocean through the merging of models and data. In addition to simply combining model estimates with observations, data assimilation provides a means to systematically compare theoretical
models with reality, leading to potential improvements in modelling and observing systems. It is therefore likely that applied ocean research will benefit strongly from operational progress, and vice versa. The question concerning the dominant energetic activity of the mesoscale ocean, its non-deterministic nature and the interactions with the large-scale circulation make the challenge unique, requiring sophisticated numerical models and assimilation methods that make the best use of sparse observations. To produce reliable forecasts, the models must be initialized with conditions that represent as accurately as possible the actual state of the ocean at eddy-resolving resolution. Due to the chaotic properties of ocean dynamics, the forecast range cannot be extended beyond the limit of predictability of the system and the model has to be reinitialized intermittently by correcting the forecast with the most recent observations. Fortunately, the arrival of satellite observations, in particular satellite altimetry, has provided the observational basis needed to respond appropriately to the “high-resolution challenge”. In order to extract the best possible information from the new data, it is necessary to assess how reliable the model forecast and the observations are. Therefore, error estimates on the measurements and the model prediction are inherent in the assimilation process.

Data assimilation is traditionally formulated as a least-squares estimation problem. Among the various methodological approaches, the theory of optimal statistical estimation, and more specifically the Kalman filtering approach, is well suited to provide a solution to the Best Linear Unbiased Estimation. Since Kalman in 1960, sequential filtering methods have been thoroughly explored and applied to state estimation. An extended version of the Kalman Filter (KF) has been derived for non-linear models, known as the Extended Kalman Filter (EKF) [Jazwinski, 1970; Gelb, 1974]. In spite of a fairly simple theoretical framework, the question of its applicability in assimilating observations into high-resolution, non-linear numerical models of the ocean circulation is far from trivial. As stated by Courtier [1997], the scientific difficulty associated with data assimilation is in finding algorithms which simplify the search for an affordable solution in terms of computer resources, while preserving some of the essential characteristics. A hierarchy of approximations to the Kalman filter has been defined to make the methodology suitable for solving large-dimension problems. These developments represent a substantial part of the research effort devoted to oceanic and atmospheric data assimilation over the past 10 years.

Compared with variational approaches such as the 4D-VAR, statistical algorithms require less initial investment in terms of coding and are naturally designed to incorporate gradual developments. This is one of
several reasons why most assimilation methods used today in operational forecasting systems are inspired by the statistical approach. In the context of GODAE (Global Ocean Data Assimilation Experiment), one objective will be to test and compare the bunch of algorithms implemented in operational systems, in order to better understand the importance of the various possible approximations made in each of these.

The objectives of this chapter are to review the fundamentals of sequential data assimilation for ocean state estimation and to expose the basic ingredients of practical assimilation algorithms developed for applied ocean research and operational systems, focusing mainly on high-resolution applications. Section 2 is dedicated to the fundamentals of applied estimation methods leading to the KF equations. Numerical representations of the mathematical objects introduced by the theory will then be illustrated using oceanographic examples in Section 3. Section 4 will provide a brief description of traditional simplifications of the KF. In Section 5, we discuss various approaches to reduce the size of the estimation problem and, in Section 6, we derive the framework of low-rank Kalman filters. The important question of the verification of consistency will be addressed in Section 7, where the concept of adaptivity is also mentioned. Finally, a number of advanced implementation issues such as the transition to incremental/smoothing algorithms will be discussed in Section 8, before concluding the chapter.

## 2. Kalman filtering: Fundamentals

### 2.1 Problem definition

In this section, we introduce the basic assimilation problem in the state space using the conventional notations proposed by Ide et al. [1997]. The goal here is not to present a rigorous and comprehensive derivation of the Kalman filter, which can be found elsewhere in dedicated text books (e.g., Gelb [1974]), but rather to introduce a simplified framework that still contains the essential characteristics needed to illustrate the more advanced concepts and implementation issues discussed in the following sections.

To start, let us assume that some *a priori* knowledge about the state of the ocean is available at time $t_i$, represented by vector $\mathbf{x}_i^a$. A physical model is also available to describe the transition of the state vector from time $t_i$ to time $t_{i+1}$, which is represented by a numerical (matrix) operator $\mathbf{M}(t_i, t_{i+1})$. A linear model will be considered at this stage. The dimension of the state space is noted as $n$. The state vector $\mathbf{x}$ contains the minimum set of independent variables needed to perfectly characterize the state of the system at any time. The model can be used
to simulate the transition of the state vector up to time $t_{i+1}$,

$$ \mathbf{x}_{i+1}^f = \mathbf{M}(t_i, t_{i+1})\mathbf{x}_i^a \quad (1) $$

where $\mathbf{x}_{i+1}^f$ is the vector describing the “forecast” state of the system. At instant $t_{i+1}$, another piece of useful information is available about the state, collected in the observation vector $\mathbf{y}_{i+1}$ of dimension $p$. From the two independent pieces of information $\mathbf{x}_{i+1}^f$ and $\mathbf{y}_{i+1}$, how can the true state of the system $\mathbf{x}_{i+1}^t$ be best estimated at time $t_{i+1}$? To answer this question, we need to know more about the precision of the different pieces of information.

## 2.2 Uncertainties and PDFs

The precision of the forecast $\mathbf{x}_{i+1}^f$ can be quantified in terms of errors on the initial guess and numerical model errors. The difference between the initial guess $\mathbf{x}_i^a$ and the true state at time $t_i$ is the error vector noted $\epsilon_i^a = \mathbf{x}_i^a - \mathbf{x}_i^t$. Of course, its value is unknown but we can make a number of assumptions about its statistical properties: we will assume that the estimation $\mathbf{x}_i^a$ is unbiased ($\overline{\epsilon}_i^a = 0$, where the overbar represents the expected value), and its error $\epsilon_i^a$ is distributed as a gaussian, multivariate random variable. The corresponding probability density function (pdf) is

$$ \epsilon_i^a \rightarrow N(0, \mathbf{P}_i^a) \sim \exp \left[ -\frac{1}{2} \epsilon_i^{aT} \mathbf{P}_i^{a-1} \epsilon_i^a \right] \quad (2) $$

where $\mathbf{P}_i^a = \overline{\epsilon_i^a \epsilon_i^a T}$ is the $n \times n$ error covariance matrix associated with $\mathbf{x}_i^a$ and $T$ denotes the transpose. Error covariances are formally obtained by multiplying an error vector by its transpose and averaging over many realizations, leading to symmetric and positive definite matrices. Morrison [1988] contains excellent background information on multivariate statistical methods. Similarly, the model operator $\mathbf{M}(t_i, t_{i+1})$ is imperfect and the simulation error is noted as follows:

$$ \eta = \mathbf{M}(t_i, t_{i+1})\mathbf{x}_i^t - \mathbf{x}_{i+1}^t \quad (3) $$

Again, the individual realization of this error is unknown (otherwise a perfect model operator could be run) but its statistical distribution is assumed to be gaussian and centered ($\overline{\eta} = 0$):

$$ \eta \rightarrow N(0, \mathbf{Q}) \sim \exp \left[ -\frac{1}{2} \eta^T \mathbf{Q}^{-1} \eta \right] \quad (4) $$

where $\mathbf{Q} = \overline{\eta \eta^T}$ is the $n \times n$ model error covariance matrix. We assume in addition that $\epsilon_i^a$ and $\eta$ are uncorrelated: $\overline{\epsilon_i^a \eta^T} = 0$. In general, these
statistical assumptions are quite crude approximations of the actual error distributions (in particular, a bias in the model is very common), but they are very convenient in deriving a baseline of optimal estimation. A schematic description of the error diagram in the state space is illustrated in figure 1.

<figure><img src="image_2.png" /><figcaption>Figure 1. Vector diagram of analysis and forecast errors in the state space.</figcaption></figure>

With the definitions introduced above, the forecast error $\epsilon_{i+1}^f$ can be broken down as:

$$ \epsilon_{i+1}^f = \mathbf{x}_{i+1}^f - \mathbf{x}_{i+1}^t = \mathbf{M}\mathbf{x}_i^a - (\mathbf{M}\mathbf{x}_i^t - \eta) = \mathbf{M}\mathbf{\epsilon}_i^a + \eta \quad (5) $$

and the statistical properties of the forecast error can be determined easily if the model is linear. Indeed, Eqs. (5), (2) and (4) implies that the forecast state is unbiased ($\overline{\epsilon_{i+1}^f} = \overline{\mathbf{M}\mathbf{\epsilon}_i^a} + \overline{\eta} = 0$) and normally distributed, i.e.

$$ \epsilon_{i+1}^f \rightarrow N(0, \mathbf{P}_{i+1}^f) \sim \exp \left[ -\frac{1}{2} \epsilon_{i+1}^{fT} \mathbf{P}_{i+1}^{f-1} \epsilon_{i+1}^f \right] \quad (6) $$

with the forecast error covariance matrix given by:

$$ \mathbf{P}_{i+1}^f = \overline{\epsilon_{i+1}^f \epsilon_{i+1}^{fT}} = \overline{\mathbf{M}\mathbf{\epsilon}_i^a \epsilon_i^{aT}} \mathbf{M}^T + \overline{\eta\eta^T} = \mathbf{M}\mathbf{P}_i^a \mathbf{M}^T + \mathbf{Q} \quad (7) $$

This equation is the first fundamental equation of the KF which can be interpreted as follows: the error on the initial state is transformed during the forecast step by the model dynamics (the error being amplified by unstable modes, while it is damped by stable modes) and by the model
imperfections which increase the forecast error covariance. Understanding the actual benefits and the practical limitations of error covariance propagation through model dynamics has been a major research issue over recent years. Although simple in its algebraic form, this equation contains several major difficulties that prevent an explicit computation in the context of realistic models, as discussed in the following sections.

In order to optimally combine the forecast with new data, the precision of the observations arising at time $t_{i+1}$ must be quantified. The vector of observations is related to the true state as follows:

$$ \mathbf{y}_{i+1} = \mathbf{H}\mathbf{x}_{i+1}^t + \epsilon_{i+1}^o \quad (8) $$

where $\mathbf{H}$ is the observation operator which computes the equivalent of the observations from the model state. The observational errors $\epsilon_{i+1}^o$ are assumed to be centered ($\overline{\epsilon_{i+1}^o} = 0$), uncorrelated with the forecast error and having a covariance matrix $\mathbf{R} = \overline{\epsilon_{i+1}^o \epsilon_{i+1}^T}$. Observation errors measure the misfit between the data and the equivalent of the observations in the true state, i.e. $\mathbf{H}\mathbf{x}_{i+1}^t$. They include not only the errors of the observational system but also the errors associated with the operator $\mathbf{H}$ arising, for example, from the numerical interpolation of the data. Again, a gaussian pdf can be assumed for the statistical distribution of the observation errors

$$ \epsilon_{i+1}^o \rightarrow N(0, \mathbf{R}) \sim \exp \left[ -\frac{1}{2} \epsilon_{i+1}^{oT} \mathbf{R}^{-1} \epsilon_{i+1}^o \right] \quad (9) $$

to make the rest of the development easier and derive the KF equations.

## 2.3 Optimal analysis

The pdf given by (6) determines the a priori statistical distribution of the true state $P(\mathbf{x}_{i+1}^t)$, while (9) provides the probability of getting measurements $\mathbf{y}_{i+1}$ given the true state, i.e. $P(\mathbf{y}_{i+1} | \mathbf{x}_{i+1}^t)$. It is then straightforward to deduce the a posteriori probability of the truth, given the observations, by using the Bayes formula:

$$ P(\mathbf{x}_{i+1}^t | \mathbf{y}_{i+1}) = \frac{P(\mathbf{y}_{i+1} | \mathbf{x}_{i+1}^t) \cdot P(\mathbf{x}_{i+1}^t)}{P(\mathbf{y}_{i+1})} \quad (10) $$

The state maximizing the posterior probability distribution is the maximum likelihood solution of this inverse problem. A comprehensive formulation of data assimilation and inverse methods using the bayesian approach has been proposed by van Leeuwen and Evensen [1996]. In Eq. (10), the denominator is just a scaling factor (the integral of the
numerator for all possible states) which can be ignored in determining
the maximum of the posterior pdf. The gaussian distributions (6) and
(9) imply that

$$
P(\mathbf{y}_{i+1} | \mathbf{x}_{i+1}^t) \cdot P(\mathbf{x}_{i+1}^t) \sim \exp \left[ -\frac{1}{2} \left\{ \epsilon_{i+1}^{oT} \mathbf{R}^{-1} \epsilon_{i+1}^o + \epsilon_{i+1}^{fT} \mathbf{P}_{i+1}^{f-1} \epsilon_{i+1}^f \right\} \right] \quad (11)
$$

and the optimal estimation of $\mathbf{x}_{i+1}^t$ is the state vector maximizing (11)
or, equivalently, minimizing

$$
J = \left\{ \epsilon_{i+1}^{oT} \mathbf{R}^{-1} \epsilon_{i+1}^{o} + \epsilon_{i+1}^{fT} \mathbf{P}_{i+1}^{f-1} \epsilon_{i+1}^{f} \right\} \quad (12)
$$

As a result of error definitions, the optimal combination of the fore-
cast and observed information corresponds to the minimum of the cost
function

$$
J(\mathbf{x}) = (\mathbf{y}_{i+1} - \mathbf{H}\mathbf{x})^T \mathbf{R}^{-1} (\mathbf{y}_{i+1} - \mathbf{H}\mathbf{x}) + (\mathbf{x}_{i+1}^f - \mathbf{x})^T \mathbf{P}_{i+1}^{f-1} (\mathbf{x}_{i+1}^f - \mathbf{x}). \quad (13)
$$

This quadratic form contains two terms measuring the misfit with the
data and the misfit with the forecast, weighted by the inverse of their
respective error covariances. Using the calculus of variations, an implicit
equation for the optimal state noted $\mathbf{x}_{i+1}^a$ is obtained:

$$
\delta J(\mathbf{x}) = 0 \Rightarrow \mathbf{x}_{i+1}^a = \mathbf{x}_{i+1}^f + \mathbf{P}_{i+1}^f \mathbf{H}^T \mathbf{R}^{-1} (\mathbf{y}_{i+1} - \mathbf{H}\mathbf{x}_{i+1}^a) \quad (14)
$$

which can be solved explicitly after some algebra, yielding

$$
\mathbf{x}_{i+1}^{a} = \mathbf{x}_{i+1}^{f} + \mathbf{P}_{i+1}^{f} \mathbf{H}^{T} (\mathbf{H} \mathbf{P}_{i+1}^{f} \mathbf{H}^{T} + \mathbf{R})^{-1} (\mathbf{y}_{i+1} - \mathbf{H} \mathbf{x}_{i+1}^{f}) \quad (15)
$$

The optimal state is obtained by correcting the forecast with a weighted
measure of the misfit between the observations and the prior estimate
(i.e. the innovation vector $\mathbf{d}_{i+1} = \mathbf{y}_{i+1} - \mathbf{H}\mathbf{x}_{i+1}^f$). The analysis is thus
simply the result of the combination of two Gaussian probability density
functions. The weight matrix of dimensions $n \times p$

$$
\mathbf{K}_{i+1} = \mathbf{P}_{i+1}^{f} \mathbf{H}^{T} (\mathbf{H} \mathbf{P}_{i+1}^{f} \mathbf{H}^{T} + \mathbf{R})^{-1} \quad (16)
$$

is the so-called Kalman gain, which involves the forecast and observation
error covariance matrices. It can be interpreted as a ratio between the
error variance of the forecast and the total error variance (the sum of the
forecast and the observation error variance) projected in the observation
space: the larger the forecast errors, the larger the correction to the
forecast. In the limit of perfect observations ($\mathbf{R} \sim \mathbf{0}$) of the whole state
vector ($\mathbf{H} \sim \mathbf{I}$), the Kalman gain matrix converges towards the identity
and the optimal estimate becomes a perfect fit of the observations. By contrast, if the forecast is extremely accurate ($\mathbf{P}^f \sim 0$) compared with the observations, the correction is negligible. Interesting similarities between this equation and scalar formulations of least squares problems can be found in Kalnay [2003]. Equation (15) is the second fundamental equation of the Kalman filter.

## 2.4 The sequential assimilation cycle

The optimal state estimate (15) at time $t_{i+1}$ can be used as the initial conditions for a new forecast up to time $t_{i+2}$ when new observations become available, and the process can be repeated recursively. To sum up, the algorithm of an assimilation cycle contains two main steps: the forecast step for transitioning the model state and the associated error covariance between time $t_i$ and time $t_{i+1}$, and the analysis step for correcting the forecast using the data available at time $t_{i+1}$. We reproduce here the complete set of the KF equations extended to non-linear models $M$ and observation operators $H$.

Starting from initial conditions $\mathbf{x}_i^a$ and $\mathbf{P}_i^a$, the forecast step equations are:

$$ \mathbf{x}_{i+1}^f = M(t_i, t_{i+1}) \{\mathbf{x}_i^a\} \quad (17) $$

and

$$ \mathbf{P}_{i+1}^f = \mathbf{M}\mathbf{P}_i^a\mathbf{M}^T + \mathbf{Q} \quad (18) $$

where $\mathbf{M}$ is the tangent linear operator derived from $M(t_i, t_{i+1})$. Thus, a linearization of the model about the non-linear evolution between $t_i$ and $t_{i+1}$ is performed to propagate the error covariance.

The forecast step is followed by an analysis step in which $\mathbf{y}_{i+1}$ is used to correct $\mathbf{x}_{i+1}^f$:

$$ \mathbf{x}_{i+1}^a = \mathbf{x}_{i+1}^f + \mathbf{K}_{i+1} (\mathbf{y}_{i+1} - \mathbf{H} \{\mathbf{x}_{i+1}^f\}) \quad (19) $$

using the Kalman gain

$$ \mathbf{K}_{i+1} = \mathbf{P}_{i+1}^f \mathbf{H}^T [\mathbf{H} \mathbf{P}_{i+1}^f \mathbf{H}^T + \mathbf{R}]^{-1} \quad (20) $$

where $\mathbf{H}$ is the gradient of $H$ computed about $\mathbf{x}_{i+1}^f$. It can be demonstrated that the matrix $\mathbf{K}_{i+1}$ corresponds to the minimization of the trace of the analysis error covariance on $\mathbf{x}_{i+1}^a$ [Miller 1989], given by

$$ \mathbf{P}_{i+1}^a = \mathbf{P}_{i+1}^f - \mathbf{P}_{i+1}^f \mathbf{H}^T [\mathbf{H} \mathbf{P}_{i+1}^f \mathbf{H}^T + \mathbf{R}]^{-1} \mathbf{H} \mathbf{P}_{i+1}^f = [\mathbf{I} - \mathbf{K}_{i+1} \mathbf{H}] \mathbf{P}_{i+1}^f \quad (21) $$
This allows us to write the gain also as

$$ K_{i+1} = P_{i+1}^a H^T R^{-1}, \qquad (22) $$

provided that $R$ is invertible. Equation (21) shows that the uncertainty in the forecast is reduced during the analysis according to the amount of additional information assimilated in the system.

A sequential assimilation run is then conducted by repeating this forecast/analysis cycle in sequence. Since only data from the past influence the best estimate at a given time, the assimilation procedure belongs to a class of filtering methods. These contrast with smoothing methods (e.g. Fukumori [2001]) in which data from both the past and the future are used to estimate the optimal state of the system at a given time. The analysis error covariance reflects the competition in the Kalman filter between this accumulation of past information and the error growth due to instability mechanisms and model imperfections. Figure 2 conceptually illustrates the filtering process in sequential data assimilation.

<figure><img src="image_3.png" /><figcaption>Figure 2. Conceptual representation of filtering with sequential assimilation.</figcaption></figure>

### 3. From theory to real ocean applications

The Kalman Filter has been primarily developed in the context of ballistic applications, involving dynamical models of fairly low dimension. The more recent interest in the KF in Earth sciences (numerical weather prediction or oceanography) has raised new issues related to the huge number of degrees of freedom taken into account by the models, with consequences for the size of the discretized operators and the quantity of information to be manipulated. It is therefore important to examine at
this stage the practical representation of the mathematical objects intro-
duced in the previous section, with examples taken from high-resolution
circulation systems developed from an operational perspective.

**3.1 The state vector and model operator**

The state vector **x** is a discrete representation of the variables involved
in the description of the system state. It characterizes information about
the space variability of the physical or biological quantities, and about
the multivariate relationships between the different dynamical variables.
In a numerical model of the Primitive Equations (PE), **x** typically con-
tains the 3D discretized temperature, salinity, zonal and meridional
velocities, and the 2D sea-surface height or barotropic streamfunction
field computed prognostically at every gridpoint of a finite-difference
mesh. The typical size of a PE state vector **n** is approximately given by
N<sub>X</sub> × N<sub>Y</sub> × (N<sub>Z</sub> × 4 variables + 1 variable), where N<sub>X</sub>, N<sub>Y</sub> and N<sub>Z</sub> are
the horizontal and vertical grid dimensions. The length **n** of the state
vector is typically 10<sup>6</sup> to 10<sup>8</sup> in scientific or operational applications.
A biological state vector would contain, for instance, the concentration
distributions of nutrients, plankton, dissolved and particulate matter,
etc. [Carmillet et al., 2001].

Let us consider, for example, the 1/12° model configuration of the
North Atlantic ocean that has been developed using the HYbrid Coor-
dinate Ocean Model (HYCOM) with a horizontal grid size of approx-
imately 1400 x 1400, and 26 layers in the vertical direction. Figure 3
shows a snapshot extracted from a simulation, illustrating the space
variability described by the multivariate HYCOM state vector.

Note that the vertical hybrid coordinate used in HYCOM is a com-
bination of geometric and dynamic vertical coordinates that evolves dy-
namically with the state of the system itself [Bleck, 2002]. Due to the
occurrence of outcropping layers at the base of the mixed layer, the
number of discrete variables may be different from one timestep to
another, and the dimension of the state vector is therefore dynamically-
dependent. This feature slightly complicates the handling of the state
vector [Brankart et al., 2003; Birol et al., 2004], compared with more
conventional models based on static vertical coordinates such as OPA.

The size problems arise mainly from equations (18), (20) and (21)
which involve the manipulation of *n* × *n* matrices. With a state dimension
of 10⁸, the storage of a full *n* × *n* matrix would require a memory of
10⁵ gigabytes, representing about 1000 times the capacity of the largest
computers available today. For this reason, it is necessary to by-pass the
explicit representation of such matrices in the algorithms. For instance,
the model $M(t_i, t_j)$ will never be available as an explicit operator (even for linear dynamics) but the numerical code of the ocean model will be used as a routine to compute the time evolution of the model state and forecast error components, as required by (17) and (18).

<figure><img src="image_4.png" /><figcaption>*Figure 3*. Physical representation of a portion of the HYCOM state vector in a North Atlantic model. The panel on the right shows a snapshot of the SSH field; the panels on the left show a vertical section in the 3D temperature and salinity fields at 30°W, superimposed to the vertical grid oriented along the layer interfaces. The full state vector also includes the velocity components (not shown here).</figcaption></figure>

## 3.2 The observation vector

The dimension $p$ of the observation vector $y$ depends on the capacity of the observation system and the frequency of data assimilation, but in general it is much smaller than the dimension of the state vector. The data sets available to control ocean circulation models in scientific or operational exercises can be categorized into measurements from space, which primarily reflect the surface signature of the ocean circulation and ocean-atmosphere interactions, and *in situ* measurements devoted to the monitoring of the ocean's interior.

Oceanic quantities measured from space include essentially Sea-Surface Temperature (SST), Sea-Level Anomalies (SLA) and ocean colour (which can be used to estimate the chlorophyll concentration in the upper ocean). Note that other important satellite data types will become available in the near future, such as surface salinity measurements from the SMOS (Soil Moisture and Ocean Salinity) mission and sea-ice observations from CRYOSAT. Unlike conventional measurements from field campaigns, satellite-based instruments are operated in routine over long
periods of time and provide considerable information about the ocean's
horizontal and temporal variability. Figure 4 illustrates a typical SST
picture obtained from composite AVHRR (Advanced Very High Resolution Radiometer) images, and a network of altimetric ground tracks covered by two satellites flying simultaneously over a 7-day period. Along the altimeter tracks, the measurements are available every 7 km and the separation between two tracks of Topex-Poseidon is about 300 km at the Equator. These data can be assimilated “along track” as shown in figure 4, providing in principle the best possible utilization of the observed information if the assimilation scheme is optimally tuned. Otherwise, the raw measurements can be transformed into gridded products using sub-optimal interpolation methods before assimilation.

The surface nature of satellite data poses specific challenges with
regard to data assimilation because the surface information has to be
projected downward to reconstruct the 3D-content of the ocean signal.
As illustrated in the next section, the extrapolation process is achieved
through the use of 3D, multivariate error covariances in the assimilation
scheme. These surface data remain insufficient, however, for describing
aspects of sub-surface variability, and other data available in the form of
vertical temperature and salinity profiles from hydrographic casts, moor-
ings or expandable bathythermographs (XBT) and profiling floats from
the ARGO international program provide valuable information about
the vertical stratification. It is therefore essential to assimilate both data
types in a consistent manner.

Two additional operators related to the data must be introduced: the observation error covariance matrix **R** and the observation operator **H** (or a non-linear **H**). The observation operator converts the forecast state into “first guesses” of the observations. This conversion is needed because the observed variables may not be located on the model grid points so that horizontal or vertical interpolations are necessary. In addition, relationships of varying complexity may exist between the observed quantity and the model variables. One example is the ocean colour, which is related to the phytoplankton concentration through a fairly complex, and sometimes approximate, relationship. Another example is the relationship between sea-level (measured by altimetry) and the prognostic variables of a rigid-lid ocean model [Pinardi et al., 1995]. In order to avoid too complex observation operators in the algorithms, it is sometimes advisable to augment the state vector with diagnostic variables (such as the surface pressure of a rigid-lid model which can be diagnosed as a function of the prognostic variables) that can be linked to the observed variables more easily.
<figure><img src="image_5.png" /></figure>

Figure 4. Multi-satellite data set collected during the week August 19-26, 1993, illustrating the different sampling properties. On the left: altimeter measurements (in cm) obtained by merging the Topex/Poseidon and ERS ground tracks; on the right: composite AVHRR picture of the sea-surface temperature.

The observation error accounted for by **R** has different possible sources. One source is the instrumental error which can often be considered as spatially uncorrelated. Another is the so-called “representativeness error”, associated with variability described by the data at scales that cannot be faithfully represented by the model grid. A third source is the error associated with the mapping **H** between the model and the observation space. The spectrum of these errors is mainly concentrated on the short scales, and it is often a reasonable approximation to represent **R** by a diagonal matrix of error variances. Note, however, that the reduction in the dimensionality of the estimation problem introduced in Section 5 will possibly be a source of correlated observation error.

## 3.3 Error covariance matrices

As mentioned above, the assimilation sequence must be initialized with some initial guess for the state $\mathbf{x}_0$ and the associated error covariance $\mathbf{P}_0$. The initialization of $\mathbf{P}$ critically determines the functioning of the filter beyond the first assimilation stages through the process of error dynamics [Ballabrera et al., 2001]. In order to understand the role played by error covariance matrices, let us consider the idealized case of an analysis step with one single observation, $\eta_0$. The observed variable, noted $\eta$, is assumed to be one of the discrete elements of the state vector
for simplicity, so that $\mathbf{y}$ is a scalar ($p = 1$) and $\mathbf{H}$ is a single-row vector of the form

$$ \mathbf{H} = [0, \dots, 0, 1, 0, \dots, 0] \qquad (23) $$

The Kalman gain (20) is then a single-column vector which simplifies as

$$ \mathbf{K} = \frac{1}{p_{\eta\eta} + \epsilon^2} \{\mathbf{P}_0\}_{\eta} \qquad (24) $$

where $\epsilon^2$ is the error variance of the single observation, $p_{\eta\eta}$ is the error variance of the observed variable (the $\eta$-diagonal element of $\mathbf{P}_0$) and $\{\mathbf{P}_0\}_{\eta}$ is the $\eta$-column of the error covariance matrix. The correction to the initial guess is then proportional to $\{\mathbf{P}_0\}_{\eta}$ weighted by the prior model-data misfit:

$$ \mathbf{x}^a - \mathbf{x}_0 = \frac{1}{p_{\eta\eta} + \epsilon^2} \{\mathbf{P}_0\}_{\eta} (\eta - \eta_0) \qquad (25) $$

Thus, any row or column of the error covariance matrix can be interpreted as a multivariate influence function associated with the observed state variable. This explains the crucial role of the error structures specified in $\mathbf{P}_0$ and shows the importance of considering dynamically-balanced error covariance matrices.

The last operator to be prescribed is the model error covariance matrix $\mathbf{Q}$. This includes all the errors associated with the various physical parameterizations necessary in the model (mixing, diffusion, turbulent closure, hydrostatic approximation, etc.), the errors in the atmospheric forcings and more generally in the boundary conditions, and the errors due to the numerical discretization on the horizontal and vertical dimensions. Note that $\mathbf{Q}$ represents the model error statistics accumulated during an assimilation interval, and should not be confused with errors generated at every time step. Those errors are clearly distributed over a wide spectrum of space scales, and it would be extremely difficult to prescribe a full $\mathbf{Q}$ matrix. By considering the prior misfit between model simulations and observations, it is possible to derive some general properties of the model errors as manifested in the observation space. These difficulties promote the adoption of simplified parameterizations for $\mathbf{Q}$ following, for instance, the approach proposed by Dee [1991], and to adjust those parameterizations by sensitivity experiments. Note that systematic biases in the models often make the statistical assumption (4) inappropriate: such biases can be detected by examining the statistics of the innovation sequence, and in principle they should be removed from the forecast to preserve the essential properties of optimal analysis.
## 4. Simplified schemes based on static background errors

In the KF algorithm, the forecast error covariance is updated before every analysis step by using the model dynamics. Equation (18) can be rewritten as

$$ \mathbf{P}_{i+1}^{f} = \mathbf{M}(\mathbf{M}\mathbf{P}_{i}^{a})^{T} + \mathbf{Q} \quad (26) $$

showing that the model code has to be used to propagate the $n$ columns of $\mathbf{P}_i^a$ (in addition to propagate the model state itself). Given the huge size of realistic ocean models, this leads to computing requirements that even the largest computers in the world will not meet in a foreseeable future. In addition, the many imperfections inherent in the representation of the error covariance matrices $\mathbf{P}_i^a$ and $\mathbf{Q}$ make the explicit computation of this equation questionable.

If the flow of assimilated observations is fairly regular over time, it makes sense to assume that the errors in an assimilation sequence tend to fluctuate around some average level after a couple of cycles. This asymptotic behaviour reflects a balance between the increase of uncertainty during the forecast step and the error reduction during the analysis step. The existence of such an asymptotic limit for $\mathbf{P}_{i+1}^f$ provides justification for simply using a static error covariance, noted $\mathbf{B}$ for “background error” (a term usually adopted when the error statistics are not propagated from one assimilation cycle to the next), instead of explicitly computing the forecast error according to (18). Different techniques have been developed to prescribe the background error covariances in ocean and atmospheric assimilation schemes. These are discussed below.

### 4.1 Optimal Interpolation

Optimal Interpolation (OI) designates a wide range of statistical assimilation schemes in which the $\mathbf{B}$ matrix is pre-determined empirically. The main advantages of OI methods are their cost, their ease of use and the possibility of conducting sensitivity studies to test many different models of error covariances. Today, the vast majority of operational systems involved in GODAE (e.g., FOAM, MFS, HYCOM, MERCATOR) are based on OI schemes. In contradiction with its name however, OI is a sub-optimal assimilation process and it would actually be more correct to designate this class of methods as “statistical interpolation” [Daley, 1991].

Any covariance matrix can be normalized into a correlation matrix $\mathbf{C}$, dividing each component by the product of the error standard deviations.
Thus, **B** can still be written as

$$ \mathbf{B} = \mathbf{D}^{1/2} \mathbf{C} \mathbf{D}^{1/2} \qquad (27) $$

where **D** is the diagonal matrix of the variances. OI schemes commonly adopt further simplifications for modelling the background error covariances. Assuming that the correlations are functions of the distance only, analytical functions with open parameters such as correlation scales have been used for a long time to build **C** (e.g. Thiebaux [1985]). A significant property of well-posed correlation functions is their compact support which reflects the negligible influence of observations at large distances. It is also generally assumed that the correlations can be separated into a product of horizontal and vertical correlations. The concept of separability is related to the predominant role of stratification and the very different scales involved horizontally and vertically in the open ocean. However, different behaviours are expected near boundaries and coastal regions [Echevin et al., 2000]. The SOFA (System for Ocean Forecasting and Analysis) scheme used in the first operational MERCATOR prototype performs OI in the horizontal direction to combine altimeter data with the model predictions [De Mey and Benkiran, 2002]. It offers the possibility of using different error covariance models, such as :

$$ \mathbf{C}(\mathbf{x}_i, \mathbf{x}_j) = (1 + a l + \frac{1}{3} a^2 l^2) e^{-al} \quad \text{with} \quad l = \text{dist}(\mathbf{x}_i, \mathbf{x}_j) \qquad (28) $$

where the parameter *a* determines the horizontal correlation scale [Gavart et al., 1999]. Fine tuning of these correlation parameters can be achieved objectively via sensitivity studies or Monte-Carlo experiments [Auclair et al., 2003]. In the MERCATOR system, the value of *a* is defined geographically, reflecting longer horizontal correlation scales in the tropics than at high latitudes.

It is worth noting that some applications of (27) allow for the error variances to be updated during the assimilation sequence using empirical prediction schemes, while still keeping the correlation matrix unchanged [Rienecker and Miller, 1991].

## 4.2 Asymptotic approximation

Other interesting approaches such as the time-asymptotic filter approximation have been proposed to determine a background error covariance matrix that is consistent with the error dynamics of the KF [Fu et al., 1993; Fukumori et al., 1993; Fukumori and Malanotte-Rizzoli, 1995]. In equations (18), (20) and (21) above, operators **M**, **H**, **Q**, and **R** have been assumed time-invariant for simplicity. Satellite altimeters,
for instance, can be considered as a time-invariant observation system given that the mission tracks are repeated in an exact manner. By combining (18) and (21), the so-called Riccati equation is obtained:

$$ \mathbf{P}_{i+1}^f = \mathbf{M} \left[ \mathbf{P}_i^f - \mathbf{P}_i^f \mathbf{H}^T [\mathbf{H}\mathbf{P}_i^f \mathbf{H}^T + \mathbf{R}]^{-1} \mathbf{H}\mathbf{P}_i^f \right] \mathbf{M}^T + \mathbf{Q} \quad (29) $$

which can be iterated prior to the assimilation sequence. It can be shown that the solution of the Riccati equation converges towards a steady-state solution, say $\mathbf{B}^f$, provided the $\mathbf{M}$, $\mathbf{H}$ and $\mathbf{Q}$ matrices have desirable properties. The details of these conditions are explained in Fukumori et al. [1993]. In addition, a fast convergence can be obtained by using appropriate recursion methods. The $\mathbf{B}^f$ solution results from a balance between three effects: evolution due to model dynamics, increase due to model errors and decrease due to the new information from assimilated data. Interestingly, the utility of a stationary $\mathbf{B}^f$ matrix in the Kalman gain was demonstrated even in the presence of evolving properties of the observation system [Fukumori, 1995].

The interest of steady-state filters has been illustrated in a number of oceanographic case studies, although additional approximations such as order reduction and model simplifications are generally needed to make the assimilation effective with realistic models.

## 5. Reduced-order Kalman filters

The previous section has shown that OI-based assimilation schemes over-simplify the propagation of errors by neglecting dynamical principles and statistical information. An alternative way to make the KF tractable with large-size models has been explored with the concept of reduced order, which aims at decreasing the computational burden of the algorithm while preserving the essential characteristics of error dynamics. Experiences from atmospheric re-analyses indicate the importance of the “errors of the day” which can be dominated by short time scales, and which are ignored when the forecast error is approximated by a constant as with OI [Kalnay et al., 1997]. Similar behaviour is observed in the ocean, especially at scales dominated by instability mechanisms such as the scale of eddies, western boundary currents, etc [Ballabrera et al., 2001].

Several arguments support the concept of order reduction. Firstly, the ocean can be considered as a driven/dissipative dynamical system governed by an “attractor” of finite dimension. The existence of a global attractor has been proved for the Navier-Stokes equations, with a dimension bounded by a function of the Reynolds number [Lions et al., 1997]. Geostrophy, for instance, is one of the dominant properties of
the ocean attractor which cannot be ignored during the assimilation process. Therefore, it makes sense to renounce making corrections in the directions where the errors eventually die away because of the attractive nature of the system. Secondly, the state space as defined in Section 3 is determined by the number of degrees of freedom implied by the model discretization, which is much larger than the actual number of dynamical features of interest in the system. Thirdly, less is probably known about the errors than about the dynamics, and the lack of statistical information will make a full KF superfluous anyway [Cane et al., 1996], especially in the ocean where the flow of observations accumulates at a fairly slow rate.

Reducing the dimensionality of the problem can be formulated in terms of state space or error space with quite different implications. A variety of approaches widely used in oceanography are discussed below.

## 5.1 Reducing the dimensionality of the state space

The reduced state can be defined explicitly using a transformation operator $T$ to convert the full model state $x$ (of dimension $n$) into a reduced state $w$ of dimension $r < n$:

$$w = Tx \qquad (30)$$

The statistical properties of $x$ such as error covariance matrices can be easily transformed in the low-dimension space using (30). Dynamical equations can be derived for $w$ using a pseudo-inverse of $T$ and the KF can be entirely reformulated in the low-dimension space, with the condition that the null space associated with (30) be dynamically uncoupled from the reduced space [Fukumori and Malanotte-Rizzoli, 1995]. By introducing this transformation, the $r$ elements of $w$ become the actual degrees of freedom of the estimation problem.

Many possibilities exist for defining the transformation, such as a truncation of the model spectrum or a selection of multivariate modes of system variability [Cane et al., 1996]. Another approach explored by Dee [1991] takes advantage of physical relationships between certain model variables. The reduction of the state space dimension can also be achieved simply by building the estimation vector from a selection of model state variables on a coarser grid [Fukumori, 1995] or with dynamical variables closely correlated with the observed signal. It is worth noting that, for sub-spaces which only preserve the scales larger than those of the observed signal, an extra term must be added to the observation error covariance matrix $R$ to account for the truncated modes
[Cane et al., 1996]. This kind of representativeness error typically cor-
responds to spatially correlated signals, and in this case off-diagonal
elements must be included in $\mathbf{R}$.

In the open ocean, dynamical considerations can justify the omission
of some of the PE variables in the state vector, such as the horizon-
tal velocity components which are expected to adjust very quickly to
the density properties. Different choices of estimation space, motivated
also by the faith in the error statistics to be prescribed for the selected
variables, are discussed by Brankart et al. [2003] for primitive equation
circulation models.

The question of assimilating satellite altimetric data into ocean circu-
lation models has been approached many times by reducing the statis-
tical estimation problem to the surface variables. The downward pene-
tration of assimilated information is then achieved empirically, using a
variety of statistical (e.g. Hurlburt [1986]; Mellor and Ezer [1991]; De
Mey and Benkiran [2002]) or physical (e.g. Cooper and Haines [1996];
Oschlies and Willebrand [1996]) extrapolation schemes on the vertical
direction. These different vertical projection methods, as well as the
dynamical adjustment process described by Brankart et al. [2003] in the
context of HYCOM, can be considered as approximations of the pseudo-
inverse of $\mathbf{T}$ needed to convert the reduced state back to the full model
space after statistical estimation.

## 5.2 Reducing the dimensionality of the error space

A number of other methods are based on approximations of the error
with fewer degrees of freedom than the model itself. Error sub-spaces
are built with the aim of selectively correcting the model state along
the most representative directions of the forecast error. The concept fits
well with statistical assimilation schemes because the analysis increment
computed by (20) can only take place within the sub-space spanned by
the forecast error covariance. The unknown are the coefficients of the
correction projected along the multivariate error directions. The concept
of Error Subspace Statistical Estimation (ESSE) introduced by Lermu-
siaux [1999] has been developed on the basis of this principle. Note that
recent oceanographic studies have also started to explore the potential
benefit of order reduction using variational assimilation methods [Robert
et al., 2005].

The sub-space can be prescribed as a time-invariant set of error direc-
tions. The success of assimilation depends essentially on the capacity of
the sub-space to capture the observed variability of the system. A num-
ber of applications have been conducted successfully along these lines in the tropical oceans, where the dynamics is “slow” and nearly linear [Verron et al., 1999]. Alternatively, the error sub-space can be allowed to evolve using deterministic or stochastic approaches, or a mixture of both. Evolving sub-spaces are in general more suitable to track the non-linear evolution of “fast” energetic error modes, or to dynamically adjust the unbalanced components of the sub-space [Lermusiaux, 2001]. Related concepts have been introduced in the context of atmospheric data assimilation by Cohn and Todling [1996], who proposed approximate schemes for error covariance propagation in the case of stable and unstable dynamics.

## 5.3 Low-rank error covariance matrix

By construction, a covariance matrix is symmetric and positive definite and can always be decomposed as $\mathbf{P} = \mathbf{N}\mathbf{\Lambda}\mathbf{N}^T$. The columns of $\mathbf{N}$ are formed by the orthonormalized eigenvectors $\mathbf{n}_k$ of $\mathbf{P}$, and $\mathbf{\Lambda}$ is a diagonal matrix formed with the corresponding eigenvalues $\lambda_k$. The inverse of $\mathbf{P}$ is then given by $\mathbf{P}^{-1} = \mathbf{N}\mathbf{\Lambda}^{-1}\mathbf{N}^T$. The pdf defined by (9) can be reformulated as follows:

$$ \epsilon \rightarrow N(0, \mathbf{P}) \sim \prod_{k=1}^{n} \exp \left[ -\frac{1}{2} \lambda_k^{-1} \epsilon_k^2 \right] \quad (31) $$

where $\epsilon_k$ is the component of the error vector $\epsilon$ in the $\mathbf{n}_k$ direction.

A reduced-order Kalman filter can be implemented by approximating the error covariance matrix $\mathbf{P}$ with only the “leading” columns of $\mathbf{N}$ associated with the $r$ largest eigenvalues. The pdf defined by such a low-rank matrix is obtained by taking the limit of (31) when $\lambda_k \to 0$ for $k > r$ (assuming that the eigenvalues have been sorted in decreasing order). Equation (31) shows that the probability tends to zero if $\epsilon_k \neq 0$ and the error vectors are confined in a sub-space of dimension $r$. From a stochastic point of view, the leading vectors describe the principal axes of the probability ellipsoid oriented along the dominant directions of uncertainty; from an algebraic point of view, they define the basis of a sub-space where the error is expected to lie.

## 5.4 Specification of error sub-spaces using EOFs

Several strategies can be adopted to determine the leading directions of the error sub-space. One of them is the Singular Value Decomposition (SVD) of a covariance matrix constructed with prescribed analytical functions. Other methods have been proposed in the literature which utilize singular, Lyapunov or breeding vectors of the transition matrix
(e.g., Miller and Erhet [2002]). An approach based on EOFs was put forward by Cane et al. [1996] to elaborate a reduced state Kalman filter, and by Pham et al. [1998] in the context of the Singular Evolutive Extended Kalman (SEEK) filter.

A practical way to estimate a low-rank error covariance matrix is to perform an EOF analysis of state vectors generated by a prior model simulation. The EOF analysis provides a compact description of the spatial and temporal variability of the model in terms of orthogonal functions. Usually, most of the variance of the time sequence is described by the first few orthogonal functions whose patterns may then be linked to dynamical mechanisms [Emery and Thomson, 1998]. In order to compute an EOF basis, a series of s model state vectors $\mathbf{x}(t_i)$ is extracted at regular intervals from a free model simulation, and the vectors

$$ \frac{\mathbf{x}(t_i) - \overline{\mathbf{x}(t_i)}}{\sqrt{s-1}} \qquad (32) $$

form the columns of a “scatter matrix” $\mathbf{X}$ of dimensions $n \times s$. The normalization factor $\frac{1}{\sqrt{s-1}}$ is introduced so that the unbiased estimation of the covariance matrix is given by the product $\mathbf{X}\mathbf{X}^T$. The size $s$ of the sample is always much smaller than the dimension $n$ of state vectors involved in realistic ocean models. The EOFs correspond to the orthonormalized eigenvectors $\mathbf{N}$ of the $n \times n$ matrix $\mathbf{X}\mathbf{X}^T$ which has a rank necessarily smaller than or equal to $s$. The explicit calculation of the matrix $\mathbf{X}\mathbf{X}^T$ is not required to compute the eigenmodes. Indeed, the matrix $\mathbf{X}^T\mathbf{X}$ of much lower dimensions ($s \times s$) has the same eigenvalues $\mathbf{\Lambda}$, and its eigenvectors $\mathbf{V}$ allow the computation of the “large” eigenmodes $\mathbf{N}$ by simple multiplication, $\mathbf{N} = \mathbf{X}\mathbf{V}$. This provides a first practical method to calculate the EOFs. A second possible approach is to perform a singular value decomposition of the matrix $\mathbf{X}$ directly, using standard SVD algorithms [Kelly, 1988; Emery and Thomson, 1998]. This provides a decomposition of the form:

$$ \mathbf{X} = \mathbf{N}\sqrt{\mathbf{\Lambda}}\mathbf{V}^T \qquad (33) $$

The EOFs determined by the two methods are identical, but the SVD has the advantage of greater computational stability. One should keep in mind that the EOF decomposition of multivariate state vectors also depends on the units adopted to represent the different physical quantities, and the choice of a particular metric may be more critical than computational stability in determining the structure and ordering of the dominant EOFs.

The $r$ dominant eigenvectors form the columns of an array noted $\mathbf{N}_0$ with dimensions $n \times r$, and the diagonal matrix of the associated
eigenvalues is noted $\Lambda_0$. The covariance can be approximated by the matrix

$$ \mathbf{X}\mathbf{X}^T \approx \mathbf{N}_0\mathbf{\Lambda}_0\mathbf{N}_0^T = \mathbf{S}_0\mathbf{S}_0^T \quad (34) $$

of rank $r$, where $\mathbf{S}_0 = \mathbf{N}_0\sqrt{\mathbf{\Lambda}_0}$ contains dimensional quantities. Assimilation experiments can be initialized using $\mathbf{S}_0\mathbf{S}_0^T$ as a guess of the initial error covariance.

The minimum dimension of the error sub-space has to be determined by practical considerations. A condition of stability for a reduced-rank KF in idealized conditions (linear model and autonomous system) is shown to be that $r$ should be larger than $r^*$, where $r^*$ is the number of eigenvalues of $\mathbf{M}$ having an absolute value larger than or equal to 1 (Pham et al. [1998]; Carme et al., [2001]). Since this number cannot be systematically evaluated with large ocean models, more pragmatic criteria must be considered to define the truncation threshold. For example, the statistical representativity of the $r$ retained EOFs can be measured by the percentage

$$ I = \frac{\sum_{k=1}^{r} \lambda_k}{\sum_{k=1}^{s} \lambda_k} $$

where the $\lambda_k$ are the eigenvalues. The choice of an acceptable threshold of explained variance provides a means to determine the minimum dimension of the EOF basis for practical assimilation problems. This number typically ranges from a few tens to a few hundreds.

The assumption underlying the procedure described above is that the model is unbiased (which is a necessary condition for the Kalman filter to yield optimal estimations) and sufficiently “good” for its intrinsic variability to be statistically representative of real world variability. The hypothesis that a model provides an unbiased simulation implies the use of the mean vector $\mathbf{x}(t_i)$ as a first guess. Sometimes, these ideal conditions are not verified with a free model simulation and more complex strategies must be developed. For example, an assimilation sequence may be carried out in a first stage with a simplified method which does not require EOFs (such as an OI or nudging scheme) so as to bring the model and observations together; an EOF analysis of the resulting fields is then performed, and a new assimilation sequence can be obtained in a second stage with a reduced-rank KF. The process can be iterated several times to provide an improved EOF basis which combines the variabilities of the model dynamics and of the observations.
## 5.5 Examples of sub-spaces based on EOFs

The EOFs calculated with primitive-equation models are fully multi-variate and three-dimensional, and cover the whole model domain: all the variables of the state vector (SSH, temperature, salinity, zonal and meridional velocities) are considered together in a dynamically consistent manner. The extrapolation of the data from observed to non-observed variables is performed along the directions represented by these EOFs which connect all grid points of the numerical domain.

<figure><img src="image_6.png" /><figcaption>Figure 5. Sea-surface height component of the first EOF (in cm) of the interannual variability of the circulation in the Tropical Pacific (Parent [2000]).</figcaption></figure>

The physical nature of an EOF basis is discussed by Parent *et al.* [2003], who studied the variability of the circulation in the Tropical Pacific ocean during the period 1994-1998 using a primitive-equation model of the Equatorial basin between 20°N and 20°S and a SEEK filter to assimilate satellite altimeter data. They computed an EOF decomposition of the simulated variability over the 5-year integration period and selected the first 15 EOFs to build a reduced basis for assimilation. The first dominant EOF illustrates the well-known west-east seesaw, which is the most important feature of the El Niño and La Niña phases of the ENSO phenomenon. In 1997, the warm pool of the western basin migrated towards the eastern basin and produced a positive SLA, whereas during the La Niña phase (second part of 1998), the opposite movement was observed (figure 5).

An issue of practical interest is the estimation of small correlations associated with distant variables, which is a well-known difficulty of finite
sampling procedures such as that described above. As pointed out by Houtekamer and Mitchell [1998], the accurate estimation of small correlation coefficients would necessitate a very large number of independent model samples to guarantee the statistical convergence of the computations. In order to keep the number of samples within tractable limits and prevent the data from exerting a spurious influence at remote distances through large-scale signatures in the EOFs, a technique based on EOFs with compact support has been setup by Testut et al. [2003] in which regional sub-domains of adjustable size are considered to characterize the error sub-space. From a theoretical point of view this approximation can be justified by the argument that, when the ocean surface is divided into local regions of moderate size, the background error in such regions tends to lie in a sub-space of much smaller dimension than the full ocean state. Similar hypotheses have been put forward to develop local ensemble Kalman filters for atmospheric data assimilation [Ott et al., 2004].

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Multivariate local EOFs of the mesoscale ocean variability in the Gulf Stream region simulated by an eddy-permitting model: surface signature of temperature and sea-surface height of the first three dominant modes, and associated vertical extensions of the temperature and salinity structures (reproduced from Testut [2000]).</figcaption></figure>
This technique has been developed more specifically for eddy-resolving
ocean assimilation models in which the signal dominates in the small
scales. Figure 6 illustrates the surface signature and vertical extension
of the three dominant modes calculated with an eddy-permitting model
in the Gulf Stream region. Theses modes reflect to some extent the
anisotropic nature of the multivariate covariances associated with the
local dynamics, with smaller scales represented by higher modes.

The local representation of the error sub-space has been shown to be
particularly effective for capturing the mesoscale features of the turbu-
lent ocean. In the example described by Penduff et al. [2002], an assim-
ilative system based on local EOFs was implemented in a 1/3° South
Atlantic OPA model to perform hindcast experiments for the 1993-1996
period. The assimilated data are similar to those shown in figure 4, con-
sisting of composite AVHRR observations of SST and altimetric mea-
surements of the sea-level anomalies. The results show that the assim-
ilation is able to successfully reproduce the Agulhas Rings present at
that time in the real ocean. In addition to correcting the variables ob-
served at the surface, the three-dimensional multivariate properties of
the EOFs also permitted a correction of the non-observed variables in
the ocean's interior. The beneficial impact of the assimilation was par-
ticularly impressive on the mean salinity in the Confluence region down
to about 1500 meters depth.

6. Low-rank Kalman filters

Several low-rank filters based on static or evolving error sub-spaces
have been developed over the past ten years, such as the SEEK filter
introduced by Pham et al. [1998] and the Reduced-Rank-SSquare-RooT
(RRSQRT) formulation explored by Verlaan and Heemink [1997]. Many
features of the Ensemble Kalman Filter (EnKF) put forward by Evensen
[1994] can be discussed using a similar framework.

The Ensemble OI scheme [Evensen, 2003] and the SEEK filter with
static error sub-space are two sub-optimal schemes which preserve a
number of important properties of statistical estimation but require only
a small fraction of the computer resources needed by the model. In con-
trast to those simplified schemes, the EnKF or the SEEK filter with
evolutive sub-space propagate the error statistics according to the model
dynamics, but they need the simultaneous integration of model states
perturbed along each direction of the error sub-space.

The consequences of using a low-rank error covariance matrix to com-
pute the forecast and analysis steps of the KF are examined in the
following sections.
## 6.1 Forecast error with a reduced basis

The reduced basis concept allows drastic simplifications to compute the evolution of the error statistics over the assimilation window. A hierarchy of algorithms of increasing sophistication and computer requirements has been proposed to compute the forecast step. Assuming that the analysis error covariance is represented as a low-rank matrix at time $t_i$

$$ \mathbf{P}_i^a = \mathbf{S}_i^a (\mathbf{S}_i^a)^T, \qquad (35) $$

where $\mathbf{S}_i^a$ (of dimension $n \times r$) defines the error sub-space associated with $\mathbf{x}_i^a$, Eq. (18) becomes

$$ \mathbf{P}_{i+1}^f = \tilde{\mathbf{S}}_{i+1}^f (\tilde{\mathbf{S}}_{i+1}^f)^T + \mathbf{Q} \quad \text{with} \quad \tilde{\mathbf{S}}_{i+1}^f = \mathbf{M}\mathbf{S}_i^a. \qquad (36) $$

The computer load associated with Eq. (36) is primarily determined by the rank $r$ of $\mathbf{P}_i^a$ which specifies the number of model integrations needed to evaluate the forecast error covariance matrix. As originally proposed by Pham et al. [1998], this algorithmic variant known as “Extended Evolutive Basis” requires the derivation of the tangent linear model $\mathbf{M}(t_i, t_{i+1})$ to update the error directions, i.e. the $r$ columns of $\mathbf{S}_i^a$. The evolved sub-space $\tilde{\mathbf{S}}_{i+1}^f$ reflects how the model dynamics affects uncertainty during the forecast.

An alternative scheme to Eq. (36) can be used for the calculation of the time evolution of the reduced basis as follows:

$$ \{\tilde{\mathbf{S}}_{i+1}^f\}_j = \frac{1}{\alpha} [M(t_i, t_{i+1}) \{\mathbf{x}_i^a + \alpha \{\mathbf{S}_i^a\}_j\} - M(t_i, t_{i+1}) \{\mathbf{x}_i^a\}] \quad (37) $$

where $\{\}_j$ is the $j$th column of the matrix, and $\alpha$ is an adjustable parameter that determines the size of the perturbations along each error direction. Equation (37) is a finite-difference approximation of the linear error evolution if $\alpha$ is small. However, the value of $\alpha$ is usually taken to be of the order of 1 to simulate the non-linear evolution of model perturbations that have an amplitude comparable to the error covariances. This algorithm known as “Interpolated Evolutive Basis” has a two-fold benefit: firstly, it avoids the computation of the tangent linear model which numerically can be a delicate task; and secondly, it seems more robust with regard to the model non-linearities because the finite difference takes into account the amplitude of the uncertainties, while the classic linearization does not.

Due to the recursive character of the Kalman filter, $\mathbf{P}_{i+1}^f$ should have the same rank as $\mathbf{P}_i^a$ in order to preserve the advantage of a low-dimensional space for the subsequent assimilation cycles. The difficulty
arises from the fact that the rank of $\mathbf{P}_{i+1}^f$ intimately depends on the structure of $\mathbf{Q}$. A possible method is simply to project the error vector $\eta_i$ and its associated covariance on the sub-space generated by the columns of $\tilde{\mathbf{S}}_{i+1}^f$. As it is often impossible to specify the model error perfectly, for the reasons discussed above, a simple parameterization of system noise can be introduced, assuming for instance that

$$ \mathbf{Q} = \frac{1 - \rho}{\rho} \tilde{\mathbf{S}}_{i+1}^f (\tilde{\mathbf{S}}_{i+1}^f)^T \quad (38) $$

where $\rho$ is a scalar quantity called the “forgetting factor” ($0 < \rho < 1$), by analogy with the approach used in automatic control algorithms. This kind of model error parameterization leads to a forecast error covariance matrix of the form:

$$ \mathbf{P}_{i+1}^f = \frac{1}{\rho} \tilde{\mathbf{S}}_{i+1}^f (\tilde{\mathbf{S}}_{i+1}^f)^T = \mathbf{S}_{i+1}^f \mathbf{S}_{i+1}^{fT} \quad \text{with} \quad \mathbf{S}_{i+1}^f = \frac{1}{\sqrt{\rho}} \tilde{\mathbf{S}}_{i+1}^f \quad (39) $$

which is singular and has the same rank $r$ as the previous analysis error covariance. The forgetting factor is one of the many possible options for accounting for some simple form of model error in the assimilation scheme. Other approaches may be implemented, however, such as using a perturbed model $\tilde{M}(t_i, t_{i+1})$ instead of the original model to dynamically update the error modes through Eq. (37). For instance, in most EnKF implementations, stochastic perturbations are introduced in the surface forcings to update each ensemble member, accounting in this way for the uncertainty in the atmospheric fluxes.

With respect to Eq. (36), an even more drastic simplification of the forecast step can be obtained by simply neglecting the dynamical transformation of the error directions during the assimilation period $(t_i, t_{i+1})$, leading to the “Fixed Basis” algorithm:

$$ \tilde{\mathbf{S}}_{i+1}^f = \mathbf{M}\mathbf{S}_i^a \approx \mathbf{I}\mathbf{S}_i^a \quad (40) $$

where $\mathbf{I}$ is the identity matrix. As in the Ensemble OI scheme [Evensen, 2003], temporal persistence of the error sub-space basis is assumed in this variant. Static error sub-spaces have been successfully used in a variety of assimilation applications (Verron et al. [1999]; Gourdeau et al. [1999]; Parent et al. [2003]; Penduff et al. [2002]), being justified by two basic arguments. The first one is of practical interest: the cost of a Fixed Basis assimilation experiment is of the same order of magnitude as a model simulation, with only a few additional computations needed to perform the algebraic operations of the analysis step. This makes the algorithm extremely useful in evaluating the overall performance of the
system, testing new parameterizations, examining the impact of different observation configurations, etc. The second argument constitutes a more objective justification in terms of statistical estimation: it is often impossible to properly characterize the statistical structure of the model error originating from imperfections in the model forcings, discretization schemes, space/time resolution, etc. As a result, the approximation made by assuming persistence of the error sub-space directions is often negligible compared to the mis-specification of the systematic errors. When the amplitude of the model error increment $Q$ dominates the difference between $P^a$ and $MP^aM^T$ in the error propagation Eq. (18), it is a waste of CPU resources to explicitly compute the error propagation with the model dynamics, because the result will eventually be polluted by mis-specified model error estimates. In such situations, it is sensible to by-pass the error propagation equation and concentrate on identifying the dominant forecast error directions using, for instance, adaptivity mechanisms as described in Section 7.

## 6.2 Analysis step with a reduced basis

The linear variance-minimizing analysis (19) can be re-formulated nicely when the forecast error covariance used to compute the Kalman gain is of low rank. Using equations (20), (39) and the matrix equality

$$[\mathbf{X}_1 + \mathbf{X}_{12}\mathbf{X}_2^{-1}\mathbf{X}_{21}]^{-1} = \mathbf{X}_1^{-1} - \mathbf{X}_1^{-1}\mathbf{X}_{12} [\mathbf{X}_2 + \mathbf{X}_{21}\mathbf{X}_1^{-1}\mathbf{X}_{12}]^{-1}\mathbf{X}_{21}\mathbf{X}_1^{-1}, \quad (41)$$

the expression for the Kalman gain, after some mathematical manipulations, can be transformed as follows:

$$\mathbf{K}_{i+1} = \mathbf{S}_{i+1}^{f} [\mathbf{I} + (\mathbf{H}\mathbf{S}_{i+1}^{f})^{T}\mathbf{R}^{-1}(\mathbf{H}\mathbf{S}_{i+1}^{f})]^{-1}(\mathbf{H}\mathbf{S}_{i+1}^{f})^{T}\mathbf{R}^{-1}. \quad (42)$$

Equation (42) shows that the size of this inversion problem is determined by the error sub-space dimension, while the original Kalman gain (20) requires an inversion in the observation space. As the number of observations is usually much larger than the rank of the error sub-space, the inversion step becomes less expensive than the corresponding computation of the original Kalman gain. By combining equations (19) and (42), the correction of the forecast state can be written as a linear combination of the error modes

$$\mathbf{x}_{i+1}^{a} - \mathbf{x}_{i+1}^{f} = \mathbf{S}_{i+1}^{f} \mathbf{c}_{i+1} \quad (43)$$

with

$$\mathbf{c}_{i+1} = [\mathbf{I} + (\mathbf{H}\mathbf{S}_{i+1}^{f})^{T}\mathbf{R}^{-1}(\mathbf{H}\mathbf{S}_{i+1}^{f})]^{-1}(\mathbf{H}\mathbf{S}_{i+1}^{f})^{T}\mathbf{R}^{-1}(\mathbf{y}_{i+1}^{o} - H\mathbf{x}_{i+1}^{f}) \quad (44)$$
The vector $c_{i+1}$ contains the $r$ amplitudes of the error modes that need to be estimated at each analysis step. Finally, the scheme evaluates the analysis error as follows

$$ \mathbf{P}_{i+1}^{a} = [\mathbf{I} - \mathbf{K}_{i+1}\mathbf{H}]\mathbf{P}_{i+1}^{f} = \mathbf{S}_{i+1}^{a}\mathbf{S}_{i+1}^{aT} \quad (45) $$

and updates the vectors of the reduced basis according to

$$ \mathbf{S}_{i+1}^{a} = \mathbf{S}_{i+1}^{f}[\mathbf{I} + (\mathbf{H}\mathbf{S}_{i+1}^{f})^{T}\mathbf{R}^{-1}(\mathbf{H}\mathbf{S}_{i+1}^{f})]^{-1/2}. \quad (46) $$

This shows that, if the rank of $\mathbf{P}_{i+1}^{f}$ is $r$, then the rank of $\mathbf{P}_{i+1}^{a}$ is equal to $r$ also; therefore, recursivity in the forecast/analysis cycles is allowed. For better numerical conditioning, it is possible to re-orthonormalize the reduced basis by recomputing a SVD decomposition of the analysis error covariance matrix.

There is one additional consideration that should be mentioned concerning the practical computation of the Kalman gain in the reduced space (42). The robust estimation of small correlations associated with remote observations is a common difficulty that can be avoided by computing EOFs with compact support, as mentioned above. Instead of computing local EOFs, a simplification of the analysis scheme can be designed by enforcing to zero the error covariances between distant variables which are believed to be uncorrelated in the real ocean [Houtekamer and Mitchell, 1998].

This simplification is implemented by assuming that distant observations have negligible influence on the analysis. The global system is split into sub-systems, and for each of these the traditional analysis is computed. Only data points located within individual regions, centered on a sub-domain of one or several grid points to be updated, actually contribute to the gain in (42). This approach can be understood as a tuning of the observation operator according to the sub-domain in question. Intuitively, this approximation makes sense because only data points located in the neighborhood of a model grid point should objectively have an impact on the analysis for that grid point. The size of the regions is determined in such a way that the distribution of the observations available on the model domain always provides at least a few data points within each region of influence (if there were no data available in the region of influence, no correction would apply). This algorithmic simplification also improves the analysis because it enables a larger part of the estimation space to be spanned over a particular subdomain (see, for instance, Brusdal et al. [2003]).

Other interesting procedures exist for addressing the locality issue. One of these is the partitioning of the large estimation problem into
a series of separate smaller calculations using the partitioned Kalman filter and smoother proposed by Fukumori [2002]. The objective of the partition is to make global eddy-resolving data assimilation problems computationally viable. In their example, the reduced state consists of perturbations of the barotropic mode and the first five baroclinic modes defined in eight overlapping cells covering the globe.

## 6.3 Stochastic vs. deterministic filters

In a low-rank filter like SEEK, the error directions are determined at the initialization step and their evolution can be anticipated in a deterministic way. As these filters do not require any further randomization, they are termed deterministic filters in contrast to the EnKF which can be considered as a stochastic filter because randomization has to be repeated at each assimilation cycle.

Evensen [2003] reviews the theoretical formulation and practical implementation of the EnKF. The EnKF was introduced to avoid the occurrence of instabilities found with the Extended KF due to the non-linear evolution of the probability density functions [Evensen 1994]. However, the EnKF only solves half of the non-linearity problems because it still combines the model prediction and the data by using only the first two moments of the pdfs, assuming that the distributions are nearly Gaussian. The non-linear analysis equations would be more difficult to use in practical applications, as discussed by Evensen and van Leeuwen [2000].

The EnKF was initially proposed by Evensen [1994] based on the following general procedure: a sampling of the state space is achieved using Monte-Carlo methods to generate an ensemble of $r$ model states $\mathbf{x}_i^{a,j}$ representing the spread of possible initial conditions at time $t_i$ around the mean $\overline{\mathbf{x}}_i^{a,j}$. Each member is then propagated individually as

$$ \mathbf{x}_{i+1}^{f,j} = M(t_i, t_{i+1}) \{\mathbf{x}_i^{a,j}\} + \eta^j, \quad j = 1, \dots, r \quad (47) $$

using the non-linear model with stochastic perturbations. The vectors $\eta^j$ have a covariance matrix $\mathbf{Q}$ and are generally introduced through the perturbation of the atmospheric forcings. The new ensemble $\mathbf{x}_{i+1}^{f,j}$ provides implicitly the forecast error covariance matrix

$$ \mathbf{P}_{i+1}^f = \frac{1}{r-1} \sum_{j=1}^{r} \left( \mathbf{x}_{i+1}^{f,j} - \overline{\mathbf{x}}_{i+1}^{f,j} \right) \left( \mathbf{x}_{i+1}^{f,j} - \overline{\mathbf{x}}_{i+1}^{f,j} \right)^T . \quad (48) $$
The statistical analysis equation of the KF is then applied to each individual member,

$$ \mathbf{x}_{i+1}^{a,j} = \mathbf{x}_{i+1}^{f,j} + \mathbf{K}_{i+1} (\tilde{\mathbf{y}}_{i+1} - \mathbf{H} \{\mathbf{x}_{i+1}^{f,j}\}) , \quad j = 1, \dots, r \quad (49) $$

providing implicitly the analysis error covariance matrix

$$ \mathbf{P}_{i+1}^a = \frac{1}{r-1} \sum_{j=1}^{r} \left( \mathbf{x}_{i+1}^{a,j} - \overline{\mathbf{x}_{i+1}^{a,j}} \right) \left( \mathbf{x}_{i+1}^{a,j} - \overline{\mathbf{x}_{i+1}^{a,j}} \right)^T . \quad (50) $$

In order to avoid the problem of systematic underestimation of the analysis covariance that occurs when the same data and the same gain are used in the set of analysis equations, an ensemble of perturbed observations $\tilde{\mathbf{y}}_{i+1}$ is considered in (49) instead of the original $\mathbf{y}_{i+1}$ [Houtekamer and Mitchell, 1998; Burgers et al., 1998]. The stochastic nature of the EnKF filter arises as a consequence of using perturbed observations [Lawson and Hansen, 2004]. Deterministic variants of the EnKF, which do not require perturbed observations, have been proposed recently involving square-root analysis schemes [e.g. Whitaker and Hamill, 2002].

An interpretation of the error propagation scheme in the SEEK filter has been proposed by Ballabrera et al. [2001] in terms of ensemble model integrations. Indeed, Eq. (37) depicts the natural dispersion of an ensemble of different model trajectories initialized in the vicinity of the initial guess; it represents the amplification of unstable error modes (or the damping of stable modes) inherent in the system's dynamics. A more exhaustive exploration of the similarities and differences between the SEEK and Ensemble Kalman Filter (EnKF) is discussed in Brusdal et al. [2003] and Nerger [2004].

## 7. Statistical consistency of assimilation schemes

A key issue concerning statistical assimilation systems relates to their capacity to produce reliable error statistics about the ocean state estimates, and to propagate those error statistics properly from one assimilation cycle to the next. In the linear KF, the specification of the system noise $\mathbf{Q}$, the observation error $\mathbf{R}$ and the background error covariance at the initial time $\mathbf{P}_0$ perfectly determines the subsequent evolution of the error statistics throughout the assimilation sequence. This is because the observations actually control the trajectory of the model state, but they have no impact on the evolution of the error statistics themselves (except through the observation network $\mathbf{H}$). For a KF to yield optimal performances, it is necessary to provide the correct a priori description of these error covariance matrices. As only guesses of these quantities
are available in real cases, it is necessary to verify the consistency of the underlying assumptions.

Comparisons with independent information provide, of course, the ideal way to assess the quality of the assimilation and to detect what has possibly gone wrong in the system. Testing the ability to produce a forecast from an analyzed state (i.e., using data that have yet to be assimilated) is also very helpful in checking whether the dynamical model and assimilation scheme are consistent with one another. The most common test is to compare the forecast to persistence [De Mey et al., 2002]. Additionally, information acquired during system operation can be used to verify if the prior error statistics have been prescribed in a way which represents the actual errors in the model and the observations. The difference between the observations y used in the assimilation system and the forecast or analysis fields (called innovations d and residuals r, respectively) provide two series of data that contain a wealth of important information about the consistency of the prior error estimates. It can be shown that examining diagnostics on the innovation vector is essentially equivalent to examining them on the residual vector, but in practice some features of the analysis may be easier to diagnose from one or the other vector [Talagrand, 1999].

In the section below, we will examine a number of useful criteria that can be diagnosed from the sequence of innovations and residuals to detect imperfections in the specification of the error statistics. Since such imperfections can be a source of drift in the assimilation system or, even more dramatically, can affect the stability of the filter, we will introduce the concept of adaptivity, the aim of which is to enforce consistency between the error statistics predicted by the filter and the observed misfits.

## 7.1 Verification of statistical consistency

Simple diagnostics can be implemented quite easily in assimilation systems, providing interesting tools to monitor overall performance and detect anomalies in system operations. A first-order check can be made by computing the mean innovation, which is expected to vanish over a sufficiently long assimilation sequence:

$$ \overline{d}_i = \overline{y}_i - \overline{Hx}_i^f = \overline{Hx}_i^t + \epsilon_i^o - \overline{H}(\overline{x}_i^t + \epsilon_i^f) = \overline{\epsilon}_i^o - \overline{H\epsilon}_i^f = 0 \quad (51) $$

A non-centered innovation sequence is the obvious indication of biases in the model and/or observations which in principle should be removed from the assimilation system. If the source of bias is in the model, $\overline{d}_i$ only provides information about the projection of this bias in the observation
space, and the difficulty remains of inversing this information back to the model space.

Similarly, the expectation of the residual $\overline{\mathbf{r}}_i = \overline{\mathbf{y}}_i - \overline{\mathbf{H}\mathbf{x}}_i^a$ should be zero for an optimal system as well as the mean increment $\overline{\mathbf{H}}(\mathbf{x}_i^a - \mathbf{x}_i^f) = \overline{\mathbf{d}}_i - \overline{\mathbf{r}}_i$ [Talagrand, 1999]. Figure 7 illustrates the mean assimilation increment of sea-surface height computed from a 10-year analysis sequence of the MERCATOR global prototype (at 2° × 2° horizontal resolution) assimilating altimeter data during the 1993-2003 period [Ferry et al., 2005]. Over large portions of the ocean, the amplitude of the increment is indeed very small, but in some regions the assimilation system systematically corrects the predicted sea-surface pressure towards higher (e.g. in the southern oceans) or lower values (e.g. in the sub-polar gyre). Ferry et al. [2005] discuss how the detection of such biases can be used to improve the modelling or assimilation components of the MERCATOR system.

<figure><img src="image_8.png" /><figcaption>Figure 7. Mean SLA increment (in cm) diagnosed from the global MERCATOR prototype assimilating along-track altimeter data (according to Ferry et al. [2005]).</figcaption></figure>

By taking the covariance of the innovation and remembering the assumptions made in Section 2 about observation and forecast errors, we obtain

$$ \overline{\mathbf{d}_i \mathbf{d}_i^T} = (\epsilon_i^o - \mathbf{H}\epsilon_i^f)(\epsilon_i^o - \mathbf{H}\epsilon_i^f)^T = \mathbf{R} + \mathbf{H}\mathbf{P}_i^f \mathbf{H}^T \quad (52) $$

if the error covariances are correctly specified. The comparison between the matrix $\overline{\mathbf{d}_i \mathbf{d}_i^T}$ and the sum of the observation and forecast error covariances used by the assimilation scheme indicates whether the forecast misfits “seen” by the filter are compatible with the prior information.
As all these quantities are full covariance matrices in the data space, the comparison is often made by looking only at a few selected diagnostics such as the trace of the matrix or its diagonal elements. After simple manipulations, it can be shown using the residuals that

$$ \overline{\mathbf{r}_i \mathbf{r}_i^T} = \mathbf{R} - \mathbf{H} \mathbf{P}_i^a \mathbf{H}^T \quad (53) $$

This formula accounts for the fact that the same observation error affects both the data and the analysis. It shows that, for an asymptotically perfect estimation system (with $\mathbf{P}_i^a \sim 0$), the residual error covariance converges towards the observation error.

Another more synthetic diagnostic can be implemented in both statistical and variational assimilation systems. If the covariances are correctly estimated, the scalar quantity

$$ J_i = \mathbf{d}_i^T (\mathbf{H} \mathbf{P}_i^f \mathbf{H}^T + \mathbf{R})^{-1} \mathbf{d}_i \quad (54) $$

behaves as a chi-squared variable with as many degrees of freedom as there are independent data (say, $p$) [Bennett, 1992]. It can be regarded as a particular norm (the so-called Mahalanobis norm) of the innovation vector. Significance tests may therefore be used to accept or reject the covariance models. In the context of low-rank KFs based on the decomposition (39), the norm (54) can be written as

$$ J_i = \mathbf{d}_i^T \left( (\mathbf{H}\mathbf{S}_i^f)(\mathbf{H}\mathbf{S}_i^f)^T + \mathbf{R} \right)^{-1} \mathbf{d}_i \quad (55) $$

and can be used to objectively evaluate the suitability of the error subspace $\mathbf{S}_i^f$. Instead of testing the complete $\chi_p^2$ behaviour of $J_i$, only the first and second statistical moments need to be examined. These should show that

$$ \overline{J_i} = p \quad \text{and} \quad \text{Var}(J_i) = p \quad (56) $$

In particular, a too small (resp. large) value of $\overline{J_i}$ is symptomatic of an overestimated (resp. underestimated) innovation amplitude, which may be the consequence of too large (resp. small) observation or forecast error covariances.

## 7.2 Adaptivity

In most applied assimilation systems, large deviations may be observed with respect to the theoretical criteria (51), (52), (53) or (56), reflecting some flaws in the prior statistical assumptions. Feedback mechanisms can be implemented on-line or off-line to adjust the prior error statistics in such a way as to restore consistency between the errors diagnosed by the filter and the innovation or residual information. Inherent
in the family of techniques used for this purpose is the concept of adaptivity. A comparison between different sophisticated adaptive KFs is provided by Blanchet et al. [1997]. In this section, we examine more closely different examples of adaptive methods that have been explored in the literature to identify and correct model biases, to tune the parameterization of model errors and to build the error sub-space of low-rank KFs.

Research efforts aimed at improving error covariance modelling in assimilation systems are of limited interest if biases are left in the models or the observations, i.e. if criterion (51) is not verified. The correction for biases in operational systems is expected to have a very strong impact on assimilation performances. In the context of sequential assimilation, Dee and Da Silva [1998] proposed a rigorous method for estimating the forecast bias and correcting the forecast prior to the analysis, assuming unbiased observations. The algorithm is designed to perform on-line and its implementation does not require substantial modifications to the assimilation system. The basic idea consists in running a simplified KF-type algorithm to estimate $\overline{d}_i$ in addition to the KF for the state itself. Assuming temporal persistence of the bias, the extra cost of the method is equivalent to one additional computation of the statistical analysis step.

The question of model error parameterization can been addressed by means of adaptive methods too. In many practical studies, the model error covariance $Q$ is probably the least well known statistical quantity impacting the forecast error. Mitchell and Houtekamer [2000] developed an adaptive EnKF using a maximum likelihood method to estimate the parameterization of model errors from the innovation sequence. The approach recently developed by Brankart et al. [2003] can also be viewed as an adaptive parameterization of the model error. By adjusting the $\rho$ coefficient of Eq. (38) according to the local innovation variance, the method is able to account for regional properties of the ocean dynamics. The estimation of the innovation variance is based on a weighted average of the latest innovations, using a weight that decreases exponentially with the age of the innovation.

This mechanism was explored in the context of hindcast experiments conducted for the 1993-1996 period. The SEEK filter was implemented in two different models : the 1/3° North Atlantic OPA model of the MERCATOR prototype system, and the Atlantic/Arctic MICOM model of the European DIADEM system. Sea-surface temperature from the NASA Pathfinder project and altimetric data from the Topex/Poseidon and ERS missions were assimilated in both systems every 10 days. Testut [2000] studied the distribution of the 10-day forecast error in the
OPA system once the assimilation experiment had reached an asymp-
totic regime (figure 8). The largest forecast errors are found along the
Gulf Stream path between Cape Hatteras and 40°N, where the standard
deviation exceeds 20 cm in some places. Local maxima are also detected
in the North Atlantic Current extension and along the Azores Current
at 35°N. Using the MICOM system, Brankart et al. [2003] compared
the estimated innovation variance with the sum of the observation and
forecast error variance diagnosed by the filter in a zonal section across
the Gulf Stream region at 30°N. The comparison indicates a fairly good
agreement between the 2 estimates in most regions, reflecting a success-
ful adaptive procedure.

<figure><img src="image_9.png" /><figcaption>*Figure 8.* Distribution of the 10-day forecast error diagnosed with the adaptive SEEK filter implemented in the OPA model (reproduced from Testut [2000]).</figcaption></figure>

Another interesting application of the adaptivity concept is the up-
date of the initial structure of the error sub-spaces. Using a low-rank
KF with an idealized primitive-equation model of the Gulf Stream, the
assimilation error dynamics was investigated by *Ballabrera et al.* [2001].
They examined how the filter performances were affected by imperfect
specifications of the initial error statistics. They observed that the trun-
cation error does not impact the control of the solution when at least one
of the following three conditions is met: (i) the initial error is perfectly
described by the reduced basis; (ii) the truncation error is dynamically
uncoupled with the error components of the low-dimension space; (iii)
the sub-space must contain all the components of the growing error
during the forecast time period. As these three conditions are never
perfectly verified in realistic assimilation systems, an adaptive proce-
dure was developed in order to introduce some feedback between the
data and the error sub-space used by the filter. The algorithm proposed
by Brasseur et al. [1999] updates the error sub-space of the SEEK fil-
ter along the geostrophic attractor by extracting information left in the
residual vector after each analysis step. This update of the reduced basis
is performed in such a way as to attenuate the truncation error and to
improve the projection of the next innovation onto the error sub-space:
this leads to the concept of adaptive sub-space. An advantage of this
variant of the SEEK filter is that it allows the evolution of the error sub-
space without incurring the cost of propagating the whole set of error
directions dynamically.

8. Intermittent vs. time-continuous filtering

In the basic assimilation problem introduced in Section 1, two major
simplifications were considered: (i) the observations were available at
discrete time intervals, and (ii) the analysis was performed at the exact
time of the measurements. In oceanographic and atmospheric applica-
tions, the situation is quite different since the flow of observations can
be considered as almost continuous in time (for instance, the acquisition
of altimeter data). Therefore, ocean data assimilation with intermittent
sequential filters necessarily involves approximations.

In principle, sequential filters could perform an analysis step as of-
ten as a new piece of information arrives. Time-continuous formulations
of the KF exist [Gelb 1974] and have been applied to analogic signal
processing. However, their application to oceanographic or atmospheric
models is inappropriate: for practical reasons, it would be impossible
to incorporate the data at their exact time of acquisition. Experience
shows that it is necessary to accumulate a certain number of observa-
tions between two successive analysis steps to correct the ocean state
with sufficient impact. Besides, operational assimilation systems must
be scheduled on a regular temporal basis so as to avoid unnecessary
algorithmic complications, to account for human intervention, delay in
data delivery, etc.
## 8.1 Computing innovations using FGAT

Typical lengths of assimilation cycles are 3 to 7 days for mesoscale ocean current predictions, and 10 to 30 days for initialization of the oceanic component of seasonal climate predictions. In spite of the fact that the ocean cannot be considered as static over these time scales, intermittent sequential filters incorporate at one single instant the set of observations collected during the whole assimilation interval. This is a major difference with 4D-VAR algorithms, which can take full advantage of the temporal distribution of the data within an assimilation window.

Fairly simple solutions can be set up to alleviate these problems in the context of statistical filters. For example, the FGAT (First Guess at Appropriate Time) method initially introduced in meteorology can be used to evaluate the innovation vector more correctly: instead of computing the difference between the time-distributed data set and the model forecast at time $t_{i+1}$ as in Eq. (19), the innovation is evaluated “on the flight” by cumulating the differences between each piece of observation and the corresponding element of the model forecast at the measurement time. This approach has also been used with 3D-VAR assimilation systems [Weaver et al., 2003] to benefit from the temporal dimension. Due to the fast propagation of equatorial waves, the FGAT feature may be particularly important in the tropical oceans.

## 8.2 Incremental Analysis Update

A direct consequence of intermittency is the discontinuity of the forecast/analysis estimates, which is recognized as a major drawback of sequential methods. Two related problems, - shocks to the model and data rejection -, arise with intermittent corrections. It is found that observations assimilated into models may introduce transient waves excited by the impulsive insertion. These waves are often the result of imperfections in the corrected state associated with physically unbalanced error covariances. In the example illustrated in figure 9, six vertical profiles of temperature and salinity measurements are assimilated into the HYCOM model using the SEEK filter. The profiles are inserted in distant regions (Labrador Basin, Irminger Sea, Gulf Stream, Azores and North Brazil currents and Caribbean Sea) in order to avoid mutual interference. The corrected state is then integrated using the model with realistic forcings for one month after analysis time. Figure 9 depicts the SSH increment after 3 days of simulation, showing the occurrence of spurious transients (particularly in the Gulf Stream region) which the model generates to dynamically adjust the new state.
<figure><img src="image_10.png" /><figcaption>Figure 9. SSH increment (in m) obtained with the HYCOM model 3 days after assimilation of vertical T/S profiles.</figcaption></figure>

In order to incorporate analysis increments in a more gradual manner, a new algorithm based on Incremental Analysis Updates (IAU) has been proposed by Bloom et al. [1996], which combines aspects of intermittent and continuous assimilation schemes. Using the regular KF equations, the IAU algorithm first computes the analysis correction ($x_{i+1}^a - x_{i+1}^f$). This correction is then distributed (uniformly or not) over the assimilation window and inserted gradually into the model. The state obtained at the end of the assimilation window can be used as the initial conditions for the next assimilation cycle, leading to time-continuous filtered trajectories. The concept of IAU can also be implemented by computing corrections to the atmospheric forcing or by introducing dynamical relaxation terms in the prognostic equations, according to the analysis increments. The IAU temporal strategy can be complemented by the FGAT scheme which computes the innovation “on the flight”.

A possible implementation of the IAU scheme in primitive equation models may include the following steps: (i) first-guess model integration and evaluation of the innovation vector using FGAT; (ii) computation of the analysis increment using the Kalman gain; (iii) second integration on ($t_i, t_{i+1}$) using primitive equations modified by the increment. The additional cost incurred by this process will be one model integration or less, depending on the details of the algorithmic implementation. Noting by $\delta T$ the increment of the temperature field, a possible modification of
the PE equation for the temperature may be written as follows:

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla_h T + w \frac{\partial T}{\partial z} = D_h(T) + \frac{\partial}{\partial z} (\tilde{\lambda} \frac{\partial T}{\partial z}) + \frac{\delta T}{(t_{i+1} - t_i)} \quad (57) $$

where $\delta T/(t_{i+1} - t_i)$ is a new term acting as a body force, in a similar way to the nudging assimilation technique [Verron and Holland, 1989; Blayo et al., 1994]. The extra term is, however, different from newtonian relaxation since it does not explicitly involve the current value of T and has the advantage of being multivariate in general and weighted consistently with error statistics. A linear analysis of the IAU procedure shows that it has the properties of a low-pass temporal filter. Variants of the IAU procedure are currently being explored with the MERCATOR assimilation systems.

## 9. Conclusions

The Kalman filter provides a theoretical framework from which a hierarchy of algorithms of increasing sophistication and increasing computer requirements can be derived. These algorithms range from Optimal Interpolation schemes, which require only a few percent of the computer resources allocated to run the model, to non-linear Kalman filters such as the SEEK or the EnKF, which require the simultaneous integration of perturbed model trajectories in equal number to the dimension of the error sub-space. The research conducted in support of operational oceanography has shown that error statistics is at the heart of applied assimilation systems and remains a major challenge for ocean data assimilation. In this respect, the derivation of a hierarchy of simplified filters offers a unique possibility to test different approaches for specifying and calibrating the background, systematic and observation error statistics.

A salient feature of statistical estimation methods is their multivariate nature: observations related to several different model variables (e.g., SST, SLA, SSS, in situ temperature and salinity data) are used to correct the whole model state in a consistent manner. This opens promising avenues for effectively envisaging the simultaneous assimilation of in situ data with satellite data from different sensors. Error bars on state estimates can be computed by the algorithms on a rigorous basis, making the schemes useful in qualifying the reliability of different forecasts, or comparing the relevance of different observation systems. For practical implementations, the adjoint of the model code is not necessarily needed as with 4D-VAR, and the basic architecture of the algorithms is modular. The transition from one model version to another, or from one model code to another can be made smoothly without much recoding.
The large variety of sequential statistical methods developed in the context of scientific or operational applications is an indication that a given assimilation technique cannot be considered as a “plug-and-play” system, capable of solving universal problems. To design the best possible assimilation system, it is necessary to clearly define the purpose of data assimilation (forecast initialization, reanalysis, model error detection, etc.), the physical characteristics of the processes of interest, the sampling properties of the observation systems, and the limitations of the assimilation techniques. Moving towards hybrid sequential/variational methods is probably an astute way of taking the best from both approaches. Due to the complexity of models and algorithms, the success of an assimilation system also depends on the community of people working within a common framework. In the future, the sharing of generic assimilation tools between operational and research teams should accelerate the progress of the methods and provide feedback from more intensive utilizations.

A number of important issues for ocean data assimilation have yet to be fully resolved, such as the incorporation of inequality constraints in statistical estimation algorithms. Such constraints are needed, e.g. in the context of tracer data assimilation to maintain the positiveness of tracer concentrations, or for the proper handling of the physical properties of the water column such as static stability. Another important challenge will be to further develop assimilation systems for coupled physical-biological models, with the aim of demonstrating the capacity to estimate and forecast marine ecosystems and biogeochemical variables in the ocean on a routine basis. As a first step, the incorporation of biogeochemical models into the physical assimilative systems developed within GODAE will provide new “bio-metrics” to evaluate the representation of the physical processes of critical importance to biology. A further challenge will be to implement suitable methods to assimilate ocean colour data in the coupled assimilative systems.

## Acknowledgements

The author gratefully acknowledges his colleagues from LEGI and MERCATOR Ocean for numerous discussions and valuable suggestions. He is supported by the Centre National de la Recherche Scientifique. Computations were carried out at the IDRIS/CNRS computing centre. The contribution of students from the GODAE school who commented an earlier version of these lecture notes has been much appreciated. This work has been partly supported by the MERSEA project of the European Commission under Contract SIP3-CT-2003-502885.
References

Auclair, F., P. Marsaleix and P. De Mey, 2003 : Space-time structure and dynamics of the forecast error in a coastal circulation model of the Gulf of Lions. *Dyn. Atmos. Oceans*, **36**, 309-346.

Ballabrera J., P. Brasseur and J. Verron, 2001: Dynamical evolution of the error statistics with the SEEK filter to assimilate altimetric data in eddy-resolving ocean models. *Quart. J. Roy. Meteor. Soc.*, **127**, 233-253.

Bennett A., 1992: *Inverse methods in physical oceanography*. Cambridge Monographs on Mechanics and Applied Mathematics, Cambridge University Press, 346 pp.

Birol F., J.M. Brankart, F. Castruccio, P. Brasseur and J. Verron, 2004: Impact of the Ocean Mean Dynamic Topography on Satellite Data Assimilation, *Mar. Geodesy*, **27**, 59-78.

Blanchet I., C. Frankignoul and M.A. Cane, 1997: A comparison of adaptative Kalman filters for a tropical Pacific ocean model. *Mon. Wea. Rev.*, **125**, 40-58.

Blayo E., J. Verron and J.M. Molines, 1994: Assimilation of Topex/Poseidon altimeter data into a circulation model of the North Atlantic. *J. Geophys. Res.*, **99** (C12), 24691-24705.

Bleck R., 2002: An oceanic general circulation model framed in hybrid isopycnic-cartesian coordinates. *Ocean Modelling*, **4**, 55-88.

Bloom S.C., L.L. Takacs, A.M. Da Silva and D. Ledvina, 1996: Data Assimilation Using Incremental Analysis Updates, *Mon. Wea. Rev.*, **124**, 1256-1271.

Brankart J.-M., C.E. Testut, P. Brasseur and J. Verron, 2003: Implementation of a multivariate data assimilation scheme for isopycnic coordinate ocean models: Application to a 1993-96 hindcast of the North Atlantic Ocean circulation, *J. Geophys. Res.*, **108**, 3074, doi: 10.1029/2001JC001198.

Brasseur, P., J. Ballabrera and J. Verron, 1999: Assimilation of altimetric data in the mid-latitude oceans using the Singular Evolutive Extended Kalman filter with an eddy-resolving, primitive equation model. *J. Mar. Syst.*, **22** (4), 269-294.

Brusdal K., J.M. Brankart, G. Halberstadt, G. Evensen, P. Brasseur, P.J. van Leeuwen, E. Dombrowsky and J. Verron, 2003: A demonstration of ensemble-based assimilation methods with a layered OCGM from the perspective of operational ocean forecasting systems, *J. Mar. Syst.*, **40-41**, 253-289.

Burgers G., P. van Leeuwen and G. Evensen, 1998: Analysis scheme in the ensemble Kalman filter, *Mon. Wea. Rev.*, **126**, 1719-1724.

Cane M.A., A. Kaplan, R.N. Miller, B. Tang, E.C. Hackert and A.J. Busalacchi, 1996: Mapping tropical Pacific sea level: Data assimilation via a reduced state Kalman filter. *J. Geophys. Res.*, **101**, C10, 22,599-22,617.

Carme S., Pham D. T., et J. Verron, 2001: Improving the Singular Evolutive Extended Kalman filter for strongly non-linear models for use in ocean data assimilation. *Inverse Problems*, **17**, 1535-1559.

Carmillet V., J. M. Brankart, P. Brasseur, H. Drange, G. Evensen et J. Verron, 2001: Assimilation of ocean color data into an ecosystem model of the North Atlantic ocean. *Ocean Modelling*, **3**, 167-192.

Cohn S.E. and R. Todling, 1996: Approximate data assimilation schemes for stable and unstable dynamics, *J. Meteo. Soc. Japan*, **74**, 63-75.

Cooper M. and K. Haines, 1996: Altimetric assimilation with water property conservation. *J. Geophys. Res.*, **101** (C1), 1059-1077.

Courtier P., 1997: Variational Methods, *J. Met. Soc. Japan*, **75**, 211-218.
Daley R., 1991: *Atmospheric Data Analysis*, Cambridge Atmospheric and Space Science Series, Cambridge University Press, 457 pp.

Dee D.P., 1991: Simplification of the Kalman filter for meteorological data assimilation, *Quart. J.Roy. Meteor. Soc.*, **117**, 365-384.

Dee D.P and Da Silva A., 1998: Data assimilation in the presence of forecast bias, *Quart. J. Roy. Meteor. Soc.*, **124**, 269-295.

De Mey P. and M. Benkiran, 2002: A multivariate reduced-order optimal interpolation method and its application to the Mediterranean basin-scale circulation. In: *Ocean Forecasting: Conceptual basis and applications*, N. Pinardi and J.D. Woods, Eds, Springer-Verlag, Belin, Heidelberg, New York, 472 pp.

De Mey P., T. Awaji, M. Bell, A.F. Bennett, P. Brasseur, G. Evensen, K. Haines, I. Fukumori, O.M. Smedstad, D. Stammer and A. Weaver, 2002: Approaches to Data Assimilation within GODAE. *Proceedings of the international symposium "En route to GODAE"*.

Echevin V., P. De Mey and G. Evensen, 2000: Horizontal and Vertical Structure of the Representer Functions for Sea Surface Measurements in a Coastal Circulation Model. *J. Phys. Oceanogr.*, **30**, 2627-2635.

Emery W.J. and R.E. Thomson, 1998: *Data analysis methods in physical oceanography*, Pergamon, 634 pp.

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

Evensen G. and P.J. van Leeuwen, 2000: An Ensemble Kalman Smoother for nonlinear Dynamics, *Mon. Wea. Rev.*, **128**, 1852-1867.

Evensen G., 2003: The Ensemble Kalman Filter: theoretical formulation and practical implementation. *Ocean Dynamics*, **53**(4), 343-367.

Ferry N., E. Remy, P. Brasseur and C. Maes, 2005: The MERCATOR global ocean operational analysis/forecast system: assessment and validation of a 11-year reanalysis, submitted to *J. Mar. Syst*.

Fu L.L., I. Fukumori and R.N. Miller, 1993: Fitting dynamic models to Geosat sea-level observations in the tropical Pacific ocean. Part II: A linear, wind-driven model. *J. Phys. Oceanogr.*, **23**, 2162-2181.

Fukumori I., J. Benveniste, C. Wunsch and D.B. Haidvogel, 1993: Assimilation of sea surface topography into an ocean circulation model using a steady state smoother. *J. Phys. Oceanogr.*, **23**, 1831-1855.

Fukumori I., 1995: Assimilation of Topex sea level measurements with a reduced-gravity, shallow water model of the tropical Pacific ocean. *J. Geophys. Res.*, **100** (C12), 25,027-25,039.

Fukumori I. and P. Malanotte-Rizzoli, 1995: An approximate Kalman filter for ocean data assimilation: an example with an idealized Gulf Stream model. *J. Geophys. Res.*, **100** (C4), 6,777-6,793.

Fukumori I., 2001: Data Assimilation by Models. In: *Satellite Altimetry and Earth Sciences*, L.L. Fu and A. Cazenave Eds., International Geophysics Series, Vol. **69**, Academic Press, 463 pp.

Fukumori I., 2002: A partitioned Kalman Filter and Smoother, *Mon. Wea. Rev.*, **130**, 1370-1383.

Gavart M., P. De Mey and G. Caniaux, 1999: Assimilation of satellite altimeter data in a primitive-equation model of the Azores-Madeira region, *Dyn. Atmos. Oceans*, **29**, 217-254.

Gelb A., 1974: *Applied optimal estimation*. MIT Press, Cambridge, MA, 374 pp.
Gourdeau L., J. Verron, T. Delcroix, A.J. Busalacchi et R. Murtugudde, 1999 : Assimilation of Topex/Poseidon altimetric data in a primitive equation model of the tropical Pacific ocean, 1992-1996. *J. Geophys. Res.*, **105** (C4), 8473-8488.

Hurlburt H.E., 1986: Dynamic transfer of simulated altimetric data into subsurface information by a numerical ocean model. *J. Geophys. Res.*, **91** (C2), 2372-2400.

Houtekamer P.L. and H.L. Mitchell, 1998: Data Assimilation Using an Ensemble Kalman Filter Technique, *Mon. Wea. Rev.*, **126**, 796-811.

Ide K., A.F. Bennett, P. Courtier, M. Ghil and A.C. Lorenc, 1997: Unified notation for data assimilation: operational, sequential and variational. *J. Meteor. Soc. Japan*, **75**, 181-189.

Jazwinski A.H., 1970: *Stochastic Processes and Filtering Theory*. Academic Press, 376 pp.

Kalnay E., D.L.T. Anderson, A.F. Bennett, A.J. Busalacchi, S.E. Cohn, P. Courtier, J. Derber, A.C. Lorenc, D. Parrish, J. Purser, N. Sato and T. Schlatter, 1997: Data assimilation in the ocean and in the atmosphere: what should be next ? *J. Meteor. Soc. Japan*, **75**, 489-496.

Kalnay E., 2003: *Atmospheric Modeling, Data Assimilation and Predictability*, Cambridge University Press, 341 pp.

Kelly K., 1988: Comments on "Empirical Orthogonal Function Analysis of Advanced Very High Resolution Radiometer Surface Temperature Patterns in Santa Barbara Channel" by G.S.E. Lagerloef and R.L. Bernstein, *J. Geophys. Res.*, **93**, 15753-15754.

Lawson W.G. and J.A. Hansen, 2004: Implications of Stochastic and Deterministic Filters as Ensemble-Based Data Assimilation Methods in Varying Regimes of Error Growth, *Mon. Wea. Rev.*, **132**, 1966-1981.

Lermusiaux P.F.J., 1999: Data assimilation via error subspace statistical estimation: Part II. Middle Atlantic Bight shelfbreak front simulations and ESSE validation. *Mon. Wea. Rev.*, **127**, 1408-1432.

Lermusiaux P.F.J., 2001: Evolving the subspace of the three-dimensional multiscale ocean variability: Massachusetts Bay. *J. Mar. Syst.*, **29**, 385-422.

Lions J.L., O.P. Manley, R. Temam and S. Wang, 1997: Physical Interpretation of the Attractor Dimension for the Primitive Equations of Atmospheric Circulation, *J. Atmos. Sci.*, **54**, 1137-1143.

Mellor G.L. and T. Ezer, 1991: A Gulf Stream model and an Assimilation Scheme, *J. Geophys. Res.*, **96**, 8779-8795.

Miller R.N., 1989: A Kalman filter analysis of sea level height in the tropical Pacific, *J. Phys. Oceanogr.*, **19**, 773-790.

Miller R.N. and Ehret L., 2002: Ensemble Generation for Models of Multimodal Systems, *Mon. Wea. Rev.*, **130**, 2313-2333.

Mitchell H.L. and P.L. Houtekamer, 2000: An Adaptive Ensemble Kalman Filter, *Mon. Wea. Rev.*, **128**, 416-433.

Morrison D.F., 1988: *Multivariate Statistical Methods, Probability and Statistics Series*, McGraw-Hill, 415 pp.

Nerger L., 2004: Parallel Filter Algorithms for Data Assimilation in Oceanography, *Reports on Polar and Marine Research*, Alfred Wegener Institute for Polar and Marine Research, Germany, Volume **487**, 175 pp.

Oschlies A. and J. Willebrand, 1996: Assimilation of Geosat altimeter data into an eddy-resolving primitive equation model of the North Atlantic ocean. *J. Geophys. Res.*, **101**, 14175-14190.
Ott E., B.R. Hunt, I. Szunyogh, A.V. Zimin, E.J. Kostelich, M. Corazza, E. Kalnay, D.J. Patil and J.A. Yorke, 2004: A local ensemble Kalman filter for atmospheric data assimilation, *Tellus*, **56A**, 415-428.

Parent L., 2000 : *Assimilation de données dans l'océan Pacifique Tropical sur la période 1994-1998*, Ph. D. dissertation, Univ. J. Fourier, Grenoble, 222 pp.

Parent L., C.E. Testut, J.M. Brankart, J. Verron, P. Brasseur and L. Gourdeau, 2003 : Comparative assimilation of Topex/poseidon and ERS altimetric data and of TAO temperature data in the tropical Pacific Ocean during 1994-1998, and the mean sea-surface height issue, *J. Mar. Syst.*, **40-41**, 381-401.

Penduff Th., P. Brasseur, C.E. Testut, B. Barnier and J. Verron, 2002: Assimilation of sea-surface temperature and altimetric data in the South Atlantic Ocean : impact on basin-scale properties, *J. Mar. Res.*, **60**, 805-833.

Pham D. T., J. Verron and M.C. Roubaud, 1998 : A singular evolutive extended Kalman filter for data assimilation in oceanography. *J. Mar. Syst.*, **16**, 323-340.

Pinardi N., A. Rosati and R. Pacanowski, 1995: The sea-surface pressure formulation of rigid lid models. Implications for altimetric data assimilation studies. *J. Mar. Syst.*, **6**, 109-120.

Rienecker M. and R. Miller, 1991: Ocean Data Assimilation Using Optimal Interpolation with a Quasi-Geostrophic Model, *J. Geophys. Res.*, **96** (C8), 15093-15103.

Robert C., S. Durbiano, E. Blayo and J. Verron, 2005: A reduced order strategy for 4D-Var ocean data assimilation, *J. Mar. Syst.*, accepted.

Talagrand O., 1999: A posteriori verification of analysis and assimilation algorithms, *Diagnosis of Data Assimilation Systems Workshop Proceedings*, ECMWF, Reading U.K., November 1998, 17-28.

Testut C.E., 2000 : *Assimilation de données satellitales avec un filtre de Kalman de rang réduit dans un modèle aux Equations Primitives de l'océan Atlantique*, Ph. D. dissertation, Univ. J. Fourier, Grenoble, 228 pp.

Testut C.-E., P. Brasseur, J.M. Brankart and J. Verron, 2003 : Assimilation of sea-surface temperature and altimetric observations during 1992-1993 into an eddy-permitting primitive equation model of the North Atlantic Ocean, *J. Mar. Syst.*, **40-41**, 291-316.

Thiebaux H.J., 1985: On approximations to geopotential and wind-field correlation structures, *Tellus*, **37A**, 126-131.

van Leeuwen P.J. and G. Evensen, 1996: Data Assimilation and Inverse Methods in Terms of a Probabilistic Formulation, *Mon. Wea. Rev.*, **124**, 2898-2913.

Verlaan M. and A.W. Heemink, 1997: Tidal flow forecasting using reduced-rank square root filter, *Stochastic Hydro. Hydraul.*, **11**, 349-368.

Verron J. and W.R. Holland, 1989: Impact de données d'altimétrie satellitaire sur les simulations numériques des circulations générales océaniques aux latitudes moyennes, *Ann. Geophys.*, **7**, 31-46.

Verron J., L. Gourdeau, D.T. Pham, R. Murtugudde and A.J. Busalacchi, 1999: An extended Kalman filter to assimilate satellite altimeter data into a nonlinear numerical model of the tropical Pacific Ocean: Method and validation, *J. Geophys. Res.* **104** (C3), 5441-5458.

Weaver A., J. Vialard and D.L.T. Anderson, 2003: Three- and Four-Dimensional Variational Assimilation with a General Circulation Model of the Tropical Pacific Ocean. Part I: Formulation, Internal Diagnostics, and Consistency Checks, *Mon. Wea. Rev.*, **131**, 1360-1378.

Whitaker J. and T.M. Hamill, 2002: Ensemble Data Assimilation without Perturbed Observations, *Mon. Wea. Rev.*, **130**, 1913-1924.