# Improved global sea surface height and currents maps from remote sensing and in situ observations

Maxime Ballarotta, Clément Ubelmann, Pierre Veillard, Pierre Prandi, Hélène Etienne, Sandrine Mulet, Yannice Faugère, Gérald Dibarboure, Rosemary Morrow, Nicolas Picot

► To cite this version:

Maxime Ballarotta, Clément Ubelmann, Pierre Veillard, Pierre Prandi, Hélène Etienne, et al. Im-
proved global sea surface height and currents maps from remote sensing and in situ observations.
Earth System Science Data: Papers in open discussion, 2022, 10.5194/essd-2022-181 . hal-04631712

HAL Id: hal-04631712

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

Submitted on 2 Jul 2024

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
<figure><img src="image_2.png" /></figure>

# Improved global sea surface height and currents maps from remote sensing and in situ observations

Maxime Ballarotta¹, Clément Ubelmann², Pierre Veillard¹, Pierre Prandi¹, Hélène Etienne¹, Sandrine Mulet¹, Yannice Faugère¹, Gérald Dibarboure³, Rosemary Morrow⁴ & Nicolas Picot³

¹Collecte Localisation Satellites, 31520 Ramonville-Saint-Agne, France

²Datlas, 38400 Saint Martin d'Hères, France

³Centre National d'Études Spatiales, 31400 Toulouse, France

Corresponding author: M.Ballarotta (mballarotta@groupcls.com)

## Abstract.

We present a new gridded sea surface height and current dataset produced by combining observations from nadir altimeters and drifting buoys. This product is based on a multiscale & multivariate mapping approach that offers the possibility to improve the physical content of gridded products by combining the data from various platforms and in resolving a broader spectrum of ocean surface dynamic than in the current operational mapping system. The dataset covers the entire global ocean and spans from 2016-07-01 to 2020-06-30. The multiscale approach decomposes the observed signal into different physical contributions. In the present study, we simultaneously estimate the mesoscale ocean circulations as well as part of the equatorial wave dynamics (e.g., tropical instability and Poincaré waves). The multivariate approach is able to exploit the geostrophic signature resulting from the synergy of altimetry and drifter observations. Sea level observations in Arctic leads are also used in the merging to improve the surface circulation in this poorly mapped region. A quality assessment of this new product is proposed against the DUACS operational product distributed in the Copernicus Marine Service. We show that the multiscale & multivariate mapping approach offers promising perspectives for reconstructing the ocean surface circulation: leads observations contribute to improve the coverage in delivering gap free maps in the Arctic; drifters observations help to refine the mapping in regions of intense dynamics where the temporal sampling must be accurate enough to properly map the rapid mesoscale dynamics; overall, the geostrophic circulation is better mapped in the new product, with mapping errors significantly reduced in regions of high variability and in the equatorial band; the effective resolution of this new product is hence between 5% and 10% finer than the Copernicus product.

## 1 Introduction

Several oceanographic applications (e.g., operational oceanography, marine weather, climate monitoring...) rely on high-quality observational datasets. The European Union (E.U.) Copernicus Marine & Climate Change Services provide
<figure><img src="image_3.png" /></figure>

operational services and indicators on the observed state of the climate. Sea level and surface currents are, among others, key variables distributed by the services. There are also listed as Essential Climate Variables (ECVs) for the detection of climate change and the characterization of climate system variability (Bojinski et al., 2014).

As part of the Copernicus Services, the Thematic Assembly Centres (SL-TAC & MOB-TAC) deliver near-real time and delayed time sea level and surface currents products (satellite and in situ level 3 and 4 products) that are used by the ocean science community to study, understand and monitor the evolution of the ocean system. These products do not resolve the entire spectrum of the ocean surface variability; they have resolution limits of about 60 km for the along-track products (Dufau et al., 2016) and >200 km x 20 days for the Level 4 gridded products (Ballarotta et al., 2019), but recent nadir altimetry instruments, such as the new Sentinel-3A and 3B SAR missions, or future missions based on large swath technologies (e.g., the upcoming SWOT mission) offer, for example, the possibility of observing finer ocean structures (Morrow et al., 2019) which could be used to provide better gridded product resolution.

In addition, the growing needs to develop observing systems or methods with finer spatial scales / higher frequencies have been identified by the ocean scientific community and the Copernicus Services as R&D priorities to serve Copernicus marine users and decision-makers (see, e.g., Abdalla et al., 2021, or the “Copernicus Marine Service Evolution Strategy: R&D priorities - Version 5 June 30, 2021” document, https://marine.copernicus.eu/sites/default/files/media/pdf/2021-09/CMEMS%20Service_evolution_strategy RD priorities v5-June-2021.pdf, last-access: 20211009). Therefore, with the support of the French Space Agency (CNES), the development of new experimental products has been undertaken, aiming at improving the resolution of the current Level-3 and Level-4 Sea level products (Mulet et al. 2019, Ballarotta et al., 2020, Ubelmann et al., 2021a, Prandi et al., 2021, Pujol et al., in prep.) and preparing operational systems for the SWOT era (Ubelmann et al., 2015, Ubelmann et al., 2021b, Le Guillou et al., 2021, Beauchamp et al., 2020).

The present study focuses on the development and assessment of experimental global gridded products based on a recent multiscale & multivariate mapping approach (Ubelmann et al., 2021a) and applied to real Earth observations. We here investigate the possibility of improving the content of gridded products in combining the data from various platforms (in situ and satellite) and in resolving a larger spectrum of the ocean surface dynamic than in current operational products.

The paper is structured as follows: the data sources and merging methods used in this study are described in section 2. Section 3 presents the experiments and validation metrics. The quality assessment of the new products is proposed in section 4. The key results are then summarized in section 5.

# 2 Data & Methods

## 2.1 Data sources

The mapping method used in this study takes input data from remote sensing or in situ observations, which are summarized in Table 1 and described below.
<figure><img src="image_4.png" /><figcaption>Table 1: List of observation datasets used in this study</figcaption></figure>

<table><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>Product type</td><td>Global Altimeter SLA products</td><td>Arctic leads Altimeter SLA products</td><td>Drifters' geostrophic velocity product</td></tr><tr><td>Product ref.</td><td>SEALEVEL_GLO PHY L3_REP OBSERVATIONS_008_062</td><td>Experimental</td><td>AOML</td></tr><tr><td>Spatial coverage</td><td>[0:360] [-90:90]</td><td>&gt;60°N</td><td>[0:360] [-90:90]</td></tr><tr><td>period</td><td>20160115-20200630</td><td>20160115-20200630</td><td>20160115-20200630</td></tr></tbody></table>

## 2.1.1 Sea level anomaly products

The global ocean SSH observations are from the (Delayed-Time DT) Level-3 altimeter satellite along-track data, reprocessed in 2021 and distributed by the E.U. Copernicus Marine Service (product reference SEALEVEL_GLO PHY L3_MY_008_062, https://doi.org/10.48670/moi-00146). These data cover the period 19930101-20201231 over the world ocean (excluding ice-covered areas, e.g., Figure 1) and are available at a sampling rate of 1Hz (~7km spatial spacing). Homogenisation and cross-validation are applied to the dataset to remove any residual orbit error, long-wavelength error (lwe), large-scale biases and discrepancies between different data streams. The list of geophysical and environmental corrections applied to the data sets is described in the Quality Information Document (Taburet et al., 2021) and summarized below in Equation (1). In this study, unfiltered sea level anomalies (SLA) corrected with dynamic atmospheric correction (dac), ocean tide and lwe corrections are considered in the multi-scale/multivariate mapping.

$$SLA = \text{Orbit} - \text{Range} - \sum (\text{Environmental Corrections}) - \sum (\text{Geophysical Corrections}) - \text{Mean Sea Surface} \quad (1)$$

with $\sum (\text{Environmental Corrections}) = \text{wet tropospheric} + \text{dry tropospheric} + \text{ionospheric} + \text{sea-state-bias}$, $\sum (\text{Geophysical Corrections}) = \text{solid earth tide} + \text{load tide} + \text{ocean tide} + \text{pole tide} + \text{dynamic atmospheric correction}$. The Mean Sea Surface used here is the CNES-CLS18 (Mulet et al., 2021).
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 1: Example of sea level altimetry coverage for a 7-day period (20190701-20190707). Colour scale represents the sea level anomaly amplitude in meters</figcaption></figure>

### 2.1.2 Sea level anomaly products in artic leads

To fill polar gaps and improve SSH maps in ice-covered areas, specific Arctic Level-3 SSH products are also used in the mapping. These products cover the Arctic region (up to 88°N) and are available at a sampling rate of 20Hz (~350 metres) for 3 altimetry missions: SARAL/AltiKa, Sentinel-3A and CryoSat-2 (Figure 2 & Table 2). The Level-2 data processing is described in Prandi et al, (2021). To ensure continuity with the open ocean, the corrections are derived from the global Level-3 along-track processing (Taburet et al., 2021) when possible. The noticeable exceptions concern 1) the wet tropospheric correction that comes from ECMWF model since on-board radiometer estimates are not reliable over ice, 2) the sea state bias correction is not applied since waves and winds are considered small over leads, 3) orbit error corrections are not applied as they are difficult to compute over this small region. Then a constant bias of ~8cm is applied for each mission to ensure continuity with DUACS DT2021 open ocean SLA.

