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

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.


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 TOPEXtiPoseidon (TIP) radar altimeter measurements in the tropical Pacifi c Ocean demonstrates an accuracy of 2 em root-mean square (r.m.s.) for the instantaneous signal, and 3-4 em r. m.s. for the low-frequency signal (Picaut et al. 1995). In addition to TIP, the European Remote Sensing satellites currently provide complementary SSH measurements. In the near future, several projects (Geosat Follow On, Envisat/RA-2 and Jasonl) 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 defi ning the system at a given time. The evolution of the state vector can be written as where M represents the propagation operator, computing the evolution of the system from time tk-1 to time tk. The combination between observations, y 0 , 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 0 = H(x).
(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 assim ilation 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 xf , is corrected every time a set of observations y 0 is available. The updated field, xa, 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 197 4; 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 x 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 = (!) (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 simplifi cations. In brief, the SEEK filter reduces the cost of the KF by assuming that the n x n error covariance matrix, P, is singular, i.e. the eigenvector decomposition P=ST AS, A =diag(.A,, ... , .An), is such that nr eigenvalues Ai 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 (4) where S T is now the n x r matrix of the eigenvectors associated with non-zero eigenval ues. Matrix S T is also known as the pseudo-inverse of the projection that maps x onto JL.
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, 11-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, ax, is computed as follows: (5) where 8xk = 8x(tk), ilk= il(tk), and M = M� is the linear tangent model associated with the transition matrix Eq. (1 ). Because of the orthogonality of the eigenvectors, SST = I, a linear transition model for the principal components was derived: where A is an r x 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 ST, i.e. 8xk � MSJ_ 1 ilkl = sJ ilk -! , with sJ = MSJ_ 1 , 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 x n covariance matrix may be computed by integrating the r elements of the subspace basis defined by the columns of ST.
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 obser vations 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 ex periments. 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 fi lter 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 fi lter 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 confi guration is presented in section 4. Despite its simplicity the confi guration 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 fi nal discussion is given in section 6.

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: Ak+l = {Ai: 1 + (HkSJ)TRi: 1 HkSJ} -1 , Kk = SJ Ak+l (HkSJ) T Ri: 1 , M a x k+l = x k, sJ +1 = MSJ, Pk +l = SJ +I Ak +1Sk+l· In these equations pf, and pa represent the covariance matrices associated with the errors included on the forecast xf and the analysis xa states, respectively, i.e. (15) where a 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 xf at time tb i.e. the analysis step. In particular, Eqs. (8) and (9) express the gain matrix of the classical KF fi lter 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 ST.

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 Ak+ 1 is an r x r, symmetric and defi nite positive matrix. In this case, there is a unique triangular matrix f, verifying Ak+ 1 = ff T (Cholesky Theorem). Such a triangular matrix r is used to define an n X r matrix, sJ r. Let s;' m = 1, . .. , r, denote the m-th vector column of r 1 12sJr. Then, because of Eq. (11), the ensemble covariance of the r n-dimensional vectors sJ: coincides with the analysis error covariance P�.
If vectors s'J: are used as a set of arbitrary perturbations around the analysis state xk, a set of r + 1 model integrations is computed: xk + l = M (xk + s'J:), m = 1, ... , r. (17) Now we assume that the spread of these vectors around its mean value is represen tative of the forecast error P k + 1 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 x (r + 1) scatter matrix s�:;_ 1 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. s�:;_ 1 = UDV T , where U is an n x (r + 1) matrix verifying U T U = I, D is a square diagonal matrix, and V is a square matrix verifying V T V = 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 D has at least one diagonal element equal to zero. Because of the orthogonality properties of the SV D, we have By comparing Eqs. (3) and (18) it may be written where only the vectors of U associated with non-zero elements of 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 Ak+ 1 related to the scatter of the ensemble at time tk+ 1, and not from the estimation of the analysis error at time tk (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.

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 as similation 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 ST.
When the projection operator is not complete, or the description of the errors is not exact, Eq. (4) must be substituted by (20) ,-... . where e r is a vector belonging to the nulls pace of S T, called the truncation error, i.e. se r = 0.
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 A (Fukumori and Malanotte period. For each mode k, the curve gives the explained variance accounted for by modes 1 to k.
The impact of the truncation error on the forecast step of the fi lter is studied here in the simple case of a perfect, linear model. Let eo be the true error at the initial time, to, The time evolution of such an error is given by The true projection of e1 on SJ is ST r 2 e 1 = 1 Ill + e 1.
If sT is orthogonal, i.e. slsJ = I, a relation between the true projection at time tj and the values at the initial time can be found: Equation (26) states that the error projection on the reduced-order subspace, /.to, is correct only if the term S1Me� can be neglected. This can be verified if one of these three conditions is satisfi ed: (i) eo = 0, i.e. the initial error is perfectly described by the reduced-order subspace.
(ii) S1 Me0 = 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  In a real application, none of these constraints is expected to be verifi ed, and the fi lter defined by Eqs. (8)-(14) will fail. These conditions, however, should be considered as a guide to build strategies to improve the fi lter. For example, the identifi cation 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, misfi ts between the analysis state, x�, and observations at that time, yk, are inverted, using a combination of statistical inversion and dynamical integration of the model, to provide a new direction susceptible to improve the fi lter 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 fi rst 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 fi lter when it is properly initialized. f3 o = 2 x w-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-Vaisala frequency of N 2 (z) = 5.9 X w-s 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 where Ps is the surface dynamic pressure, i.e. surface pressure divided by a reference density, (u , v) is the horizontal current velocity, and p 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 rx (y) = -10-4 cos(2rry I Ly) m 2 s-2 , where Ly is the meridional extension of the domain. Bottom friction is parametrized by a linear friction law -CDubottom. with the drag coefficient CD = 2.65 x w-4 m s-1 . Lateral dissipation is provided by a biharmonic operator with a coefficient of 8 x 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 Corio lis 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) E x periment strateg y 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 fi lter 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. to= 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): 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 R is given by a diagonal matrix corresponding to a random error with standard deviation equal to 5 em. Initialization  � 1 000 ���-il�----  of the SEEK filter requires an initial fi rst guess of the state of the system, of an initial set of basis functions S;j, and of the initial reduced-rank error covariance Ao.

(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 initializa tion. 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 fi ve 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;j, the associated eigenvalues define Ao, and the time-average state, x, is x� . 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) -( 1 0) are used to correct x� , and the analysis field x0 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 Ak is updated by the sequential application of Eq. (8). Orthogonal matrices SJ, and diagonal matrices Ak 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  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, x{ 8 0 , 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 S T still come from a singular value decomposition of a cqvariance matrix, higher modes (small eigenvalues) contain smaller structures. However, the size of the 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 ST. The representativity of the new set of modes is measured by the ability of the new modes to describe the forecast error (ef = x t -xf). As the true state is known in the case of twin experiments, a principal component analysis of ef 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 IIM(x +ox) -M(x)ll/118xll . 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.

(d) Imperfect filter initialization
In order to study the degradation of the performance of the fi lter 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, xo, and same set of modes SJ. However, the initial spectrum of the error Ao is now false. This has been done by inverting the order of eigenvalues, i.e. A.�= Ar-k+l· where A. and A.' 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 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 fi lter 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 simplifi es the problem of the fi lter initialization, because the dynamical fi lter only needs a fi rst guess of the state of the system xo, and a set of modes SJ". Figure 10 shows the evolution of the maximum and minimum amplifi cation 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 fi lter is initialized with a set of modes SJ" chosen to be not representative of the initial error. This has been done by choosing the columns of SJ" 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 fi lter 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 fi lter does not provide any advantage with respect to the steady fi lter.

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 specifiqtion of a first g�ess for the system state, xo, a set of basis functions, i.e. the qplumn vectors of matrix SJ, and the description of the initial error covariance in terms of the subspace defined by SJ, Ao. Because of the possible existence of error pqmpom�nts not accounted for by the subspace, the performance of the reduced-�rder 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 SJ, than to the description of the error covariance itself, Ao. That is, the problem of the filter initialization may be reduced only to the determination of XQ and SJ.
Finally, a simple experiment has shown that the dynamic evolution of the subspace plays no role if the initial value of SJ 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.