*3.3. The Uncertainties on the Geostrophic Velocities (σu, σv)*

Another input parameter for the reconstruction method is the error on both the components of the surface geostrophic currents, previously defined as *σ*<sup>u</sup> and *σ*<sup>v</sup> for the zonal and meridional components of the current, respectively. This was estimated in a similar manner as for the forcing term, via comparison with the in-situ derived surface currents in the Mediterranean basin (data doi = 10.6092/7a8499bc-c5ee-472c-b8b5-03523d1e73e9). Using 24 years of altimeter-derived GV (introduced in Section 2) we could colocate the GV along the buoys trajectories and generate two-dimensional maps of the RMSE between the satellite and the in-situ measured currents. This was achieved selecting all the buoys evolving in presence of their drogue and binning the RMSE values in 1◦ × 1◦ boxes, as illustrated in Figure 3.

The basin-scale mean error on the geostrophic velocities is of the order of 10 cm·s−<sup>1</sup> for both the components of the motion and it is slightly larger for the meridional component than for the zonal component of the motion (10.5 cm·s−<sup>1</sup> versus 9.6 cm·s−1, respectively). The maximum RMSE values can exceed 20 cm·s−<sup>1</sup> in correspondence of the Alboran gyres for both the surface flow components and their spatial variability is also correlated with the mean positions of the major currents in the basin. Indeed, the largest RMSE values (≥15 cm·s−1) are found in correspondence of the Algerian, Liguro-Provenzal, Libyo-Egyptian and Cilician currents as well as the Atlantic Ionian Stream [3].

**Figure 3.** RMSE between the altimeter and the drifter-derived surface currents, binned on 1◦ × 1◦ boxes.

#### **4. Results**

Relying on the CMEMS altimeter-derived GV, SST and on the error fields described in the previous section, we implemented the RS18 method to generate five years of Optimal currents in the Mediterranean Sea (from 2012 to 2016), a time window compatible with all the datasets used in the study. Such currents have spatial temporal resolution of 1/24◦ and 1 day, respectively. In the following sections we discuss the Optimal currents performances with respect to the geostrophic estimates, mainly via validation against in-situ or model-derived data.

#### *4.1. Case Study 1: Vortex Dynamics*

Here we give a first example of the Optimal Currents capabilities in presence of an oceanic vortex (generally referred to as eddy in the oceanographic community). An anticyclonic eddy located halfway between the Balearic Islands and Sardinia is detected on 25 April 2016 and is recognisable by a ring-shaped SST anomaly centered at approximately 39◦45 N and 5◦50 E and is associated to a clockwise circulation (Figure 4).

**Figure 4.** Top: Comparison between the Optimal currents (red arrows) and the GV (black arrows) in presence of anticyclonic circulation in the Western Mediterranean (25th of April 2016). The surface currents fields are superimposed to the SST (in colours, ◦C); the white ellipse is referenced in Section 4.1. Bottom left: Soil moisture and ocean salinity (SMOS)-derived sea-surface salinity (SSS) salinity anomaly from Barcelona Expert Center (25th of April 2016). The black closed contour is the 15.9 ◦C level indicating the position of the vortex. Bottom right: Same as bottom left, using the LOCEAN L3 SMOS SSS data, on the 26th of April 2016 .

Cold SST anomalies in presence of surface-intensified anticyclones have already been reported in literature [36,37]. In the Mediterranean sea, they can result from SSS compensation mechanisms, as also described in [16]. Using the Barcelona Expert Center SMOS L4 SSS for the Mediterranean Sea as well as the LOCEAN L3 SSS estimates (the latter being available on 26 April 2016), we found that the eddy lies in a negative salinity anomaly area. This is confirmed by the lower panels of Figure 4. Here, the eddy is highlighted by the 15.9 ◦C closed contour and the SSS anomaly is evaluated as the SSS in the core of the eddy minus the SSS at the eddy periphery.

