Assimilation of SLA along track observations in the Mediterranean with an oceanographic model forced by atmospheric pressure

A large number of SLA observations at a high along track horizontal resolution are an important ingredient of the data assimilation in the Mediterranean Forecasting System (MFS). Recently, new higher-frequency SLA products have become available, and the atmospheric pressure forcing has been implemented in the numerical model used in the MFS data assimilation system. In a set of numerical experiments, we show that, in order to obtain the most accurate analyses, the ocean model should include the atmospheric pressure forcing and the observations should contain the atmospheric pressure signal. When the model is not forced by the atmospheric pressure, the high-frequency filtering of SLA observations, however, improves the quality of the SLA analyses. It is further shown by comparing the power density spectra of the model fields and observations that the model is able to extract the correct information from noisy observations even without their filtering during the pre-processing.


Introduction
The Mediterranean Forecasting System (MFS) (Pinardi and Coppini, 2010) each day provides 10-day-long oceanographic forecasts of the Mediterranean Sea.The forecasts are initiated by the analyses of the state of the ocean at the daily frequency.The analyses combine estimates by the ocean model (Oddo et al., 2009), satellite observations of sea surface temperature (SST) and sea level anomaly (SLA), and in situ observations of temperature and salinity by Argo floats and of temperature by expendable bathythermograph (XBT) instruments (e.g.Dobricic et al., 2007).
The SLA along track observations have the largest impact on the state estimates by the MFS system.Their important impact arises from the fact that, in addition to barotropic effects, SLA represents a vertically integrated effect of variations in temperature and salinity over the whole water column.While in situ observations provide more accurate direct observations of the temperature and salinity variability along the water column, their spatial scarcity limits their impact only to small areas of the Mediterranean.There are also many SST observations, but their impact is often limited only to the estimate of the temperature distribution close to the surface of the ocean.On the other hand, the SST observations can be combined with in situ observations of vertical profiles to get a three-dimensional estimate of the state of the ocean.Even in this case, the ability to construct three-dimensional estimates is, however, limited due to the small number of in situ observations.Another important property of SLA observations is that they provide information on geostrophic component of surface currents, a product that is used in a number of applications requiring MFS analyses and forecasts.For example, surface currents are implemented to forecast oil spill movement, assist search and rescue operations and detect sources and the spreading of the pollution.Figure 1 shows the spatial and temporal coverage by SLA and in situ observations in 2009.It can be seen that the SLA observations repeatedly observe the whole Mediterranean Sea, whilst the in situ observations sample only some limited parts of the basin.The relative importance of the large number of SLA observations in the Mediterranean has been investigated in detail in Pujol et al. (2010).
Published by Copernicus Publications on behalf of the European Geosciences Union.The daily frequency of MFS analyses and the ability of the MFS system to simulate a large number of high-frequency oceanographic processes may require SLA products that are specially prepared for the MFS system.In particular, recently the atmospheric pressure forcing has been added to MFS model equations (Oddo et al., 2012).In a semi-enclosed basin like the Mediterranean, the inverse barometer estimate cannot fully remove the impact of the atmospheric pressure gradient, because the slow water exchange between the Mediterranean and the Atlantic through the narrow Gibraltar Strait often generates large-scale oscillations that can last for several days (e.g.Candela and Lozano, 1994).In the past, these oscillations were removed from SLA observations by adjusting several SLA satellite tracks in consecutive days (Le Traon and Gauzelin, 1997).Since 2004, the ocean response to the high-frequency atmospheric pressure, winds and tidal potential has been simulated with a barotropic ocean model, and this simulated high-frequency signal has been removed from SLA observations (Carrere and Lyard, 2003).
This procedure, however, creates inconsistencies when the MFS assimilation system uses the ocean model forced by the atmospheric pressure.In order to compare observations and background model fields, which contain the atmospheric pressure effects on sea level, it would be necessary to subtract from the model background field the same SLA high-frequency correction of the SLA observations.This quantity may be difficult to assess, and, therefore, the use of a more sophisticated ocean model forced by the atmospheric pressure gradient could even degrade the quality of the analyses.
In this study we will describe the procedure implemented for the assimilation of SLA observations when the model is forced by the atmospheric pressure.We will assimilate the TAPAS observations with and without the high-frequency correction.These experiments will help us to draw the guidelines for the future use of SLA along track observations in the MFS data assimilation system in which the model is forced by atmospheric pressure gradient.In Sect. 2 we describe the TAPAS data product, and the data assimilation system methodology for the use of the model forced by atmospheric pressure.Section 3 describes the setup of numerical experiments and compares the results of the assimilation with different model implementations, with and without the atmospheric pressure forcing and with TAPAS observations containing and excluding the high-frequency correction.Section 4 provides a discussion and offers some guidelines for the future use of the SLA products in the MFS data assimilation system.

