Sequential assimilation of multi-mission dynamical topography into a global finite-element ocean model

This study focuses on an accurate estimation of ocean circulation via assimilation of satellite measurements of ocean dynamical topography into the global finite-element ocean model (FEOM). The dynamical topography data are derived from a complex analysis of multi-mission altimetry data combined with a referenced earth geoid. The assimilation is split into two parts. First, the mean dynamic topography is adjusted. To this end an adiabatic pressure correction method is used which reduces model divergence from the real evolution. Second, a sequential assimilation technique is applied to improve the representation of thermodynamical processes by assimilating the time varying dynamic topography. A method is used according to which the temperature and salinity are updated following the vertical structure of the first baroclinic mode. It is shown that the method leads to a partially successful assimilation approach reducing the rms difference between the model and data from 16 cm to 2 cm. This improvement of the mean state is accompanied by significant improvement of temporal variability in our analysis. However, it remains suboptimal, showing a tendency in the forecast phase of returning toward a free run without data assimilation. Both the mean difference and standard deviation of the difference between the forecast and observation data are reduced as the result of assimilation.


Introduction
Reliable estimation of the ocean circulation is one of the central topics in the oceanography.This problem has two aspects.The first one is the availability of an adequate ocean Correspondence to: S. Skachko (skachko.sergey@uqam.ca)model capable of representing well the modeled phenomena.The second aspect is the data assimilation approach combining a numerical model with direct observations of ocean state.Many authors have already demonstrated the suitability of satellite altimetry for the estimation of ocean circulation variability (Fu and Chelton, 2001;Fukumori, 2001;Le Traon and Morrow, 2001).A major advantage of satellite observations are their global coverage, continuity and synopticity.This is especially important for studies of the Southern Ocean where observations with conventional technique are sparse.However, it is known that the use of only temporal sea surface height anomalies obtained from altimetry is not enough to correct the mean ocean state (see, e.g., Kivman et al., 2005).To this end assimilation of the absolute dynamical topography is necessary.This presents a difficult task because ocean general circulation models commonly show systematic deviation from the mean dynamic topography.
On the model side, this systematic deviation can be caused by a variety of reasons including poor knowledge of surface forcing and details of the bottom topography representation selected during the model design.Slow equilibration processes mediated by baroclinic Rossby waves can lead to changes in barotropic circulation via changes in deep pressure across the major topographic features and, in this way, to systematic difference between the observed mean dynamic topography and model elevation.
This is a well known problem which can be partly solved by adjusting the model forcing.Ji and Leetmaa (1997) address this problem by a careful choice of the time-mean wind stress used to force the model and by calculating observed anomalies and assimilating them as anomalies for their model climatology.Yu and O'Brian (1991) tried a variational scheme in which the wind-stress field is a control variable.Derber (1989) developed a continuous variational Published by Copernicus Publications on behalf of the European Geosciences Union.
scheme specifically to address the problem of model bias.Frieland (1969) proposed augmenting the state vector by a model bias.Dee and da Silva (1998) have applied this idea to numerical weather prediction.Bell et al. (2004) built on the ideas of the latter work and applied a term of "pressure correction" combining it with the augmenting state vector concept.
A natural way to address this problem would be to use a variational assimilation scheme, where the significant systematical errors can be eliminated by slight changes in temperature and salinity profiles at ocean depth (Wenzel et al., 2001;Stammer et al., 2002).However, such a procedure is hard to set up and computationally expensive.This is why it is useful to look for alternative schemes that are numerically less expensive and easy to implement.In this paper, we will use a sequential technique based on the sequential evolutive interpolated Kalman filter (SEIK) (Pham, 2001).
Further, an open problem is how to project the assimilation update derived from the altimetry data into the interior of the ocean.In order to resolve this problem De Mey and Robinson (1987); Dombrowsky and De Mey (1992); Fischer and Latif (1995) used vertical Empirical Orthogonal Functions (EOFs) for the assimilation of altimetric sea level anomalies (SLA) into limited-area numerical ocean models.They have shown that the use of vertical EOFs is very effective and applicable to assimilate satellite products.However, an open question of these works is the nature of the space in which the statistics are calculated.Gavart and De Mey (1997) introduced the isopycnal EOF method for data assimilation.The authors used isopycnal EOF profiles calculated from hydrographic data in the Azores region.The use of isopycnal EOFs proved efficient in capturing the vertical structure of dominant processes in ocean, a quasi-homogeneous vertical displacement of isopycnals with quasi-conservation of water masses and potential vorticity on the isopycnal grid.Besides, isopycnal EOFs appeared to be more observable from altimetry than vertical EOFs.
Another approach of vertical update projection is based on a water property conservation principle and exploits the idea of water column vertical shift (Cooper and Haines, 1996).One more approach departs from the fact that extracting an efficient and coherent description of the large-scale oceanic thermohaline field is useful for data assimilation (Sun and Watts, 2001).The authors proposed a gravest empirical mode approach by projecting hydrographic profiles onto a geostrophic stream function for a Southern Ocean application.Finally, Fukumori et al. (1999) proposed the method of vertical modes decomposition where the temperature and salinity are updated according to the vertical structure of the first baroclinic mode.In this work we rely on this latter method because the empirical approaches require additional data profiles to represent the variability.This fact makes the use of such methods in global ocean models problematic at present time.Besides, the method of Fukumori et al. (1999) is easy to implement and computationally effective.
The goal of this paper is exploring the feasibility of assimilation of the global altimetric signal based on sequential assimilation technique.As follows from the discussion above this requires different technologies with respect to the systematic (mean) part of the difference between the model and observational data, and its variable part.The success of the approach suggested in this paper is based on two key elements.
First, recognizing that systematic drift of mean surface elevation in numerical ocean circulation models is associated with systematic changes in their thermohaline structure we suggest to use the methods of adiabatic correction of the model (Sheng et al., 2001;Eden et al., 2004) to effectively reduce systematic difference between mean state of the model and the mean state derived from the dynamic topography.This method of correction consists in representing the dynamically active density field as a linear combination of the density given by the model ρ m and the climatological in-situ density ρ c : On average, Sheng et al. (2001) find α=0.5 to be the most appropriate value in Eq. ( 1).Second, the proper functioning of the sequential assimilation part requires of an appropriate algorithm of mapping the statistically derived surface elevation updates into updates of the temperature and salinity fields.This paper uses the method proposed by Fukumori et al. (1999) according to which the temperature and salinity updates follow the first baroclinic mode in the vertical direction.
We demonstrate that combining these two key elements allows one to noticeably reduce the difference between the model state and observations.This reduction, however, remains suboptimal partly because the vertical modes used in the method by Fukumori et al. (1999) deviate from real modes of variability.The latter are affected by thermal wind and variable bottom topography and are sensitive to the horizontal size of perturbations.The other reason is the remaining systematic bias which cannot be fully removed by the adiabatic correction method.
The structure of the paper is as follows.Sections 2 and 4 present the ocean model and the assimilation methodology.The observational data are compared with the model in Sect.3; Sect. 5 describes the experimental setup and Sect.6 presents the results of the experiments.The final section concludes.

