Towards high resolution mapping of 3-D mesoscale dynamics from observations

The MyOcean R&D project MESCLA (MEsoSCaLe dynamical Analysis through combined model, satellite and in situ data) was devoted to the high resolution 3-D retrieval of tracer and velocity fields in the oceans, based on the combination of in situ and satellite observations and quasi-geostrophic dynamical models. The retrieval techniques were also tested and compared with the output of a primitive equation model, with particular attention to the accuracy of the vertical velocity field as estimated through theQ vector formulation of the omega equation. The project focused on a test case, covering the region where the Gulf Stream separates from the US East Coast. This work demonstrated that innovative methods for the high resolution mapping of 3-D mesoscale dynamics from observations can be used to build the next generations of operational observation-based products.


Introduction
Ocean science and operational oceanography are based on the analysis of the space and time distribution of the parameters that characterize the state of the sea under observation.In principle, these state variables could be estimated directly from measurements, at least for the physical component of the system, usually including velocity, pressure, temperature and salinity (density).The system evolution could then be forecasted through the fundamental laws of oceanic physics, provided the external forcings are also known (un-less a fully coupled ocean-atmosphere model is considered, in which case the state variables include part of the forcings).In practice, the determination of the distribution of the ocean state variables is restrained by sampling, instrumental and resource limitations, and their forecasting is complicated by the non-linear nature of the equations that drive the system, generally requiring numerical solution, and by the huge number of degrees of freedom involved in corresponding models.All these elements make the combination and analysis of the few available observations a challenge for physical oceanographers, especially if the phenomena of interest are ubiquitous in the global oceans and involve relatively small space and time scales, such as mesoscale processes.
The assimilation of in situ and satellite observations in high resolution (eddy resolving) numerical models is an efficient strategy to simulate mesoscale ocean dynamics (e.g.Stammer et al., 2010, and references therein).Data assimilation is based on statistical or variational approaches that combine the observational data with the underlying (approximated) dynamical principles that govern the evolution of the system, taking into account some estimates of the errors and uncertainties associated with both components.As a consequence, while data assimilation is clearly crucial to obtain accurate analyses and realistic forecasts, its results are clearly influenced by specific model configurations, which represent important sources of uncertainties (e.g.external forcings, parameterization of processes acting at scales smaller than the model grid resolution, choice of the grid and numerical schemes, assimilation technique, etc.).
Published by Copernicus Publications on behalf of the European Geosciences Union.

