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

Maxime Ballarotta, Clement Ubelmann, Pierre Veillard, Pierre Prandi, Helene Etienne, Sandrine Mulet, Yannice Faugere, Gerald Dibarboure, Rosemary Morrow, Nicolas Picot

► To cite this version:

Maxime Ballarotta, Clement Ubelmann, Pierre Veillard, Pierre Prandi, Helene Etienne, et al. Improved global sea surface height and current maps from remote sensing and in situ observations. Earth System Science Data, 2023, 10.5194/essd-15-295-2023. hal-04720455

HAL Id: hal-04720455

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

Submitted on 4 Oct 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 current 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⁴, and 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

⁴Centre de Topographie des Océans et de l'Hydrosphère, Laboratoire d'Etudes en Géophysique et Océanographie Spatiale, CNRS, CNES, IRD, Université Toulouse III, Toulouse, France

**Correspondence:** Maxime Ballarotta (mballarotta@groupcls.com)

Received: 25 May 2022 – Discussion started: 5 July 2022

Revised: 15 December 2022 – Accepted: 17 December 2022 – Published: 17 January 2023

**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 and multivariate mapping approach that offers the possibility to improve the physical content of gridded products by combining the data from various platforms and 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 1 July 2016 to 30 June 2020. 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 with regard to an operational product distributed in the Copernicus Marine Service. We show that the multiscale and multivariate mapping approach offers promising perspectives for reconstructing the ocean surface circulation: observations of leads contribute to improvement of the coverage in delivering gap-free maps in the Arctic and observations of drifters 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 resolved scales of this new product are therefore between 5% and 10% finer than the Copernicus product (https://doi.org/10.48670/moi-00148, Pujol et al., 2022b).

# 1 Introduction

Several oceanographic applications (e.g. operational oceanography, marine weather, and climate monitoring) rely on high-quality observational datasets. The European Union (EU) Copernicus Marine and Climate Change Services provide 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.

They 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 Sea Level The- matic Assembly Centre (SL-TAC) delivers near-real-time and delayed-time sea-level and surface current products (along-track Level-3 and gridded Level-4 products) that are used by the ocean science community to study, un-
derstand, 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 × 20 d for the 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 Surface Water and Ocean Topography (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 and higher frequencies have been identified by the ocean scientific community and the Copernicus Services as research and development priorities to serve Copernicus marine users and decision-makers (see, for example, International Altimetry Team, 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/CMEMSService_evolution_strategy_RD_priorities_v5-June-2021.pdf, last access: 1 December 2022). 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., 2021a; Ballarotta et al., 2020; Ubelmann et al., 2022; Prandi et al., 2021) and preparing operational systems for the SWOT era (Ubelmann et al., 2015, 2021; 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 and multivariate mapping approach (Ubelmann et al., 2021, 2022) and applied to real Earth observations. Here we 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 Sect. 2. Section 3 presents the experiments and validation metrics. The quality assessment of the new products is proposed in Sect. 4. The key results are then summarized in Sect. 5.

# 2 Data and methods

## 2.1 Data sources

The mapping method used in this study takes input data from remote sensing and in situ observations, which are summarized in Table 1 and described below.

### 2.1.1 Sea-level anomaly products

The global ocean sea surface height (SSH) observations are from the (delayed time, DT) Level-3 altimeter satellite along-track data, reprocessed in 2021 and distributed by the EU Copernicus Marine Service (product reference SEALEVEL_GLO PHY_L3_MY_008_062, https://doi.org/10.48670/moi-00146, Pujol et al., 2022a). These data cover the period from 1 January 1993 to 31 December 2020 over the world ocean (excluding ice-covered areas; see, for example, Fig. 1) and are available at a sampling rate of 1 Hz (~7 km spatial spacing). Homogenization 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 datasets is described in the quality information document (Taburet et al., 2021) and summarized below in Eq. (1). In this study, unfiltered sea-level anomalies (SLAs) corrected with dynamic atmospheric correction (dac), ocean tide, and lwe corrections are considered in the multiscale and multivariate mapping.

$$SLA = \text{orbit} - \text{range} - \sum (\text{environmental corrections}) - \sum (\text{geophysical corrections}) - \text{mean sea surface}, \quad (1)$$

where $\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}$ (see Taburet et al., 2021, for the references associated with each mission correction). The mean sea surface used here is the CNES-CLS18 (Mulet et al., 2021b).

### 2.1.2 Sea-level anomaly products in Arctic leads

In the polar regions, satellite sea-level observations are limited by the sea ice. Thanks to dedicated processing, sea level can however be estimated within fractures in the ice (leads). The echoes from the altimeters over the ice-covered region are classified to identify peaky waveforms corresponding to lead echoes. Range estimation is then made with specific retracking methods, and it is corrected from instrumental and geophysical corrections to obtain a sea-level anomaly (Prandi et al., 2021). To ensure continuity with the open ocean, the corrections are derived from the global ocean Level-3 along-track processing (Taburet et al., 2021) when possible. The noticeable exceptions concern (1) the wet tropospheric correction that comes from the European Centre for Medium-Range Weather Forecasts (ECMWF) model since onboard radiometer estimates are not reliable over ice, (2) the sea-state bias correction which is not applied since waves and winds are considered small over leads, and (3) orbit error corrections which are not applied as they are difficult to compute over this small region. Then a constant bias of
<table><caption>**Table 1.** List of observation datasets used in this study.</caption><thead><tr><th>Product type</th><th>Global altimeter SLA products</th><th>Arctic leads' altimeter SLA products</th><th>Drifters' geostrophic velocity product</th></tr></thead><tbody><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>[90° S–90° N], [0° E–360° E]</td><td>>60° N</td><td>[90° S–90° N], [0° E–360° E]</td></tr><tr><td>Period</td><td>From 15 January 2016 to 30 June 2020</td><td>From 15 January 2016 to 30 June 2020</td><td>From 15 January 2016 to 30 June 2020</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of sea-level altimetry coverage for a 7 d period (from 1 July 2019 to 7 July 2019). Colour scale represents the sea-level anomaly amplitude in metres. For this time interval, data originate from six altimeters: Jason-3, Sentinel-3A, Sentinel-3B, SARAL/AltiKa, Cryosat-2, and Haiyang-2A.</figcaption></figure>

~ 8 cm is applied for each mission to ensure continuity with the SEALEVEL_GLO PHY L3_MY_008_062 open-ocean SLA previously described. These products cover the Arctic region (up to 88° N) at a sampling rate of 20 Hz (~ 350 m) for three altimetry missions: SARAL/AltiKa, Sentinel-3A, and CryoSat-2 (Fig. 2 and Table 2).

### 2.1.3 Geostrophic current anomaly products

To further constrain the surface circulation, we used delayed-time horizontal surface velocities from the NOAA's Atlantic Oceanographic and Meteorological Laboratory (AOML) Surface Velocity Program (SVP; Lumpkin and Centurioni, 2019). The data cover the entire world ocean and are available at a 6 h frequency. The SVP drifters are designed to follow the 15 m 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 IN SITU The- matic Assembly Centre of the EU Copernicus Marine Ser- vice (see Product User Manual; http://marine.copernicus.eu/ documents/PUM/CMEMS-INS-PUM-013-044.pdf, last ac- cess: 9 January 2023) with an additional wind slippage cor- rection for undrogued buoys derived from the Rio (2012) methodology. For the study, the undrogued and drogued drifters are selected over the global ocean and the period from 1 June 2016 to 31 July 2020. Note that for specific

experiments described hereafter, we excluded drifters' tra- jectories between –10° S and 10° N (e.g. Fig. 3) to isolate and evaluate only the impact of the equatorial wave's mode in this region. As in Mulet et al. (2021a), 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} \qquad (2) $$

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

where $U_{buoy}$ ($V_{buoy}$) is the drifter's zonal (meridional) veloc- ity. Each component is corrected as follows.

- The wind-driven component $U_{ekman}$ ($V_{ekman}$) uses an update of the model used in Mulet et al. (2021a) and described in Etienne (2021a). The Ekman component is not available in the Mediterranean basin, so there is no drifter used in this region for the study. In this re- cent version, ERA5 wind stress (Hersbach et al., 2018) replaces the ERA Interim data, and the equatorial sym- metry 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 15 m depth velocity, as this compo- nent is supposed to mostly vanish in the first 2–4 m.
<table><caption>**Table 2.** Arctic leads' product characteristics. LRM: low-resolution mode. TFMRA: threshold first-maximum retracker algorithm.</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° N</td><td>81.5° N</td><td>88° N</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_4.png" /><figcaption>**Figure 2.** Example of Arctic leads' sea-level altimetry coverage for a 7 d period (from 1 July 2019 to 7 July 2019). Colour scale represents the sea-level anomaly amplitude in metres.</figcaption></figure>

- The wind slippage 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 are 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., 2021b) $U_{mdt}$ ($V_{mdt}$) is subtracted to obtain the geostrophic velocity anomaly.

