# Dynamical evolution of the error statistics with the SEEK filter to assimilate altimetric data in eddy-resolving ocean models

J. Ballabrera-Poy, Pierre Brasseur, Jacques Verron

► To cite this version:

J. Ballabrera-Poy, Pierre Brasseur, Jacques Verron. Dynamical evolution of the error statistics with the SEEK filter to assimilate altimetric data in eddy-resolving ocean models. Quarterly Journal of the Royal Meteorological Society, 2001, 127, pp.233-253. 10.1002/qj.49712757113. hal-00211729

HAL Id: hal-00211729

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

Submitted on 14 Jan 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
Dynamical evolution of the error statistics with the SEEK filter to assimilate altimetric data in eddy-resolving ocean models

By J. BALLABRERA-POY¹, P. BRASSEUR² and J. VERRON²*

¹Universities Space Research Association, USA

²LEGI, France

## SUMMARY

The Singular Evolutive Extended Kalman (SEEK) filter introduced by Pham et al. is applied to a primitive-equation model in order to reconstruct the mesoscale circulation typical of the mid-latitude ocean from altimetric data. The SEEK filter is a variant of the Kalman-filter algorithm based on two concepts: the order reduction of the initial-error covariance matrix, and the dynamical evolution of the reduced-order basis. This makes the method potentially suitable for problems with a high number of degrees of freedom.

Previous work has shown the ability of a steady version of the filter to improve the vertical structure of the ocean thermocline in the case of the quasi-linear dynamics associated with the equatorial tropical Pacific Ocean, and the need to combine the dynamical evolution of the basis with an adaptive scheme in a mid-latitude ocean model of the Gulf Stream region.

This work examines the potential advantages of the dynamical evolution of the basis functions with simple assimilation experiments. It demonstrates the ability of the method to propagate in time the statistical properties of the system when the filter is initialized properly. However, the lack of robustness of the filter is investigated theoretically and experimentally, showing the need to consider variants of the method when the filter is not properly initialized.

KEYWORDS: Altimetry Data assimilation Kalman filter Primitive-equation models

## 1. INTRODUCTION

During the last decade, the observation of the ocean from space has led to an unprecedented amount of oceanographic data. Satellites are the only observing system so far able to give a global and quasi-synoptic picture of the state of the ocean. Amongst the many types of data measured by satellites, altimetry plays a particular role because of the spatial coverage and accuracy of sea surface height (SSH) observations. For example, validation of TOPEX†/Poseidon (T/P) radar altimeter measurements in the tropical Pacific Ocean demonstrates an accuracy of 2 cm root-mean square (r.m.s.) for the instantaneous signal, and 3–4 cm r.m.s. for the low-frequency signal (Picaut et al. 1995). In addition to T/P, the European Remote Sensing satellites currently provide complementary SSH measurements. In the near future, several projects (Geosat Follow-On, Envisat/RA-2 and Jason1) guarantee the continuity of the altimetric observation of the ocean over the next decade. However, as the ocean is opaque to electromagnetic radiation, remote sensing gives information about the state of the ocean surface only. Data assimilation can therefore be used as a tool for vertical extrapolation, as well as temporal or horizontal extrapolation between satellite ground tracks if necessary.

More generally, data assimilation designates a concept that includes any statistical or deterministic method combining observations and numerical models. Formally, these methods need a representation of the system in terms of a state vector x. The number of components, n, is the number of state variables defining the system at a given time. The evolution of the state vector can be written as

$$x(t_k) = M\{x(t_{k-1})\}, \quad (1)$$

<footnote>* Corresponding author: LEGI, UMR 5519 CNRS, BP 53 X, 38041 Grenoble, France.</footnote>

<footnote>e-mail: verron@hmg.inpg.fr</footnote>

<footnote>† Typhoon Operational Experiment.</footnote>
where M represents the propagation operator, computing the evolution of the system from time $t_{k-1}$ to time $t_k$. The combination between observations, $y^o$, and model solutions, $x$, is possible only if a relation can be established between them. In other words, data assimilation is possible only if there is an application, H, for which,

$$y^o = H(x). \qquad (2)$$

Review papers on data assimilation methods and their applications to meteorology and oceanography are, for example, Ghil and Malanotte-Rizzoli (1991), Bennett (1992), Miller and Cane (1996), Talagrand (1997) and De Mey (1997). In brief, simple assimilation methods, like nudging, direct insertion, or vertical statistical correlation, have demonstrated their utility to modify the time evolution of simple models (Verron 1990; Hurlburt et al. 1990; Haines et al. 1993; Smedstad and Fox 1994). They have also been applied to the more complex ocean general-circulation models (OGCM) with relative success (e.g. Oschlies and Willebrand 1996). However, more sophisticated approaches have been investigated also, because they offer more guarantee in terms of optimality and provide quantitative error estimates.

The Kalman-filter (KF) algorithm (Kalman and Bucy 1960) is based on a sequence of forecast/correction cycles in which an estimate of the state of the system, symbolized by $x^f$, is corrected every time a set of observations $y^o$ is available. The updated field, $x^a$, traditionally called the analysis field, is then taken as the initial condition for a forecast step in which the numerical model is used to carry the information forward in time.

As the mathematical derivation, advantages and drawbacks of the KF algorithm have been widely reported (e.g. Gelb 1974; Ghil and Malanotte-Rizzoli 1991; Bouttier 1996), we only mention here that the application of the KF on meteorology and oceanography is limited by its numerical cost, associated with the size of the $n \times n$ error covariance matrices, and the number of model integrations (proportional to $n$) needed to compute the time evolution of these covariance matrices. As a consequence, the practical application of the full KF has been limited to systems with a small number of components (Miller and Cane 1989; Fu et al. 1993). The current state of the art of computers does not allow the application of the full KF to OGCM simulations with $n = O(10^6)$. On the other hand, the KF provides the optimal blend between observations and model solutions only in the case of linear systems. Several simplifications have been developed accordingly for high-dimensional, nonlinear OGCM simulations.

The Singular Evolutive Extended Kalman (SEEK) filter introduced by Pham et al. (1998) is one of these simplifications. In brief, the SEEK filter reduces the cost of the KF by assuming that the $n \times n$ error covariance matrix, $P$, is singular, i.e. the eigenvector decomposition