<table><caption>Table 2: Arctic leads product characteristics</caption><thead><tr><th>Altimeter</th><th>SARAL/AltiKa</th><th>Sentinel-3A</th><th>CryoSat-2</th></tr></thead><tbody><tr><td>Latitude max.</td><td>81,5°</td><td>81,5°</td><td>88°</td></tr><tr><td>Retracking</td><td>Adaptive (LRM)</td><td>TFMRA 50% (SAR)</td><td>TFMRA 50% (SAR)</td></tr></tbody></table>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 2: Example of arctic leads sea level altimetry coverage for a 7-day period (20190701-20190707). Colour scale represents the sea level anomaly amplitude in meters</figcaption></figure>

### 2.1.3 Geostrophic current anomaly products

We used delayed-time horizontal surface velocities from the AOML Surface Velocity Program (SVP, Lumpkin and Centurioni, 2019). The data cover the entire world ocean and are available at a 6-hour frequency. SVP are designed to follow the 15m depth circulation, which is the centre depth of their drogues. When the drogue is lost, they follow the surface current, but are also under the direct influence of the wind.

AOML distributes a flag to indicate whether the drogue is lost or not (Lumpkin et al, 2013). These data are also distributed by the INSITU Thematic Assembly Centre of the E.U. Copernicus Marine Service (see Product User Manual, http://marine.copernicus.eu/documents/PUM/CMEMS-INS-PUM-013-044.pdf) with an additional wind slippage correction for undrogued buoys derived from the Rio (2012) methodology. For the study, the undrogged and drogued drifters are selected over the global ocean and the period from 20160601 to 20200731. Note that for specific experiments described hereafter, we excluded drifters' trajectories between -10°S and 10°N (e.g., Figure 3) to isolate and evaluate only the impact of the equatorial wave's mode in this region. In this study, we computed the geostrophic velocity anomaly components, which are defined as:

$$U_{anom} = U_{buoy} - U_{ekman} - U_{stokes} - U_{inertial} - U_{tidal} - U_{ahf} - U_{slip} - U_{mdt} \quad (2)$$

$$V_{anom} = V_{buoy} - V_{ekman} - V_{stokes} - V_{inertial} - V_{tidal} - V_{ahf} - V_{slip} - V_{mdt} \quad (3)$$

120
<figure><img src="image_9.png" /></figure>

With $U_{buoy}$ ($V_{buoy}$) is the drifter's zonal (meridional) velocity. Each component is corrected from:

- the wind-driven component $U_{ekman}$ ($V_{ekman}$) using an update of the model used in Mulet et al (2021) and described in Etienne (2021). In this recent version, ERA5 wind stress (Hersbach et al, 2018) replaces the ERAinterim data and the equatorial symmetry of the wind driven parameters is removed.

- The Stokes drift $U_{stokes}$ ($V_{stokes}$) from ERA5 reanalysis (Hersbach et al, 2018) is also removed from the surface drifter velocity (undrogued drifters). No Stokes drift is removed from the 15m depth velocity, as this component is supposed to mostly vanish in the first 2-4 m.

- The wind slippage which is the direct effect of the wind on the buoy $U_{slip}$ ($V_{slip}$). This correction is significant only in the case of drogue loss (Etienne et al, 2021), when the drifters are advected by the surface current.

Then the data is filtered from the tidal and inertial velocities $U_{inertial}+U_{tidal}$ ($V_{inertial}+V_{tidal}$) as well as the residual high frequency ageostrophic signal $U_{ahf}$ ($V_{ahf}$). Finally, the mean geostrophic velocity (CNES-CLS2018, Mulet et al., 2021) $U_{mdt}$ ($V_{mdt}$) is subtracted to obtain the geostrophic velocity anomaly.

Drifters velocity

<figure><img src="image_10.png" /><figcaption>Figure 3: Example of drifter's trajectories coverage for the 20190101 to 20191231 period. Colour scale represents the velocity amplitude in m.s⁻¹</figcaption></figure>

## 2.2 Methods

Two mapping methods are compared in this study: the operational DUACS mapping approach and the experimental MIOST multiscale & multivariate mapping approach. Each method is briefly described in the following sections.
<figure><img src="image_11.png" /></figure>

## 2.2.1 The DUACS mapping approach

The DUACS mapping approach constructs a SSH field on a regular grid by combining measurements from various altimeters. It is based on a global suboptimal space-time objective analysis that considers along-track correlated errors as described for instance in Ducet et al., (2000) or Le Traon et al. (2003). The mathematical formulation, known as Optimal Interpolation, is described hereafter.

We assume a state to estimate, denoted x, and partial observations, denoted y, which can be related to the state by a linear operator H such as:

$$y = Hx + \epsilon \quad (4)$$

where $\epsilon$ is an independent signal (e.g., observation error) not related to the state. If we define B the covariance matrix of x and R the covariance matrix of $\epsilon$, both variables being assumed Gaussian, then the linear estimate is written:

$$x^a = BH^T(HBH^T - R)^{-1}y \quad (5)$$

The observation vector y represents the SLA observations. The state vector x is the gridded SLA. The operator H (formally a tri-linear interpolator transforming the gridded state SLA to the equivalent along-track SLA) is not considered explicitly. The matrices $BH^T$ and $HBH^T$, representing the covariance of the signal in the (grid, obs) and (obs, obs) spaces, are directly written with the analytical formula of the Arhan and Colin de Verdière (1985) covariance model as described in Ducet et al., (2000), Le Traon et al. (2003) or Pujol et al. (2016). This covariance model is mainly optimized for mesoscale signal reconstruction. The R matrix represents the representativity and instrumental errors. Since the covariance of mesoscale SLA is assumed to vanish beyond a few hundreds of kilometres in space and beyond 10–20 days in time (Le Traon & Dibarboure, 2002), separate inversions are performed locally selecting observations over time and space windows adjusted to these values. In practice, since the number of observations is limited to less than 1000, the inversion in observation space is computationally manageable. More details on the map production are given in Pujol et al. (2016).

In DUACS, the geostrophic current ($U_g, V_g$) is then directly derived from the mapped SSH:

$$U_g(x, y) = -\frac{g \partial SSH(x, y)}{f_c \partial y} \quad (6)$$

$$V_g(x, y) = \frac{g \partial SSH(x, y)}{f_c \partial x} \quad (7)$$

where g is the gravity, $f_c$ is the Coriolis frequency, which is a function of latitude.

## 2.2.2 A multiscale & multivariate mapping approach

The Optimal Interpolation requires the inversion of a matrix of the same size as the observation vector y. When the number of observations exceeds the size of the state to resolve, it can be interesting to use an equivalent formulation given by the
<figure><img src="image_12.png" /></figure>

Sherman-Morrison-Woodbury transformation, allowing an inversion in state space, with a matrix of the size of the state vector x,

$$x^a = (H^T \llbracket R^{-1}H + B^{-1} \rrbracket)^{-1} H^T R^{-1} y \quad (8)$$

The formulation of the multiscale & multivariate mapping algorithm is detailed in Ubelmann et al. (2021). We here recall the main principle. We consider an extended state vector x composed by N physical components that will be later assumed independent (in this study N=3 for 1) geostrophy and equatorial waves: 2) Tropical Instability Waves (TIW) and 3) Poincaré waves):

$$x = (x_1^T, \dots, x_N^T)^T \quad (9)$$

Each component $x_k$ represents the state of the surface topography and surface current to be resolved in the grid space, noted $x_k = (h_k^T, u_k^T, v_k^T)^T$. The key aspect of the method is a rank reduction of the state vector, through a subcomponent decomposition, such as $x_k$ can be written as:

$$x_k = \begin{bmatrix} \Gamma_{k,h} \\ \Gamma_{k,u} \\ \Gamma_{k,v} \end{bmatrix} \eta_k = \Gamma_k \eta_k \quad (10)$$

where $\eta_k$ is the reduced state vector for component k, $\Gamma_{k,h}$, $\Gamma_{k,u}$, and $\Gamma_{k,v}$ are the subcomponent matrices expressed in topography and currents, respectively. Note that for some components, one of the blocks can be set to zeros (e.g., if geostrophy component is considered with zero contribution on SSH, which is the case for the equatorial wave components). Their concatenation is called $\Gamma_k$ which is the matrix transforming the reduced state vector in the grid space for topography and currents. In practice, $\Gamma_k$ will be a wavelet decomposition of the time-space domain, with elements of appropriate temporal and spatial scales to represent the component k. These wavelet scales, and their specified variance set with a diagonal matrix noted $Q_k$, will define the equivalent covariance model $B_k$ in the grid space for component k:

$$B_k = \Gamma_k Q_k \Gamma_k^T \quad (11)$$

The observation vector y is also extended to the observed surface topography and surface current noted y = ($h^o^T, u_r^o^T)^T$. Then, if $H_k$ is the observation operator for component k (from grid space to observation space), we note $G_k = H_k \Gamma_k$ the subcomponent matrix expressed in observation space. In these conditions, the observation vector y is the sum of all component contributions plus the unexplained signal $\epsilon$ (instrument error and representativity),