Ocean model
The ocean model used to perform this study is the Finite-Element Ocean circulation Model (FEOM) (Wang et al., 2008;Danilov et al., 2004).The model is configured on a global almost regular triangular mesh with the spatial resolution of 1.5 • .There are 24 unevenly spaced levels in the vertical direction.The model uses continuous linear representation for the horizontal velocity, surface elevation, temperature and salinity, and solves the standard set of hydrostatic ocean dynamic primitive equations.It uses a finite-element flux corrected transport algorithm for tracer advection (Löhner et al., 1987).
The model is forced at the surface with momentum fluxes derived from the ERS scatterometer wind stresses complemented by TAO derived stresses (Menkes et al., 1998).Vertical mixing is parameterized by Pakanowsky-Philander scheme (Pakanowski and Philander, 1981).The thermodynamic forcing is replaced by restoring of surface temperature and salinity to monthly mean surface climatology of WOA01 (Stephens et al., 2002).The model is initialized by mean climatological temperature and salinity of Gouretski and Koltermann (2004).
This configuration of the model is further referred to as V 1 .It is used only for comparison.The assimilation experiment is based on a configuration where the adiabatical pressure correction method of Sheng et al. (2001); Eden et al. (2004) (see Sect. 1) is applied.This configuration is denoted as V 2 (we have found that introducing this correction strongly reduces the drift of model state with time).Climatologic fields used for this correction are monthly mean temperature and salinity of Gouretski and Koltermann (2004).The success of this correction is discussed in the Sects.3 and 5.Both model variants are spun up from the state of rest over a 10year time period before the 1-year experiment is run.