B. Buongiorno Nardelli et al.: Towards high resolution mapping of 3-D mesoscale dynamics
Conversely, a purely observation-based approach, aiming to retrieve the three-dimensional (3-D) mesoscale dynamics from a statistical and/or empirical combination of in situ and satellite observations, has received limited recognition until now.This is partly related to the limits of present observation-based products, namely to their relatively low resolution, and to the difficulties in providing any indirect estimate of the ocean currents going beyond the simple geostrophic velocities (e.g. as usually obtained from satellite altimetry or dynamic heights).Moreover, especially for model validation purposes, there is a general tendency to look at each measured parameter separately (i.e. through univariate approaches), instead of considering the combinations of available observations as incomplete realizations of the system state.Conversely, our experience of the way the oceans behave indicates that the effective degrees of freedom in the system are significantly less than the number of variables involved (i.e.covariance/autocorrelation reduces the effective number of independent data), so that multivariate reduced-space analyses can provide a more efficient description of the system state than univariate approaches.
In fact, even if different multivariate techniques have been proposed until now to retrieve 3-D fields of temperature and salinity from combined observations, only a few of these have been translated into operational products (e.g.Fox et al., 2002;Guinehut et al., 2004).One of these products, named ARMOR3D, has been included in the GMES My-Ocean project catalogue as a "core" product (see also Guinehut et al., 2012), but its spatial resolution can be presently classified only as "eddy permitting", not exceeding 1 / 3 • .
The development of higher resolution observation-based 3-D fields might thus contribute to a more efficient analysis of observations and model validation through more advanced comparisons than those based on climatologic fields or sparse observations of temperature, salinity or velocities.
In this context, the MyOcean project (through its first open call for research and development) funded a small research initiative (MESCLA -MEsoSCaLe dynamical Analysis through combined model, satellite and in situ data, 2010-2012) devoted to the high resolution 3-D retrieval and analysis of tracer fields, horizontal and vertical velocities in the oceans based on the combination of in situ and satellite observations and simplified diagnostic models, and on their comparison with primitive equation model output.MESCLA rationale lies in the fundamental role played by the mesoscale in modulating the ocean circulation and the fluxes of heat, freshwater and biogeochemical tracers between the surface and the deeper layers.
In fact, in order to correctly retrieve the vertical velocities associated with mesoscale features, one should be able to resolve very small scales, i.e. down to less than 10 km (for a complete review on the topic, see Klein and Lapeyre, 2009).As said, present 3-D observation-based systems are far enough from being able to correctly reproduce the global variability at these scales (Willis et al., 2003;Roemmich and Gilson, 2009;Von Schuckmann et al., 2009;Guinehut et al., 2004Guinehut et al., , 2012;;Larnicol et al., 2006).Moreover, while the vertical component of ocean currents can be diagnosed in primitive equation numerical models by solving the continuity equation, the same technique is not applicable to direct observations.This is due, on one hand, to the few current measurements available, and, on the other hand, to the high error that would result from the computation of the divergence from measured horizontal velocities, which may include significant instrumental errors.It is also clearly impossible to use the continuity equation to estimate the vertical velocities from dynamic heights (computed from temperature and salinity profiles), as the geostrophic velocities are non-divergent by definition.On the other hand, simplified diagnostic models can be applied to retrieve the vertical velocity field from both 2-D and/or 3-D estimates of geostrophic currents and density fields (e.g.Tintoré et al., 1991;Allen and Smeed, 1996;Buongiorno Nardelli et al., 2001;Pascual et al., 2004;Klein et al., 2009;Isern-Fontanet et al., 2008;Ruiz et al., 2009).
In this framework, the first step in our work consisted of improving the existing MyOcean observation-based products (ARMOR3D, Guinehut et al., 2012) and in the development and testing of new high resolution horizontal interpolation and vertical extrapolation techniques (Buongiorno Nardelli, 2012;Buongiorno Nardelli et al., 2006;Buongiorno Nardelli andSantoleri, 2004, 2005), analysing the scales they are effectively resolving.As a second step, a quasi-geostrophic (QG) diagnostic numerical model (the Q vector formulation of the omega equation) has been used to estimate the vertical velocities (Pascual et al., 2004;Ruiz et al., 2009).The omega equation was applied to different MyOcean products (both model and observation-based) in order to quantify the differences/limitations in the diagnostic tools used and the impact of the spatial resolution on the retrieved velocity.The results of these two steps are the subject of the present paper.
It is worth noting that this work represents the first attempt to apply purely observation-based 3-D retrieval techniques at high resolution (also resolving mesoscale quasi-geostrophic dynamics) to obtain data that could be produced routinely within an operational system (namely from near real-time, freely available data, and potentially with global coverage).However, given the lack of independent (direct) measurements of the vertical velocities, a full validation of the new products is clearly not possible.Consequently, the approach followed here relies on the comparison between all the different products, concentrating on a test case.
The area selected for our test case lies in the North Atlantic between 32 • N-44 • N and 75 • W-40 • W, and covers the region where the Gulf Stream separates from the US East Coast, downstream of Cape Hatteras.The Gulf Stream is an extensive western boundary current that plays a fundamental role in the poleward transfer of heat and salt, and is one of the world's most intensely studied current systems.Several studies were conducted to evaluate the effect of the strong mesoscale activity associated with its flow on the vertical and cross-front exchanges (e.g.Bower, 1989;Bower and Rossby, 1989;Lindstrom et al., 1997;Joyce et al., 2009;Thomas and Joyce, 2010).It is also one of the few areas where direct measurements of the vertical exchanges associated with frontal meanders and mesoscale instabilities have been collected and analysed (Bower and Rossby, 1989).
Our test case focused on a specific day, 17 October 2007, when three large and well developed Gulf Stream meanders were observed between 65 • W and 50 • W. The core of the current is well identified by the comparison of the surface absolute dynamic topography (ADT) field with corresponding SST (sea surface temperature) and SSS (sea surface salinity) patterns, shown in Fig. 1.Upstream is a thin warm water tongue which develops into a steep meandering thermal and salinity front between 38 • N and 42 • N. The strongest gradients are observed at the first meander, around 62 • W/41 • N.
To summarize, after presenting the dataset used (Sects.2 and 3), this paper will focus on: -the strategies adopted to increase the effective resolution of the observation-based products (Sect.4), namely: -the integration of different high resolution Sea Surface Temperature level 4 products (SST L4, i.e. interpolated data) in the ARMOR3D processing (Sect.4.1); -the integration of the new high resolution Sea Surface Salinity (SSS L4) product (Buongiorno Nardelli, 2012) in the ARMOR3D processing (Sect.4.1); -the implementation of additional extrapolation methodologies to obtain high resolution 3-D re-analyses based on the high resolution Sea Surface Salinity product, on one selected high resolution SST L4 product and standard altimeter products (Sect.4.2); -the comparison of the various 3-D reconstructed fields, mainly focusing on the spatial scales that they are effectively able to resolve (Sect.4.3); -the diagnostic model used to retrieve the quasigeostrophic vertical velocity field from the improved observation-based density and geostrophic velocity fields and its validation by comparison with the My-Ocean Mercator Océan 1 / 4 • and 1 / 12 • resolution model vertical velocities (Sects.5.1 and 5.2); -the impact of the effective product resolution on the estimation of the vertical velocity field from the 3-D observation-based products (Sect.5.3) and an example of the dynamical interpretation of the vertical velocity field concentrating on a specific mesoscale feature (frontal meander).

Observations
In this study, two 3-D T /S retrieval techniques are considered: ARMOR3D (Guinehut at al., 2012) and the multivariate EOF reconstruction (mEOF-r) (Buongiorno Nardelli and Santoleri, 2005).Both methods require, on one hand, a historical in situ dataset to estimate the correlations between the parameters of interest (namely, temperature, salinity and steric height) or to identify their main statistical or empirical modes of variability, and, on the other hand, surface measurements of at least some of these parameters to be able to extrapolate their vertical profiles.In the following, the datasets used for this work are briefly described.

In situ data
The ARMOR3D systems (Sect.Nardelli, 2012).These observations are pre-processed according to Argo recommendations for data quality control (Wong et al., 2012).