SLA observations
An experimental product called "Tailored Altimetry Products for Assimilation Systems" (TAPAS) was developed that provides SLA observations without the high-frequency signal corrections, i.e. the inverse barometer corrections and others.The TAPAS data set further addresses issues encountered with advanced modeling systems that may contain most of the high-frequency signal in the sea level.A first TAPAS along-track product based on delayed-time (DT) observations was delivered in 2010 for the global ocean.These SLA observations are without along-track filtering and subsampling (raw 1Hz altimeter data at the horizontal resolution of 7.5 km).In addition, a new subset of TAPAS observations is produced over the Mediterranean Sea for this study.For each SLA observation, it gives the information on the modifications made in different stages of the processing: the longwave error correction, the low-frequency inverse barometer, the correction for ocean tides, and the high-frequency correction obtained after the application of a barotropic model forced by high-frequency winds and atmospheric pressure.This allows us to put the high-frequency signal back into SLA observations, and to estimate the impact of SLA high-frequency components on the MFS assimilation system.Therefore, with the TAPAS data set, we could test the effect of unfiltered SLA data in the MFS assimilation system that includes the atmospheric pressure forcing.The MFS data assimilation system uses the Mediterranean setup of the NEMO ocean model (Oddo et al., 2009).The model has the horizontal resolution of 1/16 • and 72 unevenly distributed layers in the vertical direction.In the Atlantic it is nested with the global 1/4 • model by Mercator-Ocean (Barnier et al., 2006;Oddo et al., 2009).The atmospheric turbulence fluxes are calculated in bulk formulae by using the atmospheric fields of wind, temperature, humidity and surface pressure from the ECMWF operational analyses, while short-and longwave radiation are parameterized with the use of the cloud cover operational analyses by the ECMWF (Pettenuzzo et al., 2010).Recently, the forcing by atmospheric pressure has been introduced into the MFS model.It is applied by adding an additional pressure gradient term into horizontal momentum equations (Oddo et al., 2012):

Ocean
where v is the horizontal velocity, ρ 0 the reference density, and p A the atmospheric pressure.