## 2.2 Methods

Two mapping methods are compared in this study: the operational DUACS (Data Unification and Altimeter Combination System) mapping approach and the multiscale and multivariate MIOST (Multiscale Inversion of Ocean Surface Topography) mapping approach. Each method is described in detail in reference articles, such as Le Traon et al. (1998, 2003), Ducet et al. (2000), or Pujol et al. (2016) for the DUACS method and Ubelmann et al. (2021, 2022) for the

MIOST method. A description of the methods is given in Appendix A, and we propose hereafter to focus on the specific developments and processes that are considered in this study.

It is important to mention that DUACS maps are constrained by a single-scale covariance function (Arhan and Colin de Verdière, 1985; Le Traon et al., 1998) and focus mainly on the geostrophic circulation (i.e. processes with typical space and timescales >100 km, 10 d). Consequently, they do not resolve the full spectrum of ocean surface variability. It is, for example, the case for the equatorial surface dynamics (see, for example, Fig. 7). While slow Rossby waves are already resolved within geostrophy in DUACS maps, faster equatorial waves such as Poincaré waves are filtered out, even though the space-time coverage of altimetry data allows for sampling of large-scale waves with periods of 4–10 d and more (Farrar and Durland, 2022). The multiscale approach proposed by the MIOST method offers the possibility to solve some of the missing surface variabilities in DUACS, accounting for the covariances of various surface processes in a single inversion. The covariance functions in
<figure><img src="image_5.png" /><figcaption>Figure 3. Example of the drifter’s trajectory coverage for the period 1 January 2019 to 31 December 2019. Colour scale represents the velocity amplitude in metres per second (m s⁻¹).</figcaption></figure>

the MIOST system are expressed as wavelet modes, and the inversion is performed in this space using a variational approach (Ubelmann et al., 2021). In the following, we focus on the main components that have been tested in this study with the MIOST method: the geostrophy component already investigated in Ubelmann et al. (2021) and two new components associated with the equatorial wave’s dynamic.

The geostrophy component follows the same formulation provided in Ubelmann et al. (2021) (see their Sect. 2.3.2.1, where the analytical formula of the ensemble of wavelet elements is given) and is also reported in Appendix A. The covariance function associated with the geostrophy component is plotted for a given point (5° N, 210° E) in Fig. 4a and c, shown as a function of space (bottom left panel) and as a function of time (top left panel). This covariance function is similar to what is currently used for altimetry mapping with DUACS.

In the present study, we simultaneously estimate the surface signatures of the geostrophy and equatorial tropical instability waves (TIWs) and Poincaré waves. As for the geostrophy component, the equatorial wave covariances are expressed as a reduced wavelet basis, with typical wavelength and propagation speed given in the literature (e.g. Shinoda et al., 2009; Farrar, 2008, 2011; Farrar and Durland, 2012; Tanaka and Hibiya, 2019). For Poincaré waves, we built an ensemble of wavelets between 10° S and 10° N which follow the dispersion relation (Matsuno, 1966):

$$ \omega = \sqrt{k^2 \cdot c^2 + \beta \cdot c \cdot (2 \cdot n + 1)}, \quad (4) $$

where $\omega$ is the time frequency, $c = \pm 2.8 \text{ m s}^{-1}$ is the Poincaré wave propagation speed (considered a constant here), $k$ is the spatial wavenumber, and $n$ is a positive integer defining the wave mode. The wavelets are localized with a Hamming window with half-widths of 1000 km in the zonal direction, 300 km in the meridional direction, and 5 d in the temporal direction. For the TIW component, we also built an

ensemble of wavelets between 10° S and 10° N which follow the dispersion relation (Matsuno, 1966):

$$ \omega = c \cdot k, \quad (5) $$

where $\omega$ is the time frequency, $c = -0.5 \text{ m s}^{-1}$ is the TIW propagation speed (considered a constant here), and $k$ is the spatial wavenumber. The wavelets are localized here with a Hamming window with half-widths of 500 km in the zonal direction, 300 km in the meridional direction, and 20 d in the temporal direction. The covariance function for a westward-propagation-wave-like TIW is illustrated in Fig. 4b and d for a given point (5° N, 210° E), shown as a function of space (bottom right panel) and as a function of time (top right panel). Note that for Poincaré waves, both eastward and westward propagation is considered. A more detailed description of the equatorial wave’s components implemented in MIOST is provided in Appendix A.

# 3 Experiments and validation metrics

## 3.1 Experiments