SSS
The Sea Surface Salinity L4 data used as input to the 3-D reconstructions (section 4) has been developed in the framework of the MESCLA project (Buongiorno Nardelli, 2012).Its space-time resolution is 1 / 10 • , daily.The method used to retrieve this high resolution SSS field is based on an optimal interpolation (OI) algorithm that interpolates in situ salinity observations including satellite high-pass filtered SST in the determination of the weights used to interpolate SSS observations (using 1 / 10 • ODYSSEA SST L4 as input).SSS is represented as a function of space, time and SST (it is thus defined in a four-dimensional space), and a "generalized" distance is used to define a new covariance model that thus includes a thermal decorrelation term.In practice, this covariance function associates a higher weight to the SSS observations that lie on the isothermal of the interpolation point with respect to observations taken at the same temporal and spatial separation but characterized by different SST values.As satellite SST coverage and resolution is significantly higher with respect to in situ SSS observations, this method improves the interpolated field resolution not just in terms of grid spacing, but also in terms of the space-time features that are effectively retrieved.The covariance function parameters (i.e.spatial, temporal and thermal decorrelation scales) and the noise to signal ratio have been determined empirically, as fully described in Buongiorno Nardelli (2012).Hereafter, this SSS product will be called MESCLA HR SSS.

SLA/ADT
The altimeter Sea Level Anomalies (SLA) and Absolute Dynamic Topography (ADT) gridded data used for the 3-D retrieval are those produced and disseminated by the SSALTO/DUACS centre and represent the MyOcean Sea Level Thematic Assembly Centre intermediate product (AVISO, 2012).They are obtained as daily combined maps from all processed altimeters (Jason-1, Jason-2 and Envisat for the DT products used in our study) with a 1 / 3 • horizontal resolution.

SST
Different satellite SST datasets have been used for the different phases of the work, each characterized by different nominal and effective resolution (as summarized in Table 1 and discussed in Sect.4.3).In fact, as discussed by Reynolds and Chelton (2010), the true resolution of a L4 product is given by a combination of the grid spacing and of the analysis procedures and configurations applied (e.g.weighting functions and background fields).As a consequence, while combined ARMOR3D product is based on the Reynolds Optimally Interpolated L4 SST, at 1 / 4 • , corresponding to MyOcean V1 product (Larnicol et al., 2006;Guinehut et al., 2012), the higher resolution tests on ARMOR3D system have been performed on the ODYSSEA L4 data produced by Ifremer in the framework of the MERSEA project and maintained as part of MyOcean (Autret and Piollé, 2007), and on the Operational SST and Sea Ice Analysis system (OSTIA, see Donlon et al., 2011) also distributed as part of the MyOcean Sea Surface Temperature Thematic Assembly Centre.ODYSSEA provides daily SST estimates on a 1 / 10 • grid for the Global Ocean, based on both infrared and microwave measurements, while OSTIA L4 is available on a 1 / 20 • horizontal grid and includes also in situ SST measurements.

