Estimating ocean currents from the joint reconstruction of absolute dynamic topography and sea surface temperature through deep learning algorithms

Our study focuses on absolute dynamic topography (ADT) and sea surface temperature (SST) mapping from satellite observations, with the primary objective of improving the satellite-derived ADT (and derived geostrophic currents) spatial resolution. Retrieving consistent high-resolution ADT and SST information from space is challenging, due to instrument limitations, sampling constraints, and degradations introduced by the interpolation algorithms used to obtain gap-free (L4) analyses. To address these issues, we developed and tested different deep learning methodologies, specifically convolutional neural network (CNN) models that were originally proposed for single-image super resolution. Building upon recent findings, we conduct an Observing System Simulation Experiment (OSSE) relying on Copernicus numerical model outputs (with respective temporal and spatial resolutions of 1 d and 1/24°), and we present a strategy for further refinements. Previous OSSEs combined low-resolution L4 satellite equivalent ADTs with high-resolution “perfectly known” SSTs to derive high-resolution sea surface dynamical features. Here, we introduce realistic SST L4 processing errors and modify the network to concurrently predict high-resolution SST and ADT from synthetic, satellite equivalent L4 products. This modification allows us to evaluate the potential enhancement in the ADT and SST mapping while integrating dynamical constraints through tailored, physics-informed loss functions. The neural networks are thus trained using OSSE data and subsequently applied to the Copernicus Marine Service satellite-derived ADTs and SSTs, allowing us to reconstruct super-resolved ADTs and geostrophic currents at the same spatiotemporal resolution of the model outputs employed for the OSSE. A 12-year-long time series of super-resolved geostrophic currents (2008–2019) is thus presented and validated against in situ-measured currents from drogued drifting buoys and via spectral analyses. This study suggests that CNNs are beneficial for improving standard altimetry mapping: they generally sharpen the ADT gradients, with consequent correction of the surface currents direction and intensities with respect to the altimeter-derived products. Our investigation is focused on the Mediterranean Sea, quite a challenging region due to its small Rossby deformation radius (around 10 km).

1 Introduction

Oceanic currents play a pivotal role in influencing both short-term and long-term dynamics within the ocean–atmosphere system. Monitoring these currents on a large scale is essential for assessing the transport of heat and salt and enhancing our ability to predict variations and shifts in ocean dynamics and their impact on marine ecosystem. At smaller oceanic scales, such as mesoscale and submesoscale, observing and predicting ocean currents is key to understand and model Earth system dynamics. Mesoscale eddies, oceanic features with length scales from 10 to 100 km, and persisting on timescales from weeks to months, can migrate over considerable distances while carrying heat, salt, and nutrients. These eddies also introduce perturbations that drive significant vertical exchanges. Submesoscale features, including eddies, fronts, and filaments, operate at spatial and temporal scales of 0.1 to 10 km and hours to days. These transient, small-scale features introduce notable variations in both horizontal and vertical velocities, profoundly affecting local three-dimensional transport properties (Clarke and Li, 2004; Li and Clarke, 2004; Carlson and Clarke, 2009; Siokou-Frangou et al., 2010; Buongiorno Nardelli, 2013; Barbosa Aguiar et al., 2013; Ponte et al., 2013; Frenger et al., 2013; Bashmachnikov et al., 2015; Chenillat et al., 2016). The monitoring of ocean currents is also crucial for various societal and environmental purposes, including aiding ship navigation, supporting safety and rescue operations, and managing marine ecosystem services (Pisano et al., 2016, and Onink et al., 2019). These applications necessitate precise, high-resolution tracking of surface oceanic currents.

Since the early 1990s, satellites equipped with nadir-looking radar altimeters have provided indirect observations of the marine surface circulation at a global scale. This is accomplished by measuring the absolute dynamic topography (ADT) with respect to a reference surface (the geoid) along 1D tracks and inferring surface motion from interpolated 2D ADT maps using the geostrophic approximation. However, this method has inherent limitations tied to ADT sampling and the geostrophic approximation, primarily capturing only larger mesoscale geostrophic processes, O(100 km, 10 d) (Pascual et al., 2006; Pujol et al., 2012, 2016; Ballarotta et al., 2019).

The direct estimation of marine surface (or near-surface) currents relies on satellite radar interferometry techniques or in situ measurements from Lagrangian buoys, ship-mounted devices such as the acoustic Doppler current profiler (ADCP), or high-frequency radar (HFR) platforms. Lagrangian observations can serve as reference points to validate remotely sensed surface currents and, when properly organized in space and time, can yield pseudo-Eulerian surface circulation estimates. However, this approach faces constraints related to the spatial–temporal coverage of Lagrangian platforms and their tendency to become trapped in oceanic recirculation or convergence zones. Conversely, HFR systems offer comprehensive maps at fine spatial and temporal resolutions (less than 10 km and 1 h, respectively), albeit only within coastal regions (Chapron et al., 2005; Falco and Zambianchi, 2011; Lumpkin et al., 2017; Laurindo et al., 2017; Capodici et al., 2019; Ribotti et al., 2023; Fanelli et al., 2024b).

The fusion of altimeter-derived and in situ-measured currents represents a valuable strategy for enhancing surface current estimates from altimetry, in both coastal and open-ocean regions (Mulet et al., 2021; Ballarotta et al., 2023). Nonetheless, the effectiveness of this approach is contingent on the availability of in situ measurements. Alternatively, several studies proposed the merging of altimeter-derived products with independently observed satellite-derived tracers, like sea surface temperature (SST) and surface chlorophyll concentration data (González-Haro and Isern-Fontanet, 2014; Rio and Santoleri, 2018; Ciani et al., 2020; González-Haro et al., 2020; Miracca-Lage et al., 2022). Such methodologies turned out to be useful in improving the altimeter-derived geostrophic circulation with limitations related to the season and/or the geographic location, making it challenging to implement them operationally at a global scale.