We produced 4 years (from 1 July 2016 to 30 June 2020) of SSH maps using the MIOST multiscale and 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, and Sentinel-3B missions, the Level-3 Arctic lead sea-level anomaly products from SARAL/AltiKa, Sentinel-3A, and CryoSat-2 missions, and geostrophic current anomaly data from the AOML drifter database. These MIOST products are available on the AVISO + (Archivage, Validation et Interprétation des données des Satellites Océanographiques) website (see Sect. 5, “Data availability”, for more details).

Specific maps were also made to quantitatively assess the quality of these MIOST products. Table 3 summarizes the
<figure><img src="image_6.png" /><figcaption>Figure 4. Example of spatio-temporal covariance models at (5° N, 210° E) for (a, c) the geostrophy component and for (b, d) a westward-propagating wave component, e.g. TIW.</figcaption></figure>

list of experiments conducted in this study, indicating the input data used in the mapping and the physical content of the maps.

DUACS allsat-1 and MIOST allsat-1 experiments focus on the geostrophic variability. These SSH maps were produced from six altimeters (Jason-3, Cryosat-2, Sentinel-3A, Sentinel-3B, Haiyang-2A, Haiyang-2B) for the period 1 January 2019 to 31 December 2019, excluding one altimeter (Saral/AltiKa, over open-ocean region) from the mapping to perform independent assessments. The MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiment focuses on the geostrophic and equatorial wave variabilities. This experiment is based on (1) 80% of the drifter data, (2) the six altimeters previously mentioned over ocean, and (3) lead altimeter observations. The Saral/AltiKa dataset (over open-ocean region) and the remaining 20% of the drifter trajectories were here excluded from the mapping to perform independent assessments. Note that for these specific maps, drifter trajectories between -10° S and 10° N (e.g. Fig. 3) were also excluded to evaluate only the impact of the equatorial wave's mode in this region.

## 3.2 Validation 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 three main steps: (1) the SSH gridded data are 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>along-track</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 [65–500 km] 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 >500 km). In the study, the validation metric is based on the error variance scores in 1° × 1° longitude × latitude boxes (or averaged over a specific region of interest), defined as

$$ \sigma_{\text{err}}(x, y) = \frac{\sum_{t=1}^{N} (\text{SSH}_{\text{error}}(x, y, t) - \overline{\text{SSH}_{\text{error}}(x, y, t)})^2}{N}. \quad (6) $$
<table><caption>**Table 3.** List of mapping experiments with the input data and physical content considered.</caption><thead><tr><th rowspan="2">Experiment</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>DUACS allsat-1</td><td>All w/o AltiKa</td><td>No</td><td>No</td><td>Yes</td><td>No</td></tr><tr><td>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>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></tbody></table>

The similar statistical analysis can also be performed on the geostrophic velocity errors $U_{error} = U_{map} - U_{drifter}$ for the zonal component and $V_{error} = V_{map} - V_{drifter}$ for the meridional component.

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_{\text{err}}(\text{EXP2}) - \sigma_{\text{err}}(\text{EXP1})}{\sigma_{\text{err}}(\text{EXP1})} \quad (7) $$

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>along-track</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), and (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 are selected to compute the mean power spectral densities of the independent signal (SSH<sub>along-track</sub>) and the mapping error (SSH<sub>map</sub> − SSH<sub>along-track</sub>). Before the spectral calculation, the signals are detrended, and a Hanning window is applied. The signal-to-noise (SNR) ratio (Eq. 8) is then derived from the power spectral density (PSD) of the along-track SSH (SSH<sub>along-track</sub>) and the PSD of the error (SSH<sub>map</sub> − SSH<sub>along-track</sub>). As in Ballarotta et al. (2019), the effective resolution is then given by the wavelength $\lambda$s where the SNR($\lambda$s) is 2 (Eq. 9), i.e. the wavelength where the

SSH<sub>error</sub> is 2 times lower than the signal SSH<sub>along-track</sub>.

$$ \text{SNR}(\lambda) = \frac{\text{PSD}(\text{SSH}_{\text{along-track}})(\lambda))}{\text{PSD}(\text{SSH}_{\text{error}})(\lambda)} \quad (8) $$

$$ \text{SNR}(\lambda s) = 2 \quad (9) $$

# 4 Results

## 4.1 Qualitative assessment

Here we qualitatively assess the gridded products from the DUACS allsat-1 and MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiments. The SLA maps from the DUACS and MIOST mapping approaches are relatively similar in the subpolar region, as illustrated in Fig. 5 by an example of SLA reconstruction on 15 February 2019 for (a) the DUACS mapping approach and (b) the MIOST mapping approach. More significant differences take place in the Arctic basin: in contrast to the DUACS products, the use of Arctic lead observations in MIOST offers the possibility to extend sea-level mapping into ice-covered area and thus to deliver gap-free maps to end users (Fig. 5b).

From a global perspective, the MIOST maps are slightly more energetic than the DUACS maps as illustrated in Fig. 6 with the variance maps and their differences. The difference between MIOST and DUACS variance maps (Fig. 6c) 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, and coastal and polar regions. Tropical ocean regions are prone to lower SSH variability (10%) in the MIOST maps than in the DUACS maps.

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 [10° S–10° N], [180°–280° E]) for the period 2008 to 2018, for (a)
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example of sea-level anomaly maps on 15 February 2019 over the Arctic region constructed with the DUACS mapping approach (a) and with the MIOST mapping approach (b). The black line contour indicates the 15% sea-ice concentration from the OSI-SAF product.</figcaption></figure>

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

DUACS, (b) MIOST with equatorial wave modes, and (c) the GLORYS12V1 reanalysis (Lellouche et al., 2018). The rapid equatorial wave dynamics are resolved in the GLORYS12v1 ocean numerical simulation (Fig. 7c): the zonal wavenumber-frequency spectrum of the SSH in the Pacific reveals significant spectral peaks at periods close to 4, 5, and 7 d for a wavelength >20° in longitude. These peaks are associated with inertia-gravity (Poincaré) waves. These SSH variabilities for timescales smaller than 10 d are filtered in the DUACS mapping approach (Fig. 7a). In contrast, the MIOST multiscale mapping approach (MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic) resolves spectral peaks near 4, 5, and 7 d for wavelengths >20° in longitude (Fig. 7b). We show in the next section that these equatorial wave modes in MIOST also contribute to a significant reduction of the mapping error in this region. For timescales >10 d, each dataset has relatively similar spectral contents, particularly the energetic westward propagation of equatorial Rossby waves for negative wavenumbers.

## 4.2 Quantitative assessment

### 4.2.1 Mesoscale mapping assessments

The first assessment is a comparison of the DUACS allsat-1 and MIOST 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 but are used as independent validation.

#### Sea-level anomaly quality