$$P = S^\top \Lambda S, \quad \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n), \qquad (3)$$

is such that $n-r$ eigenvalues $\lambda_i$ are equal to zero. That is, only $r$ eigenvectors are needed on the assimilation algorithm instead of the full covariance matrix. This hypothesis is equivalent to the condition of order reduction, in which the state vector may be expressed as

$$x = S^\top \mu, \qquad (4)$$

where $S^\top$ is now the $n \times r$ matrix of the eigenvectors associated with non-zero eigenvalues. Matrix $S^\top$ is also known as the pseudo-inverse of the projection that maps $x$ onto $\mu$. The need for order reduction in the KF algorithm may be justified, not only to overcome the problem of the numerical cost, but also for cognitive reasons. These problems are
related to the difficulty of identifying, from a large-scale observing system, the true structure of the error for all the scales resolved by the numerical model. A more detailed discussion about this topic can be found in Cane et al. (1996).

A straightforward construction of $S^T$ may be obtained with a truncated set of multivariate empirical orthogonal functions (EOFs) derived from a reference run of the model. In that case, $\mu$ are the principal components of the state vector. Several authors have used multivariate EOFs to reduce the degrees of freedom of the assimilation problem. For example, Blayo et al. (1998) reduce the cost of a variational adjoint method by defining the control space as the principal components of the initial condition instead of the initial condition itself. In the context of a sequential scheme, Cane et al. (1996) project the KF equations onto an EOF basis in order to construct an analog of the Kalman filter for the principal components. In their application, the time evolution of the principal components of a small perturbation of the state vector, $\delta\mathbf{x}$, is computed as follows:

$$ \delta\mathbf{x}_k = \mathbf{M}(\delta\mathbf{x}_{k-1}) \approx \mathbf{M}\delta\mathbf{x}_{k-1} = \mathbf{MS}^\top\mu_{k-1}, \quad (5) $$

where $\delta\mathbf{x}_k = \delta\mathbf{x}(t_k)$, $\mu_k = \mu(t_k)$, and $\mathbf{M} = M'_x$ is the linear tangent model associated with the transition matrix Eq. (1). Because of the orthogonality of the eigenvectors, $\mathbf{S}\mathbf{S}^\top = \mathbf{I}$, a linear transition model for the principal components was derived:

$$ \mu_k = \mathbf{SMS}^\top\mu_{k-1} \equiv \mathbf{A}\mu_{k-1}, \quad \text{with} \quad \mathbf{A} = \mathbf{SMS}^\top, \quad (6) $$

where $\mathbf{A}$ is an $r \times r$ matrix. In the SEEK filter, Pham et al. (1998) take a different approach to Eq. (5), by allowing the dynamical model to modify the directions defined by the columns of $S^T$, i.e.

$$ \delta\mathbf{x}_k \approx \mathbf{MS}_{k-1}^\top\mu_{k-1} \equiv \mathbf{S}_k^\top\mu_{k-1}, \quad \text{with} \quad \mathbf{S}_k^\top = \mathbf{MS}_{k-1}^\top, \quad (7) $$

i.e. the linear evolution of the current state of the system is expressed by the evolution of each one of the members of the subspace basis. This means that, if Eq. (4) is true, the time evolution of the $n \times n$ covariance matrix may be computed by integrating the $r$ elements of the subspace basis defined by the columns of $S^T$.

The SEEK filter has been applied with success to a reduced gravity, primitive-equation model of the tropical Pacific Ocean to explore the impact of altimetric observations on the structure of the thermocline (Verron et al. 1999). The filter has also been used with a nonlinear, primitive-equation model of the Gulf Stream region (Brasseur et al. 1999). In that case, the SEEK filter was found to be able to recover the vertical structure of the ocean from simulated altimetric observations in the case of twin experiments. Because of the numerical cost associated with the computation of the time evolution of the basis functions, a common point of these two works was the comparison between the dynamic filter and a steady filter, in which the basis is kept constant in time. In the case of the quasi-linear equatorial dynamics, both versions of the filter provide similar results (Verron et al. 1999). This is not the case in the mid-latitude region studied by Brasseur et al. (1999), where the steady filter displays a progressive degradation of the results. On the other hand, the sensitivity of the assimilation with respect to the choice of the basis at the initial time was found critical in the mid-latitude experiment. Therefore, a combination of the dynamic filter and an adaptive scheme was implemented to ensure the robustness of the method.

The work presented here is a follow-up of the assimilation experiments of Brasseur et al. (1999). However, the objective here is not to provide a recipe for the construction of a reduced-order algorithm, but to investigate the advantages of computing the dynamic
evolution of the reduced-order basis, and to study the conditions for which such a filter may or may not work. The outline of the paper is as follows. The equations of the filter, and the computation of the time evolution of the basis are presented in section 2. The sensitivity of the filter to the initial conditions is examined in section 3. The ideal twin-experiment configuration is presented in section 4. Despite its simplicity the configuration allows the development of mid-latitude mesoscale turbulence. The dimension of the system is significantly high (n = 202 700). The numerical experiments are presented in section 5, and a final discussion is given in section 6.

## 2. THE EQUATIONS OF THE SEEK FILTER

A complete derivation of the SEEK filter equations can be found in Pham *et al.* (1998) or Verron *et al.* (1999). Here, we summarize the equations, using the notations recommended by Ide *et al.* (1997), of the filter in the case of a linear perfect system:

$$ \Lambda_{k+1} = \{\Lambda_k^{-1} + (\mathbf{H}_k \mathbf{S}_k^{\top})^{\top} \mathbf{R}_k^{-1} \mathbf{H}_k \mathbf{S}_k^{\top}\}^{-1}, \quad (8) $$

$$ \mathbf{K}_k = \mathbf{S}_k^{\top} \Lambda_{k+1} (\mathbf{H}_k \mathbf{S}_k^{\top})^{\top} \mathbf{R}_k^{-1}, \quad (9) $$

$$ \mathbf{x}_k^a = \mathbf{x}_k^f + \mathbf{K}_k \{\mathbf{y}_k^o - \mathbf{H}_k \mathbf{x}_k^f\}, \quad (10) $$

$$ \mathbf{P}_k^a = \mathbf{S}_k^{\top} \Lambda_{k+1} \mathbf{S}_k, \quad (11) $$