High-resolution monitoring of ocean dynamic topography using imaging sensors, which natively provide two-dimensional observations, has only recently become available thanks to the successful launch of the Surface Water and Ocean Topography (SWOT) mission in December 2022 (https://swot.jpl.nasa.gov/mission/, last access: 16 August 2024). This is currently providing 2D images of sea surface height at unprecedented spatial resolutions, allowing us to better characterize the signatures of oceanic mesoscale features from altimetry maps. The SWOT revisit time is, however, set to 11 d, limiting the high rate monitoring of the fast-evolving or persistent oceanic features (Fu et al., 2009; Fu and Ubelmann, 2014; Morrow et al., 2019, 2023).

Here, we propose an ocean surface current reconstruction methodology based on artificial intelligence, also following recommendations from the international altimetry team (Abdalla et al., 2021) and recently implemented by Beauchamp et al. (2022), Buongiorno Nardelli et al. (2022), Martin et al. (2023, 2024), Fablet et al. (2023), Moschos et al. (2023), Kugusheva et al. (2024), and Archambault et al. (2024).

In particular, we rely on a family of deep learning methods known as convolutional neural networks (CNNs) for super resolution, firstly proposed for computer vision purposes (e.g., Dong et al., 2015). CNNs for super resolution learn a direct mapping between low-resolution and high-resolution images and are here applied to the case of satellite-derived images of absolute dynamic topography (ADT) and sea surface temperature (SST). Our primary focus is to super-resolve ADT maps exploiting information from low-resolution interpolated ADTs and higher-resolution interpolated SSTs, improving both in a joint reconstruction. Then, the super-resolved ADT maps are used to derive the oceanic surface circulation via the geostrophic balance. The study builds on, but also constitutes a substantial advance with respect to, the results firstly presented by Buongiorno Nardelli et al. (2022) (BBN22 hereinafter), which consisted in training a CNN through an Observing System Simulation Experiment (OSSE) and then using the same CNN model to reconstruct ocean surface currents from true satellite products. BBN22 also presented a sensitivity study to assess the performances of different CNN architectures, mainly testing networks with different numbers of tunable parameters and different perceptive capabilities. However, BBN22 performed the OSSE only considering the limitations of the altimeter products (i.e., accounting for the interpolation of observations in order to obtain 2D ADT maps from the along-track observations), thus assuming a perfectly known SST. Our work improves BBN22 in two particular ways: (i) we generate a satellite equivalent SST time series accounting for a realistic operational level 4 (gap-free) analysis based on gapped input data, as in Buongiorno Nardelli et al. (2013), and (ii) we train our CNN relying on updated loss functions, which now include physics-based constraints, and exploit the joint reconstruction of ADT, SST, and ∂_tSST data, also imposing physics-based constraints. The study focuses on the Mediterranean Sea area, dominated by motions with length scales down to 6 km (Malanotte-Rizzoli, 2014), thus constituting a challenging test bed for the super resolution of standard altimetry products.

2 Materials and methods

The datasets involved in our study are basically those described in BBN22 and Ciani et al. (2021), with the exception of the newly generated synthetic SSTs described in Sect. 2.5.2. For the sake of clarity, we provide a brief description of all datasets in the present paper as well.

2.1 Numerical model

The Mediterranean Forecasting System (MFS) is a hydrodynamic model designed for the Mediterranean Basin and the easternmost section of the Atlantic Ocean near the Strait of Gibraltar. It provides 3D horizontal current and sea surface height (SSH) outputs ranging from monthly to 15 min intervals, as well as 3D temperature and salinity fields with monthly to hourly estimates. These data are accessible via the Copernicus Marine Service web portal (product ID: MEDSEA-ANALYSIS-FORECAST-PHY-006-013). For our current research, we utilized daily SSH and SST data extracted within the boundaries of the Mediterranean Basin (30 to 46° N and 6° W to 37° E). These datasets are provided on a 1/24° regular grid with 125 unequally spaced vertical levels. The simulations are based on the NEMO model (Nucleus for European Modelling of the Ocean) in combination with WAVEWATCH III for the wave component. The MFS simulations also incorporate data assimilation from 2D satellite-derived SST, vertical salinity profiles, and sea level anomaly observations along satellite tracks (Clementi et al., 2019).

2.2 Altimeter ADT and derived quantities

The altimeter-derived ADTs were obtained from the Copernicus Marine Data Store. The surface currents, for consistency with the results presented in Sect. 3.2, were derived from ADT using the geostrophic approximation equation, i.e., by applying a finite-central-differences operator to the gridded L4 ADTs in order to compute its spatial gradients. The Copernicus ADTs are provided as daily data with a nominal 1/8° horizontal resolution. We extracted the time series covering the years 2008 to 2019. The corresponding Copernicus Marine Service product and dataset IDs are SEALEVEL-MED-PHY-L4-REP-OBSERVATIONS-008-051 and dataset-duacs-rep-medsea-merged-allsat-phy-l4, respectively (last access: 1 March 2021 and now included as part of the SEALEVEL-EUR-PHY-L4-MY-008-068 product). To match the resolution of the numerical model outputs used in our study, the altimeter ADTs are up-sized to 1/24° via cubic interpolation.

2.3 Satellite SST

We obtained remotely sensed SST data from the Copernicus Marine Service (https://doi.org/10.48670/moi-00172). These are L4 products, ensuring gap-free estimates of the foundation temperature (i.e., at ≃ 10 m depth) provided on a regular grid, and are produced and distributed operationally in near real time. Specifically, we used 12 years of the ultra-high-spatial-resolution (UHR) Mediterranean dataset spanning from 2008 to 2019, with a nominal resolution of 1/100° (product ID: SST-MED-SST-L4-NRT-OBSERVATIONS-010-004-c-V2). This SST product is generated by combining nighttime images collected by satellite infrared sensors after rigorous quality control, removal of cloudy pixels, and application of an optimal interpolation algorithm. In order to homogenize the satellite SST L4 data with that used for training the neural network model, we evaluated the effective spatial scales captured by the model SST and filtered the UHR via a low-pass, cut-off wavelength Lanczos filter to obtain comparable effective resolutions. Successively, the satellite SSTs were remapped onto the final (model) 1/24° grid through bilinear interpolation.

2.4 In situ measurements

In situ measurements of sea surface currents were obtained from autonomous Lagrangian drifting buoys, which are transported passively by ocean surface currents. During the buoy's drifting process, positional data are interpolated at regular intervals (approximately every 30 min) using the kriging interpolation method developed by Poulain et al. (2012). Velocities are subsequently computed through a finite-differences method applied to the interpolated positions and are provided with 6-hourly temporal resolution. The data covering the period of our study were originally provided by the Italian Institute of Oceanography and Experimental Geophysics (OGS) for the ESA CIRCOL project (http://circol.artov.ismar.cnr.it/, last access: 29 September 2024). These time series are accessible via https://doi.org/10.6092/7a8499bc-c5ee-472c-b8b5-03523d1e73e9. It is worth noting that buoy-derived surface current values are retained only if the buoy is equipped with a drogue, a device that ensures the buoy's movement is primarily driven by ocean currents rather than surface winds (Menna et al., 2018).

2.5 Synthetic satellite equivalent ADT and SST

2.5.1 Satellite equivalent ADT

We generated 1 year (2017) of synthetic, satellite equivalent altimeter-derived absolute dynamic topography (SE-ADT hereinafter) maps using outputs from the Copernicus Marine Service MFS hydrodynamic simulation, employing the Data Unification and Altimeter Combination System (DUACS) mapping method. The steps involved in this process are detailed below. Initially, the sea level anomaly (SLA) was calculated from model outputs according to the following expression:

$$\begin{array}{}\text{(1)}& \mathrm{SLA}=\mathrm{SSH}-(\mathrm{MDT}-\mathrm{0.344}).\end{array}$$

Here, the mean dynamic topography (MDT) is provided as a static field along with the model outputs. A constant value of 0.344 (expressed in meters) is used to adjust the SLA values in the Mediterranean Sea to ensure that the spatiotemporal average of SLA is zero for the year 2017. To remove the large-scale, high-frequency variability typically due to dynamic atmospheric processes, we applied a Loess filter to these synthetic data. Subsequently, the SLA was sampled along the actual paths of a synthetic constellation comprising four radar altimeters: Jason-3, Sentinel-3A, SARAL/Altika, and CryoSat-2 missions. This step was executed using the SWOT simulator software, which accounts for the actual orbits, errors, and noise associated with each mission. The chosen four-satellite constellation represents the constellation used in the Copernicus Marine Service processing during 2017. These along-track synthetic measurements were then incorporated into the DUACS processing chain to generate L4 SLA maps. The optimal interpolation (OI) scheme follows the DUACS DT2018 (delayed time) configuration for the Mediterranean area, as described in Taburet et al. (2019). The reconstructed L4 maps were then combined with the filtered large-scale maps in order to obtain the ADT. These data are provided on a daily basis and are available on a regular 1/8° grid (further details can be found in Ciani et al., 2021). A subsequent up-size to 1/24° is applied before using this dataset for our application.

2.5.2 Satellite equivalent SST

The generation of the synthetic satellite equivalent SST (SE-SST) follows the same principle as for the SE-ADT. The SE-SSTs combine information from the model-derived SST and the Copernicus Marine Service merged multi-sensor L3S (Level 3, super-collated) satellite SSTs for the Mediterranean area (product ID: SST_MED_SST_L3S_NRT_OBSERVATIONS_010_012, last access: 20 October 2023). The L3S SSTs are a satellite-derived product that considers information from several infrared radiometers and provides daily ocean surface temperature on a regular latitude–longitude grid at high (1/16°) and ultra-high (1/100°) spatial resolution, representative of nighttime SST values (00:00 UTC). As such, the L3S SSTs contain gaps whenever the infrared SST retrieval is not possible (e.g., due to cloud cover) or the single-sensor satellite SSTs have been labeled as poor-quality observations. Starting from 1 year of model-derived SSTs, we remapped gappy (missing values) patches found in the 1/16° L3S product (year 2017) onto the original modeled SST, thus generating a synthetic model-derived L3S SST time series. Finally, the gap-free (Level 4) SE-SSTs, including an estimate of the uncertainty on the Level 4 analysis, are obtained through standard OI via a dedicated algorithm reproducing the main steps described in Buongiorno Nardelli et al. (2013), constituting the backbone of the present-day L4 satellite SST operational production for the Mediterranean area. An inter-comparison of the synthetic gapped SST, SE-SST generated via OI, and ground truth SST is depicted in Fig. 1, where we notice a strong smoothing and underestimation of the mesoscale SST features in correspondence of cloudy areas, as expected for present L4 SST satellite products. In the end, the ∂_t (SE-SST) and its error are computed from the SE-SST time series via central finite differences and standard error propagation, respectively. The SE-SST, ∂_t (SE-SST), and their errors are daily fields provided on 1/24° regular grid.

Figure 1 Results of the OSSE for the satellite equivalent (SE) SST generation on 16 May 2017. (a) SE L3S gapped SST; (b) SE-SST generated via OI algorithm; (c) ground truth SST.

2.6 Ocean currents and SST reconstruction methodology

In this study, ocean currents are reconstructed through a joint ADT–SST super resolution relying on a CNN approach, with a primary focus on ADT. The workflow of the reconstruction exercise is depicted in Fig. 2. The work consists of three phases: (i) the satellite equivalent data generation (detailed in Sect. 2.5.1 and 2.5.2), (ii) the training of the network by means of an Observing System Simulation Experiment (OSSE), and (iii) the application of the neural network model (optimized via the OSSE) for the reconstruction of the ocean surface currents from satellite-derived input ADTs and SSTs. The OSSE consists in the generation of super-resolved ADTs and SSTs which are then directly compared with an independent test dataset detailed in Sect. 2.6.1. In addition, in the present study we will also inter-compare our results with the ones obtained by BBN22, i.e., an earlier version of the dilated adaptive multi-scale residual super-resolution CNN, detailed in Sect. 2.6.1. The super-resolved geostrophic currents obtained from satellite-derived data are provided at 1/24° (to match the data used for training the neural network with the OSSE), and their performance assessment is carried out via direct comparison with in situ-measured currents (see Sect. 2.4) and spectral analysis.

Figure 2 Workflow of the study. Blue stands for generation of the satellite equivalent input data, orange blocks indicate the OSSE, and green blocks refer to the application of the reconstruction methodology to satellite-derived data.

2.6.1 Convolutional neural network architecture and training strategy

Convolutional neural networks used for super resolution are built to learn a mapping between low-resolution and high-resolution images. This is achieved by relying on convolution operators designed to detect and learn specific features in the input images. The network learns the end-to-end mapping from low- to high-resolution images via a minimization of the error between the output and the validation data. This is quantitatively achieved using loss functions which are iteratively minimized by relying on a validation dataset (i.e., a fraction of the dataset involved in the training procedure). The overall performance of the CNN is finally quantified by means of an independent test dataset.

The dataset preparation is carried out as follows. (i) Firstly, we extracted 40 dates from the model outputs detailed in Sect. 2.1, which served as independent test dataset. The 40 dates are distributed along the year in order to cover all dynamical regimes. (ii) Starting from the remaining model data, we performed a resampling into 76 × 100 pixels (corresponding to ≃ 300 × 400 km) tiles. The tile extraction also considered a 50 % spatial overlap, resulting in a total of 42 250 samples. All the samples then undergo a normalization procedure with respect to the maximum observed in the time series. In addition, for SE-ADT, SE-SST, and ∂_t(SE-SST), we compute anomalies filtering out signals larger than approximately 200 km. This operation aims at retaining signals related to the local SST variations due to the horizontal current advection. In this way the extraction of patterns from the SST field can help the ADT reconstruction with contributions primarily due to the ocean surface dynamics. Such samples were thus used for the network training and validation procedure. Here, compared to BBN22, we changed the splitting strategy between training and validation datasets, whose relative amount is 85 % and 15 % of the 42 250 (total) samples, respectively. In particular, we forced the validation dataset to be a time series of samples adjacent in time (during the late fall/early winter season), instead of applying a random selection from the available samples. This strategy is used to train the network, preferably selecting scenarios where enhanced small-mesoscale and submesoscale activity is expected, as further detailed in Sect. 4.

The CNN employed here is called the dilated adaptive deep residual network for super resolution (dADR-SR), depicted in Fig. 3. In previous formulations (BBN22), the network considered

• four predictors, namely the SE-ADT, SE-ADT error, and SST and its temporal derivatives (∂_tSST), and

• one target, i.e., the super-resolved ADT.

The predictors were chosen considering the (local) modulation of SST features by water mass advection as well as the present-day altimeter observation geometry and satellite observation repetitiveness. The main upgrade of the present study is the introduction of additional inputs (predictors) like the SE-SST and ∂_t(SE-SST) errors as well as two additional targets/outputs, i.e., the super-resolved SST and its temporal derivatives (∂_t(SE-SST)). This upgraded architecture, depicted in Fig. 3, allows us to perform a more realistic ADT and SST reconstruction exercise, also accounting for the limitations of the present-day satellite SST operational retrieval.

The basic architecture of the CNN is briefly recalled here. In the dADR-SR network, the low-resolution input dataset initially passes through three parallel dilated convolutional layers, each containing 10 convolutional filters characterized by a 3×3 kernel and an increasing dilation factor (1, 3, and 5, respectively). Dilated convolution enhances the CNN perceptive capability, enabling the extraction of information at increasingly larger scales and maintaining the same computational cost of a standard 3×3 convolutional filter. Following this initial stage, the data undergo a sequence of 12 multiscale adaptive residual blocks (MARBs). Each of these blocks incorporates two sets of parallel dilated convolutional layers, featuring 120 and 10 filters, respectively, also including a squeeze-and-excitation (SQ-EX) module. The SQ-EX block functions as a channel attention mechanism, adaptively scaling information coming from the several input data before summing all contributions to produce the final high-resolution output.

The algorithm is trained applying an early stopping rule, terminating the training process when the validation loss function increases for a predefined number of epochs, as determined by the patience parameter, set to 5 in our study. An adaptive learning rate (initialized at lr = 10⁻⁴) is employed, and we rely on the Adam optimizer with the same configuration used in BBN22. The dADR-SR training model ultimately utilizes nearly 1.6 million trainable parameters.

We adopt the mean squared error as the initial loss function (LF), although physics-informed constraints are introduced, as detailed in Sect. 3.

Finally, for the CNN test/prediction phase, we perform a tile-by-tile reconstruction, and we merge the output high-resolution tiles (considering a weighted average in the overlaps) in order to produce maps covering the entire study area. This is done on a tile-by-tile basis, with each tile incorporating a 50 % overlap in both the longitudinal and latitudinal directions. During this process, the central regions of the tiles are expected to achieve superior performance, whereas edge effects can lead to spurious features near the tile boundaries related to the application of convolutional kernels. To mitigate these edge-related artifacts, a pixel-wise weighting function is applied, progressively decreasing the weight assigned to pixels with increasing distance from the tile centers. This approach enables a seamless basin-wide reconstruction from the 76×100 tiles.

Figure 3 CNN architecture. The set of input tiles, from left to right, respectively, indicate SE-ADT, SE-ADT error, SE-SST, SE-SST error, ∂_t(SE-SST), and ∂_t(SE-SST) error. The three output tiles, from left to right, respectively, indicate the SR-ADT, SR-SST, and SR-∂_tSST. In the figure, Conv2D stands for 2D convolutional filter, Dil stands for dilation factor, F is the number of filters, ADD stands for aggregation, SQ-EX is the squeeze-and-excitation module, and MARB indicates a multiscale adaptive residual block (further detailed in the panel highlighted by the dashed lines). The ADT output file is emphasized by a thick edge to indicate the focus of the present study.

3 Results

3.1 The Observing System Simulation Experiment

The results of the OSSE are presented below. We firstly discuss a test case for the ADT reconstruction on 4 January 2017, depicted in Figs. 4 and 5.

Figure 4 Results of the OSSE on 4 January 2017. (a) SE-ADT; (b) super-resolved ADT; (c) model ADT (ground truth); (d) comparative spectral analysis of the ADT maps (SE-ADT, SR-ADT, and ground truth are given in blue, red, and green, respectively). Results refer to the 2D box depicted in panel (c).

Figure 5 Results of the OSSE on 4 January 2017. (a) Satellite equivalent (SE) current velocity from SE-ADT; (b) current velocity from super-resolved (SR) ADT; (c) current velocity from model ADT (ground truth); (d) comparative spectral analysis of the current velocity maps (SE-ADT, SR-ADT, and ground truth are given in blue, red, and green, respectively). Results refer to the 2D box depicted in panel (c).

Figure 6 Results of the OSSE on 4 January 2017. Kinetic energy spectra of the surface currents derived from satellite equivalent (SE) ADT (blue), super-resolved (SR) ADT (red), and model outputs (ground truth) ADT (green). Results refer to the 2D box depicted in Fig. 4c. The continuous and dashed black lines refer to the k⁻³ and $k^{-\mathrm{5}/\mathrm{3}}$ slopes, respectively. The 50 km scale is reported for reference.

The ADT maps obtained through CNN for super resolution (panel b), compared with the ones given by standard altimetry processing (in panel a), exhibit features in good agreement with the model outputs (our ground truth, panel c). Visual inspection suggests an overall sharpening of the basin-scale ADT gradients and the potential of the CNN approach to overcome the dynamical feature distortion due to standard altimetry mapping. This is observed in many cases along the Algerian coast, in proximity of the Gulf of Lion, as well as in the Adriatic Sea and the Levantine Basin. Choosing a land-free area in the Levantine Basin (depicted in Fig. 4c) and relying on fast Fourier transform (FFT) analysis, we also quantified the spectral properties of the three aforementioned ADT estimates, as in Droghei et al. (2018). At a large scale (< 200 km) the three spectra exhibit a similar power spectral density (PSD), indicating a similar description of the largest mesoscale oceanic features. Progressively approaching smaller scales, i.e., from ≃ 100 km downward (1 per degree wavenumber onward), the super-resolved ADT spectrum (SR-ADT, red line in Fig. 4c) evolves in fairly good agreement with the ground truth (green line in Fig. 4c), confirming an improved representation of smaller mesoscale features compared to standard altimetry products. The SR-ADT spectrum eventually shows the injection of noise below scales of ≃ 20 km, as confirmed by a flattening of the spectrum. The aforementioned sharpening of the super-resolved ADT gradients is further confirmed by the analyses reported in Fig. 5. In particular, we show maps of geostrophic current velocity derived from the SE-ADT, the super-resolved ADT, and the ground truth ADT, shown in panels (a)–(c), respectively. The SE geostrophic current velocity is characterized by flows seldom exceeding 0.75 m s⁻¹, with broader spatial patterns compared to the ones derived from the ground truth ADT. On the other hand, the CNN super-resolved current velocity (Fig. 5b) exhibits sharper patterns at the basin scale, with flow intensities reaching 1 m s⁻¹, in agreement with our reference field, given in Fig. 5c. Accordingly, performing the spectral analysis on the current velocity fields (Fig. 5d) leads to the same results already discussed for ADT. The spectral analyses presented in Figs. 4 and 5 serve to assess the enhancement of the mesoscale activity in the SR fields compared to present-day altimetry mapping. A more-in-depth comparison with theoretical predictions of energy and enstrophy transfer is based on the computation of kinetic energy (KE), depicted in Fig. 6. We firstly compute the KE relying on the surface currents from the SE-ADT, SR-ADT, and model outputs (ground truth). The KE spectra are then presented along with the $k^{-\mathrm{5}/\mathrm{3}}$ and k⁻³ slopes, as in Ciani et al. (2019) and following Vallis (2006). With similar outcomes with respect to previous analyses, the ground truth KE (green line in Fig. 6) is in agreement with the predictions of energy and enstrophy transfer, following the k⁻³ slope for larger wavenumbers and the $k^{-\mathrm{5}/\mathrm{3}}$ for shorter ones. An expected exception occurs for wavenumbers approaching 10 per degree, where the spectrum exhibits an energy loss, suggesting that model outputs are unable to fully resolve the small mesoscale motion. The KE spectra derived from SR- and SE-ADTs are superimposed to the ground truth case only in the large mesoscale regime (wavenumbers < 1 per degree). As soon as the mesoscale regime is reached, the SE KE spectrum experiences a significant loss of energy, while the SR KE follows theoretical predictions until scales of ≃ 50 km, which is thus identified as the scale at which processes are fully resolved.

In addition, we present the CNN performances as root mean square error (RMSE) differences, ΔRMSE, between two versions of the ADT L4 mapping, using the independent test data (i.e., the original model-derived ADTs) as benchmark. The aforementioned L4 ADT mappings are (i) the operational standard altimetry system discussed in Sect. 2.5.1 and (ii) the CNN-based reconstruction proposed in the present study. In practice, whenever the ΔRMSE assumes positive values, the CNN-based reconstruction improves the capabilities of the standard altimetry system. In particular, Fig. 7 illustrates a comparative analysis between the optimal CNN reconstruction proposed by BBN22 (panel a) and the one resulting from upgrading the CNN architecture (panel b).

Figure 7 (a) ΔRMSE (in m) according to the CNN architecture of BBN22. (b) ΔRMSE (in m) following the CNN proposed in the present study. Positive values indicate an ADT mapping improvement with respect to standard altimetry.

In the BBN22 formulation, the OSSE indicates an overall improvement of the CNN-based reconstruction, with very few spots of degradation mainly distributed along the coasts of the Adriatic Sea, the Liguro-Provencal area, and few additional isolated spots distributed in the southern section of the basin (from 30 to 35° N). Upgrading the CNN still enabled us to observe improvements in more than 95 % of the basin, with isolated spots of degradation found in the western Mediterranean Sea and an overall decrease in the more intense ΔRMSE values. However, ΔRMSE areas around 2 cm, expressing an improvement of 60 % with respect to standard altimetry (computed as in Rio and Santoleri, 2018), are ubiquitous. The slight decrease in the ΔRMSEs compared to the BBN22 formulation is not surprising, as we presently train the neural network based on a degraded, satellite equivalent SST field, which likely impacts on the overall CNN performances, with the advantage of including information on the actual, current SST observing systems, also including a new training strategy which can likely impact the results at the OSSE level, as described thoroughly at the end of this section. Upgrading the CNN architecture, beyond relying on a realistic observing system, also enables the fine tuning of the LF during training. In the present-day CNN formulation, the LF is given by the following equation:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{LF}& =\mathit{\alpha }\left[\stackrel{\mathrm{‾}}{{\left({\mathrm{SST}}_{\mathrm{pred}}-{\mathrm{SST}}_{\mathrm{ref}}\right)}^{\mathrm{2}}}\right]\\ & +\mathit{\beta }\left[\stackrel{\mathrm{‾}}{{\left({\mathrm{ADT}}_{\mathrm{pred}}-{\mathrm{ADT}}_{\mathrm{ref}}\right)}^{\mathrm{2}}}\right]\\ & +\mathit{\gamma }\mathit{\left\{}\stackrel{\mathrm{‾}}{{\left[{\left({\partial }_{\mathrm{t}}\mathrm{SST}\right)}_{\mathrm{pred}}-{\left({\partial }_{\mathrm{t}}\mathrm{SST}\right)}_{\mathrm{ref}}\right]}^{\mathrm{2}}}\mathit{\right\}}+\mathit{\delta }{\mathrm{Loss}}_{\mathrm{phy}},\end{array}\end{array}$$

where

• ${\mathrm{Loss}}_{\mathrm{Phy}}=\stackrel{\mathrm{‾}}{\left({\left(\frac{\partial \mathrm{SST}}{{\partial }_{\mathrm{t}}}\right)}_{\mathrm{pred}}-\frac{g}{f}\frac{\partial {\mathrm{ADT}}_{\mathrm{pred}}}{\partial y}\frac{\partial {\mathrm{SST}}_{\mathrm{pred}}}{\partial x}\right.}$ $\stackrel{\mathrm{‾}}{{\left.+\frac{g}{f}\frac{\partial {\mathrm{ADT}}_{\mathrm{pred}}}{\partial x}\frac{\partial {\mathrm{SST}}_{\mathrm{pred}}}{\partial y}\right)}^{\mathrm{2}}}$;

• α=1, β=0.25, γ=0.67, and δ=0.025;

• the subscripts “pred” and “ref” stand for the output of the CNN prediction and the validation dataset used during training, respectively; and

• g is the gravity acceleration; f is the Coriolis parameter; and x, y, and t are the zonal, meridional, and temporal coordinates, respectively.

As such, the LF employed in our study quantifies the discrepancies between the predicted and the ground truth ADTs, SSTs, and ∂_tSSTs, and it includes an additional physics-informed term (Loss_phy in Eq. 2). Locally, this is equivalent to imposing the evolution of the ocean surface temperature according to the horizontal geostrophic advection on the super-resolved fields predicted by the CNN, based on the following equation:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial \mathrm{SST}}{{\partial }_{\mathrm{t}}}}+u_{\mathrm{g}}{\displaystyle \frac{\partial \mathrm{SST}}{\partial x}}+v_{\mathrm{g}}{\displaystyle \frac{\partial \mathrm{SST}}{\partial y}}=\mathrm{0},\end{array}$$