Model data
The model outputs used in this study are daily means computed from the global physical ocean forecasting system delivered as intermediate products by the global Monitoring and Forecasting Centre from MyOcean, namely Mercator Océan.In order to investigate the impact of the horizontal resolution in the vertical velocities reconstruction phase, the two components which compose the global system have been used: the global 1 / 4 • called PSY3V2R2 and the North Atlantic and Mediterranean Sea 1 / 12 • called PSY2V3R1 (Dombrowsky et al., 2009;Lellouche et al., 2012).Except for the horizontal resolution, the two configurations are really close in terms of ocean model version, numerical schemes, physical parameterizations, bathymetry, atmospheric forcing, assimilation scheme and assimilated data.The model configuration is based on NEMO1.09(Madec, 2008) with vertical z coordinates including partial step parameterization and 50 vertical levels from 1 m resolution at the surface to 400 m at the bottom.The main numerical schemes used in these configurations are a TVD (total variation diminishing) advection scheme and an isopycnal laplacian diffusion for the tracers, the energy and enstrophy conserving scheme and a biharmonic diffusion for the momentum.
The vertical mixing scheme is TKE (turbulence kinetic energy) with an enhanced convection parameterization in case of instability of the water column.All these options are classical and used in the global ocean model as mentioned in Barnier et al. (2006).The atmospheric forcing for the real time production is based on daily averages of the atmospheric variables or flux provided by the ECMWF real time forecasting system.The assimilation scheme (Tranchant et al., 2008) used in both configurations is based on the singular evolutive extended Kalman (SEEK) filter which allows assimilation of the sea level along-track satellite observations delivered by the MyOcean Sea Level Thematic Assembly Centre, the temperature and salinity profiles from the MyOcean In Situ Thematic Assembly Centre and the RTG (Real-Time-Global) sea surface temperature (http://polar.ncep.noaa.gov/sst/oper/Welcome.html).The model outputs used in this study are based on the "best analysis", which is performed every week with a one week delay in time to assimilate the most of observations over a one week assimilation window.This system was the global forecasting system operated during the V0 phase of MyOcean.The vertical velocity which is used as "reference" in this study to analyse the limits of the quasigeostrophic omega equation is computed by an upward integration of the horizontal divergence from the bottom (Madec, 2008) which is the standard way to compute the vertical velocities in the NEMO model in the case of a free surface condition.

3-D reconstruction
Several dynamic, variational, statistical and empirical techniques have been developed in the past to retrieve 3-D fields from a combination of in situ and satellite data (e.g.Carnes et al., 1994;Gavart and De Mey, 1997;Pascual and Gomis, 2003;Meinen and Watts, 2000;Watts et al., 2001;Mitchell et al., 2004).In fact, many of these methods are technically similar to some assimilation schemes (optimal interpolationlike), with the difference that the first guess used, i.e. the background analysis, is given by an average over the observations instead of a numerical model forecast.The error associated with this analysis thus represents the actual system variability.In essence, most statistical methods are based on the analysis of covariance relative to a set of in situ data profiles and on the identification of the principal modes characterizing the latter.However, the accuracy of each technique depends on the choice of the variables characterizing the state of the system, as well as on the number of degrees of freedom absorbed by each method (e.g.Buongiorno Nardelli and Santoleri, 2004).
Univariate techniques such as single EOF (empirical orthogonal function) reconstruction analyze the principal components of each parameter along the water column and hypothesize a relationship between the amplitude of such components and a (not necessarily linear) combination of surface parameters (Carnes et al., 1994).Simpler methods, as the one used within ARMOR3D, assume a direct correlation between surface and deep values (Guinehut et al., 2012).The new methods considered within MESCLA are based on multivariate approaches (multivariate EOF reconstruction, mEOF-r).These methods analyze the steric height, temperature and/or salinity covariance and reconstruct the vertical profiles via a combination of a limited number of modes.Following an idea first proposed by Pascual and Gomis (2003), they include an approximation of the geopotential stream function (the steric height profile) in the status vector, thus more directly correlating physical-chemical parameter variability to dynamics.The application of these methods already yielded promising results, also compared to empirical methods such as the computation of the gravest empirical modes (Buongiorno Nardelli and Santoleri, 2005).
A double approach has thus been followed to improve the resolution of the observation-based 3-D fields.The two approaches involve different levels of complexity and might be considered as potential successive steps in a gradual improvement of operational products.As a first step, the algorithm used to obtain the ARMOR3D product has been adapted to integrate higher resolution SST and SSS data.Then mEOF-r technique has been adapted and tested on a subset of input data (considering only the highest resolution SST data).The MyOcean combined ARMOR3D product is computed every week (Wednesday fields) on a 1 / 3 • Mercator horizontal grid, which corresponds to the altimeter SLA grid and from the surface down to 1500 m depth on 24 vertical levels.ARMOR3D method, thoroughly described in Guinehut et al. (2012), improves a climatological first guess using two main steps.At first, synthetic temperature (T ) profiles are estimated by extrapolating altimeter and SST data through a multiple linear regression method and covariances computed from historical data.For synthetic salinity (S) profiles, the method uses only altimeter data.The multiple/simple linear regression methods are expressed as: and where SLA and SST denote anomalies from the ARIVO monthly climatology (Gaillard and Charraudeau, 2008), T clim and S clim denote the ARIVO monthly fields, and α, β and γ are the regression coefficients of the SLA and SST onto temperature and of SLA onto salinity, respectively.They vary with depth, time and geographical location and are expressed as covariances between the variables (only the z variable is indicated here for clarity): and Successively, the synthetic profiles (hereafter referred to as synthetic ARMOR3D fields) are combined with in situ temperature and salinity profiles using an optimal interpolation method (Bretherton et al., 1976) to create the combined AR-MOR3D product.The current paper focuses on the synthetic ARMOR3D fields.As a preliminary step, ARMOR3D performs some crucial processing of altimeter data, being able to extract the steric contribution to the sea level variations consistent with the first 1500 m depth (filtering out the eustatic component and the deep steric contribution).This pre-processing is based on regression coefficients deduced from an altimeter/in situ comparison study (Guinehut et al., 2006;Dhomps et al., 2011).
In the present work, the three SST products described in Sect.2.2.3 have been used to test the impact of SST resolution on the synthetic T field estimation (step one of the method).Additionally, the use of MESCLA HR SSS fields has also been tested for the reconstruction of the synthetic salinity.While the synthetic ARMOR3D salinity fields is obtained with a simple linear regression to altimeter SLA, the method has been thus modified to a multiple linear regression method (as for temperature) to include also the information from SSS.The salinity field is now calculated from the following equation: with again the regression coefficients expressed as: and All tests required us to interpolate the altimeter SLA onto each SST grid (i.e.interpolating the original data at 1 / 4 • , 1 / 10 • and 1 / 20 • , respectively).Actually, this interpolation has to be performed with particular care in order not to introduce a spurious signal.After having tested different methods (simple bilinear interpolation, Akima spline), a classical spline method has been chosen.

Method
The multivariate EOF reconstruction (mEOF-r) technique is based on the analysis salinity, temperature and steric height profiles through a multivariate EOF decomposition and on the availability of corresponding surface values (Buongiorno Nardelli and Santoleri, 2005).
Here, we will briefly recall how mEOF-r works.A single 3 m × n multivariate observation matrix X is obtained from the three original sets of data, each of m × n dimensions, where n is the number of measurements (stations) and m the number of vertical levels.Data are preliminarily normalized dividing each parameter by its standard deviation (computed for the whole profile).Mean profiles estimated from the whole training dataset are removed in order to obtain anomalies and estimate the covariances.The columns of this matrix consist of the three normalized profiles of temperature (T ), salinity (S) and steric height (SH) anomalies, each taken at the same location.
To compute the multivariate EOF, the singular value decomposition of this new matrix of data is performed.In that way, "multi-coupled" modes are identified, each containing the three patterns corresponding to the parameters considered.T (z, r), S(z, r) and SH(z, r) can thus be expanded in terms of these three series of patterns.The same coefficient/amplitude (a k ) is found for all parameters, L k , M k and N k being the modes: If these expansions are limited to the first three modes, the vertical profiles can be estimated from the surface values (z = 0) of the three parameters solving the system for a 1 , a 2 and a 3 and substituting them in the truncated expansions: Of course, it is also possible to truncate the expansions to the second mode (or to the first one), which actually means that only two (one) surface parameters are sufficient to retrieve the whole profiles.Similarly, the whole analysis can be performed directly on two sets of parameters at a time.
The mEOF-r method requires a training dataset of T , S and SH profiles to extract the main vertical modes of (co)variability.This training dataset might be selected differently at each grid point on the basis of different criteria: fixing a space and/or time search radius (e.g.1000 km, week, month, or year), or keeping only the nearest n profiles, etc. Depending on this choice, one may end up with different reconstruction models (i.e.different mEOFs) for each grid point or with a single set of modes.After some preliminary hindcast tests (not shown), it was decided to select all the profiles collected in the domain within a monthly window.
Given the sparse, even though regular, distribution of data in the training set, this was found to be the simplest but also most reliable way to estimate EOFs.
Similarly to the pre-processing performed by ARMOR3D, in order to retrieve the 3-D vertical fields from surface data, a preliminary step is to estimate/extract the surface steric heights from satellite altimeter data.Actually, there is no way to evaluate the deeper baroclinic and the barotropic contributions from altimeter data and surface measurements alone.As a simple approximation, this estimation is reduced here to an adjustment of the ADT to minimize the differences between the steric height computed from simultaneous (or quasi-simultaneous) in situ profiles and co-located ADT estimates (through a simple regression).In contrast to Guinehut et al. ( 2012), who compute spatially varying climatological regression coefficients between co-located dynamic height anomaly computed from Argo T /S profiles and sea level anomaly (SLA) altimeter data, this adjustment has been performed here considering weekly matchups and a single regression.The adjusted ADT will be called in the following surface steric height (SSH).The in situ T /S profiles described in Sect.2.1 were re-interpolated at 10 dbar resolution down to 1000 dbar, and steric height profiles were obtained taking 1000 dbar as reference pressure level.These data have been used as training dataset, while the same altimeter ADT/SLA data, ODYSSEA SST L4 and MESCLA HR SSS L4 data used for the tests on ARMOR3D (Sect.2.2) were used as surface input for the mEOF-r technique.

Hindcast evaluation of the mEOF-r performance in the Gulf Stream area
The accuracy of the mEOF-r technique depends on the characteristics of the system under study.In fact, statistical modes will generally reflect the variability associated with different physical processes, depending on the area under study, and the percentage of variance associated with each mode will change accordingly.As this technique was never applied before on the Gulf Stream area, different mEOF-r configurations (namely varying the number of parameters and modes considered) have been compared.A first estimate of their accuracy has been obtained through a hindcast validation.This means that the surface values of the in situ profiles used as training datasets were taken as input data for the reconstruction.The hindcast errors were thus estimated as the mean and standard deviation of the differences between the vertically reconstructed (synthetic) profiles and the original measurements.The advantage of this kind of validation is given by the large number of profiles available (more than 80), while the disadvantage clearly resides in the fact that hindcast is not an independent validation, as the same data are used to train and test the method.

B. Buongiorno Nardelli et al.: Towards high resolution mapping of 3-D mesoscale dynamics
The hindcast validation was applied to the mEOF-r configurations listed below: 1. mEOF-r(T -S-SH): The mEOFs are computed from T , S and SH profiles, and corresponding synthetic profiles are obtained using SST, SSS and SSH as input data (the amplitude of the first 3 modes is retrieved).
2. mEOF-r(T -S-SH) SST-SSH : The mEOF are computed from T , S and SH profiles, and corresponding synthetic profiles are obtained using only SST and SSH as input data (the amplitude of the first 2 modes is retrieved).

mEOF-r(T -S-SH) SSS-SSH :
The mEOF are computed from T , S and SH profiles, and corresponding synthetic profiles are obtained using only SSS and SSH as input data (the amplitude of the first 2 modes is retrieved).

mEOF-r(T -SH):
The mEOF are computed from T and SH profiles only, and corresponding synthetic profiles are obtained using SST and SSH as input data (the amplitude of the first 2 modes is retrieved).

mEOF-r(S-SH):
The mEOF are computed from S and SH profiles only, and corresponding synthetic profiles are obtained using SSS and SSH as input data (the amplitude of the first 2 modes is retrieved).
Mean bias error (MBE) and standard deviation error (STDE) profiles for both temperature and salinity are shown in Fig. 2a and b, respectively.It is interesting to observe that the synthetic mEOF-r MBE is generally quite small.The mEOF-r provides the smallest STDE errors when only two modes are considered, both in the trivariate formulation and in the bivariate one.A simple explanation for this may be found looking at the mEOF modes (Fig. 3) and corresponding explained covariance percentage, and by comparing them with the dynamical modes that can be inferred from the mean stratification profile N 2 (Brunt-Väisälä frequency), as computed from the training dataset.In fact, the first mEOF mode explains an extremely high percentage of the variance (almost 99 %).If we look at the corresponding SH profile, it displays a quite smooth shape, increasing from zero at the reference level to a maximum at surface.This first mode closely reminds the first baroclinic mode that can be estimated from the linearized quasi-geostrophic vorticity equation, (e.g.Cushman-Roisin, 1994): This equation, written here for the stream function , can be solved searching a solution of the kind: kl+ly+ωt) .
The corresponding eigenvalue problem for the vertical component (z) can be easily integrated numerically.In our estimates, free surface boundary condition and a null stream function at the bottom were assumed.Due to this zero stream function assumption, the solution only provides baroclinic modes (Fig. 4).This analysis thus confirms that instabilities associated with the first baroclinic mode can be considered the predominant source of variability in the study area.
Conversely, only about 0.3 % of the variance is explained by the second mEOF mode.However, some important information is still contained in the second mode.In fact, we also ran a single mode mEOF(T -S-SH) reconstruction (both for temperature and salinity) which gave much worse results than the mEOF-r in the two mode configurations (see Fig. 5).Actually, T and S patterns in the first mEOF mode have the same sign, meaning that the surface anomalies with respect to the mean profile driven by this mode reflect down to the deep layers.On the contrary, the second mode basically accounts for the presence of T and S anomalies only in the upper layers, which might be related to the presence/absence of waters of coastal/riverine origin.More investigations will be needed to better understand which process drives the variability of this second mode.Meanwhile, it is clear that the third mode is related to conditions that apply to an extremely low number of profiles.Adding it to the reconstruction, in fact, leads to the typical errors associated with model over-fitting.Only the best performing techniques, as evaluated with the previous hindcast validation, have been applied to the SST, SSS and adjusted ADT maps for our study in the Gulf Stream area.
Considering the results of the tests presented in this section, the selected techniques are the mEOF-r(T -S-SH) SST-SSH and mEOF-r(T -S-SH) SSS-SSH .