The largest SSH mapping error $\sigma_{\text{err}}$ in DUACS allsat-1 reaches 50–100 cm² in the western boundary surface current and over the continental plateaus (Fig. 8a and b). In the offshore low-variability region, the error variance is <10 cm². Figure 8c and d show the difference in mapping error be-
<figure><img src="image_9.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) the GLORYS12V1 reanalysis. White lines represent the theoretical dispersion relation curves for equatorial waves corresponding to the Kelvin [1], Yanai [2], Rossby [3], and Poincaré [4] waves.</figcaption></figure>

tween the MIOST allsat-1 and DUACS allsat-1 experiments for all spatial scales and the spatial scale between 65 and 500 km, respectively. A blue (red) pattern means a reduction (increase) of the mapping error in MIOST compared to DUACS. For all spatial scales considered, MIOST mapping errors are smaller than those of DUACS, especially at mid-latitudes, 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 intertropical region, MIOST and DUACS have similar scores. For spatial scales between 65 and 500 km, MIOST mapping errors are reduced by ~10% compared to DUACS in high-variability regions at mid-latitudes. 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, the Siberian plateau, and the New Zealand plateau. Table 4 summarizes 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.

## Geostrophic current quality

Figure 9a and b 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 400 cm$^2$ s$^{-2}$ in the west-

ern boundary surface current (e.g. the Gulf Stream and the Kuroshio, Mozambique, and Agulhas currents). In offshore low-variability regions, the error variance is <80 cm$^2$ s$^{-2}$. The differences in mapping error between MIOST and DUACS are shown in Fig. 9c and d 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 summarizes 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.

## 4.2.2 Contribution of equatorial wave modes and drifters' observations

The comparison of the MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiment with MIOST allsat-1 examines the impact of the equatorial waves' mode and the drifters' observations in the MIOST mapping approach.

## Sea-level anomaly quality

The differences in mapping error between MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and MIOST allsat-1 are shown in Fig. 10a and b for all spatial scales and the spatial scale between 65 and 500 km, respectively.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Variance of the difference SSH<sub>map</sub> – SSH<sub>along-track</sub> computed for the DUACS allsat-1 experiment and in considering (a) all spatial scales and (b) the spatial scale between 65 and 500 km. Gain/loss of the mapping error variance of SLA in the MIOST allsat-1 experiment relative to the DUACS allsat-1 mapping error variance for (c) all spatial scales and (d) the scale between 65 and 500 km. Blue colour means a reduction of error variance in MIOST.</figcaption></figure>

<table><caption>**Table 4.** Regionally averaged mapping error variance and gain/reduction of error variance on the SSH variable between MIOST and DUACS.</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">All spatial scales</th><th colspan="3">Spatial scale [65–500 km]</th></tr><tr><th>Error variance DUACS [cm²]</th><th>Error variance MIOST [cm²]</th><th>Gain/loss error variance MIOST vs DUACS [%]</th><th>Error variance DUACS [cm²]</th><th>Error variance MIOST [cm²]</th><th>Gain/loss error variance MIOST vs DUACS [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>23.18</td><td>23.17</td><td>-0.02</td><td>7.07</td><td>6.84</td><td>-3.31</td></tr><tr><td>Antarctic</td><td>33.07</td><td>31.13</td><td>-5.86</td><td>7.86</td><td>7.65</td><td>-2.64</td></tr><tr><td>Equatorial band</td><td>14.07</td><td>13.96</td><td>-0.80</td><td>4.66</td><td>4.67</td><td>+0.32</td></tr><tr><td>Low variability – offshore</td><td>12.54</td><td>11.81</td><td>-5.83</td><td>3.70</td><td>3.55</td><td>-4.11</td></tr><tr><td>High variability – offshore</td><td>30.87</td><td>27.71</td><td>-10.22</td><td>14.28</td><td>12.87</td><td>-9.86</td></tr></tbody></table>

For all spatial scales considered, we observe that the equatorial wave 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 parts of Africa, and South America) are prone to deterioration. This suggests 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 percent of difference in the mapping error variance), with a reduction of error variance mainly in the high-variability regions. For a spatial scale between 65 and 500 km (Fig. 10b), the equatorial wave modes deteriorate the mapping solution in the western and central equatorial Pacific Ocean and 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 summarizes the results of the comparison over different regions of interest (Arctic, Antarctic, equatorial band, low-variability region, and high-variability region).

## Geostrophic current quality