Data assimilation scheme
The OceanVar data assimilation scheme (Dobricic and Pinardi, 2008) is a variational scheme in which the slowly evolving vertical part of temperature and salinity background error covariances is represented by seasonally and regionally variable empirical orthogonal functions (EOFs), whilst their horizontal part is assumed to be Gaussian isotropic depending only on distance.The horizontal radius of correlation is set to 10 km.In the horizontal direction, we apply the isotropic covariances, because we assume that, due to the large variability of parameters at the high horizontal resolution of the model, it could become very difficult to correctly estimate the complex structures of the horizontal background error covariances by a set of climatological EOFs.The rapidly evolving part of the background error covariances, consisting of the sea level and the barotropic velocity components, is modeled in each step of the minimization algorithm by applying a barotropic model forced by the vertically integrated buoyancy force resulting from temperature and salinity variations.The baroclinic components of velocity are then estimated by applying the geostrophic relationship, modified along the coast in order to eliminate the horizontal divergence.In this way OceanVar combines a long-term three-dimensional variational scheme for the slow processes with a scheme that fully dynamically evolves the covariances by model equations for the fast processes.The scheme is minimized in iterations that stop when the gradient of the cost function becomes smaller than 1 % of the initial gradient.The background SLA estimate is formed by subtracting the mean dynamic topography (MDT) estimate from the background sea level estimate.The MDT is obtained by combining the estimate by Rio et al. (2007) with the information from the in situ observations in a procedure similar to the one applied in Dobricic (2005).
When the assimilation uses the high-frequency corrected SLA data and the ocean model is forced by the atmospheric pressure gradient, it is necessary to modify the model background SLA field in order to represent the same quantity that has been observed.It is then necessary to subtract an estimate of the sea level response to high-frequency atmospheric pressure forcing from the background SLA field.On the contrary, when the model is not forced by atmospheric pressure and the SLA observations contain the high-frequency sea level signal, this contribution must be added to the model background SLA field.The relationship between the model background field from a numerical model with and without atmospheric pressure forcing can be written as where η A is the background sea level estimate produced by the model with the atmospheric pressure forcing, η NA is the same estimate produced by the model without the atmospheric pressure forcing, p A is the atmospheric pressure field, g is the acceleration due to gravity, and A(t) is a horizontally constant field over the whole model domain, which includes the Mediterranean and a part of the Atlantic Ocean.
The second term on the right-hand side of Eq. ( 2) is the socalled inverse barometer effect that assumes the full isostatic balance between the atmospheric pressure and the sea level (e.g.Ponte, 1993).In addition, term A(t) is estimated by Term ε represents differences due to the nonisostatic response of the model forced by the atmospheric pressure in the Mediterranean Sea.Equations ( 2) and ( 3) are also applied to model background estimates of SLA fields obtained after subtracting the MDT from the background sea level.
In the assimilation scheme, the mean residual is subtracted from the residuals along each SLA satellite track.In this way the unknown steric height signal is removed together with the largest-scale oscillations that may originate from the atmospheric pressure forcing, or from other small-scale processes, like for example the local variations of wind near the Gibraltar Strait, which may be unresolved by the ECMWF analyses (Fukumori et al., 2006).Residuals are calculated at the observational time, and the atmospheric pressure available from the ECMWF analyses at the interval of 6 h is interpolated in time to the observational time.This procedure can produce aliasing of the atmospheric pressure signal due to the atmospheric tides (e.g.Ponte and Ray, 2002), but, since this part of the signal has a very large spatial scale, it is most likely also removed by the subtraction of the mean residual along the satellite track.It should be noticed that the model has an implicit numerical scheme for the adjustment of the sea level and the vertically integrated velocity.Therefore, gravity waves in deep ocean are strongly damped instead of being dispersed.

Scales of sea level differences in simulation and assimilation experiments
Six experiments are performed to estimate the impact of the atmospheric pressure forcing on the simulated sea level and in the analyses produced by different models and SLA observation processing methods.Two simulation experiments, without data assimilation, are carried out with (SIM1) and without (SIM2) atmospheric pressure forcing only in January 2009.The four assimilation experiments cover all of 2009 and include models with (ATPR) and without (CONT) atmospheric pressure forcing and with (ATPR1, CONT1) and without (ATPR2, CONT2) high-frequency corrections to SLA observations.Table 1 summarizes the differences between the experiments.In all of these experiments, the parameters of the data assimilation scheme are kept invariant, and only observations and model are changed.All effects that could appear due to the non-linear interaction between the components of the data assimilation scheme only due to the change of the model and observations are assumed to be negligible.
When the sea level estimate from SIM1 is compared with SIM2 by using Eq. ( 2), the difference between the two solutions gives the term ε.The field of ε is shown in Fig. 2. It can be seen that the differences mostly have large scales.In particular, the basin scale dominates the differences, and there are gradients between the western and the eastern basins, and further in the semi-enclosed Adriatic and Aegean Seas.In addition, there are some smaller-scale differences, concentrated mainly in the western basin, and especially along the Northern African coast.The smaller-scale differences along the Northern African coast can be explained by the propagation of Kelvin waves along the coast generated by the differences in the transport through the Gibraltar Strait.In general these differences are expected because, as mentioned above, the Mediterranean sea level response to atmospheric forcing is not well captured by the inverse barometer effect due to the slow exchange of water with the Atlantic Ocean through the Gibraltar Strait.
The simulation differences are further compared with differences between assimilation experiments ATPR1 and CONT1 in January 2009.Figure 3 shows the root-meansquare (RMS) differences of sea level estimates between the two simulations (SIM1-SIM2) and between the two analysis estimates in January 2009 (ATPR1-CONT1).It can be seen that, when averaged over the whole month, the differences between simulations have small amplitudes and they are at large scales (top panel in Fig. 3).This happens because the small-scale differences in snapshots shown in Fig. 2 are mainly due to the barotropic Kelvin waves that appear less frequently and are rapidly dispersed.Therefore, they do not leave a strong signal over a longer period of time.On the other hand, the monthly averaged RMS of differences between the analyses in ATPR1 and CONT1 contains both large and small scales (middle panel in Fig. 3).The large-scale structure of the analysis differences is very similar to the structure of the simulation differences.Clearly, the data assimilation does not change the sea level signal at the large barotropic scales which originates from the atmospheric pressure forcing.On the other hand, the combination of the information coming from the SLA observations and the differences in the background fields due to the atmospheric pressure forcing create different sea level estimates at smaller scales.In the analyses, the small-scale differences propagate slowly, because they reflect the changes in temperature and salinity in the ocean's interior introduced by the assimilation of SLA observations.Therefore, their signal is clearly visible even when the RMS of differences is averaged over the whole month (bottom panel of Fig. 3).The next subsection will evaluate the impact of these differences on the accuracy of the analyses.