Thus, the salinity compensation mechanism is probably the cause of the observed oceanographic structure. However, the detailed investigation of the eddy formation and characteristics is outside the scope of this study. Here, we want to assess the capability of using the SST information to improve the altimeter-derived currents. In the western boundary of the vortex, the GV do not resolve such a dynamical feature, indicating a circulation which crosses the north western section of the eddy with a northeastward flow direction. This effect is a result of the mapping procedure statistics in order to obtain two dimensional maps from along-track altimeter data. On the other hand, if we consider the dynamical evolution of the SST, i.e., we account for the SST spatial-temporal gradients, the Optimal reconstruction method prescribes a correction factor that deviates the flow from the altimeter estimate, resulting in a more accurate description of the eddy circulation. This is particularly evident in the area highlighted by the white ellipse in the left panel of Figure 4.

This case study shows one of the main advantages of RS18, i.e., merging a background currents field with a higher resolution SST map. Other approaches, as SQG, rely on SST images to infer the surface flow. Nevertheless, when salinity anomalies contribute significantly to the sea surface density gradients, such approach may well lead to a current estimate which is in an opposite direction with respect to the actual one, as also found by Isern-Fontanet et al. 2016 [16]. Instead, using the RS18 method this ambiguity is overcome.

#### *4.2. Case Study 2: Coastal Upwelling in the Sicily Channel*

The Sicily Channel area, whose approximate coordinates are 11◦E to 13◦E and 36◦N to 38◦N constitutes an interesting testbed for the performances of the Optimal Currents. In this region, in proximity of the southern coast of Sicily, coastal upwelling is a very frequent phenomenon and is known to generate cold SST patches associated with an offshore directed current [38,39]. In particular Piccioni et al. 1988 [39], using a statistical approach, found that such events are mostly detectable during summer and that they mostly develop with a delay of three days with respect to the southeasterly wind input (i.e., with a non-zero alongshore component with respect to the southern coast of Sicily).

We show an example of coastal upwelling detected on 23 July 2016. The result of its dynamics is clearly visible from the SST spatial distribution and is consistent with the CERSAT/IFREMER averaged surface winds during the three previous days (20–22 July 2016) (Figure 5).

**Figure 5.** Left panel: Geostrophic Velocities (black arrows, 1/8◦ horizontal resolution) and Optimal Velocities (red arrows, 1/24◦ horizontal resolution) over sea-surface temperature (SST) in ◦C—23 July 2016 (the white boxes are referenced in Section 4.2). Right panel: CERSAT/IFREMER 6 hourly averaged surface winds—20–22 July 2016).

