# Inversion of sea surface currents from satellite-derived SST-SSH synergies
with 4DVarNets

Ronan Fablet
IMT Atlantique  Lab-STICC  INRIA team Odyssey
Brest  France
<EMAIL_ADDRESS>
&Bertrand Chapron
Ifremer  LOPS  INRIA team Odyssey
Brest  France
<EMAIL_ADDRESS>
&Julien Le Sommer
CNRS  IGE
Grenoble  France
<EMAIL_ADDRESS>
&Florian Sévellec
CNRS  LOPS
Brest  France
<EMAIL_ADDRESS>
1  2  3  4

###### Abstract

Satellite altimetry is a unique way for direct observations of sea surface
dynamics. This is however limited to the surface-constrained geostrophic
component of sea surface currents. Ageostrophic dynamics are however expected
to be significant for horizontal scales below 100 km and time scale below 10
days. The assimilation of ocean general circulation models likely reveals only
a fraction of this ageostrophic component. Here, we explore a learning-based
scheme to better exploit the synergies between the observed sea surface
tracers, especially sea surface height (SSH) and sea surface temperature
(SST), to better inform sea surface currents. More specifically, we develop a
4DVarNet scheme which exploits a variational data assimilation formulation
with trainable observations and a priori terms. An Observing System Simulation
Experiment (OSSE) in a region of the Gulf Stream suggests that SST-SSH
synergies could reveal sea surface currents for time scales of 2.5-3.0 days
and horizontal scales of 0.5∘-0.7∘, including a significant fraction of the
ageostrophic dynamics ($\approx$ 47%). The analysis of the contribution of
different observation data, namely nadir along-track altimetry, wide-swath
SWOT altimetry and SST data, emphasizes the role of SST features for the
reconstruction at horizontal spatial scales ranging from
$\nicefrac{{1}}{{20}}$∘ to $\nicefrac{{1}}{{4}}$∘.

Enter journal name here

IMT Atlantique, Lab-STICC, INRIA team Odyssey, Brest, France Ifremer, LOPS,
INRIA team Odyssey, Brest, France CNRS, IGE, Grenoble, France CNRS, LOPS,
Brest, France Ronan<EMAIL_ADDRESS>
We present end-to-end deep learning schemes to improve the reconstruction of
total sea surface currents from satellite-derived observations;

Experiments for a region of the Gulf Stream support the synergistic analysis
of sea surface temperature and sea surface height data;

The strain of sea surface dynamics is a proxy of the uncertainty of the
retrieved estimation.

## Plain Language Summary

Satellite altimetry is a unique way for the direct observation of sea surface
currents. It can however only inform their geostrophic component. Combined
with a scarce sampling, this results in a limited ability of operational
products to retrieve a significant fraction of the sea surface currents in the
meso-scale range. We introduce a novel deep learning scheme based on a
variational data assimilation formulation. It offers new means to exploit
multimodal synergies between satellite-derived observations of the sea
surface, here sea surface height (SSH) and sea surface temperature (SST).
Numerical experiments using realistic simulations for a case-study area in the
Gulf Stream highlight the potential to enhance the reconstruction of sea
surface currents and to retrieve a significant fraction of the ageostrophic
dynamics. They support the relevance of physics-informed deep learning schemes
to make the most of available multimodal observation datasets and numerical
simulations to reveal partially-observed sea surface dynamics.

##  1 Introduction

Satellite altimeters provide the main source of observations to inform sea
surface dynamics on a regional and global scale (chelton_satellite_2001).
Their scarce space-time sampling of the sea surface prevents to recover
spatial scales below 100 km and time scales below 10 days in general
(ballarotta_resolutions_2019). Also altimetry can only reconstruct geostrophic
velocities. As stressed by simulation and observational studies, the
ageostrophic components are however critical features of upper ocean dynamics
regarding for instance vertical mixing properties (mahadevan_analysis_2006),
Lagrangian dynamics at sea surface (baaklini_blending_2021; sun_impacts_2022).

Retrieving sea surface currents at finer scales with both their geostrophic
and ageostrophic components naturally arises as a key challenge. This has
motivated a large research effort both in terms of simulation studies
(uchida_cloud-based_2022), observational effort (villas_boas_integrated_2019;
ardhuin_skim_2019), and data assimilation methods (moore_synthesis_2019;
storto_ocean_2019). Regarding the latter aspect, state-of-the-art approaches
mostly rely on the one hand on optimal interpolation approaches
(taburet_duacs_2019) and on the other hand on data assimilation schemes
combined with ocean general circulation models (baaklini_blending_2021;
benkiran_assessing_2021; fujii_observing_2019). As mentionned above, both
approaches still show limitations in the ability to retrieve fine-scale
patterns, whereas both observation-driven and theoretical studies evidence the
interplay between fine-scale sea surface dynamics and some observed processes
such as sea surface tracers (ciani_ocean_2021; isern-fontanet_potential_2006)
and drifters’ trajectories (sun_impacts_2022).

From a methodological point of view, data-driven and learning-based schemes
have also received a growing attention to solve inverse problems in geoscience
(alvera-azcarate_multivariate_2007; lguensat_analog_2017; barth_dincae_2020).
Especially, deep learning schemes appears as appealing schemes for the
reconstruction of sea surface dynamics from irregularly-sampled satellite-
derived observations (fablet_end--end_2021; fablet_multimodal_2022;
george_deep_2021; manucharyan_deep_2021). Interestingly, these studies open
new research avenues to make the most of available simulation and observation
datasets. They also suggest potential breakthrough through the ability to
exploit the synergies between different sea surface observational fields with
no explicitly-known relationship (fablet_multimodal_2022).

In this study, we exploit these recent methodological advances to explore
satellite-derived SST-SSH synergies to inform sea surface currents, including
their ageostrophic component. We exploit and adapt multi-modal 4DVarNet
schemes introduced in (fablet_multimodal_2022). Through an observing system
simulation experiment for a region of the Gulf Stream, our key contributions
are four-fold:

  1. 1.

We stress the potential of physics-informed deep learning schemes to enhance
the reconstruction of sea surface currents (SSC) with a relative improvement
greater than 50% in terms of resolved space-time scales and mean-square-error
metrics compared with the state-of-the-art products;

  2. 2.

Our results also support wide-swath SWOT altimetry data to improve the
reconstruction of SSC fields with a potential gain of $\approx 30\%$ for most
performance metrics, though the reported improvement is marginal for the
divergent component of the SSC.

  3. 3.

We emphasize the contribution of SST-SSH synergies to retrieve a significant
fraction of the ageostrophic component of the SSC, typically $\approx$ 47% in
terms of divergence of the SSC fields, with a major contribution of SST
features in the spatial scale range of
$\nicefrac{{1}}{{20}}$∘-$\nicefrac{{1}}{{4}}$∘.

  4. 4.

We point out that, as hypothesized from theoretical considerations
(Sevellec:22), the strain of sea surface dynamics partially explains
($\approx$ 60%) the time-averaged mean square error of the SSC.