$$y = \sum_{k=1}^{N} G_k \eta_k + \epsilon \quad (12)$$
<figure><img src="image_13.png" /></figure>

If we use the notation $\eta = (\eta_1^T, ..., \eta_k^T)^T$ for the concatenation of the subcomponent state vectors, and $G = (G_1, ..., G_N)$, then we have,

$$y = G\eta + \epsilon \qquad (13)$$

Applying the same transformation from Equation 4, Equation 5, and Equation 8 to the reduced state vector $\eta$, the global solution is written:

$$\eta^a = (G^T R^{-1} G + Q^{-1})^{-1} G^T R^{-1} y \qquad (14)$$

where $Q$ is the covariance matrix of $\eta$, expressed as the concatenation of the diagonal matrices $Q_k$ for each component. Finally, the solution in the reduced-space projects into the grid space with the following relation:

$$x_a = \Gamma \eta^a \qquad (15)$$

In practice, to solve Equation 14, each block of $G$ is directly filled from the analytical expression of the reduced-space elements constituting the columns of the matrix. Also, in many situations, the $(G^T R^{-1} G + Q^{-1})$ matrix, noted A hereafter, would be too large to be inverted (as required by Equation 14 explicitly). We use a preconditioned conjugate gradient method to solve $\eta = A^{-1}z$ where $z = G^T R^{-1} y$ is computed initially from $G$ and the observation vector $y$. The algorithm involves many iterations of $A\eta$ computations for updated $\eta$ until convergence is reached (when $A\eta$ approaches $z$). Note that if $A$ is too large to be written explicitly, the result $A\eta$ can still be computed in two steps from a matrix multiplication of $G$ then of $G^T$. Once the solution $\eta$ is obtained, the projection in physical grid space given by Equation 15 is applied sequentially, by summing the analytical expression of the ripples applied to grid coordinates (the columns of $\Gamma$), separately for each component $k$. As in any inversion based on linear analysis, the result strongly relies on the choice of covariance models, here defined by the reduced elements of each component. The choices of these elements are discussed in the following section focussing on the components assessed in this study: the geostrophy and equatorial wave components.

### Geostrophy component

Geostrophy is the component that has a signature on both topography and currents, and on which some synergy between altimetry and drifter observations can be expected. We define here the gridded variable $H_1$ to resolve, and the corresponding gridded geostrophic current field $(U_1, V_1)$ writes

$$\left\{ \begin{array}{l} U_1 = -\frac{g}{f_c} \frac{\partial H_1}{\partial y} \\ V_1 = \frac{g}{f_c} \frac{\partial H_1}{\partial x} \end{array} \right. \qquad (16)$$
<figure><img src="image_14.png" /></figure>

The proposed reduced state for geostrophy is based on an element decomposition of $H_1$, expressed by $\Gamma_{1,h}$ with wavelets of various wavelength and temporal extensions. This will allow to approximate the standard covariance models used in altimetry mapping, accounting for specific variations with wavelength and time. A given p element of the decomposition $\Gamma_{1,h}$ is expressed as follows:

$$ \Gamma_{1,h}[i, p] = \cos(k_{x,p}(x_i - x_p) + k_{y,p}(y_i - y_p) + \Phi_p) * f_{tap} \left( \frac{x_i - x_p}{L_{x_p}}, \frac{y_i - y_p}{L_{y_p}}, \frac{t_i - t_p}{L_{t_p}} \right) \quad (17) $$

235

where the ith line of the matrix stands for a given grid index of coordinates $(x_i, y_i, t_i)$. For the ensemble of p, $\Phi_p$ is alternatively 0 and $\pi/2$, such as all subcomponents are defined by pairs of sine and cosine functions to allow the phase degree of freedom. $k_{x,p}$ and $k_{y,p}$ are zonal and meridional wavenumbers respectively, set to vary in the mappable mesoscale range (between 80 km and 900 km with a spacing inversely proportional to the wavelet extensions, allowing to represent a signal of any intermediate wavelength). $(x_p, y_p, t_p)$ are the coordinates of a space-time pavement. The function $f_{tap}$ localizes the subcomponent in time and space (at scales $L_{t_p}$, $L_{x_p}$ and $L_{y_p}$, respectively) as geostrophy has local extension of covariances. It is expressed as:

240

$$ f_{tap}(\delta x, \delta y, \delta t) = \begin{cases} \cos(\frac{\pi}{2}\delta x) \cos(\frac{\pi}{2}\delta y) \cos(\frac{\pi}{2}\delta y), & \text{for } (|\delta x|, |\delta y|, |\delta y| < (1,1,1)) \\ 0, & \text{elsewhere} \end{cases} \quad (18) $$

245

In practice, $L_{x_p}$ and $L_{y_p}$ will be set to 1.5 the wavelength of element p and $L_{t_p}$ to the decorrelation time scale. Then, the same element p of the decomposition has also an expression in geostrophic current (through the geostrophic relation Equation 16) written in the $\Gamma_{1,u}$ and $\Gamma_{1,v}$ matrices:

250

$$ \left\{ \begin{aligned} \Gamma_{1,u}[i, p] &= -\frac{g \partial \Gamma_{1,h}[i, p]}{f_c \partial y_i} \\ \Gamma_{1,v}[i, p] &= \frac{g \partial \Gamma_{1,h}[i, p]}{f_c \partial x_i} \end{aligned} \right. \quad (19) $$

255

The whole time-space domain is paved with similar subcomponents, along coordinates $(x_p, y_p, t_p)$ for wavelengths between 80 km and 900 km spanning in all directions of the plan. The ensemble can be seen as a wavelet basis. Finally, each subcomponent p is assigned an expected variance in the $Q_1$ matrix, consistent with the power spectrum observed from altimetry at the corresponding wavelength with isotropy assumption. For a given point i on the time-space grid (210°E, 5°N), the representer $\Gamma_{1,h}$ [i,:]$Q_1\Gamma_{1,h}$ is plotted on Figure 4, shown as a function of space (bottom left panel) and as a function of time (top left panel). It illustrates the equivalent covariance function, which is quite like what is currently used for altimetry mapping with OI inverted in observation space.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 4: Example of spatio-temporal covariance models at (210°E, 5°N) for a), c) the geostrophy component (left panel) and for b), d) a westward propagating wave component, e.g., TIW (right panel)</figcaption></figure>

Equatorial waves component

We define here the gridded variable $H_2$ to resolve and we consider no contributions of the equatorial wave components on the geostrophic currents, therefore the corresponding gridded geostrophic current field ($U_2, V_2$) writes: $U_2=0, V_2=0$. The reduced state is represented in the time-space domain by the following $\Gamma_{2,h}$ matrix:

$$ \Gamma_{2,h}[i, p] = \cos(\omega_{t,p}(t_i - t_p) - k_{x,p}(x_i - x_p)) * f_{tap} \left( \frac{x_i - x_p}{L_{x_p}}, \frac{y_i - y_p}{L_{y_p}}, \frac{t_i - t_p}{L_{t_p}} \right) \quad (20) $$

where $k_{x,p}$ refers to the zonal wavenumber, and $\omega_{t,p}$ is the frequency which satisfies the dispersion relation of the wave component, e.g.:

$$ \left\{ \begin{array}{l} \omega_{t,p} = c_{wave} \cdot k_{x,p} \quad \text{for the TIW, } c_{wave} = -0.5 \, \text{m} \cdot \text{s}^{-1} \\ \omega_{t,p} = \sqrt{k_{k,p}^2 \cdot c_{wave}^2 + \beta \cdot c_{wave} \cdot (2 \cdot n + 1)} \quad \text{for the Poincaré waves, } c_{wave} = \pm 2.8 \, \text{m} \cdot \text{s}^{-1} \end{array} \right. \quad (21) $$
<figure><img src="image_17.png" /></figure>

With $c_{wave}$ the wave propagation speed (the sign indicating the direction of propagation, negative for westward, positive for eastward), $\beta$ the meridional gradient of the Coriolis frequency $f_c$, and $n = 1,2,3$ ....

In the present study, we chose $L_{tp} = 20$ days, $L_{xp} = 500$ km and $L_{yp} = 300$ km for the TIW component; $L_{tp} = 5$ days, $L_{xp} = 1000$ km and $L_{yp} = 300$ km for equatorial Poincaré wave component. As for the geostrophy component, the function $f_{tap}$ localizes the subcomponent in time and space (at scales $L_{tp}$, $L_{xp}$ and $L_{yp}$, respectively).

For a given point i on the time-space grid (210°E, 5°N), the representer $\Gamma_{2,h}$ [i,:]$Q_2\Gamma_{2,h}$ for a westward propagation wave like TIW is plotted on Figure 4, shown as a function of space (bottom right panel) and as a function of time (top right panel). It illustrates the equivalent covariance function for a westward propagation wave like TIW. Note that for Poincaré waves, both eastward and westward propagation are considered.

# 3 Experiments and validation metrics

## 3.1 Experiments

We produced 4 years (2016-07-01 to 2020-06-30) of SSH maps using the MIOST multiscale & multivariate approach by combining the Level-3 altimeter dataset from SARAL/Altika, Envisat, Jason-1, Jason-2, Jason-3, Cryosat-2, Haiyang-2A, Haiyang-2B, Sentinel-3A, Sentinel-3B missions, the Level-3 arctic leads sea-level anomaly products from SARAL/AltiKa, Sentinel-3A and CryoSat-2 missions and geostrophic current anomaly data from AOML drifter database.

Specific maps were also made to quantitatively assess the quality of the MIOST products: SSH maps for the period 2019-01-01 to 2019-12-31 were produced with the available altimeter (Jason-3, Cryosat-2, Sentinel-3A, Sentinel-3B, Haiyang-2A, Haiyang-2B) and drifter data, excluding one along-track altimeter (Saral/Altika, over open ocean region) and 20% of the drifter trajectories from the mapping to perform independent assessments. Note that for these specific maps, drifter trajectories between -10°S and 10°N (e.g., Figure 3) were excluded to evaluate only the impact of the equatorial wave's mode in this region. A twin experiment was also conducted with the optimal interpolation method used by the DUACS operational system, allowing the comparison between the multi-scale approach and the operational method chosen in the Copernicus Marine Service. Table 3 summarises the list of experiments conducted in this study, specifying the input data used in the mapping and the physical content of the maps.

<table><caption>Table 3: List of mapping experiments with the input data and physical content considered</caption><thead><tr><th rowspan="2">EXP</th><th colspan="3">Input data</th><th colspan="2">Physical content</th></tr><tr><th>altimeter</th><th>drifters</th><th>L3 arctic</th><th>geostrophy</th><th>equatorial waves</th></tr></thead><tbody><tr><td><strong>EXP01</strong>: DUACS allsat-1</td><td>All w/o Altika</td><td>No</td><td>No</td><td>Yes</td><td>No</td></tr></tbody></table>
<figure><img src="image_18.png" /></figure>

<table><tr><td>EXP02: MIOST allsat-1</td><td>All w/o Altika</td><td>No</td><td>No</td><td>Yes</td><td>No</td></tr><tr><td>EXP03: MIOST allsat-1 80% drifters + equatorial waves+ L3 arctic</td><td>All w/o Altika</td><td>Yes (80%)</td><td>Yes</td><td>Yes</td><td>Yes</td></tr></table>

## 3.2 Validations metrics

The validation metrics are based on statistical and spectral analysis.

One quantitative assessment is based on the comparison between SSH maps and independent SSH along-track data. This diagnostic follows 3 main steps: 1) the SSH gridded data is interpolated to the locations of the independent SSH along-track, geo-referenced by their longitude, latitude, and time; 2) the difference SSH<sub>error</sub> = SSH<sub>map</sub> - SSH<sub>alongtrack</sub> is calculated and 3) a statistical analysis on the SSH<sub>error</sub> is performed in 1° × 1° longitude × latitude boxes. Prior to the statistical analysis, a filtering operation can be applied to isolate the spatial scales of interest. For example, the analysis can be performed over the spatial range [65km:500km] typically representative of the medium mesoscale ocean signal. This excludes the noisy part of the reference signal (along-track) as well as possible large-scale biases (scale > 500km). In the study, the validation metric is based on the error variance scores in 1°x1° longitude x latitude boxes (or averaged over specific region of interest), defined as:

$$ \sigma_{err}(x, y) = \frac{\sum_{t=1}^{N} (SSH_{error}(x, y, t) - \overline{SSH_{error}(x, y, t)})^2}{N} \quad (22) $$

The similar statistical analysis can also be performed on the geostrophic velocity errors $U_{error} = U_{map} - U_{drifter}$.

The comparison of the error variance score between two experiments informs about the gain or reduction $\Delta$ of the mapping error, for example:

$$ \Delta = 100 \cdot \frac{\sigma_{err}(EXP2) - \sigma_{err}(EXP1)}{\sigma_{err}(EXP1)} \quad (23) $$

The previous diagnosis is undertaken in physical space (space/time space). For a more descriptive assessment by wavelength and to avoid spatio-temporal filtering of independent and study datasets, diagnostics can be performed in frequency space, using spectral analysis of SSH altimetry and gridded datasets. More specifically, a spectral analysis can be applied to altimetry data to estimate the effective resolution of gridded SSH products. It is described for example in Ballarotta et al. (2019). Here, we recall the main processing steps for the estimation of the effective resolution: 1) the SSH<sub>map</sub> data are interpolated to the locations of independent SSH<sub>alongtrack</sub> data, 2) the along-track and interpolated data are divided into overlapping segments of 1500 km length every 300 km, 3) each segment is stored in a database and referenced by its median coordinates (longitude, latitude), 4) finally, between latitudes 90°N-90°S and longitudes 0°-360°E, we consider 10° × 10° longitude × latitude boxes for the global products every 1° increment. All available segments referenced in the 10° × 10° box
<figure><img src="image_19.png" /></figure>

are selected to compute the mean power spectral densities of the independent signal ($SSH_{alongtrack}$) and the mapping error ($SSH_{map}$, $SSH_{alongtrack}$). Before the spectral calculation, the signals are detrended and a Hanning window is applied. The signal-to-noise ratio (Equation 24) is then derived from the power spectral density of the PSD along the trace ($SSH_{alongtrack}$) and the power spectral density of the error ($SSH_{map}$, $SSH_{alongtrack}$). The effective resolution is then given by the wavelength $\lambda$s where the SNR($\lambda$s) is 2 (Equation 25), i.e., the wavelength where the $SSH_{error}$ is two times lower than the signal $SSH_{alongtrack}$.

$$SNR(\lambda) = \frac{PSD(SSH_{along-track})(\lambda)}{PSD(SSH_{error})(\lambda)} \qquad (24)$$

$$SNR(\lambda s) = 2 \qquad (25)$$
<figure><img src="image_20.png" /></figure>

# 4 Results

## 4.1 Qualitative assessment

The SLA maps from the DUACS and MIOST mapping approaches are relatively similar in the subpolar region, as illustrated in Figure 5 by an example of SLA reconstruction on 2019-02-15 for a) the DUACS mapping approach and c) the MIOST mapping approach. More significant differences take place in the Arctic basin: in contrast to the DUACS products, the use of arctic leads observations in MIOST offers the possibility to extend sea level mapping into ice-covered area and thus to deliver gap-free maps to the end-users (Figure 5c). Differences also appear in the relative vorticity field (Figure 5b and 5d): a more turbulent surface ocean circulation associated with finer structures is found in the maps constructed by the MIOST approach (Figure 5d) than in the DUACS maps (Figure 5b).

<figure><img src="image_21.png" /><figcaption>Figure 5: Example of sea level anomaly and relative vorticity maps over the Arctic region constructed with the DUACS mapping approach a) and b) and with the MIOST mapping approach c) and d). The black line contour indicates the 15% sea-ice concentration on 2019-02-15 from the OSI-SAF product</figcaption></figure>
<figure><img src="image_22.png" /></figure>

From a global perspective, the MIOST maps are indeed slightly more energetic than the DUACS maps (Figure 6). The difference between MIOST and DUACS variance maps (Figure 6b) indicates regions of higher variability in the MIOST maps (>10%) than in the DUACS maps, such as in the equatorial band, regions of low variability at mid-latitudes, coastal and polar regions. Tropical ocean regions are prone to lower SSH variability (10%) in the MIOST maps than in the DUACS maps.

Comparison mean variance of SLA maps

<figure><img src="image_23.png" /><figcaption>Figure 6: Variance (in m²) of sea level anomaly maps constructed with a) the DUACS approach and b) difference between the MIOST and DUACS variance maps expressed in %</figcaption></figure>

The large SSH variability in the equatorial band of the MIOST maps is mainly associated with the equatorial wave components. The zonal wavenumber-frequency spectrum of SSH in the Pacific has been investigated in several studies (e.g., Shinoda et al., 2009; Farrar, 2008, 2011) to examine the SSH variability associated with Tropical and Equatorial waves. Figure 7 shows contours of the base 10 logarithm of power in the wavenumber-frequency space calculated from SSH in the equatorial Pacific (region [180°E-280°E] [10°S:10°N]) for the period 2008 to 2018, for a) DUACS, b) MIOST with equatorial wave modes and c) in the GLORYS12V1 reanalysis (Lellouche et al., 2018). The rapid equatorial wave dynamics are resolved in the GLORYS12V1 ocean numerical simulation (Figure 7c): the zonal wavenumber-frequency spectrum of the SSH in the Pacific reveals significant spectral peaks at periods close to 4 days, 5 days, and 7 days for a wavelength > 20° in longitude. These peaks are associated with inertia-gravity (Poincaré) waves. These SSH variabilities for time scales smaller than 10 days are filtered in the DUACS mapping approach (Figure 7a). In contrast, the MIOST multiscale mapping approach
<figure><img src="image_24.png" /></figure>

(EXP03) resolves spectral peaks near 4 days, 5 days, and 7 days for wavelengths > 20° in longitude (Figure 7b). We show in the next section that parameterizing these equatorial wave modes in MIOST also contributes to significantly reduce the mapping error in this region. For time scale > 10 days, each dataset has relatively similar spectral contents, particularly the energetic westward propagation of equatorial Rossby waves for negative wavenumbers.