where

• u_g and v_g indicate the zonal and meridional surface geostrophic currents, respectively, and

• x, y, and t stand for zonal, meridional, and temporal coordinates, respectively.

The physics-informed LF is applied to the super-resolved fields predicted by the CNN. The coefficients α, β, γ, and δ have been determined via preliminary CNN training over three epochs, considering separately the four terms of the LF appearing in Eq. (2). In particular, we firstly impose α=1 and perform a first training of the CNN accounting only for SST. We then repeat the same exercise for the remaining terms of the LF separately, and we estimate the mean ratio of the training and validation losses compared to the SST case. We finally assign those ratios to the coefficients β, γ, and δ. The aim of this operation is to allow the CNN to be trained considering the weighted contributions from the four terms appearing in the LF (more details on the LF optimization are provided in Appendix A).

As mentioned in Sect. 2.6.1, in the present study we modified the CNN training strategy, systematically excluding 15 % of tiles towards the end of our available time series (year 2017). This is primarily due to the potential overfitting issue recently claimed by Martin et al. (2023). In order to quantify the CNN performances throughout the year, we computed daily, basin-scale percentages of improvement (again, as in Ciani et al., 2021; Rio and Santoleri, 2018) to quantify the relative improvement of the CNN reconstruction with respect to standard altimetry along the 40 dates of our independent test dataset, as reported by Fig. 8.