We further discuss the implications of these results in the context of ongoing
research efforts towards the monitoring of upper ocean dynamics.

##  2 Problem statement

Geostrophic and ageostrophic sea surface dynamics: The horizontal momentum
equations together with the hydrostatic equilibrium for the vertical momentum
and the non-divergence leads to the following set of equations to describe sea
surface dynamics:

$D_{t}u-fv=-\frac{1}{\rho_{0}}\partial_{x}P+\mathcal{F}_{x},$ (1a)
$D_{t}v+fu=-\frac{1}{\rho_{0}}\partial_{y}P+\mathcal{F}_{y},$ (1b)
$0=-\frac{1}{\rho_{0}}\partial_{z}P-\frac{g}{\rho_{0}}\rho,$ (1c)
$\partial_{x}u+\partial_{y}v+\partial_{z}w=0,$ (1d)

where $u$, $v$ and $w$ are the zonal, meridional, and vertical velocities,
respectively, $t$ is time, $x$, $y$ and $z$ are the longitude, latitude, and
depth, respectively, $\rho_{\left(0\right)}$ is the (reference) density for
seawater, $P$ is the pressure, $f$ is the Coriolis parameter, $g$ is the
acceleration due to gravity, $\mathcal{F}_{x}$ and $\mathcal{F}_{x}$ are the
action of the zonal and meridional viscous forces, respectively, and $D_{t}$
(=$\partial_{t}$+$u\partial_{x}$+$v\partial_{y}$+$w\partial_{z}$) is the
material derivative.

These equations are often simplified to reflect the geostrophic balance which
occurs under the small Rossby number (Ro$\ll$1 i.e., inertial terms are
negligible), slow dynamics ($\partial_{t}\rightarrow$0 for the horizontal
momemtum equation), large Péclet number (Pe$\ll$1 i.e., viscous terms are
negligible), and away from direct forcing. This reads:

$-fv_{g}=-\frac{1}{\rho_{0}}\partial_{x}P,$ (2a)
$+fu_{g}=-\frac{1}{\rho_{0}}\partial_{y}P,$ (2b)

where $u_{g}$ and $v_{g}$ are the zonal and meridional geostrophic velocities,
respectively. This formulation still implies an horizontal divergence, which,
using (2b) with (1d), reads:

$\partial_{x}u_{g}+\partial_{y}v_{g}=-\frac{\beta}{f}v_{g},$ (3)

where $\beta$ (=$\partial_{y}\,f$) accounts for the meridional variation of
Earth equivalent rotation rate. Hence, the horizontal momentum equations for
the ageostrophic components become:

$D_{t}u-fv_{a}=\mathcal{F}_{x},$ (4a) $D_{t}v+fu_{a}=\mathcal{F}_{y},$ (4b)

where $u_{a}$ and $v_{a}$ are the zonal and meridional ageostrophic
velocities, respectively. This last set of equations show the complexity of
the ageostrophic components through the action of viscous terms and of the
full inertial terms in (4b), compared to the rather simple linear
relationships between pressure gradients and geostrophic velocities in (2b).
This is worth noting that boundary conditions, such as observed at the ocean
surface, occured though the actions of the viscous terms ($\mathcal{F}_{x}$
and $\mathcal{F}_{y}$). This further suggests the key role of the ageostrophic
terms when studying the surface velocities.

In the subsequent, for the sake of simplicity, we will refer to
geostrophically-derived sea surface velocities from SSH fields from (2b) using
$\mathrm{SSH}=\nicefrac{{(P|_{z=0}-P_{\mathrm{atm}})}}{{g\rho_{s}}}$ [where
$P|_{z=0}$ is the pressure at the geoide ($z=0$), $P_{\mathrm{atm}}$ is the
atmospheric pressure, and $\rho_{s}$ is the assumed-vertically-constant
density of the surface layer] as SSH-derived sea surface currents.

Data assimilation for sea surface dynamics: Classically, inverse problems in
geoscience (evensen_data_2009) are stated as data assimilation problems
through some underlying state-space formulation