<figure><img src="image_25.png" /><figcaption>Figure 7: Zonal wavenumber-frequency spectrum of SLA in the Equatorial Pacific computed for a) DUACS, b) MIOST with equatorial wave modes and c) in the GLORYS12V1 reanalysis. White lines represent the theoretical dispersion relation curves for equatorial waves corresponding to the Kelvin, Yanai, Rossby and Poincaré waves.</figcaption></figure>

## 4.2 Quantitative assessment

### 4.2.1 Mesoscale mapping assessments

The first assessment is a comparison of the EXP01 and EXP02 experiments (MIOST allsat-1 and DUACS allsat-1 experiments). Both experiments aim to map the mesoscale circulation from altimetry data only. The SARAL/Altika altimeter and drifter sensors are not included in the mapping for independent evaluation.

Sea level anomaly quality
<figure><img src="image_26.png" /></figure>

The largest SSH mapping error $\sigma_{err}$ in DUACS allsat-1 reaches 50-100cm$^2$ in the western boundary surface current and over the continental plateaus (Figures 8a and 8b). In the offshore low variability region, the error variance is < 10cm$^2$. Figures 8c and 8d show the difference in mapping error between the EXP02 and EXP01 experiments for all spatial and spatial scales between 65 and 500km, respectively. Blue (red) pattern means a reduction (increase) of the mapping error in MIOST compared to DUACS. For all spatial scale considered, MIOST mapping errors are smaller than those of DUACS, especially at mid-latitude with an average reduction in mapping error between 5% and 10%. The largest reduction in mapping error (~10%) is found in regions of high variability. In the inter-tropical region, MIOST and DUACS have similar scores. For spatial scale between 65 and 500km, MIOST mapping errors are reduced by ~10% compared to DUACS in high variability region at mid-latitude. In low variability regions, the mapping error is between 3 and 4% smaller with MIOST than with DUACS, but the mapping errors are locally larger with MIOST than with DUACS: for example, in the Argentine Sea, in the Siberian plateau and New Zealand plateau. Table 4 summarises the results of the comparison over different regions of interest (arctic, Antarctic, equatorial band, low variability region, and high variability region). Overall, the geostrophic flows in the MIOST SSH maps are closer to the independent SARAL/AltiKa observations than those in DUACS maps.

Comparison SSH maps with independent SSH along-track

<figure><img src="image_27.png" /><figcaption>Figure 8: Variance of the difference SSH<sub>map</sub>-SSH<sub>alongtrack</sub> computed for the DUACS allsat-1 experiment and in considering a) all spatial scale and b) spatial scale between 65km and 500km. Gain/loss of the mapping error variance of SLA in MIOST allsat-1 experiment relatively to the DUACS allsat-1 mapping error variance for c) all spatial scale and d) scale between 65km and 500km. Blue colour means a reduction of error variance in MIOST.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Table 4: Regionally averaged mapping error variance and gain/reduction of error variance on the SSH variable between MIOST and DUACS</figcaption></figure>

<table><thead><tr><th rowspan="2">Region</th><th colspan="3">All spatial scale</th><th colspan="3">Spatial scale [65:500km]</th></tr><tr><th>Error variance DUACS [m²]</th><th>Error variance MIOST [m²]</th><th>Gain/loss error variance MIOST vs DUACS [%]</th><th>Error variance DUACS [m²]</th><th>Error variance MIOST [m²]</th><th>Gain/loss error variance MIOST vs DUACS [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>0,0023</td><td>0,0023</td><td>0,0,0%</td><td>0,0007</td><td>0,0007</td><td>-3,3%</td></tr><tr><td>Antarctic</td><td>0,0033</td><td>0,0031</td><td>-5,9%</td><td>0,0008</td><td>0,0008</td><td>-2,6%</td></tr><tr><td>Equatorial band</td><td>0,0014</td><td>0,0014</td><td>-0,8%</td><td>0,0005</td><td>0,0005</td><td>0,3%</td></tr><tr><td>Low variability - offshore</td><td>0,0013</td><td>0,0012</td><td>-5,8%</td><td>0,0004</td><td>0,0004</td><td>-4,1%</td></tr><tr><td>High variability - offshore</td><td>0,0031</td><td>0,0028</td><td>-10,2%</td><td>0,0014</td><td>0,0013</td><td>-9,9%</td></tr></tbody></table>

Geostrophic current quality

Figures 9a and 9b show the validation against the independent drifter velocity data in terms of mapping error $\sigma_{err}$ for the zonal and meridional velocities. The largest mapping error $\sigma_{err}$ in DUACS reaches 300 to 400cm².s⁻² in the western boundary surface current (e.g., the Gulfstream, the Kuroshio, Mozambique, and Agulhas currents). In offshore low variability region, the error variance is < 80cm².s⁻². The differences in mapping error between MIOST and DUACS are shown in Figures 9c and 9d for zonal and meridional velocities, respectively. Mapping errors are smaller in MIOST than in DUACS mainly in the core of the ocean gyres. In the intertropical region, the DUACS maps appear to be closer to the independent drifter velocities than MIOST. Table 5 summarises the results of the comparison over different regions of interest (arctic, Antarctic, equatorial band, low variability region, and high variability region). Overall, MIOST surface velocities are slightly closer to drifter velocities than the DUACS surface velocities.
<figure><img src="image_29.png" /></figure>

Comparison $U_{geos}$, $V_{geos}$ maps with independent $U_{geos}$, $V_{geos}$ from drifters

<figure><img src="image_30.png" /><figcaption>Figure 9: Variance of the difference $U_{map} - U_{drifter}$ computed for the DUACS allsat-1 experiment and in considering a) the zonal velocity component and b) the meridional velocity component. Gain/loss of the mapping error variance of currents in MIOST allsat-1 experiment relatively to the DUACS allsat-1 mapping error variance for c) the zonal velocity component and d) the meridional velocity component. Blue colour means a reduction of error variance in MIOST.</figcaption></figure>

<table><caption>Table 5: Regionally averaged mapping error variance and gain/reduction of error variance on the surface currents between MIOST and DUACS</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">Zonal velocity</th><th colspan="3">Meridional velocity</th></tr><tr><th>Error variance DUACS [m².s⁻²]</th><th>Error variance MIOST [m².s⁻²]</th><th>Gain/loss error variance MIOST vs DUACS</th><th>Error variance DUACS [m².s⁻²]</th><th>Error variance MIOST [m².s⁻²]</th><th>Gain/loss error variance MIOST vs DUACS</th></tr></thead><tbody><tr><td>Arctic</td><td>0,0153</td><td>0,0149</td><td>-2,9%</td><td>0,0134</td><td>0,0131</td><td>-1,6%</td></tr><tr><td>Antarctic</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Equatorial band</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Low variability - offshore</td><td>0,0130</td><td>0,0129</td><td>-1,4%</td><td>0,0124</td><td>0,0123</td><td>-0,9%</td></tr><tr><td>High variability - offshore</td><td>0,0386</td><td>0,0372</td><td>-3,5%</td><td>0,0410</td><td>0,0404</td><td>-1,5%</td></tr></tbody></table>

### 4.2.2 Contribution of equatorial waves modes, and drifters' observations

The comparison of experiment EXP03 with EXP02 (MIOST allsat-1 80% drifters + equatorial waves with experiment MIOST allsat-1) examines the impact of the equatorial waves' mode and the drifters' observations in the MIOST mapping approach.
<figure><img src="image_31.png" /></figure>

Sea level anomaly quality

The difference in mapping error between EXP03 and EXP02 are shown in Figures 10a and 10b for all spatial and spatial scales between 65 and 500km, respectively. For all spatial scale considered, we observe that the equatorial waves modes locally reduce the mapping error in the equatorial band by more than 10%. However, coastal equatorial regions (e.g., Indonesian Archipelago, western and Eastern part of Africa and South America) are prone to deterioration. This suggest that the equatorial wave mapping is not adapted in these coastal regions where different ocean processes are at play. In extra-equatorial regions, we evaluate the impact of drifter observations in MIOST. This impact is moderate on the SLA mapping (a few % of difference in the mapping error variance), with a reduction of error variance mainly in the high variability regions. For spatial scale between 65 and 500km (Figure 10b), the equatorial waves modes deteriorate the mapping solution in the western and central Equatorial Pacific Ocean, in the Indian Ocean, while a reduced mapping error is found in the eastern Equatorial Pacific and the Equatorial Atlantic. In the extra-equatorial region, the impact of drifter observations remains moderate (with 1.5% error variance reduction in the high variability region). Overall, the drifters reduce the mapping errors primarily in regions of intense dynamics where the temporal sampling must be sufficiently accurate to properly map the rapid mesoscale dynamics. Table 6 summarises the results of the comparison over different regions of interest (arctic, Antarctic, equatorial band, low variability region, and high variability region).

Comparison SSH maps with independent SSH along-track

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