Figure 8 Basin-scale percentage of improvement of the CNN reconstruction with respect to standard altimetry. The improvement is evaluated against the independent test dataset. The green and yellow dots stand for Julian day 4 and 324, respectively, and refer to the case studies depicted by Figs. 5 and 9, respectively.

Figure 9 Results of the OSSE on 20 November 2017. (a) Satellite equivalent (SE) current velocity from SE-ADT; (b) current velocity from super-resolved ADT; (c) current velocity from model ADT (ground truth); (d) comparative spectral analysis of the current velocity maps. Results refer to the 2D box depicted in panel (c).

As expected, the percentage of improvement drops after Julian day 300, in correspondence of the time indices for which tiles had been excluded during training and validation. However, the improvement is still exhibiting satisfying basin-scale improvements, around 30 %, indicating that the CNN prediction is still efficient in reconstructing features unseen during the training phase. Performing the spectral analysis, as depicted in Fig. 5, for Julian day 324 (20 November 2017), revealed differences in the CNN performances compared to the temporal window bounded by Julian days 1 and 310. In general, the SE-ADT spectrum is closer to the ground truth characteristics over the entire wavenumber range. However, within the range of 1 to 4 per degree (approximately 100 to 30 km), the SR-ADT spectrum demonstrates better agreement with the ground truth in comparison to the satellite equivalent spectrum. This reaffirms the CNN ability to enhance the characterization of mesoscale dynamics (Fig. 9d). The potential overfitting issue will be further addressed in future studies, in which we plan to extend the time series of the OSSE. Ideally, we aim at exploiting 1 full year for each of the following operations: training, validation, and test.