Synthetic ARMOR3D and mEOF-r comparison
The new T and S synthetic fields described in Sects.4.1 and 4.2 have been compared: four are based on the ARMOR3D algorithm (three using only Reynolds L4, OSTIA L4 and ODYSSEA L4 SST as input, and one using ODYSSEA L4 SST and MESCLA HR SSS), and one is based on mEOF-R (using ODYSSEA L4 SST and MESCLA HR SSS as input).The choice of using only ODYSSEA L4 in combination with MESCLA HR SSS for the highest resolution tests was driven, on one hand, by the fact that ODYSSEA is among the highest resolution Global SST L4 available and quality controlled operationally (see also Dash et al., 2010;Maturi et al., 2010) and, on the other hand, by the fact that it is the same product used to retrieve the high resolution SSS by Buongiorno Nardelli (2012).
The comparison was performed considering the temperature and salinity patterns at the surface, namely the input fields, and at 100 m depth, which is a particularly interesting level as it corresponds approximately to the base of the euphotic layer in the Gulf Stream area (Oschlies and Garc ¸on, 1998).
The analysis has been performed both qualitatively, i.e. looking at spatial gradients and at the way the different features are resolved/projected at depth, and quantitatively, by estimating and comparing the zonal power spectra, defined similarly to Reynolds and Chelton (2010), for each of the products.To compute the spatial spectra at each latitude reducing the spectral leaking, a Blackman-Harris windowing has been applied prior before computing the fast Fourier transform.Moreover, as the windowing attenuates the overall signal power, the spectra have been scaled with corresponding coherent power gains to make them comparable (Harris, 1978).The spectra obtained for each latitudinal band have been finally averaged in a single spectrum for each product.
Starting from the analysis of the temperature fields at the surface, Ostia SST has the lowest energy levels at all spatial frequencies (namely the lowest variance), and displays a 'red spectrum' behaviour until reaching an almost flat response at approximately 1 / 4 • resolution, though its nominal resolution is the highest among the products considered (1 / 20 • ) (Fig. 8).This is coherent with the smooth appearance of the SST field as evidenced in Fig. 7. Conversely, Odyssea SST and Reynolds SST display much higher power (more than one order of magnitude) at frequencies higher than 10 −1 (deg −1 ) and up to their Nyquist frequency (i.e. at characteristic lengths between approximately 1000 km and 25 km/50 km, respectively), even if both of them also reach an almost flat plateau at scales smaller than 1 / 4 • (though if keeping a higher energy content than Ostia).Odyssea SST very clearly displays the highest energy in the spatial range between 100 km and 25 km, and the smallest spatial features: the Gulf Stream core is visible as a thin warm tongue offshore Cape Hatteras and several small scale features are observed along the main flow.Even if some of these might be considered quite "noisy", most of them are clearly related to mesoscale features (e.g. the small cold core cyclonic eddy centered at 40 The deep (100 m) synthetic temperature fields generally reflect the spatial characteristics of the surface input data, but the ARMOR3D and mEOF-r also have a different impact on the way the surface scales are projected in the deep fields.In particular, ARMOR3D synthetic field displays smoother structures than those visible in the corresponding mEOF-r field, even when using Odyssea SST as input.The corresponding power spectrum gets closer to the Ostia/Reynolds ARMOR3D spectrum, while mEOF-r spectrum at 100 m displays basically the same features as the corresponding SST, except for a slightly higher energy at the scales smaller than 1 / 4 • .
Concerning the salinity field, it has to be stressed that MESCLA HR-SSS and the ARMOR3D SSS fields are very different, as the latter is computed by combining a climatological first guess and the covariance between altimeter SLA and salinity (Fig. 9).More in detail, the MESCLA HR-SSS field shows a very sharp gradient and several small scale meanders all along the front of the Gulf Stream, with much fresher SSS than those present in the ARMOR3D field, which appears quite uniform and smooth.Corresponding SSS spectra display very different energy content below approximately 300 km, down to the Nyquist spatial frequency.Both the mEOF-r and the ARMOR3D salinity reconstruction methods propagate the sharp gradient and the small scale features retrieved by MESCLA HR SSS down to 100 m, while the standard ARMOR3D field displays a much smoother pattern.Similar to what was found for the synthetic temperature fields, ARMOR3D vertical projection, however, significantly reduces the energy in the spatial frequency range below 10 −1 (deg −1 ) (Fig. 10).Again, mEOF-r keeps basically the same energy down to 100 m, even with a slight increase at spatial scales below 10 −1 (deg −1 ), and a slight decrease at the longer wavelengths.The higher energy level at scales below 1 / 3 • (where a plateau is again present) might anyway be considered as an indication of a too noisy reconstruction.
A preliminary evaluation of the vertical accuracy of the synthetic mEOF-r and ARMOR3D products has also been roughly performed via a weekly matchup comparison, which means that the in situ profiles collected in a temporal range of ±3 days (weekly matchups) have been taken aside as an independent test dataset, and have been compared with the co-located profiles obtained with the different reconstruction methods.The matchup comparison has the advantage of starting from fully independent surface input data, but the number of in situ profiles available within this weekly window was too low to provide a significant estimate (actually, only 13 matchups were found).In fact, all error profiles estimated from the high resolution surface input produced equivalent results (within a 1 sigma confidence level, estimated through bootstrapping).The estimated MBE and STDE are shown in Fig. 6.However, a longer test period should be used to get a real validation, and the estimated difference can possibly also be affected by the temporal variability at scales shorter than 3 days, which might not be negligible in rapidly evolving frontal areas.

QG Vertical velocity estimation
Relatively intense vertical exchanges in the oceans are associated with the mesoscale structures.Nevertheless, vertical velocities are generally lower by a factor of up to 10 4  than horizontal ones, and consequently they are not easily measured through direct observations (e.g.Klein and Lapeyre, 2009;Frajka-Williams et al., 2011).Various indirect methodologies have thus been proposed to estimate vertical velocity from observed density and geostrophic ve- locity fields.Though more complicated techniques such as the semi-geostrophic omega equation (Viúdez and Dritchel, 2004) have been proposed, the most used technique is based on the solution of the quasi-geostrophic (QG) omega equation (e.g.Tintoré et al., 1991;Buongiorno Nardelli et al., 2001;Pascual et al., 2004;Ruiz et al., 2009), which has already been shown to give reasonable estimates of the vertical velocities compared to primitive equation models (Pinot et al., 1996).

Q vector formulation of the omega equation
In this work, the algorithm for the solution of the Q vector formulation of the omega equation (as applied by Pascual et al., 2004;Ruiz et al., 2009) was adapted to the specific products considered.
Actually, the omega equation requires as input both the geostrophic field and the density stratification.The geostrophic currents have been estimated referencing the thermal wind estimates to the absolute surface altimeter velocities (when applied to observation-based products, see also Mulet et al., 2012).Two reference levels for dynamic height computation are considered for models: the surface and 1000 m depth.The code is derived from the QG vorticity and thermodynamic equation (Hoskins et al., 1978): where (U, V ) are the geostrophic velocity components, N is the Brunt-Väisälä frequency and f the Coriolis parameter.In this implementation, N only depends on depth.Different boundary conditions have been tested (i.e.Dirichelet and Neumann conditions), however, given the elliptic nature of the omega equation, no significant differences were found a few grid points away from the boundaries in the two cases.

Comparison of the MyOcean model PE and QG vertical velocities
To evaluate the accuracy of the quasi-geostrophic approximation, the QG vertical velocities The omega estimates show a high sensitivity, both in terms of shape and intensity, to the spatial resolution (Fig. 11).Vertical velocities obtained from the 1 / 12 • (PSY2w) model (Fig. 11c and d) are a factor of 2-3 larger than the 1 / 4 • (PSY3w) (Fig. 11a and b) version (maximum upward and downward velocities of the order of 40-60 m day −1 vs. 20-30 m day −1 , respectively).On the other hand, the comparison between PE (Pearson coefficient) and QG vertical velocity patterns shows a reasonable agreement in both model simulations, even if model QG vertical velocity values underestimate the PE velocities at low resolution.These differences can be better quantified looking at the scatter plot and computing correlation coefficients, as displayed in Fig. 12 (similar results were also found for deeper layers, not shown).
In fact, QG velocities displayed maximum upward and downward values of the order of 40-60 m day −1 and 10-15 m day −1 in the high resolution and low resolution tests, with correlation with PE reaching almost 0.8 and 0.7, respectively.The test of re-interpolating the PSY2 data (using bilinear interpolation) on the 1 / 4 • grid showed that both w fields (QG and PE) are dramatically reduced by the interpolation, so that the final results are similar to what was obtained with PSY3 (Fig. 12).
The results thus indicate that resolution is a key factor in the estimate of the vertical component of the ageostrophic velocity through the quasi-geostrophic omega (compared to PE estimates).If sufficient resolution is kept, velocities are more correctly reproduced even in the quasi-geostrophic approximation, and even though the patterns in the QG solution may be slightly smoother, the estimated values compare reasonably well.In the following, the QG method will thus be applied to the observation-based synthetic 3-D products, providing a fully observational estimate of the vertical velocities for the test case.

Vertical velocities estimates from observation-based 3-D products
As expected, the quasi-geostrophic vertical velocities obtained from the synthetic fields increase significantly when increasing the grid resolution.A qualitative comparison of the qgw retrieved by applying the omega equation to the synthetic ARMOR3D products and mEOF-r fields is thus illustrated in Fig. 13.The lowest velocities were estimated from the synthetic ARMOR3D using the Reynolds L4 SST (peak qgw ∼ 32 m day −1 , horizontal velocities peak ∼ 1.5 m s −1 ).
On the other hand, the synthetic ARMOR3D qgw computed from OSTIA, namely at the highest nominal resolution, do not display the highest average and peak velocities.This is not surprising considering the effective resolution of the tracer fields is lower (see Sect. 4.3).However, coherent with the findings of the Sect.5.2, they display higher values than those obtained from Reynolds (peak qgw ∼ 52 m day −1 , peak horizontal velocities ∼ 1.7 m s −1 ).Slightly higher values were computed from synthetic AR-MOR3D using Odyssea SST as input, with peak qgw values of ∼ 54 m day −1 and peak horizontal velocities of ∼ 1.7 m s −1 .The most intense velocities, however, were estimated from the two products using both Odyssea SST and MESCLA SSS as input, with qgw peak values of 58 m day −1 and 66 m day −1 in the synthetic ARMOR3D and mEOFr, respectively.More pronounced differences in both the geostrophic velocities and vertical velocity structures were found along the main Gulf Stream jet, where the mEOFr field displays higher values than synthetic ARMOR3D.Though a full dynamical description of the Gulf Stream dynamics is clearly beyond the scope of the present paper, it can be observed that the alternating upwelling/downwelling patterns found in our calculations along its main flow are coherent with the evolution of the Gulf Stream meanders as described in literature.For example, the zoom in Fig. 14 shows an alternating pattern of troughs and crests along one of the main meanders.Lindstrom et al. (1997) found upwelling and downwelling regions (at the main thermocline level) characterized by values of up to 2 mm s −1 (namely more than 150 m day −1 ) at horizontal scales of about 100 km.These scales were found to be comparable both along and across stream.As for our estimates, downwelling/upwelling generally occurs entering/exiting meander troughs in Lindstrom's et al. (1997) observations, even though the intensity of the vertical circulation is clearly dependent on the phase of the meander evolution/propagation.Similar observations were also described by Bower (1989) through the direct analysis of RAFOS floats' trajectories.It might be worth noting that in our calculations, independent of the input data and  technique used, all the synthetic QG estimates gave the same upwelling/downwelling patterns, even if the strongest velocities were retrieved from mEOF-r fields.

Conclusions
Within the MyOcean R&D project MESCLA, a step towards a more efficient combination and more complex analysis of existing observations has been made.MESCLA tested innovative methods for the high resolution mapping of 3-D mesoscale dynamics from a combination of in situ and satellite data (as described in Sect.4), developing new products that might be used as prototypes to gradually build the next generation of operational observation-based products.In order to demonstrate the new techniques' potentials, different estimates of the vertical velocities derived from different 3-D synthetic fields through a quasi-geostrophic diagnostic model have been compared (Sect.5).Resolution confirmed to be an important factor for the retrieval of the currents.However, even within the limits of a simplified dynamical framework, and knowing that most of the analysis could not necessarily go beyond a simple qualitative comparison (vertical velocities cannot be measured directly at sea), realistic estimates of the vertical field could be retrieved, at least as compared to those diagnosed through primitive equation numerical models.The ocean observation-based products tested within MESCLA might thus open a wide range of possible applications for both operational oceanography and ocean climate monitoring studies.

Fig. 6 .
Fig. 6.Weekly comparison of temperature and salinity profiles from mEOF-r and and synthetic ARMOR3D using Odyssea SST (a) and MESCLA SSS (b).

Fig. 8 .
Fig. 8. Zonal wavenumber spectra computed from the reconstructed temperature fields at the surface (a) and at 100 m depth (b).

Fig. 10 .
Fig. 10.Zonal wavenumber spectra computed from the reconstructed salinity fields at the surface (a) and at 100 m depth (b).

Fig. 13 .
Fig. 13.QG vertical velocity fields at 100 m as retrieved from the temperature and salinity fields reconstructed through the selected methods: (a) ARMOR3D using Reynolds L4 SST as input, (b) ARMOR3D using OSTIA L4 SST as input, (c) AR-MOR3D using ODYSSEA L4 SST as input, (d) ARMOR3D using ODYSSEA L4 SST and MESCLA SSS as input, and (e) mEOF-r(T -S-SH) SSS-SSH for the salinity field and mEOF-r(T -S-SH) SST-SSH for the temperature.

Fig. 14 .
Fig. 14.Zoom of the QG vertical velocity fields at 100 m as retrieved from the temperature and salinity fields reconstructed through the selected methods: (a) ARMOR3D using Reynolds L4 SST as input, (b) ARMOR3D using OSTIA L4 SST as input, (c) ARMOR3D using ODYSSEA L4 SST as input, (d) ARMOR3D using ODYSSEA L4 SST and MESCLA SSS as input, and (e) mEOF-r(T -S-SH) SSS-SSH for the salinity field and mEOF-r(T -S-SH) SST-SSH for the temperature.SSH isolines are superimposed to allow the identification of meanders' crests and troughs.

Table 1 .
Description of the L4 SST products used.