Analysis accuracy
The accuracy of the analyses is estimated by the evaluation of the RMS of residuals with respect to the SLA observations and the in situ observations of temperature and salinity by Argo floats.Assuming that the time between the two subsequent observations by the satellites or the Argo floats is only several days long, we may assume that the background errors are mainly due to the errors in the initial state represented by the analyses, and not due to the errors in the model integration or in the atmospheric forcing.Figure 4 shows the evolution of the RMS of SLA residuals in the year 2009 obtained by experiments ATPR1 and CONT1.It is shown in two panels, because each panel shows the comparison of background fields with a different SLA data set that is either filtered or unfiltered.It can be seen that the RMS of SLA residuals is very similar in all experiments.When averaged over the whole year, it was 4.00 cm for experiment ATPR1 and 4.05 cm for experiment CONT1.The difference is only about 1 % of the total RMS of residuals, but throughout the year 2009 the RMS of SLA residuals is lower when the model is forced by the atmospheric pressure gradient.In the second set of twin experiments, the analyses assimilated TAPAS observations of SLA with the high-frequency correction.This time the RMS of SLA residuals is 4.15 cm in experiment ATPR2 and 4.10 cm in experiment CONT2.On average, the CONT2 experiment gives 1 % lower RMS of SLA residuals than the ATPR2 experiment.The horizontal distribution of differences between the RMS of SLA residuals (Fig. 5) shows that, when the observations are not filtered, the forcing of the ocean model by the atmospheric pressure especially reduces the RMS of residuals in the central part of the Western Mediterranean.This area also shows large differences in January between SIM1 and SIM2 experiments (Fig. 2).The central part of the Western Mediterranean is characterized by a strong dynamical activity that appears to be sensitive to the perturbations in the atmospheric forcing (e.g.Pinardi et al., 2011).On the other hand, there is no particular area with large differences between the RMS of SLA residuals in experiments CONT2 and ATPR2, but in CONT2 the RMS of residuals has lower values over most of the Mediterranean.
The RMS of the temperature and salinity residuals (Fig. 6) is the lowest in experiment ATPR1 when the model is forced by the atmospheric pressure gradient and the observations are not corrected for the high-frequency dynamical signal.When the model is not forced by atmospheric pressure, experiment CONT1 on average gives a lower RMS of temperature and salinity residuals than CONT2, in which the observations are not corrected for the high-frequency signal (Fig. 6).It can be seen in Fig. 6 that the differences in the RMS of temperature and salinity residuals are most pronounced in the layer between 50 and 400 m.This layer roughly represents the depths of the first baroclinic mode in the Mediterranean (e.g.Molcard et al., 2002) that should be mainly affected by the differences in the assimilation of the SLA observations.There the RMS of temperature residuals in ATPR1 is 3-4 % lower than in CONT1, 3-6 % lower than in CONT2, and 3-9 % lower than in ATPR2.In ATPR1 the RMS of salinity residuals is 3-5 % lower than in CONT1, 4-8 % lower than in CONT2, and 6-10 % lower than in ATPR2.
The data assimilation combines the information from the observations with the information from the background fields and it produces analyses that contain the dynamical constraints given by the model dynamics.This feature of the MFS data assimilation system is demonstrated in Fig. 7.The power spectral density of sea surface height observations follows approximately the k −2.5 slope at scales larger than 50 km.This slope is flatter than the slope of k −3 that Fig. 3.The RMS of the differences between the sea level (cm) in January 2009 from simulations with and without the atmospheric pressure gradient forcing (top panel).The RMS of the differences between the sea level (cm) in January 2009 from the analyses with and without the atmospheric pressure gradient forcing (ATPR1-CONT1, middle panel).The RMS of the differences between the sea level differences (cm) in January 2009 obtained from the simulations and from the analyses with and without the atmospheric pressure gradient forcing (bottom panel).could be predicted theoretically for larger synoptic scales (e.g.Nastrom and Gage, 1985).In the Mediterranean, the wind is strongly modified by the surrounding topography, and therefore the topography may significantly impact the slope of the power spectrum at these scales.On the other hand at the shorter scales, we could expect that the power spectrum of the sea level depends less on the external forcing, and it should have a slope that is closer to the theoretical estimate for the geostrophic turbulence.At scales shorter than 50 km SLA observations show a power spectrum slope that is quite flatter than the theoretical one of k −5/3 , indicating the presence of the spatially correlated observational noise at these scales.Contrary to the observations, the power spectrum of the analyses (ATPR1) has a slope that is very close to the theoretical, since the model equations provide a dynamical constraint for the distribution of energy among the scales following the assimilation of observations.It be- comes flatter only at the scales comparable to or shorter than the model 3-grid spacing.It may be a consequence of a false energy aliasing due to the inaccurate simulation of the nonlinear energy cascade, or simply a consequence of an insufficiently selective horizontal diffusion operator.It should be noted that all experiments had similar power density spectra, and differences in forcing the model and processing the observations were not able to change significantly the energy distributions among the spatial scales.The comparison of the power spectrum densities from the model and the observations partly explains the ability of the model to incorporate most of the information from the observations when they contain the smallest amount of the preprocessing like in experiment ATPR1.When observations introduce variability that may partly come from errors, the model may significantly reduce those errors by imposing the dynamical constraints present in model equations.The comparison of analysis and observational power spectrum densities shows that the numerical model forced by the atmospheric pressure is capable of filtering correctly the spectrum of the noisy observational signal and that it is not necessary to filter observations in advance in order to remove the high-frequency atmospheric signal from them.