Figure 10: Gain/loss of the mapping error variance of SLA in MIOST allsat-1 80% drifters + equatorial waves experiment relatively to the MIOST allsat-1 mapping error variance for a) all spatial scale and b) scale between 65km and 500km. Blue colour means a reduction of error variance in MIOST when drifters are included in the mapping and with equatorial waves parametrization
<figure><img src="image_33.png" /><figcaption>**Table 6: Regionally averaged mapping error variance and gain/reduction of error variance on the SSH variable between EXP03 and EXP02**</figcaption></figure>

<table><thead><tr><th rowspan="2">Region</th><th colspan="3">All spatial scale</th><th colspan="3">Spatial scale [65:500km]</th></tr><tr><th>Error variance EXP02 [m²]</th><th>Error variance EXP03 [m²]</th><th>Gain/loss error variance EXP03 vs EXP02 [%]</th><th>Error variance EXP03 [m²]</th><th>Error variance EXP02 [m²]</th><th>Gain/loss error variance EXP03 vs EXP02 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>0,0023</td><td>0,0023</td><td>0,0%</td><td>0,0007</td><td>0,0007</td><td>0,0%</td></tr><tr><td>Antarctic</td><td>0,0031</td><td>0,0031</td><td>0,0%</td><td>0,0008</td><td>0,0008</td><td>0,0%</td></tr><tr><td>Equatorial band</td><td>0,0014</td><td>0,0014</td><td>-3,0%</td><td>0,0005</td><td>0,0005</td><td>0,3%</td></tr><tr><td>Low variability - offshore</td><td>0,0012</td><td>0,0012</td><td>-0,8%</td><td>0,0004</td><td>0,0004</td><td>-0,1%</td></tr><tr><td>High variability - offshore</td><td>0,0028</td><td>0,0027</td><td>-1,1%</td><td>0,0013</td><td>0,0013</td><td>-1,5%</td></tr></tbody></table>

## Geostrophic current quality

The difference in mapping error of surface geostrophic currents between EXP03 and EXP02 are shown in Figures 11a and 11b for the zonal component and the meridional component of the velocity, respectively. It is difficult to draw conclusions from this diagnosis: the mapping errors are reduced with MIOST in some regions in the tropics (such as the Bay of Bengal), in the Kuroshio extension. Overall, the contribution of drifters remains moderate for the restitution of geostrophic currents (only a few % improvement in the open ocean) as summarized in Table 7.

<figure><img src="image_34.png" /><figcaption>**Figure 11: Gain/loss of the mapping error variance of currents in MIOST allsat-1 80% drifters + equatorial waves experiment relatively to the MIOST allsat-1 mapping error variance for c) the zonal velocity component and d) the meridional velocity component. Blue colour means a reduction of error in MIOST when drifters are included in the mapping and with equatorial waves parametrization**</figcaption></figure>

<table><caption>**Table 7: Regionally averaged mapping error variance and gain/reduction of error variance on the surface currents between EXP02 and EXP03**</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">Zonal velocity</th><th colspan="3">Meridional velocity</th></tr><tr><th>Error variance EXP02 [m².s⁻²]</th><th>Error variance EXP03 [m².s⁻²]</th><th>Gain/loss error variance EXP03 vs EXP02 [%]</th><th>Error variance EXP02 [m².s⁻²]</th><th>Error variance EXP03 [m².s⁻²]</th><th>Gain/loss error variance EXP03 vs EXP02 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>0,0149</td><td>0,0145</td><td>-2,5%</td><td>0,0131</td><td>0,0128</td><td>-2,7%</td></tr><tr><td>Antarctic</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<figure><img src="image_35.png" /></figure>

<table><thead><tr><th>Equatorial band</th><th>-</th><th>-</th><th>-</th><th>-</th><th>-</th><th>-</th></tr></thead><tbody><tr><td>Low variability - offshore</td><td>0,0129</td><td>0,0128</td><td>-0,6%</td><td>0,0123</td><td>0,0122</td><td>-0,9%</td></tr><tr><td>High variability - offshore</td><td>0,0372</td><td>0,0367</td><td>-1,5%</td><td>0,0404</td><td>0,0401</td><td>-0,7%</td></tr></tbody></table>

### 4.2.3 Overall assessment

The comparison of the EXP03 and EXP01 experiments (MIOST allsat-1 80% drifters + equatorial waves with the DUACS allsat-1 experiment) allows to evaluate the complete MIOST product distributed to users against the DUACS method.

#### *Sea level anomaly quality*

The difference in mapping error between EXP03 and EXP01 are shown in Figures 12a and 12b for all spatial and spatial scales between 65 and 500km, respectively. We have the same pattern as found in the previous sections: for all spatial scale considered (Figure 12a), the equatorial waves modes help to reduce the mapping error variance in the equatorial band by more than 20% locally. At mid-latitude, the mapping error are between 5% and 10% smaller with MIOST than with DUACS. For spatial scales between 65 and 500km, MIOST and DUACS solutions are globally equivalent, except in the high variability region where the mapping error is between 10% and 20% smaller with MIOST than with DUACS. The mapping errors are locally larger with MIOST than with DUACS in regions where the circulation interact with bathymetry feature such as in the Argentine Sea, near the Siberian plateau and New Zealand plateau. Table 8 summarises the results of the comparison over different regions of interest: mapping errors are ~11% smaller in high variability region in MIOST than in DUACS. In other regions, the errors are ~3-6% smaller.

#### Comparison SSH maps with independent SSH along-track

<figure><img src="image_36.png" /><figcaption>Figure 12: Gain/loss of the mapping error variance of SLA in MIOST allsat-1 80% drifters + equatorial waves experiment relatively to the DUACS allsat-1 mapping error variance for a) all spatial scale and b) scale between 65km and 500km. Blue colour means a reduction of error variance in MIOST.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Table 8: Regionally averaged mapping error variance and gain/reduction of error variance on the SSH variable between EXP03 and EXP01</figcaption></figure>

<table><thead><tr><th rowspan="2">Region</th><th colspan="3">All spatial scale</th><th colspan="3">Spatial scale [65:500km]</th></tr><tr><th>Error variance EXP01 [m²]</th><th>Error variance EXP03 [m²]</th><th>Gain/loss error variance EXP03 vs EXP01 [%]</th><th>Error variance EXP01 [m²]</th><th>Error variance EXP03 [m²]</th><th>Gain/loss error variance EXP03 vs EXP01 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>0,0023</td><td>0,0023</td><td>0,0%</td><td>0,0007</td><td>0,0007</td><td>-3,3%</td></tr><tr><td>Antarctic</td><td>0,0033</td><td>0,0031</td><td>-5,8%</td><td>0,0008</td><td>0,0008</td><td>-2,6%</td></tr><tr><td>Equatorial band</td><td>0,0014</td><td>0,0014</td><td>-3,8%</td><td>0,0005</td><td>0,0005</td><td>0,6%</td></tr><tr><td>Low variability - offshore</td><td>0,0013</td><td>0,0012</td><td>-6,6%</td><td>0,0004</td><td>0,0004</td><td>-4,2%</td></tr><tr><td>High variability - offshore</td><td>0,0031</td><td>0,0027</td><td>-11,2%</td><td>0,0014</td><td>0,0013</td><td>-11,2%</td></tr></tbody></table>

Geostrophic current quality

The difference in mapping error of surface geostrophic currents between EXP03 and EXP01 are shown in Figures 13a and 13b for the zonal component and the meridional component of the velocity, respectively. The mapping errors are globally smaller in MIOST than in DUACS, particularly in the high variability regions. In the tropical regions, DUACS outperforms MIOST for reconstructing the surface geostrophic velocities. Overall, the mapping errors are on average between ~2% and 5% smaller with MIOST than with DUACS (Table 9).

Comparison U<sub>geos</sub>, V<sub>geos</sub> maps with independent U<sub>geos</sub>, V<sub>geos</sub> from drifters

<figure><img src="image_38.png" /><figcaption>Figure 13: Gain/loss of the mapping error variance of currents in MIOST allsat-1 80% drifters + equatorial waves experiment relatively to the DUACS allsat-1 mapping error variance for c) the zonal velocity component and d) the meridional velocity component. Blue colour means a reduction of error in MIOST</figcaption></figure>

<table><caption>Table 9: Regionally averaged mapping error variance and gain/reduction of error variance on the surface currents between EXP03 and EXP01</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">Zonal velocity</th><th colspan="3">Meridional velocity</th></tr><tr><th>Error variance EXP01 [m².s⁻²]</th><th>Error variance EXP03 [m².s⁻²]</th><th>Gain/loss error variance EXP03 vs EXP01 [%]</th><th>Error variance EXP01 [m².s⁻²]</th><th>Error variance EXP03 [m².s⁻²]</th><th>Gain/loss error variance EXP03 vs EXP01 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>0,0153</td><td>0,0145</td><td>-5,3%</td><td>0,0134</td><td>0,0128</td><td>-4,2%</td></tr><tr><td>Antarctic</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
<figure><img src="image_39.png" /></figure>