$$ \mathbf{x}_{k+1}^f = \mathbf{M} \mathbf{x}_k^a, \quad (12) $$

$$ \mathbf{S}_{k+1}^{\top} = \mathbf{M} \mathbf{S}_k^{\top}, \quad (13) $$

$$ \mathbf{P}_{k+1}^{f} = \mathbf{S}_{k+1}^{\top} \Lambda_{k+1} \mathbf{S}_{k+1}. \quad (14) $$

In these equations $\mathbf{P}^f$, and $\mathbf{P}^a$ represent the covariance matrices associated with the errors included on the forecast $\mathbf{x}^f$ and the analysis $\mathbf{x}^a$ states, respectively, i.e.

$$ \mathbf{P}^{\alpha} = \mathrm{E}\{(\mathbf{x}^{\alpha} - \mathbf{x}^{\mathrm{t}})(\mathbf{x}^{\alpha} - \mathbf{x}^{\mathrm{t}})^{\top}\}, \quad (15) $$

where $\alpha$ represents either the forecast, f, or the analysis, a, and E is the expectation operator; superscript t indicates the true state of the system.

Equations (8)–(11) correspond to the correction of the first-guess $\mathbf{x}^f$ at time $t_k$, i.e. the analysis step. In particular, Eqs. (8) and (9) express the gain matrix of the classical KF filter in terms of the reduced-order subspace, i.e. its validity is related to the validity of Eq. (3). Equations (12)–(14) correspond to the temporal evolution of the state of the system and its expected error. The time evolution of the error, discussed in more detail in the following section, is expressed by the time evolution of the basis functions $\mathbf{S}^{\top}$.

## 3. TIME EVOLUTION OF THE ERROR SUBSPACE

Verron *et al.* (1999) discuss several methods of generalizing Eq. (13) to a nonlinear model. The method used here is based on the spread of a set of *r* states (where *r* is the rank of the covariance matrix) having the same covariance as the analysis error. This is done as follows. After the analysis step Eq. (10), the estimation of the analysis error is given by Eq. (11), where $\Lambda_{k+1}$ is an $r \times r$, symmetric and definite positive matrix. In this case, there is a unique triangular matrix $\mathbf{\Gamma}$, verifying $\Lambda_{k+1} = \mathbf{\Gamma}\mathbf{\Gamma}^{\top}$ (Cholesky Theorem). Such a triangular matrix $\mathbf{\Gamma}$ is used to define an $n \times r$ matrix, $\mathbf{S}_k^{\top}\mathbf{\Gamma}$. Let $\mathbf{s}_k^m$, $m = 1, \dots, r$, denote the *m*-th vector column of $r^{1/2}\mathbf{S}_k^{\top}\mathbf{\Gamma}$. Then, because of Eq. (11),
the ensemble covariance of the $r$ $n$-dimensional vectors $\mathbf{s}_k^m$ coincides with the analysis-error covariance $\mathbf{P}_k^a$.

If vectors $\mathbf{s}_k^m$ are used as a set of arbitrary perturbations around the analysis state $\mathbf{x}_k^a$, a set of $r+1$ model integrations is computed:

$$ \mathbf{x}_{k+1}^0 = \mathbf{M}(\mathbf{x}_k^a), \qquad (16) $$

$$ \mathbf{x}_{k+1}^m = \mathbf{M}(\mathbf{x}_k^a + \mathbf{s}_k^m), \quad m = 1, \dots, r. \qquad (17) $$

Now we assume that the spread of these vectors around its mean value is representative of the forecast error $\mathbf{P}_{k+1}^f$ as it is commonly assumed in Monte Carlo methods (Epstein 1996; Evensen 1994). Then, states Eqs. (16) and (17) are used to construct an $n \times (r+1)$ scatter matrix $\mathbf{S}_{k+1}^{\mathrm{T}}$ accounting for the dispersion around the mean state.

An orthogonal basis is obtained by applying the singular value decomposition (SVD) on such a scatter matrix, i.e. $\mathbf{S}_{k+1}^{\mathrm{T}} = \mathbf{UDV}^{\mathrm{T}}$, where $\mathbf{U}$ is an $n \times (r+1)$ matrix verifying $\mathbf{U}^{\mathrm{T}}\mathbf{U} = \mathbf{I}$, $\mathbf{D}$ is a square diagonal matrix, and $\mathbf{V}$ is a square matrix verifying $\mathbf{V}^{\mathrm{T}}\mathbf{V} = \mathbf{I}$. As the ensemble mean of the columns of the scatter matrix is equal to zero, at least one column can be expressed as a linear combination of the other columns, and matrix $\mathbf{D}$ has at least one diagonal element equal to zero. Because of the orthogonality properties of the SVD, we have

$$ \mathbf{P}_{k+1}^f = \frac{1}{r+1} \mathbf{S}_{k+1}^{\mathrm{T}} \mathbf{S}_{k+1} = \frac{1}{r+1} \mathbf{U}\mathbf{D}^2\mathbf{U}^{\mathrm{T}}. \qquad (18) $$

By comparing Eqs. (3) and (18) it may be written

$$ \mathbf{S}_{k+1}^{\mathrm{T}} \equiv \mathbf{U}, \quad \text{and} \quad \mathbf{\Lambda}_{k+1} \equiv \frac{1}{r+1} \mathbf{D}^2, \qquad (19) $$

where only the vectors of $\mathbf{U}$ associated with non-zero elements of $\mathbf{D}$ are retained. For an OGCM, the additional cost of the SVD of the scatter matrix is only a small fraction of the cost of $r+1$ model integrations. However, the advantages of Eqs. (18) and (19) are the orthogonality of the basis, and the automatic detection of any reduction of the rank of the system. More important is the fact that Eq. (19) provides an estimate of $\mathbf{\Lambda}_{k+1}$ related to the scatter of the ensemble at time $t_{k+1}$, and not from the estimation of the analysis error at time $t_k$ (Eqs. (8), (9) and (11)). This is important because of the inherent deficiencies of Eqs. (8) and (9) if Eq. (4) is not exact. These deficiencies appear because Eqs. (8) and (9) neglect the possible correlations between the reduced space and its complementary subspace.

### 4. FILTER INITIALIZATION