Comparison of model results with observations
The observational data assimilated in the present work are a combination of the simultaneously measuring altimeter missions ENVISAT, GFO, Jason-1 and TOPEX/Poseidon missions interpolated onto the model grid so that the observations are available at every point of the model grid every ten days.The altimeter data of these missions were homogenized (with respect to the same reference ellipsoid and time scale), upgraded with the most recent corrections and models, and cross-calibrated.Most important, the inverted barometer correction for all missions was replaced by the dynamic atmospheric corrections (DAC) produced by CLS Space Oceanography Division using the MOG2D model from Carrère and Lyard (2003).Ocean and loading tides were computed using the FES2004 model (Lyard et al., 2006).The cross-calibration was carried out by means of the discrete multi-missions crossover-analysis (Bosch, 2007;Bosch and Savchenko, 2007) in order to remove the mission specific relative biases and to correct radial errors.The altimetric sea surface heights were finally reduced by an equipotential surface (geoid) to construct the absolute signal, the dynamical topography of the sea surface.The geoid used for this purpose, EIGEN-GL04S1 gravity field model (Förste et al., 2008) is obtained from the Geo-Forschung Zentrum (GFZ), Potsdam, Germany.The data cover the period between January 2004 and January 2005.
First of all, the sea surface topography was computed on the ground tracks of the missions similar to the approach described in (Albertella et al., 2008).The geoid heights were spectrally smoothed by a Gaussian filter (Wahr et al., 1998) with the filter length of 175 km applied on the spherical harmonic coefficients in order to remove the striping patterns plaguing the gravity fields derived exclusively by data of the GRACE gravity mission.To achieve spectral consistency the same filter length was used for the smoothing of the altimetric sea surface heights, which are smoothed along the ground track of satellite missions.The systematic difference between one-dimensional filtering of sea surface heights and two-dimensional (spectral) filtering of geoid heights was corrected by means of so called filter corrections, which were derived as differences between the mean sea surface CLS01 (Hernandez and Schaeffer, 2000) smoothed along the profiles and two-dimensionally.The geodetic sea surface topography was defined as difference between smoothed sea surface heights corrected by means of filter corrections and smoothed geoid heights.As for the assimilation the geodetic sea surface topography in the form of 10 day grids were required the block mean values of the sea surface topography available along track were computed for each time span of ten days.
The ice contamination makes the altimetry data unusable in the polar areas as far as the Sea of Okhotsk.This is the main reason of the lack of data in the Southern ocean where the ice covered area drifts zonally with the seasonal cycle so that the significant part of the surface appears to be icecovered at least for some time during the year.Hence, the variability contained in such data can not be reliable to use it in data assimilation problems.
The Indonesian region is characterized by complex bottom topography where neither geoid measurements nor model results appear to be accurate enough, and where huge discrepancies between them are observed.This is also true for the Mediterranean Sea.Hence, the observational data in these areas were substituted by the values of the RIO05 mean dynamical topography (MDT) (Rio and Hernandez, 2004;Rio et al., 2005).In order to prevent possible discontinuity appearing at the boundary of the two areas we applied a filter which smoothly extrapolates our measured dynamic topography towards that of RIO05 in the transition zone.These areas are shown in the top right panel of Fig. 1 as the deepblue rectangular areas.
The first necessary test is the comparison of the ocean model state with the dynamical topography data to be assimilated.The top left panel of the Fig. 1 shows the mean dynamical topography obtained from observations covering the period from January 2004 till January 2005.The top right panel of Fig. 1 shows the standard deviation for the same period.The mean dynamical topography corresponds to the relatively broad Gulf Stream and Kuroshio and smooth Antarctic Circumpolar Current.It qualitatively agrees with the model (not shown) and does not contain finer details than the model state.In this respect the model resolution is adequate to capture the main features contained in the data.
After a 10-year spin up from the state of rest, versions V 1 and V 2 are run for one additional year and the output is stored every 10 days.It is used to compute the mean and standard deviation presented in the middle and bottom rows of Fig. 1.The middle left panel of the figure depicts the difference between the mean dynamical topography obtained from observations and the mean calculated from V 1 model run.The difference is significant in many places over the world ocean, reaching ±0.5 m in some areas.The pattern of these discrepancies follows in many places the significant elevation of bottom topography.Many of them develop as the result of model adaptation to the bottom topography and present a systematic model bias.Probably, many other issues can be responsible for this bias.An apparent issue is the lack of realistic forcing and the absence of coupling with an ice model.In this way we cannot exclude that much finer tuning of the model is required to minimize the discrepancy between the model and data mean SSH fields, and that bottom topography is an important, but only one of many key features.Clearly, so large mean differences cannot be compensated by local changes in the steric height which makes further application of sequential data assimilation questionable.
The middle right panel of Fig. 1 shows standard deviation of the sample for the model version V 1 .As seen from the figure, the V 1 model variability differs from the variability of the observations but remains on the same scale reaching 12 cm in the most energetic areas.Compared to the data, the model overestimates variability in the Gulf Stream and Agulhas currents regions, but underestimates it in the area of Kuroshio current.The Tropical belt variability is smaller too than in the observational data.In the Antarctic Circumpolar Current (ACC) region there are regions of overestimated variability in V 1 compared to observations.
The strong systematic difference between the observations and the model demands seeking for the methods capable to reduce the model drift.The adiabatic pressure correction suggested by (Sheng et al., 2001;Eden et al., 2004) works through modifying the velocity fields leaving consistent tracer fields.In this way it suppresses the temperature and salinity advection by the erroneous velocity field which happens in the model which does not use such a correction.The utility of this correction method was consistently demonstrated in Eden et al. (2004), and is also exploited here.
The circulation in configuration V 2 is adiabatically corrected toward the velocity field which would correspond to the climatological temperature and salinity.The bottom panels of Fig. 1 show mean difference (left) and standard deviation (right) calculated using this model configuration (see Sect. 1).The difference in the mean fields is reduced compared to the V 1 version in all regions, especially in North Atlantic.Although it still remains relatively high it is much closer to the level comparable to typical variability.Thus applying the adiabatic correction by modeling the density as a linear combination of the climatology and the model density helps to reduce the deviation of the model mean compared to the mean obtained from the satellite observation.
The price paid for this reduction is that it simultaneously reduces the variability and the standard deviation of V 2 is much smaller than that of the observations (top right) and V 1 model run (middle right).Although the main regions of high variability are still visible, the amplitude is significantly smaller.The physical reason of this effect is also clear -by correcting pressure one affects (reduces) the amplitude and phase speed of baroclinic Rossby waves and in this way a  Eden et al., 2004), but it is not attempted in this paper.This fact limits the use of the proposed method for seasonal numerical experiments.To prevent the reduction of inter-annual variability it is also of interest to try more sofisticated correction methods such as (Bell et al., 2004).
Let us also consider the differences in the states obtained after ten years of the spin up with both versions of the model.The left plot of Fig. 2 shows the difference between the initial state of the SSH (here and below we use the abbreviation SSH, sea surface height, to refer to the dynamical topography) obtained via the 10-year V 1 model spin-up and the corresponding observed initial state.This plot reveals large discrepancies between the V 1 model state and the observed dynamical topography mainly in the Southern Ocean, North Atlantic and Pacific.
The 10-year V 2 model spin-up (right panel in Fig. 2) with the adiabatic pressure field correction results in a state which is much closer to the observed state.Overall, we see that the benefits of the original V 1 model are carried over to the model V 2 .However, the discrepancies compared to the observations are still present, although their magnitudes are significantly reduced almost everywhere.Hence, the V 2 model state will be used as the first input to the model run in our assimilation experiment.