3.1.1 Reconstruction of ocean surface temperature

Although our primary scope is the computation of super-resolved surface currents from super-resolved ADTs, we briefly illustrate the CNN performances in reconstructing the super-resolved SSTs, which is among the novelties of the new CNN architecture. This is summarized by Fig. 10.

Figure 10 ΔRMSE (in °C) based on the CNN proposed in the present study. Positive values indicate an SST estimate improvement with respect to standard optimal interpolation.

Overall, the CNN is indicating an improvement through the Mediterranean Basin, as illustrated by the dominance of positive ΔRMSE values. The local ΔRMSE maxima reached 0.1 °C, corresponding to ≃ 60 % improvement with respect to standard optimal interpolation mapping, while the mean basin-scale improvement is ≃ 14 %. Unlike ADT, the appearance of local degradation areas looks ubiquitous. These small-scale degradation spots are likely due to the assumptions behind our study. As claimed in Sect. 2.6.1, we build the CNN architecture considering the modulation of SST features by the surface geostrophic advection and, simultaneously, feeding the SST mapping with the dynamical information contained in ADT. This could limit the description of SST signals related to non-geostrophic phenomena or to vertical advection. Such behavior is under investigation and will be the subject of future studies.

3.2 Surface currents predictions from satellite-derived data

The neural network model trained via the OSSE is now tested to predict super-resolved ADTs, using state-of-the-art L4 ADT derived from satellite altimetry and high-resolution L4 satellite SSTs. The satellite-derived input data are produced within the Copernicus Marine Service and cover the 2008 to 2019 time frame (as detailed in Sect. 2.2 and 2.3). The evaluation of CNN performance involves deriving super-resolved (SR) geostrophic currents from SR-ADTs, through the geostrophic approximation. This is achieved applying a finite-central-differences operator to the SR-ADTs (to estimate its spatial gradients) and using drifting buoy measurements as a validation benchmark (detailed in Sect. 2.4). The validation relies on root mean square error inter-comparisons, carried out by interpolating the gridded altimeter-derived and SR currents along the 6-hourly drifter acquisitions (for both components of the circulation) at the drifter time and location (as in Rio and Santoleri, 2018). In order to maximize the number of matchups between gridded and in situ-measured currents, statistics are provided in 2° × 2° boxes, choosing a 12-year-long time series (2008–2019). This guarantees a coverage of in situ measurements as depicted in Fig. 11, with an approximate minimum number of 100 observations per box.