Brasseur et al. (1999) initialize the SEEK filter with two different sets of EOFs. The first set of EOFs accounts for the variance of the model during the time period corresponding to the synthetic observations to be assimilated. The second set of EOFs corresponds to the variability during a period different to the observed one. The assimilation experiments showed the need to improve the filter in order to increase the robustness of the filter in the second case. In this section, the reasons of such a lack of robustness is studied in terms of the dynamical evolution of the truncation error, i.e. the components of the error not accounted for by the basis $\mathbf{S}^{\mathrm{T}}$.

When the projection operator is not complete, or the description of the errors is not exact, Eq. (4) must be substituted by

$$ \mathbf{x} = \mathbf{S}^{\mathrm{T}} \mu + \mathbf{e}^{\mathrm{T}}, \qquad (20) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Typical patterns of instantaneous circulation. (a) Surface dynamical pressure (m²s⁻²). Dashed lines correspond to negative values. The thick solid line corresponds to zero values. Contour interval is 1 m²s⁻². (b) Surface horizontal velocity (m s⁻¹). Maximum velocity (longest arrow) is 1.6 m s⁻¹.</figcaption></figure>

where e^r is a vector belonging to the nullspace of S^T, called the truncation error, i.e.

$$Se^r = 0. \qquad (21)$$

The validity of Eqs. (8)–(14), as well as for any reduced-order simplification of the KF, is affected by the presence of such a truncation error. The filter is affected at both the analysis and the forecast steps. The impact on the analysis step has been discussed by Cane *et al.* (1996). In summary, Eq. (9) only extracts observational information about the retained modes, and Eq. (8) should contain additional terms accounting for the covariance of the truncated errors as well as for the correlations between the errors in the reduced space and the truncation error. Because of the extreme difficulty in determining such covariance matrices, it is a common strategy to redefine the observational error covariance **R**, or the reduced-order forecast covariance **Λ** (Fukumori and Malanotte-Rizzoli 1995; Cane *et al.* 1996; Pham *et al.* 1998). Following Pham *et al.* (1998) a compensation technique is used here. That is, the forecast covariance matrix **P**<sub>k</sub> is
<figure><img src="image_3.png" /><figcaption>Figure 2. Eigenvalue spectrum corresponding to the variability of the multivariate state vector for a 15-day period. For each mode k, the curve gives the explained variance accounted for by modes 1 to k.</figcaption></figure>

substituted by $\rho^{-1}\mathbf{P}_k^f$, where $0 < \rho < 1$. Therefore, Eq. (8) is substituted by

$$ \mathbf{\Lambda}_{k+1} = \{\rho \mathbf{\Lambda}_k^{-1} + (\mathbf{H}_k \mathbf{S}_k^\top)^\top \mathbf{R}_k^{-1} \mathbf{H}_k \mathbf{S}_k^\top\}^{-1}. \quad (22) $$

The impact of the truncation error on the forecast step of the filter is studied here in the simple case of a perfect, linear model. Let $\mathbf{e}_0$ be the true error at the initial time, $t_0$, i.e.

$$ \mathbf{e}_0 = \mathbf{S}_0^\top \mu_0 + \mathbf{e}_0^r. \quad (23) $$

The time evolution of such an error is given by

$$ \mathbf{e}_1 = \mathbf{M} \mathbf{S}_0^\top \mu_0 + \mathbf{M} \mathbf{e}_0^r = \mathbf{S}_1^\top \mu_0 + \mathbf{M} \mathbf{e}_0^r. \quad (24) $$

The true projection of $\mathbf{e}_1$ on $\mathbf{S}_1^\top$ is

$$ \mathbf{e}_1 = \mathbf{S}_1^\top \mu_1 + \mathbf{e}_1^r. \quad (25) $$

If $\mathbf{S}_1^\top$ is orthogonal, i.e. $\mathbf{S}_1 \mathbf{S}_1^\top = \mathbf{I}$, a relation between the true projection at time $t_1$ and the values at the initial time can be found:

$$ \mu_1 = \mu_0 + \mathbf{S}_1 \mathbf{M} \mathbf{e}_0^r, \quad (26) $$

$$ \mathbf{e}_1^r = (\mathbf{I} - \mathbf{S}_1^\top \mathbf{S}_1) \mathbf{M} \mathbf{e}_0^r. \quad (27) $$

Equation (26) states that the error projection on the reduced-order subspace, $\mu_0$, is correct only if the term $\mathbf{S}_1 \mathbf{M} \mathbf{e}_0^r$ can be neglected. This can be verified if one of these three conditions is satisfied:

(i) $\mathbf{e}_0^r = 0$, i.e. the initial error is perfectly described by the reduced-order subspace.

(ii) $\mathbf{S}_1 \mathbf{M} \mathbf{e}_0^r = 0$ is equivalent to the condition that the truncation error is dynamically uncoupled with the reduced-order subspace. Therefore, the validity of the reduced-order
<figure><img src="image_4.png" /></figure>

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

<figure><img src="image_6.png" /><figcaption>Figure 3. Surface dynamical pressure (m²s⁻²) associated with (a) the mean state, and (b)–(d) the first three modes. Contouring convention as in Fig. 1. Contour interval is 1 m²s⁻² in (a), and 0.2 m²s⁻² in (b)–(d).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3. Continued.</figcaption></figure>

approximation depends on the spectrum components of the truncation error and on the duration of the prediction step.

(iii) $Me_0^r = 0$ means that $\mathbf{M}$ is a singular matrix, and there is a family of states which are projected to zero by $\mathbf{M}$. If $q$ is the dimension of the nullspace of $\mathbf{M}$, the rank of $\mathbf{M}$ is $n-q$, where $n$ is the number of components of the state vector. For a nonlinear model, this condition can be written as $\mathbf{M}(\mathbf{x} + \mathbf{e}^f) \approx \mathbf{M}(\mathbf{x})$. That is, not all the fluctuations will develop, and an attractor exists inhibiting those fluctuations. The dimension of the reduced-order subspace should be given by the dimension of the attractor, and $S_0^T$ must contain all the components of the growing errors during the forecast time period.