Assimilation scheme
In this study we use the sequential evolutive interpolated Kalman filter (SEIK) introduced by Pham (2001).This method has been used in a number of studies (Hoteit et al., 2007(Hoteit et al., , 2005;;Triantafyllou et al., 2003;Nerger et al., 2007;Nerger, 2004).In the SEIK algorithm (see Pham, 2001), the forecast field is computed as an average over the ensemble members, and the forecast error covariance matrix is obtained as the corresponding covariance matrix from the ensemble.Since the forecast error covariance matrix has a rank that depends on the number of ensemble members, it is represented in the algorithm in its reduced form.This allows that the analysis covariance matrix P a k be calculated in its reduced form too.The calculated analysis error covariance P a k is then used to obtain the analysis using formulas: Here, η f k , η a k and η o k respectively denote forecast, analysis and observations of SSH field at time t k ; H k in our case is identity since the observations are interpolated prior to the assimilation onto the model grid.The matrix R k is the observation error covariance matrix that is diagonal, with the diagonal values equal to 25 cm 2 .Considering the model performance, the largest contribution to matrix R is associated with model (representativeness) and probably not data error in satelite altimetry or the geoid used.The diagonal nature of R is by convenience only.Note however, that the assimilation increment (change in model state) is correlated due to the presence of the forecasted error covariance matrix.Once the analysis is completed, the second-order accurate sampling technique is used for generation of new ensemble members that have the mean and covariance equal to η a k and P a k .This analysis ensemble is propagated with the full nonlinear model to the next assimilation time step.In this study, we apply the local version of the SEIK filter (Nerger et al., 2006) so that the analysis for each water column of the model depends only on observations within a specified influence region.Here, the influence region is a circle with a radius of 200 km.This value was chosen due to preliminary experiments with different subdomains.The experiment with the radius value of 200 km showed better results.In practice, except for high latitudes, it indeed corresponds to taking into account only all the nearest neighbours.
The strategy of our assimilation experiment is as follows.In order to minimize the deviation of the model mean from the observed mean we use the V 2 version of the model.At each time the observations are available, the analysis of the SSH field is carried out applying the local SEIK filter.Using this information, the vertical profiles of temperature and salinity are updated according to the vertical structure of the first baroclinic mode (Fukumori et al., 1999) with the amplitude computed from the elevation update, i.e. the temperature and salinity fields are updated using the following formulas: T a (x, y, z) = T f (x, y, z) + δη(x, y) gρ 0 ĥ(z) p(0) S a (x, y, z) = S f (x, y, z) + δη(x, y) gρ 0 ĥ(z) p(0) Here, overbars denote the reference state calculated as a mean from the one year free model run of the V 2 model, p and ĥ are locally defined vertical structures of the first baroclinic modes of velocity and displacement calculated using the local vertical profiles of the Brunt Väisälä frequency and density from the V 2 model (Gill, 1982).The function δη(x, y) is the analysis increment, i.e. the difference between the analysis of SSH, η a k and its forecast η f k , and g is the acceleration due to gravity.The velocity field is left unchanged so that it is simply the result of the model evolution.
In order to generate the initial error covariance matrix a model run was performed to produce a set of 10-day forecasts from a series of initial conditions distributed at 10-day intervals over a year.This procedure produces an ensemble of 37 model states.
The 10-day interval is chosen to correspond to the frequency of observational updates in the assimilation experiments, and the series of the initial conditions is used so that the covariance of the ensemble will be representative of the  1).The green and yellow solid lines show the errors corresponding to the V 1 and V 2 free simulations (without assimilation), respectively.The blue lines with bullets represent the V f 3 10-day model forecasts, while the dotted red lines correspond to the V a 3 analysis.V 1 model which was spun up was re-initialized at the initial time with the initial condition of V 2 .full period of the experiments.The initial error covariance matrix is then approximated with a lower rank matrix using the first eight empirical orthogonal functions (EOFs) of the ensemble.The first eight EOFs represent more than 90 percent of the variability.This covariance matrix is a consistent estimator of the 10-day model error covariance, and is adequate for parameterizing the filter background error covariance.The initial field is taken as a result of a 10-year model spin up of the V 2 as described in the previous section.Eight ensemble members are used in the implementation of the local SEIK algorithm.