$\left\\{\begin{array}[]{ccl}\displaystyle\frac{\partial\mathbf{x}(t)}{\partial
t}&=&{\cal{M}}\left(\mathbf{x}(t)\right)+\eta(t)\\\ {}\hfil\\\
\mathbf{y}_{k}(t)&=&{\cal{H}}_{k}\left(\mathbf{x}(t)\right)+\epsilon_{k}(t),\forall
t,k\\\ \end{array}\right.$ (5)

where $\mathbf{x}$ is the space-time process to be reconstructed and
$\mathbf{y}_{k}$ is an observation process which relates to state $\mathbf{x}$
through observation operator ${\cal{H}}_{k}$ for observation modality $k$.
When dealing with irregularly-sampled observations, operator ${\cal{H}}_{k}$
accounts for sampling masks. Processes $\eta$ and $\epsilon$ refer to random
processes to account for modeling uncertainties and observation noise,
respectively. ${\cal{M}}$ refers to the dynamical prior on state $x$. Given
this state-space formulation, the data assimilation problem for the
reconstruction of state $x$ given observation data $\\{y_{k}\\}_{k}$ comes to
the resolution of a minimization problem. Within a variational data
assimilation framework, it generally writes as:

$\widehat{\mathbf{x}}=\arg\min_{\mathbf{x}}\sum_{k}\lambda_{k}\left\|\mathbf{y}_{k}-{\cal{H}}_{k}(\mathbf{x})\right\|^{2}+\gamma\left\|\mathbf{x}-\Phi_{\cal{M}}(\mathbf{x})\right\|^{2}$
(6)

with $\\{\lambda_{k}\\}$ and $\gamma$ Lagrangian multipliers,
$\Phi_{\cal{M}}(\mathbf{x})$ the time-stepping operator to propagate one-step-
ahead state $\mathbf{x}$ at time $t$ to time $t+\Delta t$ based on dynamical
prior ${\cal{M}}$. In the above formulation, we consider a matrix form and
drop the time variable such that the norms are evaluated as a sum over a given
time interval $[0,T]$ according to time step $\Delta t$. Given dynamical prior
${\cal{M}}$ and observation operators $\\{{\cal{H}}_{k}\\}_{k}$, data
assimilation methods provide different algorithms (carrassi_data_2018;
evensen_data_2009) to solve this minimization problem especially using
adjoint-based gradient descent schemes and Kalman methods. Formulation (6)
also relates to Optimal interpolation (cressie_statistics_2015) when
considering a single observation term with a masking operator operator and a
prior given by a Gaussian process. Under these hypothesis, one can derive the
analytical solution of the resulting linear-quadratic variational cost.

The state-of-the-art methods for the reconstruction of sea surface dynamics
from satellite-derived observations rely on such data assimilation schemes.
While the reconstruction of geostrophic sea surface dynamics may be stated as
a space-time interpolation of satellite altimetry data (taburet_duacs_2019),
the inference of the total sea surface currents from satellite-derived
observations using the set of equations introduced above implies to
reconstruct not only these velocities, but the whole ocean state. The
assimilation of satellite altimetry and satellite-derived SST observations,
possibly complemented by other data sources, in ocean general circulation
models (benkiran_assessing_2021; fujii_observing_2019) follows this strategy.
An alternative approach relies on assimilation schemes with a simplified prior
on sea surface dynamics, which only depend on sea surface velocities. Quasi-
geostrophic (QG) dynamics are examples of such priors
(le_guillou_mapping_2020; ubelmann_dynamic_2014). Similarly to the above-
mentionned OI schemes, the latter only applies to geostrophic velocities and
cannot recover ageostrophic components.

Recently, a rich literature has emerged to bridge data assimilation and deep
learning (abdalla_altimetry_2021; barthelemy_super-resolution_2021;
boudier_dan_2020; bocquet_bayesian_2020; fablet_learning_2021;
nonnenmacher_deep_2021). It provides new minimization schemes as well as new
means to explore data assimilation problems when the observation operators
and/or the dynamical priors are not explicitly known. As such, it opens new
avenues to balance the complexity of the inversion problem and the genericity
of the underlying variational formulation in (6). Here, as detailed in the
subsequent, we benefit from the generic 4DVarNet framework introduced in
(fablet_learning_2021) and explore a learning-based data assimilation schemes
for the reconstruction of sea surface currents from multimodal satellite-
derived observations.

##  3 Data

###  3.1 Case-study region and NATL60 data

The considered case study focuses on a $10^{\circ}$$\times$$10^{\circ}$ region
between (33∘N, 65∘W) and (43∘N, 55∘W). As illustrated in Fig. 1, this region
involves the main meander of the Gulf Stream as well as a variety of mesoscale
eddies and finer-scale sub-mesoscale filaments. It also comprises clear
divergence features associated with the ageostrophic flow component, which
makes it suitable for the current study.

The considered simulation dataset relies on a nature run of the NATL60
configuration (ajayi_spatial_2020) of the NEMO (Nucleus for European Modeling
of the Ocean) model (madec_nemo_2022). This simulation delivers a realistic
hindcast simulation of ocean dynamics, including mesoscale-to-submesoscale
ocean dynamics, (ajayi_spatial_2020) over one year from October 2012 to
September 2013. It has been used in numerous studies regarding the
reconstruction and observability of sea surface dynamics (see below for the
related OSSE data challenge stated in (le_guillou_mapping_2020)).

|  |  |
---|---|---|---
SSH |  SST |  Altimetry data |  SSC
|  |  |
|  |  |
SSC vorticity |  SSH-derived vorticity |  SSC divergence |  SSC strain
Figure 1: Illustration of the considered case-study using NATL60 simulation
dataset: first row, from left to right, NATL60 SSH field on November 12th,
associated SST field and SSH observation data from nadir and wide-swath SWOT
altimeters, SSC field. The second row depicts the vorticity fields for the
total currents and SSH-geostrophically-derived ones along with the divergence
and strain of the total currents.

###  3.2 OSSE setting

Our numerical experiments exploit an Observing System Simulation Experiment
(OSSE). We rely on the OSSE setting 111https://github.com/ocean-data-
challenges/2020a_SSH_mapping_NATL60 proposed in (le_guillou_mapping_2020) for
the benchmarking of SSH mapping methods in the context of upcoming wide-swath
altimetry mission SWOT (gaultier_challenge_2015). For a
$\nicefrac{{1}}{{20}}$∘ spatial resolution, this OSSE dataset comprises daily-
averaged SSH fields and simulated altimetry data. The latter combine simulated
nadir along-track data according to real 4-altimeter configuration and SWOT
altimetry data using SWOT simulator (gaultier_challenge_2015). We may point
out that simulated altimetry data are created from hourly SSH fields.

In this study, we complement this initial OSSE dataset with two additional
data sources with the same $\nicefrac{{1}}{{20}}$∘ spatial resolution: the
series of daily-averaged sea surface currents (SSC) and the series of daily-
averaged sea surface temperature (SST). We also generate SST observations with
resolutions of $\nicefrac{{1}}{{10}}$∘, $\nicefrac{{1}}{{5}}$∘,
$\nicefrac{{1}}{{4}}$∘ and $\nicefrac{{1}}{{2}}$∘ using coarsening and
subsampling operations to provide a more realistic simulations of the
diversity of operational L4 SST products (donlon_operational_2012;
ocarroll_observational_2019). We report in Fig. 1 an illustration of the
considered dataset.

For the training configuration, we split the one-year time series into
training, validation, and test datasets according to the following time
periods: from February 4th 2013 to September 30th 2013, from January 1st 2013
to February 4th 2013, and from October 20th 2012 to December 4th 2012,
respectively.

##  4 Methods

This Section presents the proposed 4DVarNet scheme for the reconstruction of
sea surface currents from satellite-derived SST-SSH synergies. We first detail
the considered 4DVarNet architecture before describing our learning scheme.

###  4.1 Proposed 4DVarNet scheme

4DVarNet schemes refer to end-to-end neural architectures introduced in
(fablet_learning_2021) to solve data assimilation problems and further
extended to multimodal inversion schemes in (fablet_multimodal_2022). As
sketched in Fig. 2, the proposed 4DVarNet scheme exploits as inputs time
series of satellite altimetry data and SST fields and outputs gap-free SSH
fields and SSC fields. Importantly, it combines two main components: the
definition of the underlying variational formulation (7), the iterative update
rule of the associated trainable gradient-based solver. The latter relies at
each iteration on the gradient of the variational cost w.r.t. the state using
automatic differentiation tools embedded in deep learning frameworks such as
PyTorch. Overall, the 4DVarNet scheme implements a predefined number,
typically from 5 to 15, of the iterative update rule to map the input partial
observations to the reconstructed sea surface dynamics. We refer the reader to
(fablet_learning_2021) for a detailed presentation of the 4DVarNet framework
along with its link to variational data assimilation formulations.

Figure 2: Sketch of 4DVarNet scheme for the reconstruction of sea surface
current fields: 4DVarNet provides an end-to-end neural architecture to map
multimodal observation data, here satellite-derived SSH and SST observations,
to SSH and SSC fields. It relies on a variational data assimilation and on the
resolution of the associated minimization issue with a trainable gradient-
based descent. We refer the reader to (fablet_learning_2021) for a detailed
presentation of 4DVarNet scheme.

Let us formally introduce the considered variational formulation. Let
$\mathbf{x}$ denote hereafter the space-time state to be reconstructed,
$\mathbf{y}$ the observed altimetry data and $\mathbf{z}$ the SST
observations. The considered variational cost, denoted by
$U_{\Phi}\left(\mathbf{x},\mathbf{y},\mathbf{z}\right)$, is given by:

$\displaystyle
U_{\Phi}\left(\mathbf{x},\mathbf{y},\mathbf{z}\right)=\displaystyle\lambda_{1}\left\|\mathbf{x}-\mathbf{y}\right\|^{2}_{\Omega}+\displaystyle\lambda_{2}\left\|g\left(\mathbf{x}\right)-h\left(\mathbf{z}\right)\right\|^{2}+\displaystyle\gamma\left\|\mathbf{x}-\Phi(\mathbf{x})\right\|^{2},$
(7)

where $\lambda_{1}$, $\lambda_{2}$, and $\gamma$ are Lagrangian multipliers.
State $\mathbf{x}$ combines SSH and SSC components referred to as
$x_{\mathrm{SSH}}$ and $x_{\mathrm{SSC}}$. Besides, as in
(fablet_multimodal_2022), we adopt a two-scale decomposition of the SSH fields
to explicitly use both raw altimetry data and DUACS coarse-scale interpolation
as observation data222We refer the reader to (fablet_multimodal_2022) for a
detailed description of this two-scale decomposition of the SSH component of
the state $\mathbf{x}$.. $\Omega$ refers to a masking operator to account for
observation gaps in the altimetry data as well as for the fact that the SSC
component is never directly observed.

Operators $g()$ and $h()$ are trainable operators which aim at extracting
relevant features from SST observations $z$ to match features extracted from
state $x$. Following (fablet_multimodal_2022), we parameterize $g()$ and $h()$
as non-linear convolutional networks using four convolutional layers and tanh
activations. In our experiments, we exploit a 20-dimensional feature space for
multimodal term
$\left\|g\left(\mathbf{x}\right)-h\left(\mathbf{z}\right)\right\|^{2}$.

Operator $\Phi$ states the considered prior onto the space-time dynamics of
state $\mathbf{x}$. Within a classic model-driven approach, $\Phi$ would refer
to the time-stepping operator to forecast future states from an initial
condition. Following (fablet_learning_2021; fablet_multimodal_2022), we rather
consider a purely data-driven parameterization with a two-scale U-Net
architectures (cicek_3d_2016). The latter proved to be more relevant in
previous numerical experiments for both Lorenz’s systems
(fablet_learning_2021) and SSH mapping (beauchamp_end--end_2022;
fablet_multimodal_2022).

As sketched in Fig. 2, the considered 4DVarNet schemes implement the following
iterative gradient-based solver to minimize variational cost (7):

$\left\\{\begin{array}[]{l}\mathbf{g}^{(k+1)}={\cal{R}}\left[\cdot\nabla_{x}U\left(\mathbf{x}^{(k)},\\{\mathbf{y},\mathbf{z}\\}\right),\mathbf{h}(k),\mathbf{c}(k)\right]\\\
{}\\\
\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}-{\cal{L}}\left(\mathbf{g}^{(k+1)}\right)\\\
\end{array}\right.$ (8)

with ${\cal{R}}$ a LSTM cell, ${\cal{L}}$ a linear operator,
$(\mathbf{h}(k),\mathbf{c}(k))$ the internal state of the LSTM cell at
iteration $k$ and $\mathbf{x}^{(k)}$ the reconstructed state at iteration $k$.
As mentioned above, we consider 4DVarNet schemes with 5 to 15 iterations. The
LSTM cell is a 2D-convolutional LSTM cell with 150-dimensional internal
states. The considered iterative rule can be regarded as a momentum-based
gradient-based descent which has been widely explored for optimizer learning
problems (hospedales_meta-learning_2020).

Overall, the trainable components of the 4DVarNet scheme comprise: operators
$g$ and $h$, prior $\Phi$, LSTM cell ${\cal{R}}$, and linear mapping
${\cal{L}}$. This amounts to a total number of parameters to be trained from
data of about 1.4 M parameters for the considered case-study.

###  4.2 Learning setting

We benefit from the end-to-end feature of the considered 4DVarNet scheme to
run a supervised strategy. It aims at learning all trainable components with a
view to optimizing the reconstruction performance. The training loss then
combines reconstruction losses for SSH and SSC fields accounting respectively
for the SSH and its gradient, the zonal and meridional components of current
state $\mathbf{x}_{\mathrm{SSC}}=(\mathbf{x}_{u},\mathbf{x}_{v})$, and the
divergence of the SSC fields:

${\cal{L}}_{\nabla\mathrm{SSH}}=\sum_{i}\alpha_{\mathrm{SSH}}\left\|\mathbf{x}^{\mathrm{true}}_{\mathrm{SSH},i}-\widehat{\mathbf{x}}_{\mathrm{SSH},i}\right\|^{2}+\alpha_{\nabla\mathrm{SSH}}\left\|\nabla\mathbf{x}^{\mathrm{true}}_{\mathrm{SSH},i}-\nabla\widehat{\mathbf{x}}_{\mathrm{SSH},i}\right\|^{2}$
(9)
${\cal{L}}_{u,v}=\sum_{i}\left\|\mathbf{x}^{\mathrm{true}}_{u,i}-\widehat{\mathbf{x}}_{u,i}\right\|^{2}+\left\|\mathbf{\mathbf{x}}^{\mathrm{true}}_{v,i}-\widehat{\mathbf{x}}_{v,i}\right\|^{2}$
(10)
${\cal{L}}_{\mathrm{div}}=\sum_{i}\left\|\partial_{x}\mathbf{x}^{\mathrm{true}}_{u,i}+\partial_{y}\mathbf{x}^{\mathrm{true}}_{v,i}-\partial_{x}\widehat{\mathbf{\mathbf{x}}}_{u,i}-\partial_{y}\widehat{\mathbf{x}}_{v,i}\right\|^{2}$
(11)

where superscript true (resp. $\widehat{\ }$) refers to the true (resp.
reconstructed) SSH or SSC fields. As proposed in (fablet_learning_2021), we
complement these training losses with additional regularisation terms

${\cal{L}}_{\Phi}=\sum_{i}\left\|\mathbf{x}^{\mathrm{true}}_{i}-\Phi\left(\mathbf{x}^{\mathrm{true}}_{i}\right)\right\|^{2}+\sum_{i}\left\|\widehat{\mathbf{x}}_{i}-\Phi\left(\widehat{\mathbf{x}}_{i}\right)\right\|^{2}$
(12)

These regularisation terms better constrain the training of prior $\Phi$.

Overall, we implement the 4DVarNet scheme and the associated learning strategy
within Pytorch framework. We use Adam as optimizer with a fixed learning rate
of $10^{-3}$. After a first training procedure over 200 epochs, we fine-tune
the best model over the validation dataset over 200 more epochs. The Pytorch
code of our implementation is open source (febvre_pytorch_2022).

###  4.3 Evaluation framework and benchmarked approaches

Our numerical experiments involve a quantitative evaluation of the
reconstruction performance using metrics evaluated over the test dataset. We
adapt the metrics introduced in (le_guillou_mapping_2020) for mapping SSH
fields and geostrophic SSC fields 333We refer the reader to the following link
for the detailed presentation of the evaluation experiment and benchmarked
approaches https://github.com/ocean-data-challenges/2020a_SSH_mapping_NATL60,
especially:

  * •

$\lambda_{t,u}$ and $\lambda_{t,v}$, the minimum time scale resolved in days
for respectively the zonal and meridional velocities;

  * •

$\lambda_{x,u}$ and $\lambda_{x,v}$, the minimum spatial scale resolved in
degrees for respectively the zonal and meridional velocities.

As described in (le_guillou_mapping_2020), these metrics rely on a spectral
analysis. We also evaluate the reconstruction performance in terms of
explained variance, respectively:

  * •

$\tau_{u,v}$, the explained variance for the reconstructed SSC;

  * •

$\tau_{\mathrm{vort}}$, the explained variance for the vorticity of the
reconstructed SSC;

  * •

$\tau_{\mathrm{div}}$, the explained variance for the divergence of the
reconstructed SSC;

  * •

$\tau_{\mathrm{strain}}$, the explained variance for the strain of the
reconstructed SSC (balwada_vertical_2021; okubo_horizontal_1970;
weiss_dynamics_1991).

The last three metrics characterize the extent to which the reconstructed SSC
capture the local deformation tensor of the true velocity fields. Numerically,
the computation of the vorticity, divergence, and strain combines a Gaussian
filtering and finite difference approximation of the first-order derivatives.

Our numerical experiments evaluate two configurations of the proposed 4DVarNet
framework: one using only SSH observation data and the other one exploiting
SST-SSH synergies. For benchmarking purposes, we perform a quantitative
comparison with respect to the SSH-derived sea surface currents from the true
SSH fields and optimally-interpolated DUACS ones (taburet_duacs_2019) using
the geostrophic approximation introduced in Section 2. We also evaluate direct
learning-based inversion schemes based on U-Net architectures (cicek_3d_2016)
to directly map observation data to SSC fields either considering only SSH
observation data or jointly SSH and SST observation data.

##  5 Results

In this Section, we report the considered numerical experiments for the
evaluation of the proposed 4DVarNet schemes for the reconstruction of sea
surface currents.

Approach |  Data used |  $\lambda_{x,u}$ (∘) |  $\lambda_{x,v}$ (∘) |  $\lambda_{t,u}$ (d) |  $\lambda_{t,v}$ (d) |  $\tau_{u,v}$ |  $\tau_{\mathrm{vort}}$ |  $\tau_{\mathrm{div}}$ |  $\tau_{\mathrm{strain}}$
---|---|---|---|---|---|---|---|---|---
True SSH |  SSH only |  0.36 |  0.17 |  19.6 |  11.2 |  97.0% |  96.3% |  $-$1.0% |  92.1%
DUACS |  SSH only |  1.72 |  1.24 |  12.4 |  11.6 |  83.7% |  53.5% |  $-$0.5% |  24.8%
U-Net |  SSH only |  1.39 |  1.22 |  9.1 |  10.3 |  89.1% |  72.3% |  $-$3.0% |  65.0%
|  SSH-SST |  1.33 |  0.90 |  4.0 |  4.2 |  92.6% |  79.4% |  19.5% |  72.0%
4DVarNet |  SSH only |  0.9 |  0.7 |  4.3 |  5.6 |  94.0% |  86.1% |  12.1% |  81.3%
(ours) |  SSH-SST |  0.76 |  0.61 |  2.7 |  2.5 |  97.4% |  92.1% |  46.9% |  87.2%
Table 1: Synthesis of the reconstruction performance of the benchmarked
approaches: we report the performance metrics of the benchmarked approaches
for the reconstruction of image time series of sea surface currents from
satellite data. We refer the reader to the main text for the description of
the different metrics. We highlight in bold the best score.

Synthesis of the benchmarking experiments: Tab. 1 compares the performance of
the benchmarked schemes. DUACS SSH-derived geostrophic currents provide a
baseline performance with reconstruction score in line with those reported for
the metrics considered for SSH mapping (fablet_multimodal_2022;
le_guillou_mapping_2020) with resolved space-time scales above 1∘ and 10 days.
This baseline accounts for more than 80% of the variance of the total current
and about 50% of its vorticity and 25% of its strain. The two 4DVarNet schemes
clearly outperform this reconstruction performance with resolved space-time
scales below 1∘ and 7 days. The relative improvement is particularly
significant when exploiting SSH-SST synergies, especially for resolved time
scales below 3 days. Interestingly, with this 4DVarNet configuration, the
metrics for the explained variance of the different fields are better than
those of the SSH-derived currents using the true SSH. In this respect, the
poorly-resolved time scales of the latter are indicators of the relative
contribution of the ageostrophic component of the total currents. Both
4DVarNet configurations learn to correct geostrophic velocities to better
recover these temporal features. Though the configuration using only altimetry
data only recover on average 12% of the divergent component of the SSC, the
SSH-SST synergies captured by the 4DVarNet through the multimodal observation
term in variational formulation (7) greatly improve the ability to retrieve a
significant fraction of this divergent component (here, $\approx 47\%$ on
average). The reported relative improvement greater than 30% for most of the
performance metrics further highlights their contribution. This benchmarking
experiment also clearly supports the relevance of 4DVarNet schemes compared
with direct learning-based inversion schemes using off-the-shelf deep learning
architecture. For instance, we report for all metrics a better reconstruction
performance of the 4DVarNet scheme using only altimetry data than that of a
direct UNet-based inversion with altimetry and SST inputs.

Fig. 3 further illustrates these results through the reconstructed fields on
November 12th 2013\. This example stresses some strengthening of the sea
surface current along the main meander compared with the SSH-derived
geostrophic velocities. This is captured by 4DVarNet reconstruction using SSH-
SST synergies. They also illustrate the improvement regarding fine-scale
patterns which are smoothed by DUACS baseline. The comparison of the vorticity
also indicates some correction of the SSH-derived estimation which is revealed
by 4DVarNet scheme. Visually, the reconstructed divergence using 4DVarNet
schemes recover important features of the true divergence, though finer-scale
patterns are lost, which is in line with the performance metrics reported in
Tab. 1.

True SSC
---
|  |  |
True SSH-derived SSC
|  |  |
DUACS-based SSH-derived SSC
|  |  |
4DVarNet SSC using SST
|  |  |
4DVarNet SSH without SST
|  |  |
SSC norm |  Vorticity |  Strain |  Divergence
Figure 3: Reconstructed SSC fields on November 12th 2012:. From left to right, we depict the norm (i.e., velocity intensity), the vorticity, the divergence and the strain of SSC fields corresponding respectively to, from top to bottom, the true SSC field, the SSH-derived one, DUACS-derived one, and 4DVarNet-derived using SSH-SST fields and using only-SSH field. |  |
---|---|---
Reconstruction error |  True strain |  SST gradient
Figure 4: Time-averaged reconstruction error vs. mean strain field and SST gradient over the whole test dataset: from left to right, we display the time-averaged mean square reconstruction error field of the SSC using the 4DVarNet scheme combining altimetry and SST observations, the time-averaged true strain field and the time-averaged SST gradient. The coefficient of determination between the error field and the later fields is respectively of 59% and 24%. The reconstructed strain field depicted in Fig.3 leads to an explained variance of 57% for the time-averaged reconstruction error.  Data |  $\lambda_{x,u}$ (∘) |  $\lambda_{x,v}$ (∘) |  $\lambda_{t,u}$ (d) |  $\lambda_{t,v}$ (d) |  $\tau_{u,v}$ |  $\tau_{\mathrm{vort}}$ |  $\tau_{\mathrm{div}}$ |  $\tau_{\mathrm{strain}}$
---|---|---|---|---|---|---|---|---
NAlt |  1.5 |  1.0 |  6.5 |  6.4 |  91.7% |  79.0% |  8.8% |  72.7%
NAlt+SWOT |  0.9 |  0.7 |  4.3 |  5.6 |  94.0% |  86.1% |  12.1% |  81.3%
NALT+SST |  0.81 |  0.62 |  2.6 |  2.5 |  97.1% |  90.8% |  46.3% |  85.0%
NAlt+SWOT+SST |  0.76 |  0.61 |  2.7 |  2.5 |  97.4% |  92.1% |  46.9% |  87.2%
Table 2: Impact of SWOT and SST data onto the reconstruction performance
w.r.t. a baseline using only nadir altimetry data: we evaluate the performance
metrics over the tests dataset for the trained 4DVarNet schemes reported in
Tab.1 when considering nadir altimetry only (NAlt), nadir and wide-swath SWOT
(NAlt+SWOT) altimetry combined or not with SST data in the proposed multimodal
framework. Bold values show the highest skill.

Impact of SWOT data: We analyse in Tab. 2 how SWOT data contribute to the
reconstruction of sea surface currents. For the trained 4DVarNet models
considered in Tab. 1, we evaluate their reconstruction performance when
considering only nadir altimetry data rather than nadir and SWOT altimetry
data as used in the baseline configuration. These results support the
relevance of wide-swath altimetry data to improve the reconstruction of sea
surface currents. We report for instance a relative improvement greater than
25% for both $\tau_{u,v}$, $\tau_{\mathrm{vort}}$ and $\tau_{\mathrm{strain}}$
metrics and improved resolved time scales from 6.1 days vs. 8.3 days for the
zonal velocity. When exploiting SSH-SST synergies, we also observed some
improvement though smaller. We observe the greatest relative improvement
between 10% and 15% for $\tau_{u,v}$, $\tau_{\mathrm{vort}}$ and
$\tau_{\mathrm{strain}}$ metrics. This emphasizes for the considered case-
study that SST fields bring most of the relevant information to retrieve the
fine-scale patterns of sea surface currents.

Impact of the spatial resolution of the SST data: As the spatial resolution of
satellite-derived SST products depends both on the satellite sensors as well
as on the atmospheric conditions, we evaluate here the impact of the spatial
resolution of the SST fields onto the reconstruction performance. From the
trained mutimodal 4DVarNet, we fine-tune mutimodal 4DVarNet models using
coarsened and subsampled versions of the original SST fields for spatial
resolutions ranging from $\nicefrac{{1}}{{10}}$∘ to $\nicefrac{{1}}{{2}}$∘. We
compare in Tab.3 the performance metrics of these models. As expected, the
lower the resolution of the SST, the lower the performance. When considering
SST data with a $\nicefrac{{1}}{{2}}$∘ resolution, we only significantly
improve the reconstruction of the divergent component of the SSC (25.4% vs.
12% when using only altimetry data). From a $\nicefrac{{1}}{{4}}$∘ resolution,
we report an improvement for all metrics, which is more noticeable for the
resolved time scales. This is of key interest for the application to real
microwave SST data (ocarroll_observational_2019). Besides, the clear gain
issued from higher-resolution data also supports the potential of
multispectral satellite sensors. Overall, Tab.3 suggests that the key SST
features used to inform sea surface currents relate to horizontal scales
between $\nicefrac{{1}}{{20}}$∘ and $\nicefrac{{1}}{{2}}$∘.

SST resolution |  $\lambda_{x,u}$ (∘) |  $\lambda_{x,v}$ (∘) |  $\lambda_{t,u}(d)$ |  $\lambda_{t,v}$ (d) |  $\tau_{u,v}$ |  $\tau_{vort}$ |  $\tau_{div}$ |  $\tau_{strain}$
---|---|---|---|---|---|---|---|---
$\nicefrac{{1}}{{20}}$∘ |  0.76 |  0.61 |  2.7 |  2.5 |  97.4% |  92.1% |  46.9% |  87.2%
$\nicefrac{{1}}{{10}}$∘ |  0.81 |  0.62 |  3.0 |  2.6 |  96.9% |  90.8% |  41.8% |  85.1%
$\nicefrac{{1}}{{5}}$∘ |  0.87 |  0.68 |  3.0 |  4.4 |  96.2% |  89.0% |  37.3% |  83.2%
$\nicefrac{{1}}{{4}}$∘ |  0.90 |  0.70 |  3.0 |  3.5 |  95.6% |  87.4% |  30.0% |  81.5%
$\nicefrac{{1}}{{2}}$∘ |  0.92 |  0.82 |  4.3 |  6.3 |  94.3% |  85.3% |  25.4% |  80.0%
Table 3: Impact of the spatial resolution of the SST observations: we report
the performance metrics of the proposed 4DVarNet models using SSH-SST
synergies for SST observations with different spatial resolutions from
$1/20^{\circ}$ to $1/2^{\circ}$. We refer the reader to the main text for the
description of the different metrics. Bold values show the highest skill.

Analysis of the reconstruction error: We further analyse the reconstruction
error and illustrate in Fig. 4 the time-averaged mean square error of the
total sea surface current over the considered dataset from October 20th 2012
to December 4th 2012 for the trained 4DVarNet scheme exploiting SSH-SST
synergies. We also report the time-averaged strain of the true velocities
along with the time-averaged amplitude of horizontal SST gradients. We observe
a better match between error and strain patterns. Quantitatively, the time-
averaged strain explains about 60% of the time-averaged error over the test
dataset. We reach similar correlation statistics for the reconstructed strain
as expected from the visual similarity between the true and reconstructed
strain depicted in Fig. 3.

We know that the strain highlights regions of frontogenesis or gradient
strengthening (balwada_vertical_2021; okubo_horizontal_1970;
weiss_dynamics_1991). It could also be demonstrated that it leads to velocity
error growth (Appendix). This is also consistent with kinetic energy transfer
occurring due to the strain of the flow (Sevellec:22). Here we hypothesized
that regions of significant error growth (i.e., strain active regions) are
more difficult to reconstruct due to their chaotic nature (i.e, small
disturbances will become a significant signal). However we do not expect a
perfect relation between surface strain and velocity error growth (and so
velocity reconstruction error) since other terms, such as viscous terms or
sub-surface pressure gradients, also control the velocities, and can become
sources of error growth. Interestingly, this analysis is similar for the true
strain and the reconstructed one, which could provide a proxy for the
quantification of the uncertainty of the reconstruction using 4DVarNet
schemes.

##  6 Discussion

This study has introduced a 4DVarNet deep learning scheme, backed on a
variational data assimilation formulation, for the reconstruction of sea
surface currents from satellite-derived SSH and SST observations. Our
numerical experiments support the relevance of this learning-based approach
over state-of-the-art schemes to retrieve finer-scale sea surface dynamics
(typically 0.5-0.7∘ and 3 days for the resolved space-time scales), including
a significant fraction of the ageostrophic component of the total current
(about 47% of the divergence of the SSC). The good skill of 4DVarNet in the
reconstruction of the horizontal divergence of the flow makes it particularly
relevant for tracking surface/buoyant polluant concentration (for which
convergence regions act as attractors (Sevellec:17b)). It is also important
for the hypothetic reconstruction of vertical velocities (which vertical
gradient is directly related to the horizontal divergence of the flow), with
all their dynamical consequences for the export of nutrients, heat, and
carbon, for instance. We have also hypothesized, through theoretical
arguments, and shown the relation between the flow strain and the error of the
reconstructed horizontal velocities. This suggests that strain active regions
are good candidates for intense monitoring if one aims to better reconstruct
horizontal velocities. We further discuss our contributions according to the
following aspects: the monitoring of sea surface velocities, deep learning
inversions for unobserved upper ocean dynamics and multimodal synergies for
ocean monitoring and forecasting.

Reconstruction of sea surface velocities: Satellite altimetry plays a major
role in our ability to monitor sea surface dynamics on a global to local
scale. Through the direct measurement of the sea level anomaly, it delivers an
estimation of the geostrophic component of sea surface velocities
(chelton_satellite_2001). As the scarce sampling of the sea surface by nadir
altimeters limits our ability to retrieve dynamical features below $\approx
1^{\circ}$ and 10 days, research effort has been undertaken to observe and
reconstruct finer-scale dynamics. In this context, upcoming swide-swath
altimetry mission SWOT (gaultier_challenge_2015) will provide the first
snapshots of the sea surface height down to a $\nicefrac{{1}}{{10}}$∘ spatial
resolution. Similarly to nadir altimeters, it will only directly inform the
geostropic component of sea surface velocities. Numerous studies (e.g.,
(baaklini_blending_2021; mahadevan_analysis_2006)) have evidenced the key role
of ageostrophic dynamics in the mesoscale-to-submesoscale range. This has
motivated a large research effort towards the exploitation of other
observation sources, alone or combined with satellite altimetry, to retrieve
sea surface dynamics, including among others SST (fablet_improving_2018;
isern-fontanet_transfer_2014; rio_improving_2016), Ocean Colour
(ciani_ocean_2021), sea surface drifters (baaklini_blending_2021;
sun_impacts_2022), and SAR observations (chapron_direct_2005). From a
methodological point of view, we may distinguish three main categories of
approaches: optimal interpolation schemes (cressie_statistics_2015;
taburet_duacs_2019), data-driven approaches (fablet_data-driven_2017;
fablet_multimodal_2022; manucharyan_deep_2021), and data assimilation scheme
using OGCM (benkiran_assessing_2021; fujii_observing_2019) or QG dynamical
priors (le_guillou_mapping_2020; ubelmann_dynamic_2014). The proposed 4DVarNet
scheme benefits, on the one hand, from a variational data assimilation
formulation to make explicit observation and dynamical priors, especially the
expected though unknown relationship between SST and SSC features, and, on the
other hand, from the computational efficiency of deep learning schemes, to
learn uncalibrated terms and solvers from data in an end-to-end manner. This
seems particularly promising to make the most of available multi-source
observation datasets for the reconstruction of sea surface currents. Our study
supports the ability to retrieve finer-scale patterns than currently achieved
by operational products. While our study points out the added-value of wide-
swath altimetry data, it suggests that ageostrophic components can be revealed
by learning-based schemes from other sea surface tracers sampled at higher-
resolution as illustrated with SST fields in our study. This likely relates to
specific dynamical regimes in play in the Gulf Stream region
(mcwilliams_gulf_2019; reul_sea_2014). Key challenges for future work include
the exploration of the proposed framework for the variety of dynamical regimes
exhibited by upper ocean dynamics as well as the transfer from OSSE to real
observation dataset. The extension to other observation datasets such as
drifter trajectories (baaklini_blending_2021; sun_impacts_2022), SAR
observations, as well as multispectral satellite sensors (barnes_cross-
calibration_2021; tilstone_performance_2021; yurovskaya_ocean_2019) also
naturally arises as appealing research directions.

Deep learning inversion for unobserved upper ocean dynamics: Our study further
supports the potential for deep learning approaches to reconstruct unobserved
or partially-observed variables or processes from available satellite-derived
and/or in situ observation datasets. The reported case-study for horizontal
sea surface currents is in line with recent studies addressing ocean
processes, such as ocean eddy heat fluxes (george_deep_2021), tide-related
features (wang_deep_2022), plankton dynamics (martinez_neural_2020)… The
underlying assumption is that deep learning can disentangle the features
associated with a given process in the observations of a tracer. In this
context, classic data assimilation scheme requires defining explicitly the
observations and dynamical operator in play in the underlying state-space
formulation (5). Regarding upper ocean dynamics, this generally resorts to
considering some discretization of the primitive equations and identity
observation operators for ocean state variables and ocean-atmosphere
interactions. The complexity of the inversion problem may hinder the ability
of these schemes to retrieve fine-scale dynamics when dealing with a poorly-
constrained observation setting as with the inversion of ocean dynamics from
satellite-derived sea surface observations. By contrast, the proposed deep
learning scheme focuses only on the ocean variable of interest, here observed
and targeted sea surface variables. Through this learning paradigm, we reduce
the complexity of the inversion problem and explore the complex relationships
exhibited by sea surface tracers governed by upper ocean dynamics. We believe
this approach to be generic and to show a potential to improve our
understanding and monitoring of upper ocean dynamics which will remain
scarcely observed.

OSSEs and Deep Learning to Bridge ocean modeling and multimodal ocean remote
sensing: We may also regard the proposed deep learning schemes as alternative
means to bridge ocean modeling and ocean observation. Through the combination
of simulation dataset for ocean dynamics and observing systems, OSSE settings
(boukabara_community_2018; fujii_observing_2019) have been widely exploited as
means to assess the potential impact of new observing systems to improve the
reconstruction and/or forecasting of ocean dynamics. OSSEs in the context of
the upcoming SWOT mission fall into this category to evaluate the added-value
of wide-swath altimetry data to inform upper ocean dynamics
(benkiran_assessing_2021; ubelmann_dynamic_2014), as further supported by this
study. More recently, OSSEs have provided benchmarking frameworks to
intercompare the performance of different approaches for the same task such as
space-time interpolation problems for sea surface dynamics
(fablet_multimodal_2022; le_guillou_mapping_2020; fujii_observing_2019;
vient_data-driven_2021). As illustrated here, when dealing with unobserved or
scarcely observed processes, OSSEs also provide means to train a model from
simulation datasets with a view to applying this model to real datasets.
Previous learning-based studies dedicated to the reconstruction of ocean’s
interior processes, such as for instance the primary production
(puissant_inversion_2021) and deep geostrophic currents
(manucharyan_deep_2021), provide such examples. These studies generally learn
a mapping from observed gap-free sea surface tracers to the targeted ocean
state dynamics. Our study supports the ability to further extend this approach
to reconstruction or mapping objectives from irregularly-sampled observations.
The applicability of this general framework to real observation datasets
clearly depends on the quality of both the numerical simulations of ocean
dynamics and of the simulations of the observing systems. As such, it
advocates to pursue joint research effort into these simulation issues along
with the design and evaluation of deep learning approaches. We believe these
aspects to be critical for the development of multimodal observing systems
combining satellite-derived observations and in situ data including among
others sea surface drifters (baaklini_blending_2021; sun_impacts_2022), argo
profiles (cossarini_towards_2019; roemmich_future_2019), and underwater
acoustics data (storto_assessing_2020; storto_neural_2021)…

## Appendix: Velocity error growth

Starting from the momentum and nondivergence equations [described in (1d)] but
splitting the velocity (and pressure and viscous terms) in a targeted truth
and a small error ($X=\bar{X}+X^{\prime}$ with $|\bar{X}|\ll|X^{\prime}|$,
where $X$ is any variable) we have:

$D_{t}u^{\prime}-fv^{\prime}=-u^{\prime}\partial_{x}\bar{u}-v^{\prime}\partial_{y}\bar{u}-w^{\prime}\partial_{z}\bar{u}-\frac{1}{\rho_{0}}\partial_{x}P^{\prime}+\mathcal{F}_{x}^{\prime},$
(13a)
$D_{t}v^{\prime}+fu^{\prime}=-u^{\prime}\partial_{x}\bar{v}-v^{\prime}\partial_{y}\bar{v}-w^{\prime}\partial_{z}\bar{v}-\frac{1}{\rho_{0}}\partial_{y}P^{\prime}+\mathcal{F}_{y}^{\prime}.$
(13b)
$0=-\frac{1}{\rho_{0}}\partial_{z}P^{\prime}-\frac{g}{\rho_{0}}\rho^{\prime},$
(13c)
$\partial_{x}u^{\prime}+\partial_{y}v^{\prime}+\partial_{z}w^{\prime}=0.$
(13d)

So that error growth reads:

$D_{t}\left(u^{\prime 2}+v^{\prime
2}\right)=-\left(u^{\prime},v^{\prime}\right)\left(\begin{array}[]{cc}\partial_{x}\bar{u}&\partial_{y}\bar{u}\\\
\partial_{x}\bar{v}&\partial_{y}\bar{v}\end{array}\right)\left(\begin{array}[]{c}u^{\prime}\\\
v^{\prime}\end{array}\right)\\\
-u^{\prime}w^{\prime}\partial_{z}\bar{u}-v^{\prime}w^{\prime}\partial_{z}\bar{v}-\frac{g}{\rho_{0}}w^{\prime}\rho^{\prime}\\\
-\frac{1}{\rho_{0}}\left[\partial_{x}\left(u^{\prime}P^{\prime}\right)+\partial_{y}\left(v^{\prime}P^{\prime}\right)+\partial_{z}\left(w^{\prime}P^{\prime}\right)\right]+\left(u^{\prime}\mathcal{F}_{x}^{\prime}+v^{\prime}\mathcal{F}_{y}^{\prime}\right),$
(14)

where the left handside is the evolution of the error (or error growth), the
first line of the right handside is the horizontal transfer of momentum error,
the second line of the right handside is the baroclinic transfer of momentum
and density errors, and the third line of the right handside is the work of
the pressure and viscous force error. Here the horizontal transfer of momentum
error is a scalar product of the erroneous velocities applied to the operator
of horizontal targeted-truth-velocity gradient. This could be transformed as:

$\displaystyle\left(u^{\prime},v^{\prime}\right)\mbox{\boldmath$\mathsf{\Sigma}$}\left(\begin{array}[]{c}u^{\prime}\\\
v^{\prime}\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\partial_{x}\bar{u}&\partial_{y}\bar{u}\\\
\partial_{x}\bar{v}&\partial_{y}\bar{v}\end{array}\right)\left(\begin{array}[]{c}u^{\prime}\\\
v^{\prime}\end{array}\right),$ (21) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(u^{\prime},v^{\prime}\right)\left(\begin{array}[]{cc}\Delta+\sigma_{n}&\sigma_{s}+\zeta\\\
\sigma_{s}-\zeta&\Delta-\sigma_{n}\end{array}\right)\left(\begin{array}[]{c}u^{\prime}\\\
v^{\prime}\end{array}\right),$ (26) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(u^{\prime},v^{\prime}\right)\left(\begin{array}[]{cc}\Delta+\sigma_{n}&\sigma_{s}\\\
\sigma_{s}&\Delta-\sigma_{n}\end{array}\right)\left(\begin{array}[]{c}u^{\prime}\\\
v^{\prime}\end{array}\right),$ (31)

where $\Delta=\partial_{x}\bar{u}+\partial_{y}\bar{v}$ is the divergence,
$\zeta=\partial_{x}\bar{v}-\partial_{y}\bar{u}$ is the horizontal vorticity,
$\sigma_{n}=\partial_{x}\bar{u}-\partial_{y}\bar{v}$ is the normal strain (or
stretching), and $\sigma_{s}=\partial_{x}\bar{v}+\partial_{y}\bar{u}$ is the
shear strain (or shearing). This last expression shows the symmetric part of
the gradient deformation operator, which is virtually equivalent to the one
used for tracer gradient growth (balwada_vertical_2021) in the Okubo-Weiss
framework of frontogenesis (okubo_horizontal_1970; weiss_dynamics_1991). It is
interesting to note that here the horizontal vorticity ($\zeta$) does not
contribute to the error growth, though. $\mathsf{\Sigma}$ can be diagonalized
to show the natural growth of the error, which leads to the pair of
eigenvalues:
$\lambda_{\pm}=\nicefrac{{\Delta}}{{2}}\pm\nicefrac{{\sigma}}{{2}}$, where
$\sigma$ is the strain (such as $\sigma^{2}=\sigma_{n}^{2}+\sigma_{s}^{2}$).
Since the horizontal flow is largely strain dominated, this shows that the
error growth is controlled by the strain. This result is consistent with
energy transfer consideration suggesting the transfer from the mean state
(here the targeted-truth) to the eddy field (here the error) following the
strain of the flow (Sevellec:22). However note that the total error growth
will also be affected by baroclinic error transfers and pressure and viscous
work.

##  7 Open Research

We distribute through open source license the datasets considered in the
reported experiments, our pytorch implementation and trained pytorch models
(febvre_pytorch_2022).

###### Acknowledgements.

This work was supported by LEFE program (LEFE MANU and IMAGO projects IA-OAC
and ARVOR, respectively), CNES (OSTST DUACS-HR and SWOT ST DIEGO) and ANR
Projects Melody (ANR-19-CE46-0011) and OceaniX (ANR-19-CHIA-0016). It
benefited from HPC and GPU resources from Azure (Microsoft Azure grant) and
from GENCI-IDRIS (Grant 2021-101030).