In a real application, none of these constraints is expected to be verified, and the filter defined by Eqs. (8)–(14) will fail. These conditions, however, should be considered as a guide to build strategies to improve the filter. For example, the identification of the growing error modes over a given time interval may help to improve the error subspace. These conditions also give a better understanding about the reasons of the success of the scheme proposed by Brasseur *et al.* (1999). In their scheme, misfits between the analysis state, $\mathbf{x}_k^a$, and observations at that time, $\mathbf{y}_k^o$, are inverted, using a combination of statistical inversion and dynamical integration of the model, to provide a new direction susceptible to improve the filter basis. In essence, this adaptive scheme is not intended to identify growing error modes, nor determine the dynamically relevant modes of the flow, but to reduce, at each assimilation cycle, the truncation error, and then satisfy the first condition pointed earlier. Thus, and because of it, the numerical experiments described in the following section will focus on the study of the advantages of the dynamic filter when it is properly initialized.

# 5. NUMERICAL EXPERIMENTS

## (a) *Model and physical configuration*

The ocean model used in the numerical experiments is derived from the semi-spectral primitive-equation model (SPEM) of Haidvogel *et al.* (1991). The physical configuration consists of a rectangular box in the $\beta$-plane approximation. The domain has an extension of 2560 km zonally and 1920 km meridionally. The domain is
centred to fit the geographical position of the Gulf Stream ($f_0 = 9.3 \times 10^{-5}$ m s$^{-1}$, $\beta_0 = 2 \times 10^{-11}$ m$^{-1}$s$^{-1}$). The model has a flat bottom at a depth of 5000 m, and five levels located at 0, 400, 1100, 2500 and 5000 m depth. Vertical stratification at initial time is given by a Brunt-Väisälä frequency of $N^2(z) = 5.9 \times 10^{-5}$ exp(z/800) s$^{-2}$. The horizontal resolution of the grid is equal to 20 km in both directions, allowing the model to solve mesoscale eddies. The state vector of this configuration is defined as

$$ \mathbf{x} = \{p_s(x, y), u(x, y, z), v(x, y, z), \rho(x, y, z)\}, \qquad (28) $$

where $p_s$ is the surface dynamic pressure, i.e. surface pressure divided by a reference density, $(u, v)$ is the horizontal current velocity, and $\rho$ is the density. All the variables are normalized by the standard deviation of each one during a reference integration of the model. The number of components of the state vector is $n = 202 700$.

Lateral boundary conditions are free slip and impermeability. For simplicity, no heat or mass flux is considered across the boundaries. The system is forced by a zonal-wind stress $\tau_x(y) = -10^{-4} \cos(2\pi y/L_y)$ m$^2$s$^{-2}$, where $L_y$ is the meridional extension of the domain. Bottom friction is parametrized by a linear friction law $-C_D \vec{u}_{\text{bottom}}$, with the drag coefficient $C_D = 2.65 \times 10^{-4}$ m s$^{-1}$. Lateral dissipation is provided by a biharmonic operator with a coefficient of $8 \times 10^{10}$ m$^4$s$^{-1}$.

The model is integrated from rest over 25 years in order to reach a regime of statistically stable mesoscale turbulence. Figure 1 shows a typical circulation pattern after the spin-up. The wind forcing and the gradient of the Coriolis force are responsible for the spatial structure of the circulation: a double cell with cyclone circulation on the north and anticyclone circulation on the south. Western boundary currents associated with each cell merge and feed an eastward jet that penetrates into the domain. The size of the eddies originated by the instabilities of the central jet is governed by the first baroclinic radius of deformation (roughly 50 km, consistently with the value at the Gulf Stream region).

### (b) Experiment strategy

Three sets of assimilation experiments are presented here. The first one is intended to demonstrate the potentialities of the dynamic filter when properly initialized. Two additional sets of experiments are presented to investigate the robustness of the filter when not correctly initialized. For simplicity, a twin-experiment strategy is used to allow a full control of the true initial error statistics.

The experimental set-up is constructed from a reference trajectory of the model starting at year 25. The first state of the reference trajectory is now considered to be the origin of the time, i.e. $t_0 = 0$. A set of ‘observed’ full SSH fields is assimilated every five days, starting at time 0. The impact of the assimilation on every field $q$ of the state vector is measured in terms of the residual error (RE):

$$ \text{RE}(q) = \frac{\text{r.m.s.}(q^a - q^t)}{\text{r.m.s.}(q^r - q^t)}, \qquad (29) $$

where $q$ represents the surface dynamic pressure or a given level of the zonal velocity, the meridional velocity, or the density. Superscripts a, t, and r respectively indicate assimilation, true, and the value obtained with the assimilation routine turned off. Then, values of Eq. (29) smaller than 1 indicate a positive impact of the assimilation, and a negative impact otherwise.

In all the experiments, the observational error $\mathbf{R}$ is given by a diagonal matrix corresponding to a random error with standard deviation equal to 5 cm. Initialization
<figure><img src="image_8.png" /></figure>

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

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

Figure 4. Time evolution of the residual error during six months of assimilation of sea surface height. Perfect filter initialization. (a) Surface dynamical pressure, (b) zonal velocity at 2500 m depth, and (c) meridional velocity at 2500 m depth. Dynamic (solid line) and steady-filter (dotted line) versions of the filter are compared.
<figure><img src="image_11.png" /></figure>

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

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

Figure 5. Dynamic filter. Surface dynamical pressure (m$^2$s$^{-2}$) after six months of data assimilation. (a) Mean state $x_{180}^f$, (b) first mode of the basis $S_{180}^\top$, (c) second mode, and (d) third mode. See text for further explanation. Contouring convention as in Fig. 1. Contour interval is 1 m$^2$s$^{-2}$ in (a), and 0.001 in (b)–(d).
<figure><img src="image_14.png" /><figcaption>Figure 5. Continued.</figcaption></figure>

of the SEEK filter requires an initial first guess of the state of the system, of an initial set of basis functions $S_0^T$, and of the initial reduced-rank error covariance $\Lambda_0$.

### (c) Perfect filter initialization

In a first set of experiments, the impact of the dynamic evolution of the reduced-order basis $S^T$ on the assimilation algorithm is studied in the case of a proper initialization. Accurate description of the initial error by a small number of modes is achieved by an eigenvector analysis of the model variance around the state to be reconstructed. Here, a set of 15 multivariate EOFs is obtained from the covariance of the daily variability of the model during days 0 to 15. As usual, eigenvectors are ordered by decreasing eigenvalue. Since the decorrelation time-scale for eddy dynamics is substantially larger than 15 days, the variability of the system during such a time period is explained by a very small number of empirical modes, and the first five EOFs explain 99.3% of the variability over the 15-day period (see Fig. 2). Surface pressure of the mean value, and the first three dominant modes is presented in Fig. 3. These fields indicate that the variability during this period is associated with the meandering of the jet.