<table><thead><tr><th>Equatorial band</th><th>-</th><th>-</th><th>-</th><th>-</th><th>-</th><th>-</th></tr></thead><tbody><tr><td>Low variability - offshore</td><td>0,0130</td><td>0,0128</td><td>-2,0%</td><td>0,0124</td><td>0,0122</td><td>-1,9%</td></tr><tr><td>High variability - offshore</td><td>0,0386</td><td>0,0367</td><td>-4,9%</td><td>0,0410</td><td>0,0401</td><td>-2,2%</td></tr></tbody></table>

Effective resolution

The effective spatial resolution quantifies the minimum spatial scale resolved in the maps. Maps of the effective spatial resolution (expressed in kilometres) are presented in Figure 14a and Figure 14b for EXP01 and EXP03, respectively. For each experiment, the effective spatial resolution varies from ~500km at the equator to ~100km at high altitude, and a mean value at mid-latitude close to 200km. The difference in effective spatial resolution between the two experiments is shown in Figure 14c. The resolution of the SLA maps of the MIOST experiment is overall finer than in the SLA maps of the DUACS experiment. It is between 5% and 10% finer than the DUACS maps in regions of high variability (Gulfstream, Kuroshio, and Agulhas regions), in the Atlantic and equatorial Pacific, and in the Norwegian and Greenland seas. Some regions (e.g., tropical regions, coastal regions, the East China Sea, the New Zealand Shelf, or the Argentine Sea) are subject to a coarser effective resolution in MIOST maps than in DUACS maps. These regions will require further investigation in the near future.

515

520

525
<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /><figcaption>Figure 14: Maps of effective spatial resolution (in km) for a) the DUACS allsat-1 and b) MIOST allsat-1 80% drifters + equatorial waves experiments; and c) gain/loss of effective resolution (in %) between MIOST and DUACS. Blue means finer resolution in MIOST than in DUACS</figcaption></figure>

# 5 Summary & Conclusions

Ubelmann et al (2020, 2021) evaluated the multiscale & multivariate mapping approach in Observing System Simulation Experiment (OSSE) and Observing System Experiment (OSE) for simultaneous mapping of mesoscale circulation, coherent internal tides, surface geostrophic and ageostrophic velocities. Here, we extend the application of the MIOST solution to the simultaneous mapping of equatorial waves and mesoscale circulation from real observations. Furthermore, we investigate the levels of mapping improvement by enhancing the sampling of the ocean surface state with in-situ data and altimetry data in
<figure><img src="image_42.png" /></figure>

the Arctic sea-ice regions. We found that the Arctic leads SSH observations allow to significantly improve the monitoring coverage in this remote region. The gap-free maps, proposed with MIOST, hence offer to the end-users the opportunity to study the arctic surface circulation and its connections to the subpolar and mid-latitude regions. It is important to mention that this polar mapping will need to be validated against independent data in the near future. Drifters' observations have a moderate impact in the mapping. They mainly contribute to reduce mapping errors in regions of intense dynamics where the temporal sampling must be accurate enough to properly map the rapid mesoscale dynamics. It is important to note that drifter observations can potentially improve surface circulation in areas not or poorly sampled by altimeters. Therefore, their impact on the sea level reconstruction may be larger over period of weak altimeter sampling.

The ocean surface circulation involves a superposition of processes acting at widely different spatial and temporal scales, from the geostrophic large-scale and slow varying flow to the mesoscale turbulent eddies and at even smaller scale, the mixing generated by the internal wave field. It is also important to mention that the DUACS maps are constructed from altimetry data using an interpolation method optimized for mapping mesoscale variability. Consequently, some ocean surface variabilities are not or poorly represented in these DUACS maps: equatorial wave dynamics is thus part of the filtered ocean signals in DUACS. The multiscale approach allows to decompose the observed SSH into various physical contributions. Here, we explored and validated the possibility of improving the content of altimetry maps by simultaneously estimating the ocean mesoscale circulations as well as the equatorial wave dynamics associated to the Tropical Instability and Poincaré Waves. We show that mapping these ocean surface variabilities from altimeter observations broadens the spectrum of mappable space-time scales and reduces mapping errors by almost 20% locally relative to independent data, primarily in the equatorial Pacific and Atlantic basins. This is possible because the spatio-temporal coverage of the altimeter data allows to sample large scale waves of 4-day periods and longer. At global scale, we also found that, compared to the operational DUACS mapping approach, MIOST approach improves the surface mesoscale circulations in regions of high variability. Consequently, the effective resolution of the maps produced by the multiscale approach is finer than the DUACS maps, particularly in the western boundary currents and in the equatorial band.

To conclude, these results pave the way for the exploration of new types of ocean signals that may eventually be mapped from remote sensing and in situ observations. Future work could consist of enriching the MIOST modes in considering oceanic signals missing in the maps and yet captured by observing systems: for example, in mapping high frequency signals such as the near-inertial oscillation from drifter observations, in using SSH leads products in the Southern Ocean (Auger et al., in prep); or by enhancing the SLA maps content with dynamical model approach (Ubelmann et al, 2016) or Artificial Intelligence methods (Beauchamp et al., 2020).
<figure><img src="image_43.png" /></figure>

## 6 Data Availability

The MIOST gridded products (https://doi.org/10.24400/527896/a01-2022.009, Ballarotta et al., 2022) are hosted on the AVISO+ website at the following repository: https://data.aviso.altimetry.fr/aviso-gateway/data/SLA_MIOST_alti_drifters/. The reference DUACS maps are hosted on the E.U. Copernicus Marine Service portal (https://doi.org/10.48670/moi-00146).

570

## 7 Product description

The multiscale & multivariate products are distributed on a regular grid: the spatial grid extends from 0°E to 360°E in longitude, 80°S to 90°N in latitude, with a grid spacing of 0.1°; the temporal grid covers the period 2016-07-01 to 2020-06-30 with a time step of 1 day. The dataset is distributed in netCDF4 format. Each netCDF file contains 6 variables: sla, adt, ugosa, vgosa, ugos, and vgos.

575

## 8 Acknowledgements

The work presented here was carried out in the framework of the DUACS-R&D project funded by CNES. The authors would like to thank the AVISO+ team for their support and expertise in the distribution of the dataset.
<figure><img src="image_44.png" /></figure>

```cpp
dimensions:
    time = 1;
    latitude = 1702;
    longitude = 3600;
    bounds = 2;

variables:
    int sla(time, latitude, longitude);
        sla::FillValue = -2147483647;
        sla::coordinates = "longitude latitude";
        sla::grid_mapping = "crs";
        sla::long_name = "sea level anomaly";
        sla::standard_name = "sea surface height above sea level";
        sla::units = "m";
        sla::scale_factor = 0.0001;
    int ugosa(time, latitude, longitude);
        ugosa::FillValue = -2147483647;
        ugosa::coordinates = "longitude latitude";
        ugosa::grid_mapping = "crs";
        ugosa::long_name = "Geostrophic velocity anomalies: zonal component";
        ugosa::standard_name = "surface geostrophic eastward sea water velocity assuming sea level for geoid";
        ugosa::units = "m";
        ugosa::scale_factor = 0.0001;
    int vgosa(time, latitude, longitude);
        vgosa::FillValue = -2147483647;
        vgosa::coordinates = "longitude latitude";
        vgosa::grid_mapping = "crs";
        vgosa::long_name = "Geostrophic velocity anomalies: meridional component";
        vgosa::standard_name = "surface geostrophic northward sea water velocity assuming sea level for geoid";
        vgosa::units = "m";
        vgosa::scale_factor = 0.0001;
    int adt(time, latitude, longitude);
        adt::FillValue = -2147483647;
        adt::coordinates = "longitude latitude";
        adt::grid_mapping = "crs";
        adt::long_name = "Absolute dynamic topography";
        adt::standard_name = "sea surface height above sea level";
        adt::units = "m";
        adt::scale_factor = 0.0001;
    int ugos(time, latitude, longitude);
        ugos::FillValue = -2147483647;
        ugos::coordinates = "longitude latitude";
        ugos::grid_mapping = "crs";
        ugos::long_name = "Absolute geostrophic velocity: zonal component";
        ugos::standard_name = "surface geostrophic eastward sea water velocity";
        ugos::units = "m";
        ugos::scale_factor = 0.0001;
    int vgos(time, latitude, longitude);
        vgos::FillValue = -2147483647;
        vgos::coordinates = "longitude latitude";
        vgos::grid_mapping = "crs";
        vgos::long_name = "Absolute geostrophic velocity: meridional component";
        vgos::standard_name = "surface geostrophic northward sea water velocity";
        vgos::units = "m";
        vgos::scale_factor = 0.0001;
    float latitude(latency);
        latitude::FillValue = NaNf;
        latitude::axis = "Y";
        latitude::long_name = "Latitude";
        latitude::standard_name = "latitude";
        latitude::units = "degrees north";
        latitude::valid_max = 90.;
        latitude::valid_min = -80.1;
        latitude::bounds = "latitude bounds";
    float longitude(longitude);
        longitude::FillValue = NaNf;
        longitude::axis = "X";
        longitude::long_name = "Longitude";
        longitude::standard_name = "longitude";
        longitude::units = "degrees east";
        longitude::valid_max = 359.9;
        longitude::valid_min = 0.;
        longitude::bounds = "longitude bounds";
    double time(time);
        time::FillValue = NaNf;
        time::units = "days since 1950-01-01 00:00:00";
        time::calendar = "gregorian";
        time::axis = "T";
        time::standard_name = "time";
    float longitude bounds(longitude, bounds);
        longitude bounds::FillValue = NaNf;
    float latitude bounds(latency, bounds);
        latitude bounds::FillValue = NaNf;
```
<figure><img src="image_45.png" /></figure>

References

580

Auger M., P. Prandi, J-B. Sallée, Southern Ocean Sea Level Anomaly in the Sea-Ice Covered Sector From Multimission Satellite Observations, Nature Scientific Data, in preparation.

Ballarotta, M., Ubelmann, C., Pujol, M.-I., Taburet, G., Fournier, F., Legeais, J.-F., Faugère, Y., Delepoulle, A., Chelton, D., Dibarboure, G., and Picot, N.: On the resolutions of ocean altimetry maps, Ocean Sci., 15, 1091–1109, https://doi.org/10.5194/os-15-1091-2019, 2019.

585

Ballarotta, M., and Coauthors, 2020: Dynamic mapping of along-track ocean altimetry: Performance from real observations. J. Atmos. Oceanic Technol., 37, 1593–1601, https://doi.org/10.1175/JTECH-D-20-0030.1.

Ballarotta, M. and Coauthors, 2022: Gridded Sea Level Height and geostrophic velocities computed with Multiscale Interpolation combining altimetry and drifters, [dataset], https://doi.org/10.24400/527896/a01-2022.009

590

Beauchamp, M.; Fablet, R.; Ubelmann, C.; Ballarotta, M.; Chapron, B. Intercomparison of Data-Driven and Learning-Based Interpolations of Along-Track Nadir and Wide-Swath SWOT Altimetry Observations. Remote Sens. 2020, 12, 3806. https://doi.org/10.3390/rs12223806

Bojinski S, Verstraete M, Peterson T, Richter C, Simmons A, Zemp M. The Concept of Essential Climate Variables in Support of Climate Research, Applications, and Policy. BULLETIN OF THE AMERICAN METEOROLOGICAL SOCIETY 95 (9); 2014. p. 1431-1443. JRC87032

595

Dufau, C., Orsztynowicz, M., Dibarboure, G., Morrow, R., and Le Traon, P.-Y.: Mesoscale resolution capability of altimetry: present and future, J. Geophys. Res.-Oceans, 121, 4910–4927, https://doi.org/10.1002/2015JC010904, 2016.

Etienne, H., Verbrugge, N., Boone, C., Rubio, A., Solabarrieta, L., Corgnati, L., Mantovani, C., Reyes, E., Chifflet, M., Mader, J., Carval, T. Quality Information document for global ocean-delayed mode in-situ observations of surface (drifters and HFR) and sub-surface (vessel mounted ADCPs) water velocity, https://catalogue.marine.copernicus.eu/documents/QUID/CMEMS-INS-QUID-013-044.pdf, 2021

600

Etienne, H., Quality Information document for global ocean Multi Observation products, https://catalogue.marine.copernicus.eu/documents/QUID/CMEMS-MOB-QUID-015-004.pdf, 2021

Farrar, J. T., 2008: Observations of the dispersion characteristics and meridional sea level structure of equatorial waves in the Pacific Ocean. J. Phys. Oceanogr., 38, 1669–1689.

605

Farrar, J. T., 2011: Barotropic Rossby waves radiating from tropical instability waves in the Pacific Ocean. J. Phys. Oceanogr., 41, 1160–1181.

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thépaut, J-N.: ERA5 hourly data on single levels from 1979 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). (Accessed on 07-10-2021), 10.24381/cds.adbb2d47, 2018