The differences in mapping error of surface geostrophic currents between MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and MIOST allsat-1 are shown in Fig. 11a and b for the zonal component and the meridional component of the velocity, respectively. It is difficult to draw con-
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Variance of the difference $U_{\text{map}} - U_{\text{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 the MIOST allsat-1 experiment relative 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 [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Error variance MIOST [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Gain/loss error variance MIOST vs DUACS [%]</th><th>Error variance DUACS [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Error variance MIOST [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Gain/loss error variance MIOST vs DUACS [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>153.17</td><td>148.78</td><td>-2.87</td><td>133.50</td><td>131.34</td><td>-1.62</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>130.36</td><td>128.52</td><td>-1.42</td><td>124.36</td><td>123.20</td><td>-0.94</td></tr><tr><td>High variability – offshore</td><td>385.86</td><td>372.40</td><td>-3.49</td><td>409.75</td><td>403.54</td><td>-1.51</td></tr></tbody></table>

conclusions from this diagnosis: the mapping errors are reduced with MIOST in some regions in the tropics (such as the Bay of Bengal) and in the Kuroshio extension. Overall, the contribution of drifters remains moderate for the restitution of geostrophic currents (only a few percent improvement in the open ocean), as summarized in Table 7.

### 4.2.3 Overall assessment

The comparison of the MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and DUACS allsat-1 experiments allows the complete MIOST product distributed to users to be evaluated against the DUACS method.

#### Sea-level anomaly quality

The differences in mapping error between MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and DUACS allsat-1 are shown in Fig. 12a and b for all spatial scales and the spatial scale between 65 and 500 km, respectively. We have the same pattern as found in the previous sections: for all spatial scales considered (Fig. 12a), the equatorial wave modes help to reduce the mapping error variance in the equatorial band by more than 20% locally. At mid-latitudes, the mapping error is between 5% and 10% smaller with MIOST than with DUACS. For spatial scales between 65 and 500 km, 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
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Gain/loss of the mapping error variance of SLA in the MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiment relative to the MIOST allsat-1 mapping error variance for (a) all spatial scales and (b) the scale between 65 and 500 km. Blue colour means a reduction of error variance in MIOST when drifters are included in the mapping and with equatorial wave parametrization.</figcaption></figure>

<table><caption>**Table 6.** Regionally averaged mapping error variance and gain/reduction of error variance on the SSH variable between MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and MIOST allsat-1.</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">All spatial scales</th><th colspan="3">Spatial scale [65–500 km]</th></tr><tr><th>Error variance MIOST allsat-1 [cm²]</th><th>Error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic [cm²]</th><th>Gain/loss error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic vs MIOST allsat-1 [%]</th><th>Error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic [cm²]</th><th>Error variance MIOST allsat-1 [cm²]</th><th>Gain/loss error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic vs MIOST allsat-1 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>23.17</td><td>23.18</td><td>+0.02</td><td>6.84</td><td>6.84</td><td>+0.00</td></tr><tr><td>Antarctic</td><td>31.13</td><td>31.14</td><td>+0.02</td><td>7.65</td><td>7.65</td><td>+0.01</td></tr><tr><td>Equatorial band</td><td>13.96</td><td>13.53</td><td>-3.03</td><td>4.67</td><td>4.69</td><td>+0.32</td></tr><tr><td>Low variability – offshore</td><td>11.81</td><td>11.72</td><td>-0.77</td><td>3.55</td><td>3.54</td><td>-0.10</td></tr><tr><td>High variability – offshore</td><td>27.71</td><td>27.42</td><td>-1.06</td><td>12.87</td><td>12.67</td><td>-1.54</td></tr></tbody></table>

than with DUACS in regions where the circulation interacts with bathymetry features such as in the Argentine Sea and near the Siberian plateau and New Zealand plateau. Table 8 summarizes the results of the comparison over different regions of interest: mapping errors are ~11% smaller in high-variability regions in MIOST than in DUACS. In other regions, the errors are ~3%–6% smaller.

## Geostrophic current quality

The differences in mapping error of surface geostrophic currents between MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and DUACS allsat-1 are shown in Fig. 13a and b 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).

## Effective resolution

The effective spatial resolution quantifies the minimum spatial scale resolved in the maps (Ballarotta et al., 2019). Maps of the effective spatial resolution (expressed in kilometres) are presented in Fig. 14a and b for DUACS allsat-1 and MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic, respectively. For each experiment, the effective spatial resolution varies from ~500 km at the Equator to ~100 km at high latitudes and a mean value at mid-latitudes close to 200 km. The difference in effective spatial resolution between the two experiments is shown in Fig. 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 (the Gulf Stream, 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.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Gain/loss of the mapping error variance of currents in the MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiment relative 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 wave parametrization.</figcaption></figure>

<table><caption>**Table 7.** Regionally averaged mapping error variance and gain/reduction of error variance on the surface currents between MIOST allsat-1 and MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic.</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 MIOST allsat-1 [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Gain/loss error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic vs MIOST allsat-1 [%]</th><th>Error variance MIOST allsat-1 [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic [cm<sup>2</sup> s<sup>-2</sup>]</th><th>Gain/loss error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic vs MIOST allsat-1 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>148.78</td><td>145.04</td><td>-2.51</td><td>131.34</td><td>127.83</td><td>-2.67</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>128.52</td><td>127.80</td><td>-0.56</td><td>123.20</td><td>122.04</td><td>-0.94</td></tr><tr><td>High variability – offshore</td><td>372.40</td><td>366.81</td><td>-1.50</td><td>403.54</td><td>400.90</td><td>-0.65</td></tr></tbody></table>

## 5 Data availability

The MIOST gridded products are hosted on the AVISO+ (Archivage, Validation et Interprétation des données des Satellites Océanographiques) website (https://doi.org/10.24400/527896/a01-2022.009, Ballarotta et al., 2022).

The reference DUACS maps are hosted on the EU Copernicus Marine Service portal (https://doi.org/10.48670/moi-00148, Pujol et al., 2022b). The multiscale and multivariate products are distributed on a regular grid: the spatial grid extends from 0 to 360° E in longitude and 80° S to 90° N in latitude, with a grid spacing of 0.1°; the temporal grid covers the period 1 July 2016 to 30 June 2020 with a time step of 1 d. The dataset is distributed in netCDF4 format. Each netCDF file contains six variables: sla, adt, ugosa, vgosa, ugos, and vgos (see the list of variables available in the MIOST product in Fig. 15).

## 6 Summary and conclusions

Ubelmann et al. (2021, 2022) evaluated the multiscale and multivariate mapping approach in the Observing System Simulation Experiment (OSSE) and the Observing System

Experiment (OSE) for the simultaneous mapping of mesoscale circulation, coherent internal tides, and 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 the Arctic sea-ice regions. We found that the Arctic lead SSH observations allow the monitoring coverage in this remote region to be significantly improved. The gap-free maps, proposed with MIOST, hence offer the opportunity to end users 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 note that drifter observations can potentially improve surface circulation in areas not or poorly sampled by altimeters. Therefore,
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Gain/loss of the mapping error variance of SLA in the MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiment relative to the DUACS allsat-1 mapping error variance for (a) all spatial scales and (b) the scale between 65 and 500 km. Blue colour means a reduction of error variance in MIOST.</figcaption></figure>

<table><caption>**Table 8.** Regionally averaged mapping error variance and gain/reduction of error variance on the SSH variable between MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic and DUACS allsat-1.</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">All spatial scales</th><th colspan="3">Spatial scale [65–500 km]</th></tr><tr><th>Error variance DUACS allsat-1 [cm²]</th><th>Error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic [cm²]</th><th>Gain/loss error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic vs DUACS allsat-1 [%]</th><th>Error variance DUACS allsat-1 [cm²]</th><th>Error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic [cm²]</th><th>Gain/loss error variance MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic vs DUACS allsat-1 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>23.18</td><td>23.18</td><td>+0.01</td><td>7.07</td><td>6.84</td><td>-3.31</td></tr><tr><td>Antarctic</td><td>33.07</td><td>31.14</td><td>-5.85</td><td>7.86</td><td>7.65</td><td>-2.63</td></tr><tr><td>Equatorial band</td><td>14.07</td><td>13.53</td><td>-3.81</td><td>4.66</td><td>4.69</td><td>+0.64</td></tr><tr><td>Low variability – offshore</td><td>12.54</td><td>11.72</td><td>-6.56</td><td>3.70</td><td>3.54</td><td>-4.20</td></tr><tr><td>High variability – offshore</td><td>30.87</td><td>27.42</td><td>-11.16</td><td>14.28</td><td>12.67</td><td>-11.24</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Gain/loss of the mapping error variance of currents in the MIOST allsat-1 80% drifters + equatorial waves + L3 Arctic experiment relative 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>

their impact on the sea-level reconstruction may be larger over periods 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 scales, 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 the observed SSH to be decomposed 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 with the tropical instability waves and Poincaré waves. We show that mapping these ocean surface variabilities from altimeter observations broadens the spectrum of mappable space-timescales 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
<table><caption>**Table 9.** Regionally averaged mapping error variance and gain/reduction of error variance on the surface currents between MIOST allsat-1 80 % drifters + equatorial waves + L3 Arctic and DUACS allsat-1.</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<br>DUACS allsat-1<br>[cm<sup>2</sup> s<sup>-2</sup>]</th><th>Error variance<br>MIOST allsat-1 80 %<br>drifters + equatorial<br>waves + L3 Arctic<br>[cm<sup>2</sup> s<sup>-2</sup>]</th><th>Gain/loss error variance<br>MIOST allsat-1 80 % drifters +<br>equatorial waves + L3 Arctic vs<br>DUACS allsat-1 [%]</th><th>Error variance<br>DUACS allsat-1<br>[cm<sup>2</sup> s<sup>-2</sup>]</th><th>Error variance<br>MIOST allsat-1 80 %<br>drifters + equatorial<br>waves + L3 Arctic<br>[cm<sup>2</sup> s<sup>-2</sup>]</th><th>Gain/loss error variance<br>MIOST allsat-1 80 % drifters +<br>equatorial waves + L3 Arctic vs<br>DUACS allsat-1 [%]</th></tr></thead><tbody><tr><td>Arctic</td><td>153.17</td><td>145.04</td><td>-5.31</td><td>133.50</td><td>127.83</td><td>-4.25</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 –<br>offshore</td><td>130.36</td><td>127.80</td><td>-1.96</td><td>124.36</td><td>122.04</td><td>-1.87</td></tr><tr><td>High variability –<br>offshore</td><td>385.86</td><td>366.81</td><td>-4.94</td><td>409.75</td><td>400.90</td><td>-2.16</td></tr></tbody></table>

<figure><img src="image_16.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 + L3 Arctic experiments and (c) gain/loss of effective resolution (in %) between MIOST and DUACS. Blue means finer resolution in MIOST than in DUACS.</figcaption></figure>

spatio-temporal coverage of the altimeter data allows large-scale waves of 4 d periods and longer to be sampled. 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.

This experimental product is currently available on the AVISO + (Archivage, Validation et Interprétation des données des Satellites Océanographiques) website (see the “Data availability” section for more details), but our results suggest that the multiscale and multivariate mapping approach is very promising for use in an operational context. It is also worth mentioning that several other global gridded products exist as an alternative to the DUACS/MIOST products which provide only the geostrophic part of the surface current. Examples of these other products that provide a broader spectrum of ocean surface current variability (e.g. the total surface currents) include (1) the Copernicus GLORYS12v1 global ocean reanalysis (Lellouche et al., 2018; https://doi.org/10.48670/moi-00021, Drévillon et al., 2022), (2) the Copernicus GLOBCURRENT product (Rio et al., 2014; https://doi.org/10.48670/moi-00050, Etienne, 2021b), or (3) the OSCAR product (Dohan, 2021; https://doi.org/10.5067/OSCAR-25F20 Dohan, 2021) distributed by the NASA-JPL Distributed Physical Oceanography Active Archive Center (PO. DAAC).

To conclude, these results pave the way for the exploration of new types of ocean signals that may eventually be mapped with MIOST from remote sensing and in situ observations. Future work could consist of enriching the MIOST components 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 lead products in the Southern Ocean (Auger et al., 2022), or by enhanc-
```json
    dimensions:
        time = 1;
        latitude = 1702;
        longitude = 3600;
        bounds = 2;
    variables:
        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 service 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 service 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(latitude);
        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(latitude, bounds);
```

<caption>Figure 15. List of variables available in the multiscale and multivariate product.</caption>
ing the SLA map content with a dynamical model approach (Ubelmann et al., 2015) or artificial intelligence methods (Beauchamp et al., 2020).

# Appendix A: Mapping approaches tested in this study

## A1 The optimal interpolation (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, \qquad (A1)$$

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

$$\mathbf{x}^a = \mathbf{B}\mathbf{H}^T(\mathbf{H}\mathbf{B}\mathbf{H}^T - \mathbf{R})^{-1}\mathbf{y}. \qquad (A2)$$

The observation vector $\mathbf{y}$ represents the SLA observations. The state vector $\mathbf{x}$ is the gridded SLA. The operator $H$ (formally a trilinear interpolator transforming the gridded state SLA to the equivalent along-track SLA) is not considered explicitly. The matrices $\mathbf{B}\mathbf{H}^T$ and $\mathbf{H}\mathbf{B}\mathbf{H}^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):

$$C(x, y, t) = \left(1 + \text{ar} + \frac{1}{6}(\text{ar})^2 - \frac{1}{6}(\text{ar})^3\right) e^{-\text{ar}} e^{-\left(\frac{t}{L_t}\right)^2} \qquad (A3)$$

$$r = \sqrt{\left(\frac{x - C_{px}t}{L_x}\right)^2 + \left(\frac{y - C_{py}t}{L_y}\right)^2}, \qquad (A4)$$

where, $x$, $y$, and $t$ correspond to the zonal, meridional, and temporal position; $L_x$, $L_y$, and $L_t$ are the zonal, meridional, and temporal decorrelation scale; $C_{px}$ and $C_{py}$ denote the phase speed; and $a$ is a constant (3.337).

This covariance model is mainly optimized for mesoscale signal reconstruction. The $\mathbf{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 d in time (Le Traon

and Dibarbouré, 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 (Le Traon et al., 1998), 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}{f_c} \frac{\partial \text{SSH}(x, y)}{\partial y} \qquad (A5)$$

$$V_g(x, y) = \frac{g}{f_c} \frac{\partial \text{SSH}(x, y)}{\partial x}, \qquad (A6)$$

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

## A2 A multiscale and multivariate mapping approach

The optimal interpolation requires the inversion of a matrix of the same size as the observation vector $\mathbf{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 Sherman–Morrison–Woodbury transformation, allowing for an inversion in state space, with a matrix of the size of the state vector $\mathbf{x}$:

$$\mathbf{x}^a = (H^T[\mathbf{R}^{-1}H + \mathbf{B}^{-1}])^{-1}H^T\mathbf{R}^{-1}\mathbf{y}. \qquad (A7)$$

The formulation of the multiscale and multivariate mapping algorithm is detailed in Ubelmann et al. (2022). Here we recall the main principle. We consider an extended state vector $\mathbf{x}$ composed by $N$ physical components that will later be assumed independent. In this study $N = 3$ for (1) geostrophy and equatorial waves, (2) tropical instability waves (TIWs), and (3) Poincaré waves:

$$\mathbf{x} = (x_1^T, \dots, x_N^T)^T. \qquad (A8)$$

Each component $x_k$ represents the state of the surface topography and surface current to be resolved in the grid space, denoted $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$, which 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, \qquad (A9)$$

where $\eta_k$ is the reduced state vector for component $k$, and $\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 the geostrophy component is considered to have zero contribution to 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. \qquad (A10)$$

The observation vector $\mathbf{y}$ is also extended to the observed surface topography and surface current noted $\mathbf{y} = (h^{\mathrm{T}}, u_r^{\mathrm{T}})^{\mathrm{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 $\mathbf{y}$ is the sum of all component contributions plus the unexplained signal $\epsilon$ (instrument error and representativity):

$$\mathbf{y} = \sum_{k=1}^{N} G_k \eta_k + \epsilon. \qquad (A11)$$

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

$$y = G\eta + \epsilon. \qquad (A12)$$

Applying the same transformation from Eqs. (A1), (A2), and (A7) to the reduced state vector $\eta$, the global solution is written as

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

where $\mathbf{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 (A14)$$

In practice, to solve Eq. (A13), 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 Eq. A13 explicitly). We use a precondition 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 $\mathbf{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$ and then of $G^T$. Once the solution $\eta$ is obtained, the projection

in physical grid space given by Eq. (A14) 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.

## A2.1 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. Following the formulation provided in Ubelmann et al. (2021), here we define the gridded variable $H_1$ to resolve, and the corresponding gridded geostrophic current field ($U_1, V_1$) is written

$$\left\{ \begin{aligned} 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{aligned} \right. \qquad (A15)$$

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 the standard covariance models used in altimetry mapping to be approximated, 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) \times f_{\text{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), \qquad (A16)$$

where the $i$th 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 that all subcomponents are defined by pairs of sine and cosine functions to allow for 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 and 900 km with a spacing inversely proportional to the wavelet extensions, allowing a signal of any intermediate wavelength to be represented). ($x_p, y_p, t_p$) are the coordinates of a space-time pavement. The function $f_{\text{tap}}$ localizes the subcomponent in time and space (at scales $L_{t_p}, L_{x_p},$ and $L_{y_p}$, respectively) as geostrophy has a local extension of covariances. It is expressed as

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

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 timescale. Then, the same element p of the decomposition also has an expression in the geostrophic current (through the geostrophic relation
Eq. A15) written in the $\Gamma_{1,u}$ and $\Gamma_{1,v}$ matrices:

$$ \left\{ \begin{array}{l} \Gamma_{1,u}[i, p] = -\frac{g}{f_c} \frac{\partial \Gamma_{1,h}[i, p]}{\partial y_i} \\ \Gamma_{1,v}[i, p] = \frac{g}{f_c} \frac{\partial \Gamma_{1,h}[i, p]}{\partial x_i} \end{array} \right. \qquad (A18) $$

The whole time-space domain is paved with similar subcomponents, along coordinates ($x_p, y_p, t_p$) for wavelengths between 80 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.

## A2.2 Equatorial wave component

Here we define the gridded variables $H_2$ and $H_3$ to resolve TIW and Poincaré waves, respectively, and we consider no contributions of the equatorial wave components to the geostrophic currents; therefore the corresponding gridded geostrophic current fields ($U_2, V_2$) and ($U_3, V_3$) are written $U_2 = U_3 = 0$, $V_2 = V_3 = 0$. The reduced state is represented in the time-space domain by the following $\Gamma_{2,h}$ and $\Gamma_{3,h}$ matrix:

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

$$ \Gamma_{3,h}[i, p] = \cos(\omega_{3, t,p}(t_i - t_p) - k_{3, x,p}(x_i - x_p)) \times f_{\text{tap}} \left( \frac{x_i - x_p}{L_{3, x_p}}, \frac{y_i - y_p}{L_{3, y_p}}, \frac{t_i - t_p}{L_{3, t_p}} \right), \qquad (A20) $$

where $k_{2, x,p}$ and $k_{3, x,p}$ refer to the zonal wavenumber, and $\omega_{2, t,p}$ and $\omega_{3, t,p}$ are the frequency which satisfies the dispersion relation of the wave component (Matsuno, 1966); for example,

$$ \left\{ \begin{array}{l} \omega_{2, t,p} = c_2 \cdot k_{2, x,p} \text{ for the TIW, } c_2 = -0.5 \, \text{m s}^{-1} \\ \omega_{3, t,p} = \sqrt{k_{3, x,p}^2 \cdot c_3^2 + \beta \cdot c_3 \cdot (2 \cdot n + 1)} \end{array} \right. \qquad (A21) $$

for the Poincaré waves, $c_3 = \pm 2.8 \, \text{m s}^{-1}$, where $c_2$ and $c_3$ denote the wave propagation speed (the sign indicating the direction of propagation, negative for westward and positive for eastward), $\beta$ is the meridional gradient of the Coriolis frequency $f_c$, and $n = 1, 2, 3...$.

In the present study, we chose $L_{2, t_p} = 20 \, \text{d}$, $L_{2, x_p} = 500 \, \text{km}$, and $L_{2, y_p} = 300 \, \text{km}$ for the TIW component and $L_{3, t_p} = 5 \, \text{d}$, $L_{3, x_p} = 1000 \, \text{km}$, and $L_{3, y_p} = 300 \, \text{km}$ for the equatorial Poincaré wave component. As for the geostrophy component, the function $f_{\text{tap}}$ localizes the subcomponent in time and space (at scales $L_{t_p}$, $L_{x_p}$, and $L_{y_p}$, respectively).

**Author contributions.** CU implemented the multiscale and multivariate algorithm. PV and PP prepared the Arctic lead sea-level

dataset, and HE and SM prepared the drifter dataset. MB carried out the experiments and prepared the manuscript and figures. YF, GD, RM and NP participated in the discussion and interpretation of the results. All authors checked and provided related comments for this work.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** The authors would like to thank the AVISO + (Archivage, Validation et Interprétation des données des Satellites Océanographiques) team for their support and expertise in the distribution of the dataset. We are grateful to the three anonymous reviewers for their comments and suggestions to improve the manuscript.

**Financial support.** The work presented here was carried out in the framework of the DUACS-R&D project funded by CNES.

**Review statement.** This paper was edited by Giuseppe M. R. Manzella and reviewed by three anonymous referees.

## References

Arhan, M. and Colin de Verdière, A.: Dynamics of eddy motions in the Eastern North Atlantic, J. Phys. Oceanogr., 15, 153–170, 1985.

Auger, M., Prandi, P., and Sallée, J. B.: Southern Ocean sea level anomaly in the sea ice-covered sector from multimission satellite observations, Sci. Data, 9, 70, https://doi.org/10.1038/s41597-022-01166-z, 2022.

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-2022, 2019.

Ballarotta, M., Ubelmann, C., Rogé, M., Fournier, F., Faugère, Y., Dibarboure, G., Morrow, R., and Picot, N.: Dynamic Mapping of Along-Track Ocean Altimetry: Performance from Real Observations, J. Atmos. Ocean. Tech., 37, 1593–1601, https://doi.org/10.1175/JTECH-D-20-0030.1, 2020.

Ballarotta, M., Ubelmann, C., Veillard, P., Prandi, P., Etienne, H., Mulet, S., Faugère, Y., Dibarboure, G., Morrow, R., and Picot, N.: Gridded Sea Level Height and geostrophic velocities computed with Multiscale Interpolation combining altimetry and drifters, CNES/CLS [data set], https://doi.org/10.24400/527896/a01-2022.009, 2022.

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

Bojinski, S., Verstraete, M., Peterson, T., Richter, C., Simmons, A., and Zemp, M.: The Concept of Essential Climate Variables in Support of Climate Research, Applications, and Policy, B. Am. Meteorol. Soc., 95, 1431–1443, 2014.

Dohan, K.: Ocean Surface Current Analyses Real-time (OSCAR) Surface Currents – Final 0.25 Degree (Version 2.0), Ver. 2.0, PO.DAAC, CA, USA [data set], https://doi.org/10.5067/OSCAR-25F20, 2021.

Drévillon, M., Lellouche, J.-M., Régnier, C., Garric, G., Bricaud, C., Hernandez, O., and Bourdallé-Badie, R.: Global Ocean Reanalysis Products, Copernicus Marine [data set], https://doi.org/10.48670/moi-00021, 2022.

Ducet N., Le Traon, P. Y., and Reverdin, G.: Global high-resolution mapping of ocean circulation from the combination of T/P and ERS-1/2, J. Geophys. Res., 105, 19477–19498, 2000.

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.: Quality Information document for global ocean Multi Observation products, https://catalogue.marine.copernicus.eu/documents/QUID/CMEMS-MOB-QUID-015-004.pdf (last access: 9 January 2023), 2021a.

Etienne, H.: Global Ocean Multi Observation Products, Copernicus Marine [data set], https://doi.org/10.48670/moi-00050, 2021b.

Etienne, H., Verbrugge, N., Boone, C., Rubio, A., Solabarrieta, L., Corgnati, L., Mantovani, C., Reyes, E., Chifflet, M., Mader, J., and 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 (last access: 9 January 2023), 2021.

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

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

Farrar, J. T. and Durland, T. S.: Wavenumber–Frequency Spectra of Inertia–Gravity and Mixed Rossby–Gravity Waves in the Equatorial Pacific Ocean, J. Phys. Oceanogr., 42, 1859–1881, 2012.

Farrar, J. T. and Durland, T. S.: Equatorial waves across the Pacific (and Indian and Atlantic), 2022 Ocean Surface Topography Science Team Meeting, Venice, Italy, 31 October–4 November 2022, https://doi.org/10.24400/527896/a03-2022.3475, 2022.

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., and Thépaut, J.-N.: ERA5 hourly data on single levels from 1979 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set], https://doi.org/10.24381/cds.adbb2d47, 2018.

International Altimetry Team: Altimetry for the future: Building on 25 years of progress Advances in Space Research, 68, 319–363, https://doi.org/10.1016/j.asr.2021.01.022, 2021.

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

Le Traon, P.-Y. and Dibarboure, G.: Velocity Mapping Capabilities of Present and Future Altimeter Missions: The Role of High-Frequency Signals, J. Atmos. Ocean. Tech., 19, 2077–2088, https://doi.org/10.1175/1520-0426(2002)019<2077:VMCOPA>2.0.CO;2, 2002.

Le Traon P.-Y., Nadal, F., and Ducet, N.: An Improved Mapping Method of Multisatellite Altimeter Data, J. Atmos. Ocean. Tech., 15, 522–534, 1998.

Le Traon P.-Y., Faugere, Y. Hernamdez, F., Dorandeu, J., Mertz, F., and Abalin, M.: Can We Merge GEOSAT Follow-On with TOPEX/Poseidon and ERS-2 for an Improved Description of the Ocean Circulation?, J. Atmos. Oceanic Technol., 20, 889–895, 2003.

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

Lumpkin, R. and Centurioni, L.: Global Drifter Program quality-controlled 6-hour interpolated data from ocean surface drifting buoys, NOAA National Centers for Environmental Information [data set], 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., and Zaron, E. D.: 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, 2019.

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.

Matsuno, T.: Quasi-geostrophic motions in the equatorial area, J. Meteor. Soc. Japan, 44, 25–43, 1966.

Mulet, S., Etienne, H., Ballarotta, M., Faugere, Y., Rio, M. H., Dibarboure, G., and Picot, N.: Synergy between surface drifters and altimetry to increase the accuracy of sea level anomaly and geostrophic current maps in the Gulf of Mexico, Advances in Space Research, 68, 420–431, https://doi.org/10.1016/j.asr.2019.12.024, 2021a.

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, 2021b.

Prandi, P., Poisson, J.-C., Faugère, Y., Guillot, A., and Dibarboure, G.: Arctic sea surface height maps from multi-altimeter combination, Earth Syst. Sci. Data, 13, 5469–5482, https://doi.org/10.5194/essd-13-5469-2021, 2021.

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., Taburet, G., and SL-TAC team: Sea Level TAC – DUACS products, Copernicus Marine [data set], https://doi.org/10.48670/moi-00146, 2022a.

Pujol, M.-I., Taburet, G., and SL-TAC team: Sea Level TAC – DUACS products, Copernicus Marine [data set], https://doi.org/10.48670/moi-00148, 2022b.

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

Rio, M.-H., Mulet, S., and Picot, N.: Beyond GOCE for the ocean circulation estimate: Synergetic use of altimetry, gravimetry, and in situ data provides new insight into geostrophic and Ekman currents, Geophys. Res. Lett., 41, 8918–8925, https://doi.org/10.1002/2014GL061773, 2014.

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

Taburet, G., Pujol, M.-I., and the SL-TAC 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 (last access: 9 January 2022), 2021.

Tanaka, Y. and Hibiya, T.: Generation Mechanism of Tropical Instability Waves in the Equatorial Pacific Ocean, J. Phys. Oceanogr., 49, 2901–2915, 2019.

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

Ubelmann, C., Dibarboure, G., Gaultier, L., Ponte, A., Ardhuin, F., Ballarotta, M., and Faugère, Y.: Reconstructing ocean surface current combining altimetry and future spaceborne Doppler data, J. Geophys. Res.-Oceans, 126, e2020JC016560, https://doi.org/10.1029/2020JC016560, 2021.

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., 18, 469–481, https://doi.org/10.5194/os-18-469-2022, 2022.