Discussion
The study describes the impact of the atmospheric pressure forcing in the numerical model associated with the assimilation scheme of MFS.The analyses that use non filtered SLA observations and atmospheric pressure forcing in the numerical model (ATPR1) show differences in SLA at the small spatial scales with respect to the simple assimilation case.As demonstrated in Sect.3.1, these SLA differences are due to the temperature and salinity changes in the water column at the mesoscales, introduced by the data assimilation scheme.
Due to the semi-enclosed nature of the Mediterranean Sea, the assimilation with the model forced by the atmospheric pressure gradient required the evaluation of the level of the SLA data set processing in order to achieve the most accurate ocean state estimates.It was shown that, when the model was forced by the atmospheric pressure, the unfiltered SLA data set leads to the most accurate analyses.The improved accuracy was obtained when compared both with satellite SLA and in situ Argo observations.On the other hand, when the model was not forced by the atmospheric pressure gradient, the assimilation of the SLA data set that included the highfrequency dynamic correction on average produced more accurate analyses for SLA, while there was slight reduction in accuracy for temperature and salinity fields.It was further shown that the analyses that use a dynamical model success- fully filter the noise present in the power spectrum of the SLA observations.Analyses had the correct slope of the power spectrum at shorter scales, up to the shortest scales hardly resolved by the model dynamics.They were able to extract the useful information from the noisy observations without the need for applying additional filters on observations.
It should be noticed that the filter for the high-frequency signal removes also the impact of the high-frequency wind forcing in the Mediterranean.Unfortunately, it was not possible to have observations only without the high-frequency atmospheric pressure forcing, and we have assumed that the differences in the analysis accuracy due to the highfrequency wind forcing removal were negligible.This assumption may be justified by the final result that shows that the best result is obtained with the model forced by the full atmospheric forcing and the observations that are not filtered.Fig. 7.The power spectrum density (cm 2 km −1 ) as the function of the wavelength of the background (black line, experiment ATPR1; see Table 1) and the observed sea level estimate (green line).Both spectra are calculated at the same observational points along the SLA tracks in 2009.Only tracks longer than 650 km were taken into account in order to reduce the effects due to the variable geometry of the Mediterranean.The observed sea level estimate is obtained by adding the MDT to the SLA.The logarithmic scale is applied on both axes.Red lines indicate the −2.5 and −5/3 slopes.
One major motivation for this study was to give indications for the future design of the most efficient oceanographic analyses of SLA data in the Mediterranean.The study shows that, in order to achieve the most efficient extraction of the information from the SLA observations, the model should include as much as possible all processes influencing the sea level variability in the Mediterranean in order to assimilate the observations that contain most of the original unprocessed signal.In particular, it was found that, in the experiment that used a model forced by atmospheric pressure and observations containing the full atmospheric signal at high frequencies, the RMS of SLA, temperature and salinity residuals was consistently lowered by a few percent with respect to other experiments in which either the model was not forced by atmospheric pressure or the high-frequency signal was removed from the observations.While this is a general result, it should especially apply to closed and semi-enclosed seas like the Mediterranean.In Europe there are several other seas in which sea level responds similarly to the atmospheric forcing, like the Baltic, the Adriatic and the Black Seas.