For this experiment, eigenvectors of the 16-day variability define $S_0^T$, the associated eigenvalues define $\Lambda_0$, and the time-average state, $\bar{x}$, is $x_0^f$. In that case, a negligible truncation error is obtained by initializing the filter with the ten first modes only. Thus, $r = 10$ for these experiments. Equations (8)–(10) are used to correct $x_0^f$, and the analysis field $x_0^a$ is used as initial conditions for a five-day model run.

The evolution of the filter performance is studied for both the steady and the dynamic versions of the SEEK filter. In the case of the steady filter, matrix $S^T$ is kept constant in time, while $\Lambda_k$ is updated by the sequential application of Eq. (8). Orthogonal matrices $S_k^T$, and diagonal matrices $\Lambda_k$ of the dynamic filter are obtained with Eqs. (16)–(19).

Figure 4 illustrates the impact of the six-month assimilation of SSH fields. The curves display the temporal evolution of the residual error (RE) of the surface pressure and velocity components at 2500 m depth. Because of the perfect identification of the initial error, the state of the system is already identified at the first assimilation step. The three-dimensional EOFs instantaneously extrapolate the information from the surface
<figure><img src="image_15.png" /><figcaption>Figure 6. Truncation error ($m^2s^{-2}$) associated with the true error at day 180: (a) steady filter, and (b) dynamic filter. Contouring convention as in Fig. 1. Contour interval is $1 m^2s^{-2}$ in (a), and $0.1 m^2s^{-2}$ in (b).</figcaption></figure>

layer to all the variables and to all the grid points. Therefore, the diminution of the RE is simultaneously obtained for all the variables and for all the depth levels.

The time evolution of the RE indicates the ability of the assimilation algorithm to propagate forward in time the, initially true, error covariance. However, as the EOF basis only spans the statistical relations between the fields during a very short segment of the trajectory of the system through the phase space, the ability of the steady filter decreases as the system evolves away from the initial state. On the contrary, as the basis functions are allowed to evolve as the system moves away from the initial state, the performance of the dynamic filter is maintained over the whole assimilation period.

Figure 5 displays the forecast state, $\mathbf{x}_{180}^f$, at day 180 and the three first modes provided by Eq. (18). Comparing Figs. 3 and 5, it can be seen that at the initial time, maxima of the EOFs were located mainly over the path of the jet. Now, maxima of the modes are located around other flow structures as, for example, a large eddy located north of the jet, or regions of low variability as the north-western boundary region. As the modes of $\mathbf{S}^T$ still come from a singular value decomposition of a covariance matrix, higher modes (small eigenvalues) contain smaller structures. However, the size of the
<figure><img src="image_16.png" /><figcaption>Figure 7. Maximum and minimum amplification factor of the perturbations used to compute the time evolution of the basis $S_k^T$ in the case of perfect filter initialization (see text).</figcaption></figure>

structures of the new set of EOFs is larger than the spatial structures of the original set of EOFs. This reflects the transfer of information from small to large scales during the dynamical evolution of $S^T$. The representativity of the new set of modes is measured by the ability of the new modes to describe the forecast error ($e^f = x^t - x^f$). As the true state is known in the case of twin experiments, a principal component analysis of $e^f$ at day 180 is done, and the truncation error at that day may be computed. Figure 6 displays the surface pressure of the truncation error for both the steady and the dynamic filter. As expected, the truncation error is higher in the case of the steady filter (a rough 50% of the value of the field), and smaller in the case of the dynamic filter (2% of the value of the field). On the other hand, the error of the dynamic-filter solution does not contain large-scale patterns as in the case of the steady-filter solution error.

Finally, Fig. 7 shows the time evolution of the amplification of the perturbations, measured as the ratio $||M(x + \delta x) - M(x)||/||\delta x||$. As the components of the state vector have been normalized with respect to the standard deviation of each field, the Euclidian metric can be used to compute the norm of the perturbations. Figure 7 only shows the maximum and the minimum of the amplitude amplification for each assimilation cycle. These curves indicate that, after a while, the dynamic filter seems to discard some decaying modes, keeping only modes that are amplified. This fact indicates some analogies between the dynamic filter and the breeding method of Toth and Kalnay (1997) used to identify perturbations of rapid amplification.

As a conclusion, these experiments have illustrated the advantage of a dynamic filter to propagate the error covariance matrix in the context of mid-latitude nonlinear ocean circulations. However, the simplicity of the experimental framework used for this study places some limits on the generality of the results in the case of realistic models, for which some of the hypotheses described previously may not be fully verified.
<figure><img src="image_17.png" /><figcaption>Figure 8. Initial r.m.s. error of the surface pressure (m²s⁻²): (a) true (perfect filter initialization), and (b) false (imperfect filter initialization). Contouring convention as in Fig. 1. Contour interval is 0.5 m²s⁻² in (a), and 0.05 m²s⁻² in (b).</figcaption></figure>

(d) *Imperfect filter initialization*

In order to study the degradation of the performance of the filter when the initial error statistics are not perfectly determined, a new set of experiments is presented here. The initialization of the assimilation filter is basically the same as in the previous section: same first estimate, $x_0$, and same set of modes $S_0^T$. However, the initial spectrum of the error $\Lambda_0$ is now false. This has been done by inverting the order of eigenvalues, i.e. $\lambda'_k = \lambda_{r-k+1}$, where $\lambda$ and $\lambda'$ indicate the original eigenvalue spectrum and the new spectrum, respectively. Figure 8 compares the initial surface pressure r.m.s. error of the perfect initialization (Fig. 8(a)), and the initial r.m.s. error when the order of eigenvectors is inverted (Fig. 8(b)). Now the initial error covariance is determined by the tenth multivariate EOF, and it does not reflect either the amplitude or the spatial distribution of the initial error.

Figure 9 shows the temporal evolution of the residual error for both the steady and dynamic filter for the new initialization. The main difference from Fig. 4 is found on the first steps of the assimilation cycle: with the new initialization, more assimilation
<figure><img src="image_18.png" /></figure>

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