Figure 11 Number of in situ measurements from drifting buoys in the 2008–2019 time frame.

Figure 12 (a) Altimeter-derived surface current velocity; (b) SR current velocity. The fields refer to 19 May 2019.

As addressed by BBN22, provided clear-sky conditions, the overall effect of the CNN for super resolution is a sharpening of the mesoscale SR-ADT gradients compared to standard altimetry processing, with a significant enhancement of the mesoscale eddy activity in the derived SR currents. For clarity, a visual comparison of the altimeter-derived currents (i.e., derived from the ADT computed via standard mapping) and SR current velocity is provided in Fig. 12 and refers to 19 May 2019. Among the features of interest, we notice a significant intensification of the circulation north of the Balearic Islands, along the Algerian currents and in the Strait of Sicily, as well as a sharpening of the mesoscale features populating the southern Tyrrhenian Sea.

The SR currents, compared to standard altimeter-derived mapping, reduce the rms by about 2 to 8 cm s⁻¹ in the majority of the basin for both the zonal and meridional flow (Fig. 13). Degradations with respect to the altimeter (Alti) currents mainly occur in coastal areas and are particularly pronounced in the eastern tip of the Levantine Basin. The basin-scale rms errors of the Alti and SR currents, reported in Table 1, evidence an improvement of the SR currents presented here, with basin-scale reduction of the rms error up to ≃ 1 cm s⁻¹.

Figure 13 Differences of rms errors between the altimeter-derived (Alti) and SR currents: (a) zonal flow and (b) meridional flow. Red areas express an improvement with respect to standard altimetry.

Figure 14 Differences of rms errors between the BBN22 and the current CNN architecture: (a) zonal flow and (b) meridional flow. Red areas express an improvement of the current CNN compared to BBN22.