Fig. 1 .
Fig. 1.Number of SLA and in situ observations in the Mediterranean during 2009.Top panel shows the position of Argo (red open circles) and XBT (blue crosses) profiles in the year 2009, while the bottom panel shows the number of SLA observations along tracks from both Jason and Envisat satellites during 2009.Each Argo and XBT profile was performed only once.

Fig. 2 .
Fig.2.Snapshots of differences of sea level (cm) in January 2009 between simulations forced with and without the atmospheric pressure gradient.Simulations started on 1 January 2009, and snapshots are on days 15, 20, 25 and 30 after the start of the simulation.Isolines are drawn with the interval of 0.5 cm.

Fig. 4 .
Fig. 4. Top panel shows the RMS of SLA residuals (cm) in experiments ATPR1 (continuous red line) and CONT1 (continuous black line) during 2009, and the bottom panel shows the RMS of SLA residuals (cm) in experiments ATPR2 (dashed red line) and CONT2 (dashed black line) during 2009.The RMS of SLA residuals is calculated once a week.

Fig. 5 .
Fig.5.Spatial distribution of differences between RMS of SLA residuals (cm) in experiments ATPR1 and CONT1 (top panel), and between CONT2 and ATPR2 (bottom panel).The RMS of residuals is averaged over squares with 1/2 • resolution.The differences are displayed only in squares that had more than 100 observations during the year 2009.

Fig. 6 .
Fig. 6.The RMS of temperature ( • C) (top left) and salinity (top right) residuals calculated with respect to observations by Argo floats in ATPR1 (continuous red line), ATPR2 (dashed red line), CONT1 (continuous black line) and CONT2 (dashed black line) during all of 2009.The difference between RMS of residuals from ATPR1 and ATPR1 (continuous red line), ATPR2 (dashed red line), CONT1 (continuous black line) and CONT2 (dashed black line) normalized by the RMS of residuals in ATPR1 for temperature (bottom left) and salinity (bottom right).

Table 1 .
Summary of experiments.