<figure><img src="image_20.png" /><figcaption>Figure 9. Same as Fig. 4 but imperfect filter initialization.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 10. Same as Fig. 7 but imperfect filter initialization.</figcaption></figure>

cycles are necessary to reach the minimum of the residual error, especially in the deep ocean. Comparison of Figs. 4 and 9 also reveals a greater sensitivity of the steady filter to the false initialization of the filter. On the contrary, the dynamic filter displays, after a few assimilation cycles, a behaviour similar to the case of perfect initialization. That fact indicates that the performance of the filter is related to the determination of the error modes, but not with the precise spectral distribution of the error on such modes. This simplifies the problem of the filter initialization, because the dynamical filter only needs a first guess of the state of the system $x_0$, and a set of modes $S_0^T$.

Figure 10 shows the evolution of the maximum and minimum amplification for the new experiment. As before, the tendency of the filter to progressively discard decaying error modes is noticeable.

### (e) False error modes

In the last set of experiments, the filter is initialized with a set of modes $S_0^T$ chosen to be not representative of the initial error. This has been done by choosing the columns of $S_0^T$ as the eigenvectors of the daily variability between days 720 and 735. The spectral distribution of the new set of eigenvalues (not shown) is completely similar to the spectrum shown in Fig. 2. This is related to the nature of the forcing (steady) and the fact that the covariance matrix is also computed over a 15-day period. Results of the assimilation experiment are shown in Fig. 11. Neither the steady filter nor the dynamical filter are able to drive the system to the reference trajectory of the system. The reduction of the residual error is negligible because of the complete incompatibility between the initial error and the initial error modes. Note that under these conditions the dynamic filter does not provide any advantage with respect to the steady filter.
<figure><img src="image_22.png" /><figcaption>Figure 11. Time evolution of the residual error during two months of assimilation of sea surface height. Dynamic (solid line) and steady-filter (dotted line) versions of the filter are compared.</figcaption></figure>

## 6. SUMMARY

The singular evolutive extended Kalman (SEEK) filter is a reduced-order variant of the Kalman filter based on a set of multivariate functions to represent the error statistics. It has been applied to assimilate altimetric observations into a variety of ocean models (Verron *et al.* 1999; Brasseur *et al.* 1999). The reduction of the number of degrees of freedom by the help of multivariate EOFs has a double advantage in the context of data assimilation in meteorology and oceanography. First, it allows the use of advanced assimilation algorithms combined with complex model simulations. Second, it allows an instantaneous extrapolation of the information from observation locations to all the variables at all the model grid points. The main originality of this method is in the consistency between the dynamic evolution of such a set of multivariate functions and the current evolution of the system.

The initialization of the filter requires the specification of a first guess for the system state, $x_0$, a set of basis functions, i.e. the column vectors of matrix $S_0^T$, and the description of the initial error covariance in terms of the subspace defined by $S_0^T$, $\Lambda_0$. Because of the possible existence of error components not accounted for by the subspace, the performance of the reduced-order filter is directly related to the statistical significance of the actual, unknown initial error in that subspace. The work presented here has investigated the implications of this constraint, both theoretically and experimentally, helping in understanding the conditions for a proper initialization of the filter. There are three conditions: (i) the complete description of the initial error; (ii) the identification of all the growing error modes; or (iii) the fact that the retained modes are dynamically uncoupled from the discarded modes. Despite the difficulty in verifying these conditions in practice, they provide a guideline for further improvements of the filter.

The numerical experiments presented in section 5 compare the performances of the steady as against the dynamic filters to investigate if the dynamic evolution of the basis
plays any significant role. In these experiments, the steady filter keeps constant the structure of basis functions. On the contrary, the dynamical filter computes the time evolution of these functions from the estimated analysis error at the previous time.

The main result of these experiments is the demonstration of the ability of the dynamic filter to propagate a statistically relevant basis of multivariate functions forward in time, at least if these functions were statistically significant at the initial time. The amplification rate of the perturbations defining the time evolution of the basis seems to indicate that the dynamic filter naturally discards the decaying error modes, and only maintains the growing error modes. However, the process of mode selection is slow, and the performance of the filter might be increased if a prior selection of the growing modes was done before the assimilation cycles.

In addition, these experiments indicate that the dynamic filter is more sensitive to the initial choice of the basis functions $S_0^\top$, than to the description of the error covariance itself, $\Lambda_0$. That is, the problem of the filter initialization may be reduced only to the determination of $x_0$ and $S_0^\top$.

Finally, a simple experiment has shown that the dynamic evolution of the subspace plays no role if the initial value of $S_0^\top$ does not provide the relevant information about the true initial error, noting in that case that the steady filter fails too. That is, the dynamic filter by itself does not correct for irrelevant initialization of a reduced-order Kalman filter.

Further investigations will focus on both the improvement of the filter initialization and on the study of alternate techniques to make the filter robust with respect to the error statistics. Some options can be derived from the conditions of the proper initialization of the filter as described in section 4. For example, singular vectors of the transition operator, i.e. the family of optimal growing perturbations, may be used to capture some of the growing modes of the error. Another possibility is the combination of the SEEK filter with another assimilation method based on the adjoint model to identify the dynamically relevant modes of error.

## ACKNOWLEDGEMENTS

This work was done with the support of Centre Nationale de la Recherche Scientifique (CNRS) and Centre Nationale d'Etudes Spatiales. The calculations were performed with the support of IDRIS/CNRS and the Centre Grenoblois de Calcul Vectoriel.

## REFERENCES

