# Downscaling of ocean fields by fusion of heterogeneous observations using Deep Learning algorithms

Sylvie Thiria, Charles Sorror, Théo Archambault, Anastase Alexandre

Charantonis, Dominique Béréziat, Carlos Mejia, Jean-Marc Molines, Michel

Crépon

► To cite this version:

Sylvie Thiria, Charles Sorror, Théo Archambault, Anastase Alexandre Charantonis, Dominique Béréziat, et al.. Downscaling of ocean fields by fusion of heterogeneous observations using Deep Learning algorithms. Ocean Modelling, 2023, 182, pp.102174. 10.1016/j.ocemod.2023.102174. hal-04006736

HAL Id: hal-04006736

https://hal.science/hal-04006736v1

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

13

14

# Downscaling of ocean fields by fusion of heterogeneous observations using Deep Learning algorithms

Sylvie Thiria¹,², Charles Sorror¹, Theo Archambault¹, ⁴, Anastase Charantonis⁵, Dominique Bereziat⁴,

Carlos Mejia¹, Jean-Marc Molines³, Michel Crépon¹

¹ IPSL-LOCEAN, SU/IRD/CNRS/MNHN, Paris, France

² UVSQ, F-78035, Versailles, France

³ Univ. Grenoble Alpes, CNRS, IRD, Grenoble INP*, IGE, 38000 Grenoble, France

*Institute of Engineering and Management Univ. Grenoble Alpes

⁴ Sorbonne Université, CNRS, LIP6, Paris, France

⁵ ENSIEE, Evry, France

Corresponding author: michel.crepon@locean.ipsl.fr

## Abstract

We present a deep learning method to downscale low-resolution geophysical fields by merging them with high-resolution data. The downscaling was performed using an ensemble of convolutional neural networks (CNNs), whose prediction values are the average values of the outputs of 20 CNNs. Academic experiments were conducted on simulated ocean data in the Gulf Stream region, given by the outputs of the NATL60 model. The CNNs forced with low-resolution (120x120 km) sea surface high (SSH) data and mesoscale resolution (12x12 km) sea surface temperature (SST) data allowed us to obtain mesoscale resolution sea surface currents with good accuracy. Sensitivity experiments have shown that taking SST into account significantly increases the accuracy of the high-resolution velocity retrieval, even when noise is added to the SSH data. The velocity information embedded in the transport equation modeling the SST advection is taken into account by the CNN, which greatly increases the resolution of ocean currents provided by SSH. In the present work, we only consider spatial downscaling by assuming that SSH and SST are daily observations. The method we developed is generic and can be used to improve the resolution of a wide variety of large-scale fields by merging them with high-resolution fields.
29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

# 1 Introduction

For several decades, a wide variety of satellite sensors has allowed us to considerably improve our knowledge of the state of the planet Earth and its potential evolution. These sensors observe a multitude of geophysical parameters with various sampling strategies, both in space and in time. The objective of this study is to show how the combination of several parameters can reinforce each other to obtain geophysical fields observed with better resolution and accuracy. This problem is encountered in many aspects of geophysical sciences such as atmospheric sciences, oceanography, hydrology, fluid mechanics, geology, astrophysics. We have therefore decided to conduct an academic study that should benefit a wide range of observations of geophysical parameters. In this context, we have chosen to focus our efforts on the algorithmic aspect of downscaling these observations by combining them together.

In the present paper, we focused our research on the observation of ocean currents by satellite. The ocean has been observed for several decades by numerous satellites equipped with various sensors that have provided global coverage of dynamic parameters such as ocean circulation with sea surface topography (hereafter SSH ) from the Topex, Poseidon and then Jason altimeters, sea surface temperature (hereafter SST ) with the high-resolution AVHRR radiometers launched aboard meteorological satellites, ocean biological productivity with multispectral ocean color sensors (SeaWiFS, MERIS, MODIS, VIIRS) and ocean salinity (SMOS and Aquarius). These satellite sensors have contributed to a better understanding of the ocean and the detection of its changes in response to global warming (Cazenave et al., 2018). Ocean currents play an important role in the latitudinal transport of heat, and in advecting and structuring ocean properties such as salinity, nutrients, phytoplankton across a wide range of scales from basin to sub-mesoscale (Sasaki et al., 2014; McWilliams, 2016; Levy et al., 2018).

The estimation of ocean currents at high-resolution from satellite observations (a few kilometers) has been addressed by a large number of scientific studies that have combined SSH and SST observations (Rio et al, 2016; Rio and Santoreli, 2018, Nardelli et al, 2022) or that processed SST image sequences
(Isern-Fontanet et al, 2006; Bowen et al., 2002; Vigan et al., 2000; Kelly, 1989; Emery et al., 1986). Other studies have used intrinsic physical properties of ocean surface layers such as conservation of potential vorticity (Ubelmann et al, 2015). Recently, George et al. (2021) and Manucharyan et al. (2021) have shown that SSH interpolation can be improved by theoretical considerations based on QG turbulence and have given leads based on rational physical arguments to improve the reconstruction of mesoscale SSH from SST data.

We propose here a deep learning method (DL) to obtain high-resolution ocean currents by combining low-resolution SSH and high-resolution SST. The problem is very complicated as SSH are measured along satellite tracks which are separated of about one hundred kilometers with a repeat sampling period of about 10 days (Ballarota et al., 2019) while SST are measured as daily two-dimension images at a resolution of 4 km. The association of these two observations of different resolution implies the solving of a difficult inverse problem with time and spatial constraints as accounted in 4D Var data assimilation. In the present paper, we oversimplified the observing system by focusing our study on spatial downscaling alone, assuming that SSH and SST measurements are gridded daily observations of different resolution.