Hindcast experiment
Three simulations were performed for the period between January 2004 and January 2005.The first two are the already mentioned V 1 and V 2 model configurations described in Sect. 2. These simulations were free model runs, i.e. model integrations within the chosen time period without data assimilation.V 1 model which was spun up was re-initialized at the initial time with the initil condition of V 2 which was chosen as initial condition for all the models in our experiments.Fig. 3 shows the evolution of the RMS error of the SSH averaged over the entire ocean (except for the zones corresponding to the RIO05 MDT location in the data see Fig. 1).The green and yellow lines show the errors corresponding to the V 1 and V 2 free model simulations, respectively.Overall, we observe that in the V 2 run there is no extra model drift.The results of the V 2 are stable, and the error standard deviations are generally 5 cm lower than those of the V 1 free model run.Note that this is a global result.A similar error analysis in only the ACC area reveals a more significant improvement in the V 2 run on the order of 8 cm.
The third simulation denoted V 3 is characterized by (i) presence of the sequential data assimilation applied to the SSH field (red circles in the Fig. 3); let us denote the corre-sponding state of the sequential analyses as V a 3 ; (ii) the vertical mode update for the temperature and salinity fields as described in Sect. 4 (and horizontal structure obtained from the statistical innovation of the SSH); and (iii) adiabatic pressure correction as in the V 2 model.Hence, the initial states for every model forecast (blue circles in Fig. 3) are superpositions of statistical updates from the filtering applied on the SSH and the results of the vertical mode update applied to the temperature and salinity fields.Of course, the correction term applied to the model also influences the model forecast evolution.The model states corresponding to the forecast of the V 3 model will be denoted as V f 3 .The statistical analysis of the SSH V a 3 is stable and provides accurate estimates of the SSH during the whole period of our experiment (Fig. 3).The RMS error with respect to the observations for the V a 3 decreases monotonically from 8 to 2 cm beginning from the fourth observational update.The V f 3 model forecasts are also stable and the RMS error is generally 7 cm and 2 cm lower than those of the V 1 and V 2 free model run RMS errors, respectively.
The initial background covariance matrix is imperfect.However, its evolution in time leads to more appropriate error representation which is evidenced by monotonic decrease of the analysis error as the time goes on.On the other hand, the improved initial states (red dots) of the V 3 analyses lead to improved V 3 forecasts.Their quality is then always better than the quality of the V 2 forecasts.At the first analysis (day 10) the fit to observations is better than the fit of the following analysis steps.The price for the good fit is that innovation is propagated via the imperfect error covariance to the full model state which reacts accordingly by deviating from a balanced solution.Later the covariance becomes more "educated" and leads to increments that are better sustained.
The quality of the model forecasts is always better than the quality of the V 1 and V 2 free model simulations.the drift of forecast V f 3 towards the V 2 state is still significant which can partly be explained by the fact that the vertical mode structure does not always catch all the thermodynamic features like, for instance, the horizontal shift of the isopycnals (the first baroclinic mode used in the method by Fukumori et al., 1999 accounts for only vertical displacement).Hence, there is still an uncertainty in the resulting temperature and salinity profiles that are not always fully consistent with the updated SSH field.On the other side, the temperature profiles obtained from the modal structure are stable and realistic.
The upper plots of the Fig. 4 show spatial distributions of the RMS errors of the SSH with respect to the observations for the V 1 (left) and V 2 (right) model versions.The bottom plots in the figure correspond to the V f 3 (left) and V a 3 (right).These plots demonstrate the ability of the model to correct for the mean state, which was not subtracted in the RMS calculations.
Overall, we see that a significant part of the systematic error is already corrected by the V 2 model compared to the original V 1 model.However, the main discrepancies in the Southern Ocean, Tropical and Northern Pacific and the Gulf Stream area are still observed for the V 2 version.The V f 3 SSH is much more accurate than the V 2 free run.However, the RMS error in the Gulf Stream area is even larger than for the V 2 version.Apparently, the V a 3 SSH field is the most accurate estimate in our experiments, although the RMS error in some places of the Gulf Stream region reaches 10 cm.
Figure 5 shows spatial distribution of the SSH standard deviation with respect to the observations.Here the mean observed signal is excluded and panels show the ability of each model to correct for the temporal variability.As we see from the figure, the V 2 model does not improve significantly the temporal variability of the SSH in high-energetic areas of the world ocean.The temporal variability inferred from the V f This improvement can be explained by the fact that every 10day model forecasts start from the corrected model states V a 3 (bottom left plot in Fig. 5) which is accurately estimated in the sense of temporal variability and the mean-state error correction (cf.Fig. 4).The slight problems of the SSH analysis in the Gulf Stream area are also visible in this figure.
Finally, let us compare the differences in the mean and standard deviation between the V a 3 and V f 3 and the observations as done in Fig. 1 for the V 1 and V 2 versions of the model.In Fig. 6 on the left the differences between the mean SSH obtained from the observations and the mean calculated from the V a 3 (upper panel) and V f 3 (lower panel) samples are shown.The difference between the mean of the analysis estimates and the mean of the observations is very small throughout the world oceans.However, some discrepancies with small amplitudes can be seen in the Gulf Stream region and the ACC region.Although the difference in the mean of the forecasts sample and the observations is much larger than the analysis difference, we still see an improve-ment in the mean field compared to the V 2 version of the model.Comparing the variability of the SSH in the V 1 , V 2 (Fig. 1  which has to be solved before applying the sequential data assimilation technique is the reduction of systematic errors between the data and the model.We demonstrate that the method of adiabatic pressure correction proposed by Sheng et al. (2001); Eden et al. (2004) suppresses strong systematic deviations so that they remain on the level comparable to the variable part of the difference between the model and observations.Further improvement of the estimates of the ocean state is achieved by assimilating SSH satellite data using the local SEIK filter.The common problem encountered when applying such techniques in conjunction with the satellite altimetry data is poor covariances whereby statistical update of temperature and salinity has flaws and is also conductive to numerical instabilities.To overcome this difficulty this work uses the method proposed by Fukumori et al. (1999).According to this method, the temperature and salinity update are assumed to follow the vertical structure of the first baroclinic mode.In our experiments, we tried different sets of modes: the first, a superposition of the first two, five and all the possible modes.And we found that the first mode gives better results.
Taken together, these three techniques (adiabatic pressure correction, local SEIK filter and the method by Fukumori et al., 1999) allow a successful reduction of the errors.Having applied them together, we have managed to decrease the global mean RMS error of the reference model from 16 cm to 9 cm.
It is also clear that despite its success the approach suggested here remains suboptimal.One of the reasons why it is so is that the vertical structure of the first baroclinic mode is computed without taking into account the vertical velocity shear (the horizontal density gradient).This changes the structure of modes (compared to the true vertical structure) reducing the amplitude of temperature and salinity updates in the upper layers.Computing the true vertical structure is much more difficult and additionally it turns out to be dependent on the horizontal wavenumber of perturbations.
The other reason is applying the adiabatic correction.Although it essentially reduces the model drift and in this way is indispensable, it also reduces the sensitivity of the velocity fields to the temperature and salinity updates.In this way, even if the vertical structure of the baroclinic mode were perfectly known, the elevation which will be in equilibrium with the updated temperature and salinity is different from the updated elevation.The difference between the elevation just after the update and its equilibrium value propagate as a surface wave leading to error growth in the forecast phase.
These issues as well as several other questions remain for the future research.Namely, an impact of the ocean state update on the ocean circulation and heat and salt content changes implied by data assimilation should be analyzed and compared to results of variational data assimilation techniques.Also, one needs to find out how the use of a temporally varying geoid can help in a more precise estimation of the general ocean circulation.

Fig. 1 .
Fig. 1.The upper row: The mean dynamical topography (left) and standard deviation (right) for the period from January 2004 till January 2005.The deep-blue rectangular areas correspond to the locations where the RIO05 MDT was substituted in the data (no variability).The middle row: The mean difference between the dynamical topography obtained from the observations and V 1 model run (left) and the standard deviation of V 1 (right).The bottom row: The same as in the middle row, but for V 2 version of the model.

Fig. 2 .
Fig. 2. The difference between the first observation and initial model state obtained as a result of 10-year spin up with V 1 model (left) and V 2 model (right).

Fig. 3 .
Fig.3.Evolution of RMS error of SSH for the world ocean (except zones corresponding to RIO05 MDT location in the data see Fig.1).The green and yellow solid lines show the errors corresponding to the V 1 and V 2 free simulations (without assimilation), respectively.The blue lines with bullets represent the V

Fig. 4 .
Fig. 4. Spatial distribution of RMS errors of the SSH fields with respect to observations.The upper plots are for the V 1 (left) and V 2 (right) model versions.The bottom plots correspond to the V f 3 (left) and V a 3 (right) estimates.

Fig. 5 .
Fig. 5. Spatial distribution of standard deviation on SSH with respect to observations.The upper plots are V 1 (left) and V 2 (right).The bottom plots correspond to the V f 3 (left) and V a 3 (right).

Fig. 6 .
Fig. 6.Difference on mean SSH between: V a 3 (upper line, left), V f 3 (bottom line, left) and observations.The standard deviation (with respect to its own mean) for V a 3 (upper line, right), V f 3 (bottom line, right).