<table><tr><td>Bennett, A.</td><td>1992</td><td>Inverse methods in physical oceanography. Cambridge University Press</td></tr><tr><td>Blayo, E., Blum, J. and Verron, J.</td><td>1998</td><td>Assimilation variationnelle de données en océanographie et réduction de la dimension de l'espace de contrôle. Pp. 199-219 in Equations aux Derivées Partielles et Applications. Gauthier-Villars, Paris</td></tr><tr><td>Bouttier, F.</td><td>1996</td><td>The Kalman filter. Pp. 221-245 in Predictability, vol. I. Proceedings of seminar, European Centre for Medium-Range Weather Forecasts, Reading, UK</td></tr><tr><td>Brasseur, P., Ballabrera-Poy, J. and Verron, J.</td><td>1999</td><td>Assimilation of altimetric data in the mid-latitude oceans using the SEEK filter with an eddy-resolving primitive equation model. J. Marine Systems, 22, 269-294</td></tr><tr><td>Cane, M. A., Kaplan, A., Miller, R. N., Tang, B., Hackert, E. C. and Busalacchi, A. J.</td><td>1996</td><td>Mapping tropical Pacific sea level: Data assimilation via a reduced state space Kalman Filter. J. Geophys. Res., 101, 22599-22617</td></tr><tr><td>De Mey, P.</td><td>1997</td><td>Data assimilation of the oceanic mesoscale: a review. J. Meteorol. Soc. Jpn., 75, 415-425</td></tr></table>
<table>
  <tr>
    <td>Epstein, E. S.</td>
    <td>1996</td>
    <td>Stochastic dynamic prediction. Tellus, 21A, 739-759</td>
  </tr>
  <tr>
    <td>Evensen, G.</td>
    <td>1994</td>
    <td>Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res., 99, 10143-10162</td>
  </tr>
  <tr>
    <td>Fu, L.-L., Fukumori, I. and Miller, R. N.</td>
    <td>1993</td>
    <td>Fitting dynamic models to the Geosat sea level observations in the tropical Pacific Ocean. II: A linear, wind-driven model. J. Phys. Ocean., 23, 2162-2181</td>
  </tr>
  <tr>
    <td>Fukumori, I. and Malanotte-Rizzoli, P.</td>
    <td>1995</td>
    <td>An approximate Kalman Filter for ocean data assimilation: An example with an idealized Gulf Stream model. J. Geophys. Res., 100, 6777-6793</td>
  </tr>
  <tr>
    <td>Gelb, A.</td>
    <td>1974</td>
    <td>Applied optimal estimation. MIT Press, Cambridge, MA, USA</td>
  </tr>
  <tr>
    <td>Ghil, M. and Malanotte-Rizzoli, P.</td>
    <td>1991</td>
    <td>Data assimilation in meteorology and oceanography. Adv. Geophys., 33, 141-266</td>
  </tr>
  <tr>
    <td>Haidvogel, D. B., Wilkin, J. and Young, R. E.</td>
    <td>1991</td>
    <td>A semi-spectral primitive equation ocean circulation model using vertical sigma and horizontal orthogonal curvilinear coordinates. J. Comp. Phys., 94, 151-185</td>
  </tr>
  <tr>
    <td>Haines, K., Malanotte-Rizzoli, P., Young, R. E. and Holland, W. R.</td>
    <td>1993</td>
    <td>A comparison of two methods for the assimilation of altimeter data into a shallow-water model. Dyn. Atmos. Ocean, 17, 89-133</td>
  </tr>
  <tr>
    <td>Hurlburt, H. E., Fox, D. N. and Metzger, E. J.</td>
    <td>1990</td>
    <td>Statistical inference of weakly correlated subthermohaline fields from satellite altimeter data. J. Geophys. Res., 95, 11375-11409</td>
  </tr>
  <tr>
    <td>Ide, K., Bennett, A. F., Courtier, P., Ghil, M. and Lorenc, A. C.</td>
    <td>1997</td>
    <td>Unified notation for data assimilation: Operational, sequential and variational. J. Meteorol. Soc. Jpn., 75(1B), 181-189</td>
  </tr>
  <tr>
    <td>Kalman, R. E. and Bucy, R. S.</td>
    <td>1960</td>
    <td>A new approach to linear filtering and prediction problems. J. Basic Engineering (ASME), 82D, 35-45</td>
  </tr>
  <tr>
    <td>Miller, R. M. and Cane, M. A.</td>
    <td>1989</td>
    <td>A Kalman filter analysis of sea level height in the tropical Pacific. J. Phys. Ocean., 19, 773-790</td>
  </tr>
  <tr>
    <td>Oschlies, A. and Willebrand, J.</td>
    <td>1996</td>
    <td>Tropical data assimilation: theoretical aspects. Pp. 207-233 in Modern Approaches to data assimilation in ocean modeling. Ed. P. Malanotte-Rizzoli, Elsevier Science, Amsterdam</td>
  </tr>
  <tr>
    <td>Pham, D. T., Verron, J. and Roubaud, M.-C.</td>
    <td>1998</td>
    <td>Assimilation of Geosat altimeter data into an eddy-resolving primitive equation model of the North Atlantic Ocean. J. Geophys. Res., 101, 14175-14190</td>
  </tr>
  <tr>
    <td>Picaut, J., Busalacchi, A. J., McPhaden, M. J., Gourdeau, L., Gonzalez, F. I. and Hackert, E. C.</td>
    <td>1995</td>
    <td>Singular Evolutive Extended Kalman filter with EOF initialization for data assimilation in oceanography. J. Marine Systems, 16, 323-340</td>
  </tr>
  <tr>
    <td>Smedstad, O. M. and Fox, D. N.</td>
    <td>1994</td>
    <td>Open-ocean validation of Topex/Poseidon sea level in the western equatorial Pacific. J. Geophys. Res., 100, 25109-25127</td>
  </tr>
  <tr>
    <td>Talagrand, O.</td>
    <td>1997</td>
    <td>Assimilation of observations, an introduction. J. Meteorol. Soc. Jpn., 75, 191-209</td>
  </tr>
  <tr>
    <td>Toth, Z and Kalnay, E.</td>
    <td>1997</td>
    <td>Ensemble forecasting at NCEP and the breeding method. Mon. Weather Rev., 125, 3297-3319</td>
  </tr>
  <tr>
    <td>Verron, J.</td>
    <td>1990</td>
    <td>Altimeter data assimilation into an ocean circulation model: Sensitivity to orbital parameters. J. Geophys. Res., 95, 11443-11459</td>
  </tr>
  <tr>
    <td>Verron, J., Gourdeau, L., Pham, D. T., Murtugudde, R. and Busalacchi, A. J.</td>
    <td>1999</td>
    <td>An extended Kalman filter to assimilate satellite altimeter data into a non-linear numerical model of the tropical Pacific: Method and validation. J. Geophys. Res., 104, 5441-5458</td>
  </tr>
</table>