Table 1 The rms of the altimeter-derived (Alti) and super-resolved (SR) surface currents, computed in 2° × 2° boxes against in situ-measured currents. SR-BBN22 and SR refer to the CNN ADT mapping described in Buongiorno Nardelli et al. (2022) and the one employed in the present study, respectively. U and V stand for zonal and meridional currents, respectively.

Improvements with respect to the BBN22 results are also observed. This is expressed by the overall rms reduction of the new SR currents reported in Table 1 and is supported by the spatial distribution depicted in Fig. 14, in which red areas indicate an improvement of the current CNN formulation compared to BBN22. Such improvements cover 70 % of our study area for the two components of the flow and suggest an overall enhancement of the CNN performances mainly in the central and eastern Mediterranean, with occurrences of weak degradations (≤1 cm s⁻¹) limited to the western basin. Local occurrences of significant performance reduction with respect to BBN22 can be observed in coastal areas and indicate a worsening of the present-day SR current rms by about 5 cm s⁻¹ for the meridional flow. Nonetheless, an overall basin-scale improvement is found.

Specific spectral analyses were carried out to inter-compare the surface kinetic energies derived from standard altimeter (up-sized to the 1/24° grid) and super-resolved ADTs to have insight into the effective spatial resolution of the two datasets. This analysis was performed estimating PSD via FFT analysis over the time range 2008–2019 in two land-free areas of the Mediterranean Basin: (i) one area across the Central and Aegean Basin and (ii) one area across the Algerian Basin and Sardinian Channel, both depicted in Fig. 15c, similarly to the analyses presented for the OSSE.

Both regions are known as dynamically active areas in the Mediterranean Basin (Pujol and Larnicol, 2005). In particular, the KE spectra were inter-compared against the theoretical predictions of two-dimensional turbulence (Vallis, 2006), i.e., the k⁻³ and $k^{-\mathrm{5}/\mathrm{3}}$ slopes.

In the Central and Aegean area (Fig. 15a), the spectral analysis confirms the improvement brought by the CNN reconstruction. The SR and altimeter KE spectra are super-imposed for small wavenumbers, indicating a similar description of the large mesoscale motions, and are aligned with the predictions of energy and enstrophy transfer, following the k⁻³ slope for larger wavenumbers and the $k^{-\mathrm{5}/\mathrm{3}}$ slope for smaller ones.

The improvement of our methodology with respect to standard altimetry processing is evidenced by overall higher PSDs at larger wavenumbers and by a closer alignment with the k⁻³ slope for wavenumbers ≥4 per degree, i.e., scales ≤ 30 km (although not fully recovered through the entire range). This reflects a more efficient representation of the small mesoscale features due to our reconstruction. For the Algerian and Sardinian area, the analyses are depicted in Fig. 15b and led to similar conclusions as for the Central and Aegean area, although both KE estimates (altimeter and SR) show less agreement with the $k^{-\mathrm{5}/\mathrm{3}}$ at wavenumbers ≤3 per degree. For both areas the SR KE spectra drop at scales around 30 km (wavenumbers ≃3 per degree), which can be considered the scale at which mesoscale features are fully characterized.

Figure 15 Comparative spectral analyses of the KE maps derived from standard altimeter (blue) and SR (red) ADTs. Panels (a) and (b) refer to the Central and Aegean (Cen/Aeg) and Algerian and Sardinian (Alg/Sar) areas, respectively, depicted in panel (c).

4 Discussion and conclusions

In recent years, CNNs for super resolution have provided benefits for oceanographic applications, exploiting single-variable and multi-variate approaches (Ducournau and Fablet, 2016; Lima et al., 2017; Buongiorno Nardelli et al., 2022; Wang and Li, 2023; Fanelli et al., 2024a). Here, we present a CNN-based, multivariate approach to jointly reconstruct ocean surface currents and SST from EU Copernicus satellite data. We trained a CNN model from numerical model outputs and then applied the CNN model to satellite-derived data collected over the Mediterranean area.

With respect to previous similar exercises, we relied here on a new OSSE approach/CNN training architecture, also simulating and integrating issues/artifacts found in the present-day L4 SST satellite products (e.g., Ciani et al., 2020), as pointed out in Sect. 2. As such, this study provides a more realistic implementation of the CNN for estimating satellite-derived surface currents from ADT and SST and includes the possibility to simultaneously predict ADT and SST fields (although we mostly focused on the reconstruction of ADTs and geostrophic currents).

The OSSE study also addressed potential overfitting issues of the CNN training strategy previously adopted by BBN22. Here, the neural network validation is performed systematically, excluding tiles observed in a typical late autumn/early winter period, so that CNN prediction capabilities are tested in a worst-case scenario. In other words, the CNN is pushed to predict ocean circulation features in periods of enhanced small mesoscale and submesoscale activity (e.g., Callies et al., 2015) never seen during training. The OSSE results indicate that the CNN improves the characterization of ocean circulation patterns, intensities, and spectral properties even for periods excluded during the training phase, although showing weaker performance with respect to the rest of the year (Figs. 8 and 9). This is likely justifying the overall reduced improvement of the OSSE results compared to the BBN22 formulation (Fig. 7). However, it should also be kept in mind that the degraded synthetic SE-SST employed in our study is surely impacting these statistics as well, while previous tests were based on an excessively optimistic setup.

The application of the dADR-SR model to satellite-derived data is assessed by deriving the ocean surface currents from the CNN-derived ADTs, relying on the geostrophic approximation. Such currents are then inter-compared with standard altimeter-derived currents using (i) in situ measurements as a benchmark and relying on rms errors as validation metrics and (ii) spectral analysis to quantify the gain in the description of the mesoscale features. The CNN-derived currents presented here, compared to standard altimetry products and to previous findings of BBN22, are able to reduce the rms error by approximately 1 cm s⁻¹ for the meridional flow and 0.5 cm s⁻¹ for the zonal one. This is likely due to the adoption of an upgraded CNN, which is now trained at OSSE level, considering realistic observing system characteristics for both ADT and SST, and includes a physics-informed constraint on the LF. Degradations of the CNN-derived currents with respect to the altimetry products mostly occur in near-shore areas, potentially highlighting the limitation of employing L4 SSTs to super-resolve the surface coastal circulation. In coastal areas it is generally harder to interpret or model land contamination effects or to capture small-scale, fast coastal processes, thus affecting the quality of optimally interpolated SST products. Also, the physics-informed constraint on the LF may fail in coastal environments, where small-scale SST features do not necessarily evolve due to geostrophic advection but can result from ageostrophic processes as, e.g., coastal upwelling (particularly pronounced in the Gulf of Lion, Strait of Sicily, and south of Crete). The characterization of the mesoscale features from standard altimetry and SR-ADTs was done via spectral analysis. It confirmed an enhancement of the mesoscale activity obtained via the neural network approach. However, the representation of the 2D geostrophic motion was not fully achieved. The KE spectra obtained from the SR-ADTs depicted by Fig. 15 do not fully recover the theoretical energy KE spectrum at all scales. The use of an optimally interpolated SST is certainly a major limit for our reconstruction methodology. OI SSTs are indeed known to smooth/distort some oceanic features as a drawback of the interpolation algorithm (as pointed out recently by González-Haro et al., 2024). In addition, under prolonged cloud cover (although the neural network is fed with the information on the SST mapping error) the extraction of features becomes questionable, as the SST fields tend to a daily climatology.