610
<figure><img src="image_46.png" /></figure>

Le Guillou, F., Metref, S., Cosme, E., Ubelmann, C., Ballarotta, M. Le Sommer, J. Verron, J. Mapping Altimetry in the Forthcoming SWOT Era by Back-and-Forth Nudging a One-Layer Quasigeostrophic Model, J. Atmos. Oceanic Technol., 38, 697–710, https://doi.org/10.1175/JTECH-D-20-0104.1

Lumpkin, R., Centurioni, L.: Global Drifter Program quality-controlled 6-hour interpolated data from ocean surface drifting buoys. NOAA National Centers for Environmental Information. Dataset. https://doi.org/10.25921/7ntx-z961, 2019

Morrow, R., Fu, L.L., Ardhuin, F., Benkiran, M., Chapron, B., Cosme, E., d'Ovidio, F., Farrar, J.T., Gille, S.T., Lapeyre, G., Le Traon, P.Y., Pascual, A., Ponte, A., Qiu, B., Rascle, N., Ubelmann, C., Wang, J., Zaron, E.D., 2019. Global observations of fine-scale ocean surface topography with the Surface Water and Ocean Topography (SWOT) mission. Front. Mar. Sci. 6, 232. https://doi.org/10.3389/fmars.2019.00232.

Lumpkin, R., S. Grodsky, M.-H. Rio, L. Centurioni, J. Carton and D. Lee.: Removing spurious low-frequency variability in surface drifter velocities. J. Atmos. Oceanic Techn., 30 (2), 353–360, http://dx.doi.org/10.1175/JTECH-D-12-00139.1, 2013.

Lellouche, J.-M., Greiner, E., Le Galloudec, O., Garric, G., Regnier, C., Drevillon, M., Benkiran, M., Testut, C.-E., Bourdalle-Badie, R., Gasparin, F., Hernandez, O., Levier, B., Drillet, Y., Remy, E., and Le Traon, P.-Y.: Recent updates to the Copernicus Marine Service global ocean monitoring and forecasting real-time 1/12° high-resolution system, Ocean Sci., 14, 1093–1126, https://doi.org/10.5194/os-14-1093-2018, 2018.

Mulet, S., Rio, M.-H., Etienne, H., Artana, C., Cancet, M., Dibarboure, G., Feng, H., Husson, R., Picot, N., Provost, C., and Strub, P. T.: The new CNES-CLS18 global mean dynamic topography, Ocean Sci., 17, 789–808, https://doi.org/10.5194/os-17-789-2021, 2021.

Rio, M. H., Guinehut, S., and Larnicol, G.: New CNES-CLS09 global mean dynamic topography computed from the combination of GRACE data, altimetry, and in-situ measurements, J. Geophys. Res., 116, C07018, doi:10.1029/2010JC006505, 2011.

Pujol, M.-I., Faugère, Y., Taburet, G., Dupuy, S., Pelloquin, C., Ablain, M., and Picot, N.: DUACS DT2014: the new multi-mission altimeter data set reprocessed over 20 years, Ocean Sci., 12, 1067–1090, https://doi.org/10.5194/os-12-1067-2016, 2016.

Pujol, M-I, Dupuy, S. Vergara, O., Sánchez-Román, A. Faugere, Y., Prandi, P., Dabat, M-L., Dagneaux, Q. Lievin, M. Cadier, E., Dibarboure, G. Picot, N. (in prep.): Refining the resolution of DUACS along track (level 3) altimeter Sea Level products

Rio, M.-H.: Use of altimeter and wind data to detect the anomalous loss of SVP-type drifter's drogue, J. Atmos. Ocean. Technol., 1663–1674, doi:10.1175/JTECH-D-12-00008.1, 2012.

Shinoda, T., G. N. Kiladis, and P. E. Roundy, 2009: Statistical representation of equatorial waves and tropical instability waves in the Pacific Ocean. Atmos. Res., 94, 37–44.

Taburet, G., Pujol, M.-I. & DUACS Team, Quality information document, Sea Level TAC – DUACS Products, Copernicus Marine Service, https://catalogue.marine.copernicus.eu/documents/QUID/CMEMS-SL-QUID-008-032-068.pdf, 2021
<figure><img src="image_47.png" /></figure>

645

Ubelmann, C., P. Klein, and L.-L. Fu, 2015: Dynamic interpolation of sea surface height and potential applications for future high-resolution altimetry mapping. J. Atmos. Ocean. Technol., 32, 177–184, https://doi.org/10.1175/JTECH-D-14-00152.1.

Ubelmann, C., B. Cornuelle, and L. Fu, 2016: Dynamic mapping of along-track ocean altimetry: Method and performance from observing system simulation experiments. J. Atmos. Oceanic Technol., 33, 1691–1699, https://doi.org/10.1175/JTECH-D-15-0163.1.

650

Ubelmann, C., Dibarboure, G., Gaultier, L., Ponte, A., Ardhuin, F., Ballarotta, M., & Faugère, Y. (2021). Reconstructing ocean surface current combining altimetry and future spaceborne Doppler data. Journal of Geophysical Research: Oceans, 126, e2020JC016560. https://doi.org/10.1029/2020JC016560

Ubelmann, C., Carrere, L., Durand, C., Dibarboure, G., Faugère, Y., Ballarotta, M., Briol, F., and Lyard, F.: Simultaneous estimation of Ocean mesoscale and coherent internal tide Sea Surface Height signatures from the global Altimetry record, Ocean Sci. Discuss. [preprint], https://doi.org/10.5194/os-2021-80, in review, 2021a.

655