The SST anomaly between the interior of the upwelled waters and the surrounding areas exceeded 4 ◦C and the shape of the patch (in agreement with the mechanism described in [39]) indicates an offshore current, perpendicular to the southern coast of Sicily (between 12◦20'E and 13◦E). The offshore circulation was not entirely captured by the GV (Figure 5). Indeed, at the positions 12◦40 E–37◦30 N and 12◦20 E–37◦15 N, highlighted by the white boxes, the GV exhibit a cross thermal gradient component which is not consistent with the wind-driven circulation. A cross-thermal gradient geostrophic circulation is also found in correspondence of the lower tip of the cold SST patch. This effect is reduced when the Optimal reconstruction is implemented. The qualitative validation between the GV and the Optimal currents patterns allows to conclude that, enriching the GV with the information contained in the spatial and temporal variability of SST, it is possible to retrieve the ageostrophic circulation (as it is the case for the wind-driven upwelling dynamics) starting from the geostrophic, altimeter-derived currents. In order to further confirm this, we computed the divergence of the surface currents field. The divergence is defined as *∂*xu + *∂*yv (with u and v the generic zonal and meridional components of the motion, respectively). When applied to a geostrophic field, the divergence operator yields zero in the whole surface currents domain, while it can be O(10−<sup>5</sup> <sup>s</sup>−1) for a total surface current field, i.e., considering both the geostrophic and the ageostrophic components of the motion [40]. Indeed, choosing the same region of Figure 5, we compared the divergences obtained using the Optimal currents, the model-derived surface currents of the MFS operational model and the MERCATOR global operational model (Figure 6). In addition, a comparative analysis with the monthly outputs of the CNR-Sea Forecast operational model will be given. Such divergences were all computed using a 9-point stencil width technique, as in [41].

The GV, though not exhibiting divergence values strictly equal to zero (mostly due to the numerical scheme of the computation), were mostly O(10−<sup>8</sup> <sup>s</sup>−1) everywhere in the analyzed area, (as well as in the rest of the Mediterranean Basin (not shown). Thus, we assumed this order of magnitude to be the lower limit for the surface currents field divergence, that can be referred to as geostrophy. On the other hand, the upper limit for the divergence values is identified with the results given by the MFS operational model. In this case, the divergence could reach the value of ±<sup>1</sup> × <sup>10</sup>−5s−<sup>1</sup> (bottom panel of Figure 6). Quite interestingly, the Optimal currents exhibit divergence values at least two orders of magnitude larger than the GV, up to ±<sup>8</sup> × <sup>10</sup>−6s−1, hence comparable with the ones of a total surface current. The current field of the MERCATOR operational model also exhibits values O(10−<sup>6</sup> <sup>s</sup>−1), though its maxima only reached the 70% of the values observed in both the Optimal and the MFS-model case. Moreover, using the monthly averages of the submesoscale-permitting (resolution 2 km) operational model for the Sicily Channel and Tyrrhenian Sea (dataset labeled with "7" in Section 2), we had further confirmations that the expected value of the surface currents divergence during July 2016 must be O(10−<sup>5</sup> <sup>s</sup><sup>−</sup>1), as also shown in Figure 7.

In addition, if we compare the spatial distribution of the divergences (DIV hereinafter) of the Optimal, MFS and MERCATOR surface currents (for 23 July 2016), the following similarities can be found:


The most interesting result of this comparison is that the combination of the geostrophic current and the SST data is a promising technique for the retrieval of the total surface current field from an initial geostrophic estimate, i.e., the one provided by satellite altimetry. Indeed, the Optimal current divergence is of the same order of magnitude predicted by high-resolution numerical models and observations, also sharing common features in the spatial patterns given by the MERCATOR and MFS-model estimates (though in the numerical model outputs the structures are slightly smoothed compared to the Optimal case).

**Figure 6.** Divergence of the surface currents field in the Sicily Channel (23 July 2016). Top panel: Optimal currents; middle panel: MFS operational model for the Sicily Channel; bottom panel: MERCATOR operational model. Black contours indicate sea surface temperature, in ◦C. The light green and dark green circles are referenced in Section 4.2.

**Figure 7.** Divergence of the surface currents field in the Sicily Channel (July 2016 monthly average). Prediction of the CNR-Sea Forecast submesoscale-permitting operational model for the Sicily Channel and the Tyrrhenian Sea.

In the following sections we will describe the quantitative validation of the Optimal currents via two approaches: a first one in which the comparison will be performed with respect to the in-situ drifting buoys measurements in the Mediterranean Sea and a second in which the Optimal currents will be compared to the estimate of the High-Frequency Radar (HFR hereinafter) surface currents in the Malta–Sicily channel.

#### *4.3. Validation with the In-Situ Measured Currents*

In the period 2012 to 2016, the in-situ sea surface currents have been deduced from the trajectories of 104 15m-drogued drifting buoys at the positions indicated in Figure 8 (data doi = 10.6092/7a8499bc-c5ee-472c-b8b5-03523d1e73e9, [42]).

**Figure 8.** Positions of the sea surface currents measurements given by the surface drifting buoys in the Mediterranean Area during the period 2012–2016 [42].

The buoys have been mostly occupying the Western, the Central Mediterranean Sea, the Adriatic Sea and the area around the Island of Cyprus (eastern Mediterranean), measuring the sea surface currents with 6-hourly temporal resolution. Selecting the totality of the buoys trajectories, we compared the performances of the GV and the Optimal currents against the in-situ-derived measurements. This was achieved creating two match-up databases for the GV and the Optimal currents, which are both gridded fields. The match-up database was generated using cubic interpolation for the spatial colocalization of the gridded currents over the buoys positions while a linear interpolation was used to create the 6-hourly GV and Optimal currents (which are originally daily data) and finally perform the colocalization in time.

Similarly to RS18, the validation is based on the computation of the RMSE, correlation coefficients (CC) and biases between the GV, the Optimal currents and the drifter-derived currents, the latter being our benchmark. The computation of the biases did not evidence significant differences between the GV and the Optimal currents. For both datasets, the basin scale bias of the zonal and meridional component of the motion is around −0.3 cm/s and 0.4 cm/s, respectively. Concerning the RMS, the Optimal currents always exhibited lower values than the ones obtained with the GV, with larger RMSE reduction for the meridional component of the motion. Moreover, we found that the difference between the GV and the Optimal currents RMSE increases with the increasing magnitude of the SST spatial gradients (A,B in Equation (4)), indicating that our reconstruction method gives the best performances in large SST gradient areas. This was found selecting all the GV and Optimal current values laying in areas of progressively increasing |∇SST|, with |∇SST|=√*A*<sup>2</sup> <sup>+</sup> *<sup>B</sup>*<sup>2</sup> comprised between 0.1 and <sup>6</sup>×10−<sup>5</sup> ◦C·m−<sup>1</sup> (such interval represents more than 90% of the |∇SST| values in the basin). Based on the RMSE computation, we also defined a parameter to estimate the percentage of improvement of our reconstruction method with respect to the GV estimates. The parameter is defined as follows:

$$\text{IMPROVE}\_{\text{(U,V)}} = 100 \times \left[ 1 - \left( \frac{\text{RMSE}\_{\text{(U,V)}}^{\text{Optimal}}}{\text{RMSE}\_{\text{(U,V)}}^{\text{CV}}} \right)^2 \right] \tag{7}$$

where U,V respectively indicate the zonal and the meridional component of the sea surface currents.

In general, our reconstruction is more satisfying for the meridional component of the motion and, at the basin scale, yields improvements that can exceed 10% (while they are only about 5% for the zonal component of the motion). Such improvements increase almost monotonically with the magnitude of the |∇SST| except for the zonal currents, where a linear increasing trend is only observed until 5 × <sup>10</sup>−<sup>5</sup> ◦C·m−1. In a similar fashion, the Optimal currents show correlations with the in-situ measured currents generally higher than the GV, with a correlation improvement which is larger for increasing |∇SST| values. Such results, schematically summarized in Figure 9, indicate that the Optimal currents better represent the sea surface currents variability observed via the in-situ measurements, confirming the potentiality of bringing the spatio-temporal variability of the SST remote observations inside the large-scale geostrophic currents.

During 2012–2016 we also investigated the local variability of the Optimal reconstruction method. This has been done binning the values obtained using Equation (7) in 2◦× 2◦ boxes (Figure 10 upper panels). The available drifter-derived currents allowed to perform the analysis in the sites of validation (SOV) indicated by the black squares in Figure 10 (bottom panel). In general, the best performances are found in the Central and Western Mediterranean, where local improvements can also reach 30%. The results of the validation, both at the basin scale and in terms of local variability, are in agreement with the results found at global scale by RS18. Also, we found that in the Western Mediterranean, the summertime reconstruction can lead to improvements exceeding 40%, which was not found during winter (not shown). Considering our investigation area, the zonal component of the surface currents is improved over 60% of the domain, while this percentage raises up to 70% for the meridional component of the motion. Combining these results, we could also indicate all the SOV where the method can improve both of the surface flow components, schematically given in the bottom panel of Figure 10. As expected, such SOV are mostly located in the Western Mediterranean. Unfortunately, some of the sites of the Optimal currents validation, though never exceeding the 40% of the analyzed

domain, indicate that the method can also degrade the quality of the surface currents, when compared to the geostrophic estimates. These sites are mostly located in the eastern section of the Mediterranean basin. Quite interestingly, they are characterized by two properties: either they lay in areas of low <sup>|</sup>∇SST| ( 0.1 × <sup>10</sup>−<sup>5</sup> ◦<sup>C</sup> ·m<sup>−</sup>1), which is not the ideal condition for the method applicability, or their statistical properties are computed based on a poor number of observations, affecting their statistical significance. This can be checked comparing the maps of |∇SST| and numbers of available in-situ observations binned in 1◦× 1◦ boxes in Figure 11. In a future study we plan to compute such statistics on a longer time series of optimal currents ( 20 years) and check the effects on the local variability of the optimal reconstruction method performances.

**Figure 9.** Basin scale improvements of the Optimal reconstruction method expressed via root mean square error (RMSE), Correlation Coefficients (Corr) and percentage of improvement (% IMPROVE) in the period 2012–2016. The Root mean square errors and correlation coefficients are computed for both the Optimal currents (blue lines) and the geostrophic velocities (GV) (red lines) using the drifting buoys-derived surface currents as a benchmark. The results are presented for both the zonal (U) and the meridional component (V) of the surface motion, respectively shown in the left and right column of the figure.

**Figure 10.** Spatial variability of the percentage of improvement of the Optimal reconstruction method for the zonal (U) and meridional component (V) of the surface motion, respectively shown in the top and middle panel of the figure. The percentage of improvement, binned in 2◦× 2◦ boxes has been computed in the sites of validation (SOV) given by the black squares in the bottom panel. The crosses indicate the SOV where both the zonal and meridional components of the motion are improved.

Nevertheless, it is worth noticing that the validation against the in-situ observations confirmed the validity of the case study analyzed in Section 4.2. Indeed, in the coastal area of the Sicily Channel (close to the southern coast of Sicily), during the period 2012–2016, our surface-currents reconstruction method yielded positive improvements for both the components of the motion. Hence, we wanted to further investigate the quality of the Optimal currents via a comparison with the HFR-derived surface currents in the Malta–Sicily Channel. The results are presented in Section 4.4.

**Figure 11.** Top panel: |∇SST| binned in 1◦× 1◦ boxes, computed along the positions occupied by the drifting buoys during the period 2012–2016. Bottom panel: Number of in-situ surface-currents observations laying in each of the 1◦× 1◦ boxes.

#### *4.4. Comparative Analysis: Optimal Currents Validation in the Malta–Sicily Channel*

In this section we present an additional quality assessment for the Optimal currents via comparison with the daily surface currents derived from the Calypso High Frequency Radar platform in the Malta–Sicily Channel (during the year 2016), our benchmark. The coverage of the Calypso platform, schematically shown in Figure 12, extends from 35.7◦N to 36.8◦N and from 13.7◦E to 15.2◦E. As stated in Section 2, such platform provides hourly observations at 1/37◦ horizontal resolution, hence, our validation was realized creating daily maps of HFR currents. For each day, we computed bias and root mean square errors of the differences between the GV, the Optimal, the MERCATOR, the MFS-derived currents and the HFR estimates. Unfortunately, no comparison with the drifting buoys-derived currents was performed as, during the year 2016, the buoys have been circulating in the westernmost and easternmost sections of the Mediterranean basin only (Figure 10).

**Figure 12.** The Calypso High Frequency (HF) Radar platform: Schematics (after the Calypso Project web portal).

Unlike the case of the drifting buoys, the comparison with the High-Frequency Radar (HFR) currents involves two-dimensional gridded data. In order to carry out the comparison, all the analyzed daily currents were remapped over the HFR grid via bi-cubic interpolation. The mean results of the comparative study are summarized in Tables 1 and 2.

**Table 1.** BIAS and Root Mean Square Error (RMSE) between the Optimal-GV-MERCATOR-MFS and the High-Frequency Radar (HFR) surface currents in the Malta–Sicily Channel. Zonal Component.


**Table 2.** BIAS and Root Mean Square Error (RMSE) between the Optimal-GV-MERCATOR-MFS and the HFRsurface currents in the Malta-Sicily Channel. Meridional Component.


Focusing on both the Optimal and the GV, using the HFR estimates as benchmark, the current retrieval is more satisfying for the meridional component of the motion, exhibiting lower RMSE and BIAS (respectively 9.00 and 1 cm·s−1) compared to the zonal one, where RMSE and BIAS are respectively 10.00 and 4 to 5 cm·s<sup>−</sup>1. Moreover, if we observe the behaviour of the Optimal currents with respect to the other datasets, we find that the RS18 method brings larger improvements on the retrieval of the zonal component of the motion. Indeed, for such component the Optimal currents have the lowest RMSE of all, 9.72 cm·s−<sup>1</sup> (the MFS model yielding the larger value, 12.70 cm·s−1) and a lower BIAS (4.04 cm·s<sup>−</sup>1) than the GV and the MERCATOR estimates. The only exception is given by the BIAS of the MFS model, whose value is 3 cm·s<sup>−</sup>1.

For the meridional component of the motion, the Optimal currents have RMSE and BIAS in line with the GV and always lower than the numerical model estimates. These results may look surprising. Indeed, at the basin scale the RS18 method brought the largest improvements for the meridional component of the motion. Nevertheless, here we are restricting our analysis to a much smaller coastal

area in which our reference, the HFR currents, are known to be more accurate in retrieving the zonal component of the motion [43], probably affecting the statistics on the meridional flow.

Considering these results, we can conclude that the Optimal currents perform globally better than altimeter-derived estimates or model outputs.

#### **5. Discussion and Conclusions**

During the last decades, both physical and operational oceanography have been demanding high spatial and temporal resolution observations. Such observations are useful for the validation, optimization and run of the main operational models, especially when data assimilation is required. In the specific case of the sea surface currents retrieval, the high-resolution and repetitive observations are also necessary for several human activities in the marine context (navigation, safety and rescue activities) and for the environmental safeguard (illegal oil spills and pollutants monitoring).

The objective of our study was to improve the spatial and temporal resolution of the remotely sensed, altimeter-derived surface currents in the Mediterranean Sea. Based on the RS18 method, we could successfully introduce the information contained in a high-resolution satellite-derived SST (Δt = 1 day and Δx=1/24◦) inside the larger scale altimeter geostrophic currents. The RS18 method has the main advantage to optimize the use of SST to derive information on the surface currents field. This is achieved accounting for the source and sink terms in the SST evolution equation, unlike in maximum cross correlation techniques where only advection and diffusion can be considered. An additional advantage of the RS18 method is the possibility to rely on a background surface currents estimate, i.e., the altimeter measurements. This enables to derive the surface circulation even in areas of low tracer gradients or when the tracer alone is not a direct proxy of the ocean surface dynamics [16]. We applied the RS18 reconstruction method to the CMEMS datasets of geostrophic surface currents and sea surface temperature satellite data for the Mediterranean area. This enabled to estimate an Optimized velocity during five years.

The highlights of our study are listed here:


the satellite-derived currents (except for the bias in the retrieval of the zonal flow, where the MFS model exhibited the best performance).

The main results listed here express the capability of the RS18 method to increase the spatio-temporal variability of the altimeter-derived currents even in the Mediterranean Area. Provided the existence of non-zero sea surface temperature spatial gradients, the variability of the high-resolution SST data can be successfully introduced inside the large-scale geostrophic currents estimates. The increase in the spatio-temporal variability explains the satisfying comparison of the Optimal currents with both the in-situ measured currents and the HFR estimates. Moreover, this was confirmed by the spectral analysis shown in Figure 13. In particular, we computed the mean spatial Kinetic Energy (KE) spectrum of the Optimal, MFS and altimeter-derived surface currents (GV) as a function of the spatial wavenumber (here given as the inverse of the wavelenght.) This was done in a land-free and dynamically active area of the Western Mediterranean (37.5◦N to 38.0◦N and from 0◦E to 11◦E), delimited by the yellow box in Figure 13. Moreover, we performed this analysis done during summertime (on 23 July 2016) i.e., choosing all the favourable conditions for the Optimal reconstruction method. Looking at the behaviour of the KE power spectral density, we can see that the Optimal currents and GV power spectra are superimposed until a scale of 100 km. Afterwards the GV begin loosing energy and their spectrum falls abruptly. Such a result was expected and is intrinsic with the geostrophic large-scale estimate of the surface currents, whose spectral content is mostly in the mesoscale range. On the other hand, for scales smaller than 100 km, the Optimal Currents spectrum contains more energy than the GV one and maintains the same decreasing trend until a scale of 30 km, confirming that the resolution of the Optimal currents is actually increased with respect to altimetry (notice that in the figure we discarded the areas in which the spectra exhibit a noisy behaviour). Moreover, the Optimal Currents spectrum is in agreement with the predictions of the MFS model, that, based on a fully three-dimensional primitive equations framework, simulates the evolution of the total currents field. Quite interestingly, Figure 13 also confirms that the Optimal and MFS spectra are closer to the predictions of the two dimensional turbulence energy cascade, K−5/3 and K−<sup>3</sup> [44]. Though the K−<sup>3</sup> slope is not fully recovered at small scales, the Optimal Currents spectral separation from the K−<sup>3</sup> line is reduced compared to the altimeter estimates.

In the end, we computed the GV and Optimal Currents KE temporal spectra in the Western Mediterranean (from 6◦W to 8◦E), where the optimal reconstruction exhibited the largest improvements (Figure 14). This was done on a time window of 45 days, which is the maximum temporal correlation scale for building the L4 altimeter-derived currents from the altimetric ground-tracks [28]. In particular, the computation was performed during wintertime (Julian day 1 to 45 of the year 2016), when the basin scale averaged KE was found to be larger. Figure 14 shows the Optimal and GV currents KE temporal spectra given as a function of the temporal wavenumber (here given as the inverse of the wavelength). The Optimal currents spectrum is initially superimposed to the GV one until a timescale of about twenty days. Afterwards, due to a power spectral density drop in the geostrophic velocities, the two spectra separate, indicating that the Optimal currents contain more energy up to scales of one week, when the GV spectrum starts exhibiting a noisy behaviour. This analysis confirms the Optimal Currents enhanced temporal variability.

This work could certainly benefit from the following additional analyses. At first, we plan to extend the Optimal currents validation to the years ranging from 1993 to 2016. This would help increasing the statistical significance of the comparisons with the drifting buoys in all the areas of the basin. In addition, using a 20 years long time series, would allow to statistically determine the different contributors to the improved currents retrieval. Indeed, the improvements may come from both a better description of the non geostrophic flow and more accurate retrieval of the eddy dynamics, which were here discussed as single case studies.

**Figure 13.** Optimal (blue line), altimeter-derived (red line) and MFS (light blue line) kinetic energy spatial spectra. The mean spectra are computed in the area delimited by the yellow box (top panel). The black continuous and dashed lines represent the surface and internal predictions of the QG turbulence, respectively. The spectra are computed on the date of 23 July 2016 and are represented as a function of the spatial wavenumber.

**Figure 14.** Temporal spectral content of the Optimal currents (blue) and GV (red) kinetic energy (west Mediterranean).

New ways of estimating the forcing term in the SST evolution equation should also be exploited. In the present study, the forcing was approximated by a low pass filtered SST temporal derivative but, in principle, heat fluxes estimates given by model-derived reanalyses, e.g., the ECMWF ERA-5, or estimated using bulk formulae could be tested [45,46]. In addition, the possibility to refine the temporal resolution of the method input parameters (e.g., *σ*u,v and h), presently given by climatological two-dimensional maps, could further improve the reconstruction method capabilities. At present, the definition of *σ*u,v and h is one of the more critical and crucial points of the RS18 method, being strictly dependent on the availability of high quality in-situ measurements of SST and surface currents (as shown in Section 3).

In addition, this study could further benefit from the exploitation of the gap-free SST obtained by geostationary satellites (e.g., METEOSAT [47]). Such an approach could lead to the determination of the sub-daily currents variability, up to the hourly resolution.

In the end, this method could be further applied on other types of satellite derived tracers. For instance, the ocean colour products are available at horizontal resolutions even higher than 1 km, e.g., the observations from the OLCI sensor mounted on board Sentinel-3 or from the MSI sensor mounted on Sentinel-2 and acquiring information in coastal areas. Implementing the method with these observations could allow to further improve the retrieval of the surface and near-surface submesoscale dynamics. Nevertheless, such an approach is not straightforward as it will require to accurately assess the source and sink terms due to biological activity in the tracer conservation equation (Equation (1)). This operation will be crucial for extracting the information related to the surface currents advection.

**Author Contributions:** Conceptualization: D.C., M.-H.R. and R.S. Formal analysis: D.C. and M.-H.R. Funding acquisition: R.S. Investigation: D.C., M.-H.R. and R.S. Methodology: M.-H.R. and R.S. Supervision: M.-H.R., M.M. and R.S. Validation: D.C., M.-H.R. and M.M. Writing of original draft: D.C. Writing—review and editing: D.C., M.-H.R., M.M. and R.S.

**Funding:** This research was funded by Progetto Bandiera RITMARE (Ricerca ITaliana per il MARE), the project "Marine Strategy" of the Italian Ministry of Environment, the Consiglio Nazionale delle Ricerche-Collecte Localisation Satellites Service Agreement and the European Copernicus programme.

**Acknowledgments:** We are grateful for the helpful comments on the manuscript by the three anonymous Reviewers. Moreover, we wish to thank Antonio Olita for providing the CNR-Sea Forecast model outputs and Salvatore Marullo for the fruitful discussions. We also wish to acknowledge Sharon Fan for taking care of the reviewing process.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