Future studies on the application presented here could include the following.

• i.

  Modifying the training strategy for the CNN structure, in which the OSSE could consider ADT- and SST-derived quantities and SST as predictors and surface currents (instead of high-resolution ADTs) as targets. This would have the main advantage to go even beyond the prediction of purely geostrophic motions, as achieved in the present study, where surface currents are derived from the geostrophic approximation equations.

• ii.

  Testing a new prediction of super-resolved surface currents from satellite-derived data (as presented in Sect. 3.2), upon provision of improved L4 SSTs over long time series (i.e., at least decadal). Examples of improvements for SST L4 satellite products were recently presented by Sunder et al. (2020), Jung et al. (2022), and Fanelli et al. (2024a) and mainly propose an enhancement of feature resolution for the reconstructed L4 fields. This is thus expected to improve the description of dynamical features both in SST and, as in the present study, in L4 fields enriched with dynamical information extracted from the gap-free SST maps. The dynamical coherence of such fields could also be quantified by employing recently published metrics based on the multi-fractal theory of turbulence (González-Haro et al., 2024).

• iii.

  Reproducing the OSSE discussed in Sect. 2.6 over time series longer than 1 year. Exploiting 1 full year for training, validation, and test would make our results more robust and discard any CNN overfitting issue. Recent works on the SSH mapping based on neural networks pointed out that training based on increasingly larger time series significantly improves the network prediction performances (Archambault et al., 2024).

Finally, our approach could be adapted to directly learn an end-to-end mapping between present-day operational ADTs obtained from constellations of nadir-looking altimeters and higher-resolution observations from the SWOT wide-swath instrument. Operational L4 ADTs are presently optimally interpolated from along-track observations and achieve effective resolutions O(100 km) at mid-latitudes (Ballarotta et al., 2019; Pujol et al., 2016). With SWOT, SSH can be observed with resolutions down to 20 km (Fu and Ubelmann, 2014; Le Guillou et al., 2021), with the additional advantage of being native 2D over a ≃ 120 km wide swath. As such, potential misrepresentations of the upper-ocean dynamics by the numerical models used for the OSSE would be bypassed by directly training on high-resolution observations. In all cases, a successful implementation of dADR-SR would allow us to carry out a full reprocessing of standard ADT time series back to the start of the satellite radar altimetry era and a potential application to near-real-time processing.

Appendix A: Customization of the loss function for super resolution

A1 Determination of the loss function hyper-parameters

We show the rationale behind the choice of the hyper-parameters for the LF appearing in Eq. (2). This is shown here for the first two terms of Eq. (2), leading to the determination of the α and β parameters. Figure A1 shows the behavior of the training and validation loss curves for the CNN employing the loss function given in Eq. (2) considering only the ADT or SST term, separately. This is shown for the first three epochs, as those curves converge quite quickly afterwards.

Figure A1 Training and validation losses employing a loss function built from ADT (red and orange, respectively) and SST (blue and cyan, respectively) separately.

The ratio between the ADT and SST training and validation loss curves can be quantified by comparing the curves, and from epoch 3 onward it is around 4. As such, considering we impose α=1, β is set to 0.25 in order to equally weight ADT and SST contributions. Such an exercise has been repeated for all the other terms of the loss function, leading to the hyper-parameters appearing in the LF.

A2 Tuning of the physics-informed term

The physics-informed term Loss_Phy in Eq. (2) was shaped according to the following findings. In a first attempt, such a term was expressed by Eq. (A1),

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\mathrm{Loss}}_{\mathrm{Phy}}& ={\left({\displaystyle \frac{\partial \mathrm{SST}}{{\partial }_{\mathrm{t}}t}}\right)}_{\mathrm{pred}}-{\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial {\mathrm{ADT}}_{\mathrm{pred}}}{\partial y}}{\displaystyle \frac{\partial {\mathrm{SST}}_{\mathrm{pred}}}{\partial x}}\\ & +{\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial {\mathrm{ADT}}_{\mathrm{pred}}}{\partial x}}{\displaystyle \frac{\partial {\mathrm{SST}}_{\mathrm{pred}}}{\partial y}}.\end{array}\end{array}$$

The reader is referred to Sect. 3.1 for further details on the equation. Minimizing such Loss_Phy term equals asking the predicted SST and ADT to obey the horizontal geostrophic SST advection. During the determination of the LF hyper-parameters (as explained in Sect. A1) we obtained δ=1.38. After training the neural network and reconstructing the ocean currents with satellite-derived data, the comparison to in situ-measured currents yielded the rms errors detailed in Table A1 (linear Loss_Phy type).

Interestingly, in order to get a further rms error reduction, we had to decrease the δ factor manually, which made the approach rather empirical. We, however, noticed that the Loss_Phy term given by Eq. (A1) was linear, unlike the ADT, SST, and ∂_tSST terms appearing in Eq. (2), for which the minimization is given in a least squared sense. We thus decided to homogenize all the terms contributing to the loss function, modifying Loss_Phy as follows:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathrm{Loss}}_{\mathrm{Phy}}=& \stackrel{\mathrm{‾}}{\left({\left({\displaystyle \frac{\partial \mathrm{SST}}{{\partial }_{\mathrm{t}}}}\right)}_{\mathrm{pred}}-{\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial {\mathrm{ADT}}_{\mathrm{pred}}}{\partial y}}{\displaystyle \frac{\partial {\mathrm{SST}}_{\mathrm{pred}}}{\partial x}}\right.}\\ & \stackrel{\mathrm{‾}}{{\left.+{\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial {\mathrm{ADT}}_{\mathrm{pred}}}{\partial x}}{\displaystyle \frac{\partial {\mathrm{SST}}_{\mathrm{pred}}}{\partial y}}\right)}^{\mathrm{2}}}.\end{array}\end{array}$$

In this way, the weighting factor δ=1.38 automatically adjusted to δ=0.025, and the effect on the reconstructed surface currents was an improvement for both components of the surface motion (i.e., a lower rms error compared to in situ-measured currents), as expressed by Table A1. We thus decided to adopt the quadratic formulation for the Loss_Phy term of the LF.

Table A1 The rms error computed by means of in situ-measured currents for the two formulations of the Loss_Phy term. U and V stand for zonal and meridional currents, respectively.