Due to the lack of high-resolution ocean parameter observations and to simplify the data handling at various resolutions, we used simulated data fields. SSH and SST are the outputs of a very high-resolution ocean numerical model. SSH and SST are therefore linked to each other by well-known physical laws that are taken into account by the ocean model and that will be 'learned' by our method. We chose to conduct our study in the Gulf stream region (35°N, 55°W) where ocean currents and associated mesoscale dynamical features are important. We used the outputs given by the high-resolution ocean model NATL60 with a horizontal resolution of 1/60° in latitude. These simulated observations provide SSH, SST, and currents at very high resolution down to 2km, which is much better than operational products derived from observations such as those provided by Aviso (https://www.aviso.altimetry.fr/en/data.html).

Due to the high-resolution data fields and the large area under study, we decided to use deep learning algorithms, which are methods suitable for solving problems involving very large data sets. Recent developments in artificial intelligence, combined with the availability of large datasets such as those provided by high-resolution climate simulations, and increased computational power, have led to the
emergence of information processing capabilities to 'learn' the hidden structures of geophysical phenomena. Over the past decade, statistical and machine learning methods have greatly facilitated the study of these multidimensional data sets. As evidenced by the growing number of publications since 2018, convolutional neural networks (CNNs) which are extensions of multilayer perceptrons (Le Cun et al. 1990, 2015; Simard et al., 2003; Goodfellow et al., 2016), have proved their ability to exploit spatial and temporal data structures, find patterns, and efficiently fuse heterogeneous information sources (Reichstein et al., 2019; Rolnick et al., 2019). In addition, image analysis of complex geophysical fields using neural networks has been a very productive line of research over the past decade and has led to relevant environmental applications (Egmont-Petersen, 2002.; He et al., 2016; Wang et al., 2019; Martinez et al., 2020).

Specific neural architectures have recently been developed, which allow us to obtain numerical schemes that accurately solve partial differential equations (Ruthotto and Haber, 2019; Rousseau et al., 2019). These neural architectures can be used effectively to represent dynamic phenomena that are usually represented by numerical models. The challenge is therefore to introduce prior knowledge provided by numerical and/or statistical analyses into these architectures to improve their performance when solving specific environmental problems (Ducourneau et al., 2016; Fablet et al., 2018; de Bézenac et al., 2019). New tools are made available to the scientific community to help them implement adequate architectures to carry out their research: as an example, recent papers propose automatic differentiation tools embedded in deep learning frameworks to introduce end-to-end neural network architectures for data assimilation (Pannekoucke and Fablet, 2020).

In addition, several studies using machine learning (ML hereafter) methods have been devoted to modeling ocean surface variables with success, showing the potential power of these methods. Examples include sea surface temperature reconstruction from satellite observations using a convolutional neural network (Barth et al., 2020), interpolation of sea level anomalies using data driven methods (Lguensat et al., 2019) and CNN (Manucharyan et al., 2021; George et al., 2021), ocean ML simulations for climate modeling (Sonnewald and Lguensat, 2021).

This paper aims at showing the ability of ML to downscale SSH and associated currents by studying an academic case. Analysis of the geophysical properties of the problem, ML learning and estimation of the physical coherence of the reconstructed fields are performed using a high-resolution numerical
model output. The paper is organized as follows: in section 1, we describe the problem. In section 2, we present the data. The ML algorithm is described in section 3. In Section 4, we present the performance of our method and analyze the results. In Section 5 we discuss the results and conclude in section 6.

## 2 Data

The objective of the present study is to develop a generic machine learning prototype to perform geophysical field downscaling. To demonstrate the feasibility of such a study, we used academic datasets with physical characteristics similar to those provided by ocean satellite sensors. The simulated SSH and SST data used to construct the deep learning algorithms allow us to obtain well-documented ocean datasets with no missing data and control for the different resolutions we wish to explore. We used the SST, SSH, and horizontal velocity field (U, V) provided by the NATL60 model, which is a very high-resolution model of the North Atlantic Ocean extending from 26°N to 66°N (Figure 1). The NATL60 ocean model is based on the NEMO 3.6 code. In the version we used, the tide was not included. The horizontal grid is extracted from the global tri-polar grid of the ORCA series (Madec and Imbard, 1996) with a resolution of 1/60° (at the Equator).

131

<figure><img src="image_2.png" /></figure>

132
<caption>**Figure 1.** *The domain under study (Red Rectangle) and the NATL60 SSH for October 1, 2012, in background.*</caption>

The horizontal grid spacing ranges from 1.6 km at 26°N to 0.9 km at 65°N. The grid has been designed so that the model explicitly simulates the scales of motions that will be observed by the SWOT altimetric mission. In practice, the model's effective resolution is about 10–15 km in wavelength. The resolution increases towards the poles and is about 1x1km at the latitude of the Gulf stream. The vertical grid uses 300 levels (with a thickness of 1m near the surface, up to 50m in the deepest part), to well reproduce the vertical velocities, which condition the vertical exchanges and the interactions with the bottom topography. We used the atmospheric forcing DFS5.2 (Drakkar Forcing Set 5.2, in (Dussin et al., 2016)). Initial conditions were taken from the 1/12° MERCATOR ocean reanalysis. The output of NATL60 was analyzed in Ajayi (2020). In the present research, we studied the Gulf-stream region (26°N, 45°N; 40°W, 65°W, Figure 1) which is the most energetic region of the North Atlantic Ocean. This frontier simulation required 16McpuH, on 14000 cores, and was performed during the GENCI 2016 grand challenge.

Because we focus our study on surface currents recovered from SSH, we bypassed the Ekman surface drift by considering the horizontal current at a depth of 15 meters. We processed 366 daily images of 4 variables (SSH, SST, U, V) covering the period from October 1, 2012 to October 1, 2013 for training (dataset NATL-12 below) and four months with atmospheric forcings from the year 2008, one for validation (December, NATL-08V below) and three for testing (March, June, September, designated NATL-08T below).

Every image is composed of 1296x1377 pixels for each variable, a pixel being a square box of approximately 1.5x1.5km. We denote this dataset R01. In Figure 2, we show the mean local kinetic energy and the mean local enstrophy (see Appendix A for their definition) associated with every pixel over the 366 images of the period studied. We notice that the area can be split into two areas, the one north of 35°N, which includes the Gulf Stream, is very energetic (86% of the total energy); the second one south of 35°N, is less energetic (15% of the total kinetic energy) and corresponds to a recirculation area populated with small eddies. Characteristic statistical values of the two areas are shown in Table 1.
We note that the average current in the northern region has a fairly large eastern component (3.8cms⁻¹) corresponding to the Gulf Stream, while the other current components are far less important. Enstrophy also is much more intense in the North area than in the South one (South of 35°). This has incited us to estimate the performance of the method in the two areas separately.

Operational satellite altimeters such as Topex, Poseidon and Jason measure SSH along tracks whose spacing at the equator is of the order of 315km (Ballarota et al., 2019; Dufau et al., 2016) on a ten-day period. Ocean currents are then retrieved by geostrophy from the SSH gradients. As most altimeter observations are a combination of measurements done by different altimeters and as the studied area is situated at mid latitude (35°N), the altimeter data are available at a resolution of the order of one hundred kilometers on a ten-day period (Ballarotta et al., 2019, Figures 2 and 7). Besides, sea surface temperature can be obtained at 1.1km and 4.4km resolution and at a daily frequency from AVHRR instruments aboard NOAA polar-orbiting satellites.

174

175

<figure><img src="image_3.png" /><figcaption>**Figure 2.** NATL-12 data set, Daily mean Energy and Enstrophy over the studied period (11, 1, 2012 -11, 1, 2013).</figcaption></figure>

176

177

178

179

180

181

The combination of SSH and SST involves consideration of two constraints: a spatial and a temporal constraint, both of which present intrinsic modeling challenges. In this paper, we focus our study on spatial downscaling only, assuming that SSH and SST measurements are made at a daily frequency.
Therefore, we simulated satellite observations at different spatial resolutions from the NATL60 model outputs. We defined five distinct resolution levels for the four variables in multiples of 3X3 and formed 5 datasets (R01, R03, R09, R27, and R81) of 366 daily images each. Each level was calculated from the previous resolution, starting with the highest resolution R01 and ending with the lowest R81. At each step, the images were degraded by averaging the pixels in a 3x3 window. Table 2 shows the number of pixels for the different resolutions of each image and the corresponding size of a pixel.

These considerations have incited us to combine simulated gridded SSH observations at a resolution 120x122km (R81 hereinafter) with SST measurements at 13x14km resolution (R09 hereinafter) and to

<table><caption>**Table 1.** Characteristic statistical values of the different variables in the North and South areas of the NATL-12 dataset.</caption><thead><tr><th>Variables</th><th>Minimum</th><th>Maximum</th><th>Mean</th><th>Std</th><th>units</th></tr></thead><tbody><tr><td>SSH North</td><td>-0.67</td><td>1.36</td><td>0.36</td><td>0.30</td><td rowspan="2">m</td></tr><tr><td>SSH South</td><td>0.34</td><td>0.70</td><td>0.51</td><td>0.06</td></tr><tr><td>SST North</td><td>-2.26</td><td>28.92</td><td>20.80</td><td>4.36</td><td rowspan="2">°C</td></tr><tr><td>SST South</td><td>21.11</td><td>29.03</td><td>24.84</td><td>1.98</td></tr><tr><td>U North</td><td>-2.23</td><td>2.60</td><td>0.04</td><td>0.30</td><td rowspan="2">ms<sup>-1</sup></td></tr><tr><td>U South</td><td>-0.47</td><td>0.68</td><td>-0.01</td><td>0.08</td></tr><tr><td>V North</td><td>-2.50</td><td>2.34</td><td>-0.01</td><td>0.27</td><td rowspan="2">ms<sup>-1</sup></td></tr><tr><td>V South</td><td>-0.34</td><td>0.44</td><td>-0.00</td><td>0.05</td></tr></tbody></table>

estimate the current at the R09 resolution for which the geostrophic approximation is still valid. In the present work, we have undertaken to statistically learn the existing relationship between ocean
currents and SSTs, with dedicated machine learning algorithms, in order to improve the ocean current retrieved from SSH.

<table><caption>Table 2. Characteristic of the different resolution scales (the grid center is close to 35°N, 53°W)</caption><thead><tr><th>Dataset</th><th>Nb of pixels</th><th>Resolution in 100<sup>th</sup> degree</th><th>Resolution in km at grid center(*)</th></tr></thead><tbody><tr><td>R01</td><td>1296x1377</td><td>1.3x1.7</td><td>1.5x1.5</td></tr><tr><td>R03</td><td>432x459</td><td>4x5</td><td>4.5x4.5</td></tr><tr><td>R09</td><td>144x153</td><td>12x15</td><td>13x14</td></tr><tr><td>R27</td><td>48x51</td><td>36x45</td><td>40x41</td></tr><tr><td>R81</td><td>16x17</td><td>108x135</td><td>120x122</td></tr></tbody></table>

Figure 3 presents, for October 1, 2012, the SST and SSH simulated by the NATL60 model at the 5 different resolutions. This figure clearly shows that the degradation of the geophysical information is more and more visible as the resolution decreases. SSH can be considered as a stream function for the ocean currents with a good approximation for large scale motions (i.e. the ocean currents flow following the isolines of the stream function). We also note a similarity between SSH and SST patterns, implying that the SST field contains information on ocean currents. In the following, we label every variable by its symbol (SSH, SST, U, V) followed by its resolution (R81, R27, R09, R03, R01), the two symbols being separated by an underscore (_); as an example, SSH_R81 stands for the SSH image at the R81 resolution.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Images of the 4 variables (SSH, SST, U, V in columns) for October-1-2012. The rows correspond to the 5 resolutions, from high resolution (top) to low resolution (bottom). In parentheses we show the pixel resolution in km.</figcaption></figure>

Figure 4 presents the distribution of the four variables in the 3 datasets. The distribution of every variable in each data is similar except for the SST in the validation set (middle row). This figure allows us to estimate the range of variability of these variables, which is an important information to estimate the performances of our DL downscaling algorithm.
224

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distributions of the four variables (SSH, SST U, V in columns) for the three datasets (Training, Validation and Test, respectively in 1st, 2nd and 3rd row).</figcaption></figure>

225

226

227

228

Let us now introduce the neural methodology "RESAC" (REsolution by Stages of Altimetry and Currents) we have developed to obtain high resolution SSH and current (U, V) from low resolution SSH and high-resolution SST fields. We trained RESAC with the daily images in order to have a large variety of oceanic situations and to test how they are estimated.

229

230

231

232

233

234

235

236
## 3 The RESAC method

### 3.1 The principle of the method

From a physical point of view, the RESAC method follows the approach of (Rio et al., 2016), who showed the feasibility of extracting current information from SST images by inverting the heat conservation equation for retrieving the velocity field from simulated data. This method was applied with success on real satellite observations by (Rio and Santorelli, 2020) for improving altimeter geostrophic velocities. In the present research, we aim at improving the resolution of the ocean current field given by low resolution satellite altimeter observations (SSH) using high resolution satellite SST without applying the geostrophic approximation.

Indeed, the NEMO-NATL60 equations model the ocean surface parameters such as SST and SSH with a very good approximation, taking into account the well-defined relationship existing between them. The ocean is described by the primitive equations, which are an approximation of the Navier-Stokes equations. The surface layers are modelled by equations close to the shallow-water equation, which is of the form for momentum:

$$ \frac{D}{Dt} \mathbf{U} - \mathbf{f} \times \mathbf{U} = -g \nabla \mathbf{SSH} + \mathbf{D}^{U} + \mathbf{F}^{U} \quad (1) $$

along with the SST advection equation:

$$ \frac{\partial}{\partial t} \mathbf{SST} + \nabla (\mathbf{SST} \cdot \mathbf{U}) = \mathbf{D}^{\text{SST}} + \mathbf{F}^{\text{SST}} \quad (2) $$

where $U$ is the velocity vector, $f$ represents the Coriolis parameter vector directed upward, $D^U$ and $D^{\text{SST}}$ the parameterization of small-scale physics for momentum and temperature equations respectively, and $F^U$ and $F^{\text{SST}}$ the surface and deep ocean forcing terms of the ocean surface layer.

The interactions between the two variables (SSH and SST) are facilitated when they have the same resolution. But when they have different resolutions (low SSH, high SST resolution as in actual satellite observations), the relationship linking these two variables is much more complicated and can be
biased due to the fact that, in the Fourier space, interactions preferentially occur between close
Fourier components (Merilees and Warn, 1975; Maltrud and Vallis, 1993). The central idea is the
locality of the interactions as proposed by Kolmogorov in describing the turbulence of incompressible
fluids: only vortices of the same order of magnitude strongly interact with each other (Zakharov et al.,
2012). Besides at low resolution, Eq. (1) is approximated by geostrophy (terms 2 and 3) while at high
resolution, term (1), which has a non-time dependent component, must be added (See Appendix D for
an explanation). These considerations have led us to train a CNN architecture consisting of several
connected modules such that every module downscales the SSH field to a resolution close to that of
the SST field. We have finally obtained a SSH field and associated velocities at the high SST field
resolution. This modular CNN constitutes the framework of the RESAC model built to retrieve high-
resolution ocean currents from low-resolution SSH and high-resolution SST.
This procedure aims at learning the underlying physical laws embedded in the observations (here, the
NATL60 outputs). This study can be linked to the "physics-informed field" although our approach has
many differences. We seek an accurate approximation of the various differential equations used to
determine the components *U* and *V* of the current *U* (Eqs. (1) and (2)). The ML algorithm models the
combination of equations embedded in the NATL60 images of SSH, SST, *U* and *V*. These images take
into account various phenomena, such as momentum and SST advection, interaction with underlying
layers (see Appendix D), whose relationships are guided by the introduction of the existing coherence
between the different resolutions, as shown in Figure 5. In recent "physics-informed" studies,
dedicated numerical schemes of simulations are approximated using a neural approach (Huang et al,
2021; Sonnewald et al, 2021).

286

287

**3.2 The RESAC architeture**

Since the NATL60 ocean model represents the links between SSH and SST at every resolution with a
good accuracy, we trained the RESAC CNN on the NATL-12 dataset. In the training phase, the RESAC
strategy was to conduct the downscaling gradually, stage by stage, to avoid spurious interactions
between the Fourier components of the SSH and SST fields that are far from each other as mentioned
above. RESAC therefore consists of 3 stages, each being composed of a specific CNN providing
intermediate or final estimates of SSH, U and V. The flow chart of the RESAC training is presented in
Figure 5. The progression through the different stages requires estimating the intermediate resolutions until the final downscaling outputs are obtained. This is done by using two intermediate cost functions (J1 and J2 in Figure 5) that control the intermediate SSH resolutions (SSH_R27 and SSH_R09). The final cost functions (J3u and J3v in Figure 5) control the U and V final outputs (U_R09 and V_R09, Figure 5). The minimization considers the sum (J) of the four intermediate cost functions (J1, J2, J3u, J3v). Each one minimizes the sum of the mean square error (MSE) between the value estimated by the CNN ($y_i^{Est}$) and the NATL60 ($y_i^{Natl}$) value, the sum being computed over all the available pixels i of all images. We denote CNN_R81toR27 (resp. CNN_R27toR09) the architecture that allows us to model the downscale procedure from SSH_R81 to SSH_R27 using SST_R027 (resp. SSH_R27 to SSH_R09 using SST_R09), and we denote CNN_UV_R09 the CNN architecture that models SSH_R09 and SST_R09 to U_R09 and V_R09 (Figure 5). The first two stages model the flow of information provided by the prognostic variables, the third stage models the physics that links SSH and SST with velocities in the NATL60 ocean model.

Each CNN takes into account the existing correlations between a given resolution of SSH and the next (SSH, SST) one (stage 1 and stage 2), and between the SSH and the velocities at the same resolution. The learning process is achieved in the following manner. Once the architecture of the CNNs has been figured out, the RESAC learning phase proceeds from scratch by randomly initializing the weights of the CNNs. The weights of these CNNs are estimated at once by simultaneously minimizing J, which controls the passage from one stage to the following one. The exact specifications of the CNNs, such as the number of convolutional layers and feature maps or the size of the filters at each stage, are given in Appendix B. They were heuristically optimized at each stage, the CNN specifications being separately determined using their dedicated inputs, outputs, cost functions and learning set. In Figure 5, UP3x3 denotes the up-sampling procedure which increases the size of the image by a factor of 3x3, in order to get the higher resolution corresponding to the next stage. In RESAC, we performed a bilinear up-sampling method with a 3x3 up-sampling ratio.

We note that the CNN_R81toR27 architecture (first stage in Figure 5) is constituted of a larger number of weights than the other two (respectively the CNN_R27toR09 associated with the second stage and the CNN_UV_R09 associated with the third stage, see Figure 5). This can be explained as follows. The
322

323

324

325

326

SSH and SST NATL60 fields are weakly correlated. This correlation increases at intermediate resolution when these two fields have interacted through the CNN_R81toR27 (Figure 5), which therefore must perform a more difficult task than the other two CNNs and consequently requires a larger number of weights.

<figure><img src="image_6.png" /></figure>

327

328

329

330

331

332

333

334

335

336

337

338

339

340

<caption>**Figure 5.** Flow diagram of the RESAC training phase: the inputs of the RESAC model are in green, the observations used for controlling the RESAC weights are in blue, the double arrows indicating the four cost functions J1, J2, J3u, J3v. UP3x3 stands for the up-sampling procedure that increases the input size by a factor 3x3 in order to reach the higher resolution of the next image. During the learning phase, $\hat{H}$ is the estimated SSH at R27 and R09 respectively, $\hat{U}$ and $\hat{V}$ are the estimated ocean velocity values. The weights of all the CNNs are estimated at once by simultaneously minimizing the four cost functions. In the operational phase, RESAC is forced with SSH_R81, SST_R27 and SST_R09 (in the green frames).</caption>

The full learning process (Learning, validation, cost functions, learning steps, hyper-parameter assessments...) is detailed in Appendix B. Training, validation and test sets are completely decorrelated (see Section 2). The number of convolution layers, the size and number of filters, the activation and loss functions were chosen based on previous experiments. The learning rate schedulers and the number of epochs were empirically tuned on the validation set (see Appendix B).

341

342

343

344

At the end of the training, when all the weights of RESAC have been settled, RESAC can be used to retrieve the high-resolution ocean currents. In the retrieval phase, only the variables present in the green frames in Figure 5 (SSH_R81, SST_R27, SST_R09) are inputted in RESAC so as to estimate the current (U_R09 and V_R09). RESAC allows us to downscale SSH_R81 to SSH_R09, U_R09 and V_R09
345 gradually, the process being controlled by intermediate and high-resolution SST data (SST_R27 and SST_R09).
346

347
348 Figure 6 shows the inputs of the RESAC model and its final outputs in the operational phase.
349

<figure><img src="image_7.png" /><figcaption>350
351

**Figure 6. Retrieval phase: Inputs and outputs of RESAC** *Inputs:* SSH_R81 is the NATL-08T SSH field at low resolution (120x122km), SST_R27 and SST_R09 are the SST fields at mesoscale resolution obtained by averaging the SST_R01 (1.5km) NATL-08T observations (see Section 2). *Outputs:* SSH_R09, and (U, V)_R09.</figcaption></figure>

352
353
354
355
356
357
358
359
360

In section 4, we studied different RESAC models with different inputs and architectures. For each RESAC model, we trained 20 networks with different weights for initialization, using an adapted learning rate planner (see Appendix B). The final prediction values are the average values of the predictions of the 20 networks. This procedure allowed us to increase the performance of the final CNN.

361
362

We also tried to train a simpler architecture linking the inputs (SSH_R81, SST_R27, SST_R09) to SSH_R09 directly. But we were not able to perform the training phase correctly. The minimization was
363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

still stuck in a flat local minimum. This proves the need to process the downscaling gradually by
introducing the (SSH, SST) relationships at different resolutions.

# 4. RESAC performances

## 4.1 Qualitative comparison of the different models

To show the effectiveness of the RESAC method and the contribution of SST, we trained two separate RESAC models, called RESAC-ssh-sst and RESAC-ssh, both of which provide downscaled ocean parameters. The inputs to RESAC-ssh-sst are SSH at resolution R81 and SST observations at resolutions R27 and R09 (see Figure 5), while the inputs to RESAC-ssh are SSH_R81 data only (Appendix B, Figure B2). The architectures of these two RESAC models are very similar, but RESAC-ssh is somewhat different, due to the removal of the SST. Since the SST inputs (green in Figure 5) were not used, RESAC-ssh has a smaller number of weights (400,301) than RESAC-ssh-sst (402,797 weights).

In addition, we compared our results to those obtained by a simpler approach, a well-known unsupervised image processing algorithm (https://opg.optica.org/oe/fulltext.cfm?uri=oe-19-27-26161&id=225713). This algorithm (hereafter Bi-Cubic+geostro) consists of performing SSH interpolations with two-dimensional bicubic functions. We performed the interpolation of SSH_R81 data at R09 resolution with the bicubic algorithm and then estimated the R09 current components (U, V) using the geostrophic approximation. The SST is not used.

The training and testing of both RESAC models were performed in the same manner using the same data and conditions. As presented in Section 2, training was performed on the 366 days of the NATL-12 dataset, validation on the NATL-08V dataset (December 2008), and testing on the NATL-08T dataset (March, June, September 2008). All the statistics presented below have been computed on NATL-08T.
387

388

389

390

First, we show the downscaling of SSH and velocity values given by RESAC-ssh-sst, RESAC-ssh, and Bi-cubic+geostro (Figure 7, second row) for March 15, 2008. We also compare these SSHs with the SSH values from NATL-08T by taking their difference (third row). The fourth row shows the energy obtained from the three models.

391

<figure><img src="image_8.png" /></figure>

392
Figure 7. Comparison of the outputs of the three models (the two RESACs and the Bi-cubic+geostro) for March 15 2008. First row, true NATL-08T values: left cartoon, R81 SSH (forcing); middle, R09 SSH and velocities; right, R09 kinetic energy. Row 2, 3 and 4 outputs given by RESAC-ssh-sst (left), Resac-ssh (middle) and Bi-cubic+geostro (right) for SSH and velocities (row 2), the difference between the retrieved and the observed SSH (row 3), kinetic energy (row 4). The color bars shown for each variable are shared by all figures for the same variable, regardless of resolution.

At first glance, RESAC-ssh-sst reproduces the NATL-08T values better than the other two models. We note that the area situated North of 35°N, which is much more energetic than the South one as it is populated with a large number of energetic mesoscale eddies, is better reproduced than the South one. Clearly RESAC-ssh-sst outperforms RESAC-ssh and the worst results are provided by Bi-cubic+geostro, proving the interest of using physical information such as the SST during the downscaling procedure. We note that the geographic positions of the major oceanic structures are correctly retrieved by RESAC-ssh-sst as well as the velocity amplitudes and directions.

## 4.2 Quantitative comparison of the different models

We compared the outputs of RESAC-ssh-sst, RESAC-ssh, and Bi-cubic+geostro to the NATL-08T dataset using several statistical estimators. First, we estimated the consistency of the variable reconstruction by computing a daily RMSE for each variable and for each test image according to Eq. (3) :

$$RMSE(T_{est}) = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (y_{i}^{NATL} - y_{i}^{est})^2} \quad (3)$$

where N is the number of pixels of a R09 image (N=144x153), $y_{i}^{est}$ is the value estimated by RESAC at pixel i, $y_{i}^{Natl}$ is the NATL60 value at pixel i.

We averaged the daily RMSE of the test images. We therefore obtained a mean RMSE (denoted MRMSE in the following). Table 3 shows the MRMSE computed between the NATL-08T variables, and the RESAC-ssh-sst (first row in Table 3), RESAC-ssh (second row) and Bi-cubic+geostro (third row) outputs for SSH, U and V at the R09 resolution. It clearly appears that RESAC-ssh-sst performs better
than RESAC-ssh. The MRMSE between the NATL-08T and the RESAC-ssh-sst SSH is much lower than the NATL-08T SSH distribution whose most frequent values are of the order of 0.50 m (Figure. 4), which suggests that RESAC-ssh-sst is efficient in downscaling the SSH and the associated ocean currents.

<table><caption>**Table 3.** Performances on NATL-08T. Mean daily RMSE (MRMSE) computed by RESAC-ssh-sst (first row), by RESAC-ssh (second row) and Bi-cubic+geostrophy approximation (third row) for the three downscaled variables SSH, U and V on NATL-08T at the R09 resolution.</caption><thead><tr><th></th><th>SSH (m) - MRMSE</th><th>U (ms<sup>-1</sup>) - MRMSE</th><th>V (ms<sup>-1</sup>) - MRMSE</th></tr></thead><tbody><tr><td>RESAC-ssh-sst</td><td>0.0418</td><td>0.1107</td><td>0.1088</td></tr><tr><td>RESAC-ssh</td><td>0.0611</td><td>0.1612</td><td>0.1642</td></tr><tr><td>Bi-cubic+geostro</td><td>0.0693</td><td>0.1851</td><td>0.1797</td></tr></tbody></table>

The use of stage 3 in the RESAC training implies knowledge of the high-resolution ocean current field, which can only be given by a numerical model. We also calculated ocean currents by applying geostrophy to SSH R09 given by the two RESAC models. These outputs are labeled RESAC-ssh-sst+geostro and RESAC-ssh+geostro respectively.

In Figure 8, we present the time series of daily RMSEs given by RESAC-ssh-sst, RESAC-ssh, and Bi-Cubic+geostro for SSH (top panel) and the daily RMSEs given by these statistical models for (U, V), (middle and bottom panels). The best results for U, V are given by RESAC-ssh-sst, followed by RESAC-ssh-sst+geostro, RESAC-ssh, RESAC-ssh+geostro and the worst by Bi-cubic+geostro. The information embedded in the SST advection equation (Eq. 2) improves the R09 ocean currents given by SSH alone. Analysis of these RMSEs allows us to see the impact of the SST on the retrieval of ocean currents.
The fact that the ocean current RMSEs of RESAC-ssh and RESAC-ssh+geostro are nearly similar implies that the CNN_UV_R09 (Figure B2) included in RESAC-ssh is a function close to the geostrophy that would be difficult to improve without the use of additional parameters. The CNN_UV_R09 (Figure 5) included in RESAC-ssh-sst improves the retrieval of ocean currents over geostrophy and thus mitigates the deficiencies of the latter by introducing SST as an additional variable.

<figure><img src="image_9.png" /></figure>
**Figure 8.** SSH, U and V time series of the daily R09 RMSEs between the estimated field and the NATL-08T values (March, June and September 2008). SSH RMSEs estimated with RESAC-ssh-sst, RESAC-ssh and Bi-cubic+geostro are presented in the top panel, U and V RMSE in the middle and bottom panels. We also added the daily RMSEs for two other models, RESAC-ssh-sst+geostro and RESAC-ssh+geostro, which estimate ocean currents computed by geostrophy from the SSH provided by the two first RESAC models.

448

449

Let us now try to understand how SST acts in RESAC-ssh-sst to obtain the best SSH and ocean currents as shown in Figure 8. Schematically, the performance improvement has two main sources. First, SST helps to retrieve a better SSH, which leads to a better geostrophic approximation as we can see by comparing the ocean current RMSEs of RESAC-ssh+geostro and RESAC-ssh-sst +geostro (0.02 $ms^{-1}$ improvement on the MRMSE, Table 3). Second, SST helps the CNN_UV_R09 of RESAC-ssh-sst (Figure 5) to retrieve currents with additional performance with respect to geostrophy.

In the following, we focus our analysis on the performances of the two DL models, namely RESAC-ssh-sst and RESAC-ssh, compared to those of Bi-Cubic+geostro.

We calculated the complex correlation (see Appendix C for a definition) between the velocities given by RESAC and by Bi-cubic+geostro, and with those of NATL-08T (Table 4), in the North and South zones. We note that the correlations given by RESAC-ssh-sst are higher than those given by RESAC-ssh and Bi-cubic+geostro. The difference are of 0.11 and 0.21 respectively, values which are highly significant given

<table><caption>Table 4. Performances on NATL-08T: Mean daily Complex Correlation (CC) and mean angle between the vector modules and for RESAC-ssh-sst, RESAC-ssh and Bi-cubic+geostro for NATL-08T velocity fields at the R09 resolution.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean daily CC</th><th colspan="2">Mean daily angle in degree</th></tr><tr><th>Mean</th><th>STD</th><th>Mean</th><th>STD</th></tr></thead><tbody><tr><td rowspan="3">RESAC-ssh-sst</td><td>Total</td><td>0.92</td><td>0.01</td><td>0.40°</td><td>0.28°</td></tr><tr><td>North</td><td>0.93</td><td>0.01</td><td>0.46°</td><td>0.30°</td></tr><tr><td>South</td><td>0.82</td><td>0.02</td><td>0.98°</td><td>0.79°</td></tr></tbody></table>
<table><tbody><tr><td rowspan="3">RESAC-ssh</td><td>Total</td><td>0.81</td><td>0.03</td><td>0.53°</td><td>0.39°</td></tr><tr><td>North</td><td>0.82</td><td>0.04</td><td>0.66°</td><td>0.42°</td></tr><tr><td>South</td><td>0.73</td><td>0.03</td><td>1.39°</td><td>1.21°</td></tr><tr><td rowspan="3">Bi-cubic+geostro</td><td>Total</td><td>0.77</td><td>0.03</td><td>0.88°</td><td>0.61°</td></tr><tr><td>North</td><td>0.78</td><td>0.03</td><td>1.05°</td><td>0.74°</td></tr><tr><td>South</td><td>0.71</td><td>0.03</td><td>1.52°</td><td>1.27°</td></tr></tbody></table>

the high degree of freedom of these correlations (144x153=22,032 which corresponds to the number of pixels by image). The angle between the two vectors (NATL-08T and algorithm estimated) has no systematic bias, the mean angle between the two vector fields being close to zero (Appendix C). Once again RESAC-ssh-sst performs better in retrieving the velocities than RESAC-ssh and Bi-cubic.
We note that the standard deviation (std) of these correlations are small, which shows the good functioning of the retrieval algorithms.

## 4.3 Case of noisy SSH data

SSH measurements made by altimeters are often characterized by a high noise/signal ratio due to instrumental noise and interpolation of altimeter observations. This interpolation introduces noise since the altimeter observations are made under the satellite track and gridded by combining several tracks over a fairly long time interval depending on the satellite repeat period (10 days for TOPEX-POSEIDON, Ballarotta et al., 2019; Dufau et al., 2016). In the following, we simulated those noises by adding a Gaussian noise to the SSH_R81 of NATL-08T.

<figure><img src="image_10.png" /><figcaption>Figure 9. SSH_R81 (a) without noise, (b) with $\sigma_1=0.05$ m and (c) with $\sigma_2=0.10$ m, on 09/01/2013.</figcaption></figure>
For each pixel of each test day, we added a Gaussian random noise $N(0, \sigma^2)$ to the SSH_R81 images. The SSH signal is now SSHn_R81 = SSH_R81 + $N(0, \sigma^2)$. We tested two different noises. The standard deviation of the first one is $\sigma_1 = 0.05$m, that of the second is $\sigma_2 = 0.10$m. These two values are quite important with respect to the NATL12 SSH distribution whose most frequent values are of the order of 0.50 m (Figure 4, first row). Figure 9 shows the SSH at the R81 resolution without noise and for the two noise values on 09/01/2013. We note the strong alteration of the SSH_R81 patterns for the highest noise level.

The mean daily RMSE (MRMSE) of the two RESAC models and Bi-cubic+geostro SSH and U, V, with respect to NATL08T values are presented in Table 5 for the full area.

<table><caption>**Table 5.** Performances on NATL-08T. Mean daily RMSE (MRMSE) computed by RESAC-ssh-sst, RESAC-ssh, and Bi-cubic+geostrophy approximation for the three downscaled variables: SSH (in m), U and V (in ms⁻¹) on NATL-08T at the R09 resolution without noise and with two noise levels (0.05 and 0.10 m).</caption><thead><tr><th rowspan="2"></th><th colspan="3">RESAC-ssh-sst</th><th colspan="3">RESAC-ssh</th><th colspan="3">Bi-cubic+geostrophy</th></tr><tr><th>σ =0.0</th><th>σ =0.05</th><th>σ =0.10</th><th>σ =0.0</th><th>σ =0.05</th><th>σ =0.10</th><th>σ =0.0</th><th>σ =0.05</th><th>σ =0.10</th></tr></thead><tbody><tr><th>SSH</th><td>0.042</td><td>0.062</td><td>0.099</td><td>0.061</td><td>0.084</td><td>0.126</td><td>0.069</td><td>0.084</td><td>0.111</td></tr><tr><th>U</th><td>0.111</td><td>0.132</td><td>0.180</td><td>0.161</td><td>0.187</td><td>0.241</td><td>0.185</td><td>0.187</td><td>0.227</td></tr></tbody></table>
<table><tr><td>V</td><td>0.109</td><td>0.128</td><td>0.172</td><td>0.164</td><td>0.190</td><td>0.243</td><td>0.180</td><td>0.190</td><td>0.229</td></tr></table>

The MRMSEs increase with the noise level and are smaller when SST is considered, showing the important role of the SST in velocity retrieval. As expected, the higher the noise, the larger the MRMSE.

The complex correlation (CC) of the velocities of both RESAC and Bi-cubic+geostro models with NATL08T velocities is presented in Table 6 for the whole area and the North and South zones. RESAC-ssh and Bi-cubic+geostro show much poorer performance than RESAC-ssh-sst.

<table><caption>**Table 6.** Performances on NATL-08T. Mean daily Complex correlation CC and angle computed by RESAC-ssh-sst, RESAC-ssh, and Bi-cubic+geostro approximation for U and V on NATL-08T at the R09 resolution without noise and with two noise levels (0.05 and 0.10). Angles are given in degrees.</caption><thead><tr><th rowspan="2"></th><th colspan="3">RESAC-ssh-sst</th><th colspan="3">RESAC-ssh</th><th colspan="3">Bi-cubic+geostrophy</th></tr><tr><th>σ=0.0</th><th>σ=0.05</th><th>σ=0.10</th><th>σ=0.0</th><th>σ=0.05</th><th>σ=0.10</th><th>σ=0.0</th><th>σ=0.05</th><th>σ=0.10</th></tr></thead><tbody><tr><th>CC Total</th><td>0.92</td><td>0.88</td><td>0.80</td><td>0.81</td><td>0.73</td><td>0.60</td><td>0.77</td><td>0.71</td><td>0.59</td></tr><tr><th>Angle</th><td>0.40°</td><td>0.41°</td><td>0.47°</td><td>0.53°</td><td>0.57°</td><td>0.70°</td><td>0.88°</td><td>0.88°</td><td>1.00°</td></tr><tr><th>CC North</th><td>0.93</td><td>0.91</td><td>0.87</td><td>0.82</td><td>0.77</td><td>0.68</td><td>0.77</td><td>0.75</td><td>0.68</td></tr></tbody></table>
<table><tr><td>Angle</td><td>0.46°</td><td>0.48°</td><td>0.56°</td><td>0.66°</td><td>0.71°</td><td>0.87°</td><td>1.05°</td><td>1.06°</td><td>1.18°</td></tr><tr><td>CC South</td><td>0.82</td><td>0.68</td><td>0.49</td><td>0.73</td><td>0.53</td><td>0.34</td><td>0.71</td><td>0.53</td><td>0.35</td></tr><tr><td>Angle</td><td>0.98°</td><td>1.06°</td><td>1.47°</td><td>1.39°</td><td>1.67°</td><td>2.38°</td><td>1.52°</td><td>1.68°</td><td>2.82°</td></tr></table>

The contribution of the SST is essential as can be seen on the velocity correlations which are always better for RESAC-ssh-sst with a difference of 0.15 to 0.20 compared to the other models. We note that the downscaling procedure keeps the velocity directions with good accuracy.

We again find that the performances are better in the northern region than in the southern region. This is due to the fact that in the northern region, ocean currents are stronger and associated with large stable ocean structures corresponding to the Gulf stream, while the southern region is more turbulent with smaller currents associated with small ocean structures.

# 5 Discussion

The performance of RESAC-ssh-sst is very satisfactory and superior to that of RESAC-ssh and Bicubic+geostro, as shown by several statistics presented in section 4 (Tables 3 and 4). The daily RMSEs displayed in Figure 8 allow us to analyze the performances of the different models. The good score of RESAC-ssh-sst compared to RESAC-ssh and Bicubic+geostro is due to SST acting in two different ways on the structure of the algorithm. First, the SST input on CNN R81 to R27 and CNN R27 to R09 (Figure 5, steps 1, 2 of the flow diagram) improves SSH retrieval at the R09 resolution as can be seen by comparing the daily SSH RMSEs given by the three downscaling models (top panel of Figure 8). Second, the SST input on CNN_UV_R09 (Figure 5, step 3) provides advection information, which improves the
R09 ocean current retrieval as can be seen from the daily RMSEs of the *U* and *V* current components (Figure 8, middle and bottom panels and Table 4). Moreover, RESAC gives satisfactory performance when *SSH* is noisy. To this end, we conducted two experiments with two different noise levels on the *SSH*, one with a standard deviation (std) of 5cm, the second with a std of 10cm. These values are high compared to the dominant value of the *SSH* probability density function (pdf) which is of 50 cm according to Figure 4. In all cases (without and with noise), the complex correlations between RESAC current vectors and NATL-08T current vectors reported in Table 5 and 6 show that RESAC-ssh-sst is the best estimator of ocean currents at R09 resolution.

All these analyses show that *SST* plays a critical role in the accuracy of ocean current retrieval, especially when *SSH* measurements are noisy. The *SST* which is not considered in RESAC-ssh and Bi-cubic+geostro, adds important information on ocean currents in RESAC-ssh-sst, allowing to improve the resolution of the retrieved *SSH* field as well as that of the associated velocities. The velocity information contained in *SST* is linked to the structure of the *SST* advection equation used in the NATL60 simulation (Eq. (2)), which associates *SST* with velocity. The rationale is that *SST*, which is advected by ocean currents, is a stream function for slowly varying ocean motions (period greater than a few days) which therefore implies that the fluid velocity is almost perpendicular to the *SST* gradient (see Appendix D). In this study, Eq. (2) was implicitly modeled by RESAC so that terms 2, 3, 4 were statistically accounted for. From a theoretical point of view, we can state that RESAC-ssh-sst has approximated the physical laws included in the data fields represented by the PDEs constituting the NATL60 model (see Appendix D).

The performance of RESAC-ssh-sst is satisfactory in the northern area but deteriorates in the southern area, due to the low signal-to-noise ratio in the latter region (Tables 3 and 4 for the NATL-08T test set; Tables 5 and 6 for the noisy observations).

Furthermore, comparison between the retrieval of (*U*, *V*) by RESAC-ssh-sst and RESAC-ssh-sst+geostro (Figure 8) shows that geostrophy computed from *SSH* only is a quite weak estimator of ocean currents that can be improved with the help of additional information (*SST* in the
547

548

present case) and not by algorithmic improvements as shown by the similar *U* and *V* RMSEs provided by RESAC-ssh and Bicubic+geostro (Figure 8, second and third rows).

549

550

551

552

553

554

555

556

557

558

559

560

561

562

A major point of RESAC is to decompose the training phase into three stages as outlined in Figure 5 to facilitate the exchanges of physical properties of the fluid associated with SSH and SST, which preferentially occur in the spectral domain at wave numbers which are close, owing to the locality of the interactions in turbulence (Merilees and Warn, 1975; Maltrud and Vallis, 1993; Zakharov et al., 2012). The training phase uses an intermediate resolution of the SST observations (SST\_R27 which in the present case is simply obtained by filtering SST\_R09) and two intermediate resolutions of the SSH observations (SSH\_R27 and SSH\_R09) allowing to control the intermediate estimations of RESAC (SSH\_R27 and SSH\_R09) by using specific cost functions (see Figure 5 for the flow diagram). In the operational phase, the inputs of RESAC-ssh-sst are the SSH observations (SSH\_R81) at low resolution and the SST observations at intermediate (SST\_R27) and mesoscale (SST\_R09) resolutions. The outputs are the ocean currents at mesoscale resolution (U\_R09, V\_R09), together with the SSH at the same resolution (SSH\_R09). This decomposition in stages is an essential ingredient of the RESAC method. The performance of the algorithms without the stage decomposition was found to be much lower (not shown).

563

564

565

566

567

568

569

570

571

572

573

574

RESAC training was conducted on very high-resolution ocean model output, which allows for high resolution velocities that are not observable by in situ measurements. In a future study, we propose to train a RESAC-SAT dedicated to processing satellite observations, on SSH altimeter tracks and SST fields simulated from an operational ocean model such as Mercator models (https://www.mercator-ocean.eu/en/). As they assimilate several different observed variables, these models are close to the real ocean. It is therefore expected that the observation operator (according to data assimilation jargon) permitting to associate satellite observations with ocean model data will be close to identity as found in similar studies such as the scatterometer wind retrieval (Richaume et al, 2000). These authors calibrated a ERS1 wind retrieval neural network with ERA-interim wind fields. In the construction of RESAC-SAT, it would be necessary to take into account the complex spatio-temporal sampling along the SSH observation tracks associated with the satellite repeat period. Recent studies have proposed new efficient procedures to improve the consideration of SSH sampling; they are based on 4D Var
(Febvre et al, 2021) or on CNN (Nardelli et al, 2022). The architecture of RESAC will be adapted according to the sampling characteristics (spatial and temporal) of the satellite SSH sensors and the chosen procedure. A complementary way to improve the efficiency of such an operational algorithm would be to consider the physical properties of the links between SSH and SST. We could add new inputs to RESAC such as air-sea heat exchange (term 4 in Eq. (2)), which improves the impact of SST on SSH, as has been shown (Rio and Santorelli, 2020, Nardelli et al., 2022) for large-scale ocean motions of the order of atmospheric structures (several hundred kilometers) and should therefore improve SSH at R81 resolution. Small-scale (~10 kilometers) ocean motions are more difficult to consider because they are mainly associated with heat exchanges within the deep ocean which are related to mixed layer dynamics and horizontal eddy fluxes. Investigating the impact of these processes on the resolution of the SSH and associated ocean currents would deserve a specific study.

Besides, SSH resolution in both time and space will be significantly improved by the launch of the Surface Water and Ocean Topography (SWOT) mission, which is a joint NASA/CNES/CSA ASC/UK space agency program, which will facilitate the inclusion of SSH data.

# 6 Conclusion

We have built an efficient machine learning algorithm, the CNN-based RESAC method, to increase the resolution of ocean currents provided by low-resolution SSH (120x122km, SSH_R81) up to mesoscale resolution (13x14km, SSH_R09) by combining them with high resolution SST. We consider an oversimplified case for which both SSH and SST measurements are gridded daily observations. The essence of RESAC is to combine coarse SSH observations with high-resolution SST by introducing intermediate SST resolutions in the training and retrieval phases. To increase the accuracy of RESAC, we have constructed a CNN ensemble, whose prediction values are the average values of the outputs of 20 CNNs with the same architecture and trained on the same datasets but with different initial values. More generally RESAC can be used to merge and downscale all data fields that are observed simultaneously.
The RESAC architecture has shown the interest of decomposing the learning phase into several stages to facilitate the fusion of these data. RESAC accuracy could be improved by using a larger dataset (several years instead of one) for training. Moreover, several additional studies could be undertaken, such as training two different RESAC models, one over the northern area and the other over the southern area, with the expectation that RESAC's performance will improve especially in the southern area where ocean currents are weak and mostly turbulent. Furthermore, as designed, RESAC should be able to operate (at least in theory) in equatorial regions where geostrophy is not applicable.

The present work can be considered as a feasibility study leading to the design of a future operational DL method dedicated to the downscaling of SSH fields. In the real world, the problem is much more complicated. It is necessary to take into account the effect of SSH sampling along satellite tracks and successive interpolations of altimeter data during the repeat period. Errors in the final satellite products are not due to simple smoothing but to the complex spatiotemporal interpolation (Ballarota et al, 2019). This challenging problem implies further developments to obtain an efficient operational DL algorithm, which is beyond the scope of this paper.

The present study acts as a proof of concept regarding the merging of datasets of different resolutions to downscale low-resolution data fields up to high resolution ones. It explicitly relies on the ability of CNNs to summarize the geophysical knowledge linking SST, SSH and velocity fields (U, V) at every pixel resolution from 120 km² to 13 km². Moreover, a challenging study would be to estimate ocean currents at the R01 resolution (1.5x1.5km). This scale corresponds to ocean filaments for which the geostrophic approximation is no longer valid, which makes the velocity downscaling at high-resolution more complex. This would imply adding a new stage in the RESAC algorithm owing to the large amplitude of the downscaling range, and consequently increasing the size of the CNNs and the computational cost of the learning phase. This would require access to very powerful GPUs. These results would have been improved if a longer time series for the training had been at our disposal.

RESAC must be learned on high resolution ocean modes outputs since there is no high-resolution ocean observations. The RESAC outputs could be somewhat different from actual ocean variables but not too much because ocean models are now very efficient as they assimilate a large number of
628 observations. This constraint is inherent to any high-resolution DL downscaling operator of ocean
629 variables.

630 As mentioned in the introduction, Deep Learning (DL) is a very productive research area today. Its use
631 for solving complex problems in the environmental sciences leading to scientific discoveries is rapidly
632 increasing. Geophysical sciences and climate studies must deal with massive amounts of
633 multidimensional data. DL methodologies are well suited to solve the formidable scientific challenge of
634 multi-satellite observations and climate model outputs. High performance computing (HPC) is
635 increasingly being used to meet the exponentially growing demand for DL data processing. Multiple
636 GPUs and high-speed interconnect networks are available to support DL on HPC systems (Han J, 2019).
637 All of this research done in the field of artificial intelligence is now available to the environmental
638 research community and should enable rapid development of geophysical knowledge. This
639 development of DL should give us confidence in the applications of RESAC. Indeed, RESAC learned on
640 simulated data contains a summary of the physics that relates SSH and SST to ocean currents. Given its
641 generic nature, we expect the RESAC methodology to be applicable to a wide variety of downscaling
642 problems.

643

644

645

646

647 Data and Code Availability

648 The datasets we used (NATL-12 and NATL-08 which were extracted from the NATL60 model
649 experiments output over the Gulf stream area, in the North Atlantic Ocean) are already available in the
650 IPSL Thredds catalog at the address : https://doi.org/10.14768/c3c33afe-2a37-42e1-bb29-21428387f658

651 The Python codes of the RESAC models are hosted at GitLab:

652 https://gitlab.in2p3.fr/carlos.mejia/resacpubli_ocemod/-/releases/v1.0r

653

654 Author contributions
S. Thiria designed the research and wrote the manuscript, C. Sroror and T. Archambault developed the RESAC models, and performed the simulations and validated the results. J. M. Molines provided the NATL60 outputs (SSH, SST and ocean velocities), A. Charantonis and D. Bereziat participated in the CNNs development and manuscript writing. C. Mejia contributed to the data management and to the document preparation. M. Crepon discussed the results, wrote the manuscript, and edited it.

660

## Competing of interest

No Conflict of interest

663

## Acknowledgments

NATL60 simulations were performed on CINES Occigen supercomputer, in the frame of the 2014-2015 Great Challenge organized by GENCI. We also thank Dr Manucharyan and the anonymous reviewers who contributed to improve our document.

668

## Financial Support

The present research was supported by Sorbonne University (IPSL, SAMA, and SCAI).

671

# Appendix A

672

673

Let us denote $e_i$ the kinetic energy at a pixel i:

$$e_i = \frac{1}{2}(U_i^2 + V_i^2) \qquad A1$$

The global kinetic energy of an image for a day d is given by: $E_d = \sum_{i}^{N} e_i$

Let us denote $s_i$ the 'local enstrophy' at a pixel i:

$$s_i = \frac{1}{2} \left( \frac{\partial}{\partial x} V_i - \frac{\partial}{\partial y} U_i \right)^2 \qquad A2$$

The global enstrophy of an image for a day d is given by: $S_d = \sum_{i}^{N} s_i$

679
680 Energy characterizes the intensity of the velocity field whilst enstrophy its shear intensity.
681 Global Energy and Enstrophy are conservative quantities in the absence of forcing and dissipation for
682 two dimensional flows (Maltrud and Vallis, 1993; Merilees and Warn, 1975).
683
684

683

684

# Appendix B

685

686

687 RESACnet architecture
688 We used the Python 3.9, Keras 2.4.0 software to implement the RESAC model. RESAC was trained on
689 an Nvidia GPU with CUDA 11.0 under a Linux system.

690 The RESAC architecture consists of 3 CNNs (Figure 5).

* CNN_R81toR27 has 7 convolutional layers, each layer consisting of 36 feature maps with 6x6 filters using the *relu* activation function. Each layer is followed by a batch-normalization layer. The last layer connects the 36 characteristic maps to a 1x1 mask to fit the output image using a sigmoid activation function.

691

692

693

694

695

* CNN_R27toR09 has 3 layers consisting of 24 characteristic maps with 5x5 masks using the *relu* activation function. Each layer is followed by a batch-normalization layer. The last layer connects the 24 characteristic maps to a 1x1 mask to fit the output image using a sigmoid activation function.

696

697

698

699

700

* CNN_UV_R09 has the same structure as CNN_R27toR09 but repeated 3 times. The first one takes the SST_R09 and the output of CNN_R27toR09 to make one output UV_R09. The second and third one get that output and are dedicated to finding respectively the U and V components of the current.

701

702

703

704

705

706 RESAC_SSH has 400 301 weights and RESAC_SSH_SST 402 797 weights.

707 From a layer to the next one, the size of the input is preserved using a 'same' padding parameter. The
708 RESAC architecture weights are initialized by the 'the_normal' method (Keras documentation).
709 The cost function J is the MSE of the error for all the cost functions.

$$J = \frac{1}{N} \sum_{i=1}^{N} (y_{i}^{NATL} - y_{i}^{est})^{2} \qquad B1$$

710 The learning, validation and test sets are completely decorrelated (see Section 2). The number of 711 convolution layers, the size and number of filters, the activation and loss functions were chosen based 712 on previous experiments, where the learning rate schedulers and the number of epochs were 713 empirically tuned on the validation set.

714 For each RESAC model, we trained 20 networks with different weights for the initialization, using an 715 adapted learning rate scheduler. For RESAC-ssh-sst the learning rate starts at 2e-3 and stays constant 716 through 20 epochs, is then multiplied by exp(-0.02) at each epoch between epoch 20 and 60, and by 717 exp(-0.05) between epoch 60 and 100. For RESAC-ssh the learning rate starts at 5e-3 and stays 718 constant through 20 epochs, is then multiplied by exp(-0.02) at each epoch between epoch 20 and 40, 719 and by exp(-0.05) between epoch 40 and 100.

720 The final prediction is the average of the predictions of the 20 networks.

721 Figure B1 presents the evolutions of the cost functions during the training phase for the learning and 722 validation datasets. The x-axe stands for time evolution, y-axe represents the value of the cost 723 function. The weights are those determined by cross-validation.

<figure><img src="image_11.png" /></figure>
725

726

727

728

729

**Figure B1.** RESAC-ssh-sst cost functions (in logarithm scale) with respect to the number of iterations (Epochs). See Figure 5 for their definition and Eq. B1 for their mathematical expression.

730

Figure B2 presents the flow diagram for the training of RESAC-ssh

731

<figure><img src="image_12.png" /></figure>

732

733

734

735

736

737

738

739

740

**Figure B2.** Flow diagram of the RESAC-ssh training phase: the inputs of the RESAC-ssh model are in green, the observations used for controlling the RESAC-ssh weights are in blue, the double arrows indicating the four cost functions J1, J2, J3u, J3v. UP3x3 stands for the up-sampling procedure that increases the input size by a factor 3x3 in order to reach the higher resolution of the next image. During the learning phase, $\hat{H}$ is the estimated SSH at R27 and R09 respectively, $\hat{U}$ and $\hat{V}$ are the estimated ocean velocity values. The weights of all the CNNs are estimated at once by simultaneously minimizing the four cost functions. In the operational phase, RESAC is forced with SSH_R81, SST_R27 and SST_R09 (in green frames).

741

# Appendix C

742

743

744

## Complex correlation

745

746

747

748

749

750

751

Comparing two vector fields is a quite complicated task. We can compare the vector components with classical statistical tools (std, correlation) but their interpretation in the physical space is not easy. One of the difficulties is to take relationships existing between the vector components into account. In the present paper we propose a method leading to a straightforward geometrical interpretation for two-dimensional vector fields. The trick is to use the isomorphism existing between two-dimensional vectors and complex numbers. Let us consider two two-dimensional vector fields $(u, v)$ and $(u', v')$ that we aim at comparing. The associated complex numbers $U$ and $U'$ are written as
752 $U = u+iv = \rho e^{i\alpha}$ and $U' = u'+iv' = \rho' e^{i(\alpha+\beta)}$, where $i$ is the imaginary symbol, $\rho$ and $\rho'$ are positive numbers
753 which represent the modules of the vectors, $\alpha$ and $\alpha+\beta$ define the angles of the vectors with the real
754 axe, $\beta$ representing the angle between the two vectors.

755 Let us define the complex correlation $CC$ between $U$ and $U'$:

$$CC = \frac{\frac{1}{N} \sum_{1}^{N} U_j \bar{U}_j'}{\sqrt{\frac{1}{N} \sum_{1}^{N} U_j \bar{U}_j} \sqrt{\frac{1}{N} \sum_{1}^{N} U_j' \bar{U}_j'}} \qquad C1$$

756 where $\bar{U} = u - iv$ and $\bar{U}' = u' - iv'$ are the complex conjugates associated with $U$ and $U'$
757 respectively. The sums are taken over the $N$ vectors constituting the fields.

758 $CC$ can also be written as

$$CC = \frac{\frac{1}{N} \sum_{1}^{N} \rho_j \rho_j' e^{-i\beta_j}}{\sqrt{\frac{1}{N} \sum_{1}^{N} \rho_j \rho_j} \sqrt{\frac{1}{N} \sum_{1}^{N} \rho_j' \rho_j'}} = Ce^{-i\theta} \qquad C2$$

759 This polar representation allows us geometrical interpretations of $CC$ as shown in the following.

760 - First, let us assume that there is a systematic rotation $\beta$ between the two vector fields. $CC$
761 writes

$$C = e^{-i\beta} \frac{\sum_{1}^{N} \rho_j \rho_j'}{\sqrt{\sum_{1}^{N} \rho_j \rho_j} \sqrt{\sum_{1}^{N} \rho_j' \rho_j'}} = C e^{-i\beta} \qquad C3$$

762 where $C$ is the correlation between the vector modules and $\beta$ the systematic angle between the two
763 vector fields.

764 This case is often encountered in oceanography (Constant Angle between the wind and the
765 current in the Ekman drift; rotation of the current with depth in the Ekman spiral).
767 - Second, let us interpret the angle of the complex correlation when the two fields are quasi
768 co-linear, which is the case when we compared modeled and observed currents.

769 Eq. C2 can be written as

$$
770 \qquad CC = \frac{\sum_{1}^{N} \rho_j \rho_j' (\cos \beta_j - i \sin \beta_j)}{\sqrt{\sum_{1}^{N} \rho_j \rho_j} \sqrt{\sum_{1}^{N} \rho_j' \rho_j'}} \qquad C4
$$

771

772 Eq. C4 is somewhat difficult to interpret. Let us consider the following simple case: we now
773 assume that each vector field has a constant modulus (ρ and ρ' are constant), C4 writes

$$
774 \qquad CC = \frac{\rho \rho' \sum_{1}^{N} (\cos \beta_j - i \sin \beta_j)}{\sqrt{\sum_{1}^{N} \rho_j \rho_j} \sqrt{\sum_{1}^{N} \rho_j' \rho_j'}} = \frac{1}{N} \sum_{1}^{N} (\cos \beta_j - i \sin \beta_j) \qquad C5
$$

774 If the angles $\beta_j$ are small ($\theta < 10^\circ$) and their mean close to zero, which is a straightforward hypothesis,
775 Eq. C5 writes:

$$
776 \qquad CC = \frac{1}{N} \sum_{1}^{N} (\cos\beta_j - i \sin\beta_j) \approx 1 - \frac{1}{N} \sum_{1}^{N} (i \sin\beta_j) \approx 1 - i \varepsilon \qquad C6
$$

776 Where $\epsilon$ is a small number, due to the fact that the mean of $\beta_j$ is close to zero.

777 CC can be written as $CC = e^{-i\theta}$ where $\theta$ is a small angle.

778 This simplified example can explain the fact that the angles $\vartheta$ of the complex correlations we have
779 computed in the present study are always small. The CC angle gives an estimate of the mean angle
780 between the two vector fields which is near zero if the angle distribution is Gaussian.
Complex correlation has been used in several studies implying two-dimensional vector analysis in meteorology and in oceanography (Wang and Moers, 1977; Hardy and Walton, 1978; Legler, 1983; Horel, 1984).

Eq. C6 shows that CC gives an estimate of the mean of the angles $\beta_j$ between the two vector fields, which is close to zero and does not provide useful information on the structure of the two fields. A more informative parameter is the standard deviation of that angle which could be obtained by doing the following operation: $U_j \overline{U_j'} = \rho_j \rho_j' e^{-i\beta_j}$ which provides the $\beta_j$ and allows the computation of the standard deviation of the $\beta_j$.

In fact, the correlation between two-dimensional vector fields is a matrix $Cm$ of the form:

$$Cm = \begin{vmatrix} C_{uu'} & C_{uv'} \\ C_{vu'} & C_{vv'} \end{vmatrix}$$

Where $C_{uv'}$ is the correlation between the component $u$ of vector $U$ and $v'$ of vector $U'$ and so on

The matrix formulation of the correlation existing between two vector fields contains more information than CC, but CC is easier to interpret in the physical space. CC is mainly dedicated to compare vector fields which are related by a scaling and a systematic rotation. In the present study, we chose the simplicity of the interpretation and dealt with complex correlation.

## Appendix D

The RESAC model does not learn the time evolution of Eq. (1). But it better approximates Eq. (1) than geostrophy (balance between term 2 and 3) by taking into account a statistical modeling of term 4 and of the spatial contribution of term 1, which is of the form $(U \cdot \nabla)U$. RESAC therefore can be considered as a statistical model of the non-time dependent PDE: $(U \cdot \nabla)U - f \times U = -g \nabla \times SH + D^U$ coupled with
the slowly varying temperature equation (Eq.2) which is of the form: $\nabla(\text{SST.}U) = D^{\text{SST}}$. The operators $D^U$ and $D^{\text{SST}}$ model turbulence.

808

809

810

811

812

813

814

815

816

817

818

819

820

References

821

822

823

824

825 Ajayi, A., Le Sommer, J., Chassignet, E., Molines, J.-M., Xu, X., Albert, A. & Cosme, E. 2020. Spatial and temporal variability of the North Atlantic eddy field from two kilometric-resolution ocean models. J. Geophys. Res. Oceans, 125, e2019JC015827. https://doi.org/10.1029/2019JC015827

828 Ballarotta, M., Ubelmann, C., Pujol, M.-I., Taburet, G., Fournier, F., Legeais, J.-F., et al., 2019. On the resolutions of ocean altimetry maps. Ocean Science, 15(4), 1091-1109.

830 Barth A., Alvera-Azcárate A., Matjaz Licer M., and Beckers J.-M., 2020. DINCAE 1.0: a convolutional neural network with error estimates to reconstruct sea surface temperature satellite observations. Geosci. Model Dev., 13, 1609–1622, 2020 https://doi.org/10.5194/gmd-13-1609-2020

833 Bowen, M.M., Emery, W.J., Wilkin, J.L., Tildesley, P.C., Barton, I.J., Knewtson, R., 2002. Extracting multiyear surface currents from sequential thermal imagery using the maximum cross-correlation technique. J. Atmos. Ocean. Technol. 19 (10), 1665–1676.
836 Cazenave et al., 2018. (WCRP Global Sea Level Budget Group). Global sea level budget 1993-present.

837 Earth Syst. Sci. Data, 10, 1551–1590, 2018. https://doi.org/10.5194/essd-10-1551-

838 de Bézenac E., Pajot A., and Gallinari P., 2019. Deep learning for physical processes: incorporating

839 prior scientific knowledge. J. Statistical Mechanics: Theory and Experiment, IOP Publishing Ltd and

840 SISSA Medialab srl.

841 Ducournau A. and R. Fablet, 2018. Deep learning for ocean remote sensing: an application of

842 convolutional neural networks for super-resolution on satellite-derived SST data," 2016 9th IAPR

843 Workshop on Pattern Recognition in Remote Sensing (PRRS), 2016, pp. 1-6, doi:

844 10.1109/PRRS.2016.7867019.

845 Dufau C., Orsztynowicz M., Dibarboure G., Morrow R., Le Traon P.Y., 2016. Mesoscale resolution

846 capability of altimetry: present and future. J. Geophys. Res. Oceans, 121, 4910-4927,

847 https://doi.org/10.1002/2015JC010904.

848 Dussin R., Barnier B., Brodeau L., and Molines J. M., 2016. The Making Of the DRAKKAR FORCING SET

849 DFS5. DRAKKAR/MyOcean Report 01-04-16, (available at http://www.drakkar-ocean.eu/)

850 Egmont-Petersen, M.; de Ridder, D.; Handels, H., 2002. Image processing with neural networks—a

851 review. Pattern Recognition, 35, 2279–2301, doi:10.1016/S0031-3203(01)00178-9.

852 Emery, W.J., Thomas, A.C., Collins, M.J., Crawford, W.R., Mackas, D.L., 1986. An objective method for

853 computing advective surface velocities from sequential infrared satellite images. J. Geophys. Res.,

854 91, 12 865–12 87.

855 Fablet R., Ouala S., Herzet C., 2018. Bilinear residual Neural Network for the identification and

856 forecasting of dynamical systems. EUSIPCO 2018: European Signal Processing Conference, Sep 2018,

857 Rome, Italy. pp.1-5, 10.23919/EUSIPCO.2018.8553492. , 01686766.

858 Febvre Q, R. Fablet, J. L. Sommer and C. Ubelmann, "Joint Calibration and Mapping of Satellite

859 Altimetry Data Using Trainable Variational Models," ICASSP 2022 - 2022 IEEE International

860 Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 1536-1540, doi:

861 10.1109/ICASSP43922.2022.9746889.

862 George, T.M., Manucharyan, G.E. and Thompson, A.F., 2021. Deep learning to infer eddy heat fluxes

863 from sea surface height patterns of mesoscale turbulence. Nature communications, 12(1), pp.1-11.

864 Goodfellow I., Bengio Y., Courville A., 2016. Machine Learning Basics (chapt 4). In Machine Learning,

865 800 p, Mit Press, Cambridge, MA, USA. ISBN: 978-0262035613.

866 Han J, Xu L., Rafique MM., Butt AR, Lim SH., 2019. A quantitative study of deep learning training on

867 heterogeneous supercomputers: IEEE International Conference on Cluster Computing (CLUSTER), 1-

868 12.

869 Hardy D. M. and Walton J. J., 1978. Principal component analysis of wind measurements. J. Appl.

870 Meteor., 17, 1153-1162.
871 He, K., Zhang, X., Ren, S., and Sun, J., 2016. Deep Residual Learning for Image Recognition. IEEE
872 Conference on Computer Vision and Pattern Recognition (CVPR), IEEE, Las Vegas, NV, 27–30 June
873 2016, 770–778, https://doi.org/10.1109/cvpr.2016.90.

874 Horel J. D., 1984. Complex principal component analysis : Theory and examples. J. Clim. Appl. Meteor.
875 23, 1660-1673.

876 Huang, Y; Tang, Y; Zhuang, H; VanZwieten, J; Cherubin, L Physics-informed Tensor-train ConvLSTM
877 for Volumetric Velocity Forecasting of Loop Current, Front. Artif. Intell. 4:780271. doi:
878 10.3389/frai.2021.780271)

879 Isern-Fontanet, J., Chapron, B., Lapeyre G., Klein P., 2006. Potential use of microwave Sea Surface
880 Temperatures for the estimation of ocean currents. Geophys. Res. Lett. 33, L24608,
881 http://dx.doi.org/10.1029/2006GL027801.

882 Kelly K., Sandipa Singh, and Rui Xin Huang, 1999. Seasonal Variations of Sea Surface Height in the Gulf
883 Stream Region. J. Physc. Oceano. V29 pp 313-327.
884 DOI:https://doi.org/10.1175/1520-0485(1999)029<0313:SVOSSH>2.0.CO;2

885 Le Cun Y, Bengio Y, Hinton G., 2015. Deep learning. Nature 521, 436–444, doi:10.1038/nature14539.

886 Le Cun Y, Boser B, Denker J, Henderson D, Howard R, Hubbard W, Jackelt L., 1990. Handwritten digit
887 recognition with a back-propagation network. Advances in neural information processing systems, 2
888 ed. DS Touretzky), 396–404. San Francisco, CA: Morgan Kaufmann.

889 Legler D. M., 1983. Empirical orthogonal function analysis of wind vectors over the tropical Pacific
890 Region. Bull. Amer. Meteor. Soc., 64, 234-241.

891 Levy M., Ferrari R., Franks PJS., Martin AP, Rivière P., 2012. Bringing physics to life at the
892 submesoscale. Geophys. Res. Lett. 39, L14602.

893 Lguensat R., Huynh Viet P., Sun M., Chen G., Fenglin T., Chapron B. and Fablet R., 2019. Data-Driven
894 Interpolation of Sea Level Anomalies Using Analog Data Assimilation Remote Sens., 11(7), 858;
895 https://doi.org/10.3390/rs11070858.

896 Madec G., and Imbard M., 1996. A global ocean mesh to overcome the north pole singularity. Clim.
897 Dyn., 12, 381–388.

898 Maltrud M. E. and Vallis G.K., 1993. Energy and enstrophy transfer in numerical simulations of two-
899 dimensional turbulence. Physics of Fluids A: Fluid Dynamics 5, 1760.
900 https://doi.org/10.1063/1.858851.

901 Manucharyan, G.E., Siegelman, L. and Klein, P., 2021. A deep learning approach to spatiotemporal sea
902 surface height interpolation and estimation of deep currents in geostrophic ocean turbulence.
903 Journal of Advances in Modeling Earth Systems, 13(1), p.e2019MS001965
Martinez E., Gorgues T., Lengaigne M., Fontana C., Sauzède R., et al., 2020. Reconstructing Global Chlorophyll-a Variations Using a Non-linear Statistical Approach. Frontiers in Marine Science, 7, https://doi.org/10.3389/fmars.2020.00464.

McWilliams, J. C., 2016. Submesoscale currents in the ocean. Proc. R. Soc. A 472, 20160117 (2016). 10.

Merilees, P. and Warn, H., 1975. On energy and enstrophy exchanges in two-dimensional non-divergent flow. Journal of Fluid Mechanics, 69(4), 625-630. doi:10.1017/S0022112075001590.

Nardelli, B., D. Cavaliere, E. Charles, and D. Ciani (2022), Super-Resolving Ocean Dynamics from Space with Computer Vision Algorithms, Remote Sens., 14(5), 1159, doi:10.3390/rs14051159.

Pannekoucke O. and Fablet R., 2020. PDE-NetGen 1.0: from symbolic partial differential equation (PDE) representations of physical processes to trainable neural network representations. Geosci. Model Dev., 13, 3373-3382, https://doi.org/10.5194/gmd-13-3373-2020.

Reichstein, M., Camps-Valls, G., Stevens, B. et al., 2019. Deep learning and process understanding for data-driven Earth system science. Nature 566, 195-204. https://doi.org/10.1038/s41586-019-0912-1.

Richaume P., F. Badran, M. Crepon, C. Mejia, H. Roquet and S. Thiria, 2000. Neural network wind retrieval from ERS-1 scatterometer data. J. Geophys. Res. 105 pp 8737-8731

Rio M. H. and Santoleri R., 2018. Improved global surface currents from the merging of altimetry and Sea Surface Temperature data. Remote Sens. Envir. 216, 770-785.

Rio, M.-H., Santoleri, R., Bourdalle-Badie, R., Griffa, A., Piterbarg, L., Taburet, G, 2016. Improving the altimeter derived surface currents using high-resolution Sea Surface Temperature data: a feasibility study based on model outputs. J. Atmos. Ocean. Technol. 33. http://dx.doi.org/10.1175/JTECH-D-16-0017.1

Rolnick D., Donti PL, Kaack L., Kochanski K., Lacoste A., Sankaran K., Ross A., Milojevic-Dupont N., Jaque N., WaldmanBrown A., et al., 2019. Tackling climate change with machine learning,” arXiv preprint arXiv:1906.05433.

Rousseau, F., Drumetz, L., and Fablet, R., 2019. Residual Net-works as Flows of Diffeomorphisms, J. Math. Imaging Vis., https://doi.org/10.1007/s10851-019-00890-3.

Ruthotto, L. and Haber, E., 2019. Deep Neural Networks Motivated by Partial Differential Equations, J. Math. Imaging Vis., 62, 352-364, https://doi.org/10.1007/s10851-019-00903-1.

Sasaki, H., Klein, P., Qiu, B. & Sasai, Y., 2014. Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere. Nat. Commun. 5, 5636.

Simard P., Steinkraus, D., J. C. Platt, J. C., 2003. Best Practices for Convolutional Neural Networks Applied to Visual Document Analysis. Proceedings of the Seventh International Conference on Document Analysis and Recognition (ICDAR 2003) 0-7695-1960-1/03 $17.00 © 2003 IEEE.
938

939

940

941

942

943

944

945

946

947

948

949

950

951

952

953

954

955

956

957

Sonnewald M. Lguensat R. Jones D. Dueben P. Brajard J. and Balaji V. 2021. Bridging observations, theory and numerical simulation of the ocean using machine learning. Environ. Res. Lett. 16. 073008 https://doi.org/10.1088/1748-9326/ac0eb0

Sonnewald, M., & Lguensat, R. (2021). Revealing the impact of global heating on North Atlantic circulation using transparent machine learning. Earth and Space Science Open Archive. ESSOAr. https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021MS002496.

Ubelmann C, Klein P. and Fu L., 2015. Dynamic Interpolation of Sea Surface Height and Potential Applications for Future High-Resolution Altimetry Mapping. J. Atmos. Oceanic Technol. 32, 177-184.

Vigan, X., Provost, C., Bleck, R., Courtier, P., 2000. Sea surface velocities sea surface temperature image sequences 1. Method and validation using primitive equation model output. J. Geophys. Res. 105, 19499-19514.

"Wang, J.L.; Zhuang, H.; Chérubin, L.M.; Ibrahim, A.K.; Muhamed Ali, A. 2019. Medium-Term Forecasting of Loop Current Eddy Cameron and Eddy Darwin Formation in the Gulf of Mexico With a Divide-and-Conquer Machine Learning Approach. J. Geophys. Res. Ocean. 124, 5586-5606.

Wang D.-P. and C. N. K. Moers C. N. K., 1977. Long coast trapped waves off the West Coast of the United States, summer 1973 J. Phys. Oceanogr. 7, 856-864

Zakharov VE, L'vov VS., Falkovich G., 2012. Kolmogorov spectra of turbulence : Wave turbulence, Springer Verlag (Berlin, Heidelberg) editor, 263 p, DOI 10.1007/978-3-642-50052-7.

•