The new CNES-CLS18 global mean dynamic topography

The mean dynamic topography (MDT) is a key reference surface for altimetry. It is needed for the calculation of the ocean absolute dynamic topography, and under the geostrophic approximation, the estimation of surface currents. CNES-CLS mean dynamic topography (MDT) solutions are calculated by merging information from altimeter data, GRACE, and GOCE gravity field and oceanographic in situ measurements (drifting buoy velocities, hydrological profiles). The objective of this paper is to present the newly updated CNES-CLS18 MDT. The main improvement compared to the previous CNES-CLS13 solution is the use of updated input datasets: the GOCO05S geoid model is used based on the complete GOCE mission (November 2009–October 2013) and 10.5 years of GRACE data, together with all drifting buoy velocities (SVP-type and Argo floats) and hydrological profiles (CORA database) available from 1993 to 2017 (instead of 1993–2012). The new solution also benefits from improved data processing (in particular a new wind-driven current model has been developed to extract the geostrophic component from the buoy velocities) and methodology (in particular the computation of the medium-scale GOCE-based MDT first guess has been revised). An evaluation of the new solution compared to the previous version and to other existing MDT solutions show significant improvements in both strong currents and coastal areas.

1 Introduction

The estimation of an accurate mean dynamic topography (MDT) has been a long-standing issue for the reconstruction of the absolute dynamic topography from altimeter data (Rio, 2010). The lack of an accurate geoid at spatial scales corresponding to the along-track spatial resolution of altimeter data (7 km at 1 Hz, 300 m at 20 Hz) has led to the exploitation of the time-variable part of the sea level with respect to the sea level mean over a given reference period: the sea level anomaly (SLA). Consequently, for years, the use of altimetry for scientific ocean studies has focussed on the analysis of sea level anomalies. While providing invaluable insight into the ocean dynamics of mesoscale eddies, a large number of scientific and operational activities rely on the accurate estimate of the absolute sea level. For instance, Pegliasco et al. (2020) show that eddies are better tracked and studied using the absolute sea level instead of the sea level anomalies. Also, the MDT is the missing component for the optimal assimilation of altimeter data into operational ocean systems, such as those run under the Copernicus Marine Environment Monitoring Services (CMEMS). Most importantly, the absolute dynamic topography is directly linked, under the geostrophic assumption, to the ocean surface currents. Those are required for a wide range of applications, such as the management of fishery resources, the monitoring of potential pollution and maritime security.

Over the years, the MDT has benefitted from a number of important improvements, the major one being the measurement of the Earth geoid at 100 km resolution with centimetric accuracy thanks to the launch in 2009 of the first ESA Earth Explorer, GOCE (Gravity and Ocean Circulation Experiment). Another key element has been the accumulation and improved processing over years of altimeter data leading to significant improvements in the mean sea surface (MSS; Pujol et al., 2018). At scales larger than 100 km, the mean dynamic topography can indeed be estimated by differencing both surfaces (MSS and satellite-only geoid height) and applying the adapted filter. This method is usually referred to as the “geodetic” approach (Bingham et al., 2008). To estimate the MDT spatial scales shorter than 100 km, other information is needed. A first approach could be to add small scales to the “satellite-only” geoid model using mainly altimetry over sea to end up with a so called “combined” geoid model such as GOCO05C (Fecher et al., 2017), GECO (Gilardoni et al., 2016) or Eigen6C4 (Förste et al., 2014) among others. For instance, the geodetic DTU19 MDT is based on a combined geoid model (Knudsen et al., 2019a, b). The other approach used to calculate the CNES-CLS MDT series (Rio and Hernandez, 2004; Rio et al., 2011, 2014a) is to add small scales to a satellite-only geodetic MDT using oceanic in situ data (temperature and salinity profiles and surface drifters), processed to extract the useful information (full dynamic height from the hydrological profiles, and geostrophic component of the drifter-derived velocities).

This paper presents the latest CNES-CLS18 MDT solution. The method is reviewed in Sect. 2, while data used in the computation are presented in Sect. 3. The geodetic MDT used as first guess is described and validated in Sect. 4. The processing of surface drifters has been greatly improved: Sect. 5 explains the full processing including the new wind slippage correction and the new wind-driven current estimates, both used to remove the ageostrophic component of the drifter currents. Section 6 describes the processing of in situ measurements of dynamic heights before being merged with the first-guess and processed surface drifters through objective analysis (Sect. 7). The accumulation of in situ measurements from different international programs, such as WOCE (for the SVP drifter) and ARGO (for the hydrological profiles), allows us today for the first time to map the global MDT at $\mathrm{1}/\mathrm{8}$^∘ resolution (instead of $\mathrm{1}/\mathrm{4}$^∘ for previous MDTs), and to significantly improve its accuracy in coastal areas, as will be highlighted in the validation section (Sect. 8). Conclusions and recommendations for future work are provided in Sect. 9.

2 Method

The method used to calculate the CNES-CLS18 mean dynamic topography is similar to the one used in previous CNES-CLS MDT versions. A detailed description can be found in Rio and Hernandez (2004) and Rio et al. (2011, 2014a). It is a three-step approach.

A first MDT solution is calculated from the optimally filtered differences between an altimeter MSS and a satellite-only geoid model. The effective resolution of the obtained field depends on the level of noise of the raw difference between MSS and geoid height. It thus depends on the areas, but it is around 100–125 km (Bruinsma et al., 2014).

In the second step of the method, synthetic estimates of the MDT and mean geostrophic velocities are calculated using in situ measurements of the ocean dynamic heights and surface velocities. First the in situ measurements are processed so as to extract the geostrophic component only from the drifting buoy total velocities, and to complete the dynamic heights with the missing barotropic and deep baroclinic components. The temporal variability of the measured heights and velocities is further removed by subtracting the altimeter sea level and geostrophic velocity anomalies, respectively. The processed in situ measurements are further averaged into $\mathrm{1}/\mathrm{4}$^∘ and $\mathrm{1}/\mathrm{8}$^∘ boxes to obtain the synthetic mean heights and velocities, respectively.

These are finally used in the third step to improve the accuracy of the filtered geodetic MDT obtained in step 1 and add information at shorter scales. This is done through a multivariate objective analysis whose required inputs are the synthetic mean heights and velocities and their error, the first-guess MDT, the a priori knowledge of the MDT variance, and zonal and meridional correlation scales.

3 Data

The CNES-CLS18 MDT is calculated from a combination of altimeter and space gravity data, in situ measurements and model winds. The following datasets are used:

• MSS. The CNES-CLS15 MSS derived for the 1993–2012 reference time period by Pujol et al. (2018) is used.

• Geoid model. The satellite-only geoid model GOCO05s (Mayer-Gürr et al., 2015) is used with the CNES-CLS15 MSS in the computation of the MDT first guess.

• Altimeter sea level anomalies (SLAs). The DUACS-2018 (Taburet et al., 2019) multimission gridded sea level and derived geostrophic velocity anomaly products distributed by the Sea Level Thematic Assembly Center (SL-TAC) from the CMEMS altimeter are used.

• Dynamic heights. These are calculated from temperature and salinity ($T/S$) profiles from CORA4.2 (1993–2013), CORA5.0 (2014–2015) and CORA5.1 (2016) datasets (Cabanes et al., 2013; Szekely et al., 2019), processed by the In Situ Thematic Assembly Center (INS-TAC) of the Copernicus Marine Environment and Monitoring Service (CMEMS).

• In situ velocities. Two types of in situ drifting buoy velocities are used, the 6-hourly SVP-type drifter distributed by the Surface Drifter Data Assembly Center (SD-DAC; Lumpkin et al., 2013) and the Argo floats surface velocities from the regularly updated YOMAHA07 dataset for the period 1993–2016 (Lebedev et al., 2007). SVP-type drifters consist of a spherical buoy with a drogue attached in order to minimize the direct wind slippage and follow the ocean currents at a nominal 15 m depth. When the drogue gets lost, the drifter is then advected by the surface currents and the direct action of the wind. In this study both the 15 m drogued and the surface undrogued drifter velocities over the 1993–2016 period are considered for the MDT calculation. The year 2017 is used for independent validation of the results.

• Wind data. Wind stress data are needed for the calculation of the wind-driven velocities (Sect. 5) that is used to remove part of the ageostrophic component from drifter velocities. We use the 3-hourly, 80 km resolution wind stress fields from ERA-Interim (Dee et al., 2011) for the period 1993–2017.

• Mixed layer depth (MLD). In the computation of the wind-driven velocities, an estimation of the MLD is also needed. We use the delayed time 3D temperature and salinity fields from ARMOR3D, computed by the Multi Observation Thematic Assembly Centre (MOB-TAC) of CMEMS (Guinehut et al., 2012).

4 First-guess calculation and comparison with previous first guess

The raw difference between the CNES-CLS15 MSS and GOCO05s geoid height is filtered using an updated version of the optimal filter fully described in Rio et al. (2011). The a priori statistics (errors and variance) are estimated using the MDT CNES-CLS13 (Rio et al., 2014a). The correlation radii (x₀, y₀) are the same as those used in the computation of the first guess of MDT CNES-CLS13 (Rio et al., 2014a).

The mean geostrophic currents associated with the resulting first guess are compared with mean geostrophic currents estimated from drifters (Sect. 5) and filtered at 80 km of resolution (the smallest scales that GOCE is able to resolve considering its lowest orbit at the end of the mission). Comparison is done outside of the equatorial band [−5; 5] where the geostrophic approximation is not valid. The root mean square difference (RMSD) for the zonal component is almost everywhere lower than 30 % of the RMS of the drifters (Fig. 1a). The RMSD is higher for the meridional component (RMSDV) since the signal is smaller and thus more difficult to extract from both geodetic MDT and drifters (smaller signal-to-noise ratio). RMSDV is often less than 60 % but reaches 100 % in the middle of the Pacific (Fig. 1b). Figure 2 illustrates the improvement compared with the previous first guess estimated for the MDT CNES-CLS13. The two plots show the differences between RMSD of the MDT CNES-CLS13 first guess and RMSD of the MDT CNES-CLS18 first guess for both components of the mean geostrophic velocities. The RMSD is reduced almost everywhere (reddish colors). For instance, both components are improved along the coast in the Gulf of Maine and the meridional component is improved along the Chilean coast (the improvement in these areas is further described in Sect. 8.1.2 and 8.1.3, respectively).

Figure 1 Root mean square differences between (a) zonal and (b) meridional mean geostrophic drifter velocities and the MDT CNES-CLS18 first-guess geostrophic velocities in percentage of drifter variance (%) in $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ boxes.

Figure 2 Difference (cm s⁻¹) of the RMSD obtained using the first guess of MDT CNES-CLS13 and the first guess of MDT CNES_CLS18 for the (a) zonal and (b) meridional components in $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ boxes. Only boxes with more than 500 points and where difference is statistically significant at more than 95 % (T test) are shown. The reddish colors denote improvement while bluish colors stand for degradation (see text for more details).

5 Drifter data processing

As in previous work (Rio and Hernandez, 2004; Rio et al., 2011, 2014a), the mean synthetic velocities U_synth are calculated by removing from the drifting buoy velocities U_buoy both the ageostrophic velocity components and the temporal variability $U_{\mathrm{alti}}^{\prime }$ of the geostrophic velocity components as measured by altimetry (Eq. 1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{U}_{\mathrm{synth}}& =U_{\mathrm{buoy}}-U_{\mathrm{ekman}}-U_{\mathrm{stokes}}-U_{\mathrm{inertial}}-U_{\mathrm{tidal}}\\ & -U_{\mathrm{ahf}}-U_{\mathrm{slip}}-U_{\mathrm{alti}}^{\prime }.\end{array}\end{array}$$

In order to yield a geostrophic surface drifter velocity, the processing includes different steps. First, the wind-driven component ($U_{\mathrm{w}}=U_{\mathrm{ekman}}+U_{\mathrm{stokes}}$) is estimated and removed (Sect. 5.1). Then the drifter velocities are corrected from the direct effect of the wind on the buoy (U_slip: wind slippage, which is significant only in the case of drogue loss). This is further described in Sect. 5.2. Third, the tidal and inertial velocities (U_tidal and U_inertial) as well as the residual high frequency ageostrophic signal (U_ahf) are removed (see Sect. 5.3). The temporal variability of the drifter geostrophic component is further removed using altimeter data ($U_{\mathrm{alti}}^{\prime }$). Finally, the resulting mean synthetic velocities are averaged into $\mathrm{1}/\mathrm{8}$^∘ boxes.

5.1 Wind-driven current modeling

As described in Chapron and the GlobCurrent Team (2017), under certain conditions (spatial homogeneity of the flow and stationary temporal forcing), the equilibrium between the Coriolis force and the friction force due to wind stress leads to the classical Ekman current formulation (Eq. 2):

$$\begin{array}{ll}{\displaystyle }{U}_{\mathrm{ekman}}& {\displaystyle }={\displaystyle \frac{\sqrt{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\omega })\mathit{\delta }}}{e}^{z/\mathit{\delta }}\left[{\mathit{\tau }}_{\mathrm{e}}^x\cos \left({\displaystyle \frac{z}{\mathit{\delta }}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.-{\mathit{\tau }}_e^y\sin \left({\displaystyle \frac{z}{\mathit{\delta }}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right)\right],\end{array}$$

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{V}_{\mathrm{ekman}}& ={\displaystyle \frac{\sqrt{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}(f+\mathit{\omega })\mathit{\delta }}}{e}^{z/\mathit{\delta }}\left[{\mathit{\tau }}_{\mathrm{e}}^x\sin \left({\displaystyle \frac{z}{\mathit{\delta }}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right)\right.\\ & \left.+{\mathit{\tau }}_{\mathrm{e}}^y\cos \left({\displaystyle \frac{z}{\mathit{\delta }}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\right)\right],\end{array}\end{array}$$

where τ_e is the effective wind stress, taking into account the effect of surface-velocity dependency (depending upon oceanic and atmospheric stability correction). f+ω is the modified vorticity with $\mathrm{2}\mathit{\omega }={\partial }_{x}v-{\partial }_{y}u$, the local vorticity of the underlying flow (u, v). δ is the upper stress-driven boundary layer, i.e., the Ekman layer depth; it can be expressed as follows (Pollard et al., 1973):

$$\begin{array}{}\text{(3)}& \mathit{\delta }={\mathit{\gamma }}^{\mathrm{1}/\mathrm{4}}{w}_{*}/\sqrt{fN},\end{array}$$

where γ≈0.2, $w_{*}=\sqrt{\mathit{\tau }/{\mathit{\rho }}_{\mathrm{0}}}$, f is the Coriolis parameter and N is the Brunt–Väisälä frequency.

Following this formulation, the Ekman response to wind stress at the surface is directed at 45^∘ to the right (left) of the wind direction in the Northern (Southern) Hemisphere. Ekman currents further evolve with depth following a logarithmic profile and a further rightward (leftward) spiral. However, near the surface, the classical Ekman spiral model is modified by the influence of surface gravity waves through the inclusion of the Stokes drift. Accordingly, the vertical profile of the quasi-Eulerian current, its magnitude and its surface angle depart from the classical model of Ekman (1905). In particular, Ardhuin et al. (2009) have shown that the surface angle can be much larger than the standard 45^∘ direction. The full surface Stokes drift is estimated to be on the order of 0.5 %–1.3 % of the wind speed and the quasi-Eulerian current, U_ekman+U_stokes is then on the order of 0.6 % of the wind speed and lies at an average angle between 40 and 70^∘, to the right of the wind direction (Northern Hemisphere). Near the surface, the surface Stokes drift induced by the waves typically accounts for two-thirds of the total surface wind-induced drift. In this section, the objective is to best estimate the total wind-driven component of the current ($U_{\mathrm{w}}=U_{\mathrm{ekman}}+U_{\mathrm{stokes}}$) to further remove it from the in situ drifting buoy velocities. For that purpose, we extend the approach described in Rio et al. (2014a) for Ekman current modeling to further include the Stokes drift component.

A two-parameter (β, θ) model (Eq. 4) is used.

$$\begin{array}{}\text{(4)}& u_{\mathrm{w}}\left(z\right)=\mathit{\beta }\left(z\right)e^{i\mathit{\theta }\left(z\right)}{\mathit{\tau }}^{c\left(z\right)}\end{array}$$

The β and θ parameters are obtained through a least square fit between the wind-driven velocity u_w(z) extracted from the drifting buoy measurements and the wind stress values interpolated at the buoy times and positions.

u_w(z) is obtained for z=15 m by removing the altimeter-derived geostrophic velocities from the 15 m drogued SVP drifter velocities and further applying a low-pass filter. The filtering cutoff length is taken as the maximum between the 24 h tidal period and the inertial period. The wind stress along the buoy trajectory is also filtered consistently.

u_w(z) is obtained for z=0 m by removing the altimeter-derived geostrophic velocities from the surface Argo float velocities. Unlike drogued SVP drifters, Argo float trajectories are not filtered since there is only one velocity estimation every 10 d. Note also that Argo floats were found to be much less affected by wind slippage compared with undrogued drifting buoys, certainly thanks to their design (Rio et al., 2014a). Consequently, we do not need to correct Argo float drift from wind slippage, and thus we can use Argo floats to estimate the wind-driven model at the surface (z=0), which are unlike undrogued drifting buoys that are affected by wind slippage (see Sect. 5.2).

In Rio et al. (2014a), c(z) was set to 1 at both level (c(z=0 m$)=c(z=\mathrm{15}$ m)=1). In this study we determined c(z) at each level so as to maximize the percentage of global variance explained by the model (Eq. 4). We found c(z=0 m)=0.6 and c(z=15 m)=0.7. For comparison, in the tropical band, Ralph and Niiler (1999) found an optimal value c(z=15 m)=0.6.

From Ekman theory (Eq. 2), we expect both β and θ parameters to present regional and seasonal variabilities, in correlation with the varying upper stress-driven boundary layer δ (which depends on latitude and ocean stratification). Indeed, following Eq. (2), in summer, when Ekman depth decreases, β increases and $\left|\mathit{\theta }\right|$ increases.

In order to take into account these variations, the approach used by Rio et al. (2011, 2014a) was to fit both parameters by month and by 4^∘ by 4^∘ boxes. In Rio et al. (2014a), the wind-driven response at the surface was found to be located at around 20–40^∘ to the right of the wind direction in the Northern Hemisphere (to the left in the Southern Hemisphere), and the angle then increases to 40–60^∘ at 15 m depth. In addition, a clear seasonal cycle was obtained for both parameters and at both depths with larger angles and amplitudes in summer than in winter in good consistency with stronger summer stratification.

In this paper, in order to further resolve the temporal variability of the β and θ parameters, we fit the model from Eq. (4) as a function of latitude (1^∘ bins) and MLD (5 m bins). Weekly MLD values are calculated from the ARMOR3D $T/S$ grids (Guinehut et al., 2012; see Sect. 3) using as criteria the minimum depth between the isotherm layer depth and the isopycnal layer depth (De Boyer Montégut et al., 2004). We checked (not shown) that β and θ values are symmetrical with respect to the Equator, so that the parameters for the Northern and Southern Hemisphere are fitted together (β and θ are computed as a function of MLD and |lat|).

The obtained β and θ parameters are displayed in Fig. 3. As expected, at a given latitude, both parameters increase with increasing stratification (smaller MLD values). We also obtain a smaller angle response at the surface compared to 15 m depth, and a larger-amplitude parameter, in good agreement with an Ekman-like spiral. Finally, as expected from the inverse dependency of the Ekman amplitude response with the Coriolis parameter f (Eq. 2), the β parameter decreases with increasing latitude. We also observe a decrease with latitude of the surface angle parameter. This might be due to an increase in the Stokes drift amplitude at high latitude (due to steeper wind waves). At 15 m depth, the angle response is not linear with latitudes but shows a more complex latitude dependency with a minimum near 30^∘. The angle is dependent on the Ekman depth (Eq. 2), which depends on many different parameters all varying with latitudes: the wind stress, the Coriolis parameter and the stratification (through the Brunt–Väisälä frequency).

Figure 3 The β (a, c) and θ (b, d) parameters calculated in this study as a function of mixed layer depth and |latitude|. The (a, b) (c, d) plots show the parameters at the surface (15 m depth). Only bins with more than 50 values are shown.

Another interesting aspect from Fig. 3 is that the latitudinal dependency almost vanishes for deep-MLD values.

This is expected from theoretical considerations: from Eqs. (2) and (3) we get

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\mathit{\beta }(z=\mathrm{0}\mathrm{m})={\displaystyle \frac{\sqrt{\mathrm{2}}\sqrt{N}\sqrt{\mathit{\tau }}}{{\mathit{\gamma }}^{\mathrm{1}/\mathrm{4}}\sqrt{{\mathit{\rho }}_{\mathrm{0}}}\sqrt{f}}},\text{(6)}& {\displaystyle }& {\displaystyle }\mathit{\beta }(z=\mathrm{15}\mathrm{m})={\displaystyle \frac{\sqrt{\mathrm{2}}\sqrt{N}\sqrt{\mathit{\tau }}}{{\mathit{\gamma }}^{\mathrm{1}/\mathrm{4}}\sqrt{{\mathit{\rho }}_{\mathrm{0}}}\sqrt{f}}}{e}^{\frac{\mathrm{15}\sqrt{fN}}{{\mathit{\gamma }}^{\mathrm{1}/\mathrm{4}}{w}_{*}}}.\end{array}$$

In the case of a well mixed upper layer (strong MLD values), the Prandtl number $N/f\approx \mathrm{1}$ so that the β value does not depend any more on the f values (i.e., β becomes independent of latitude).

Table 1 Percentage of variance explained by the surface wind-driven model.

Table 1 shows the percentage of variance explained by the new surface model, compared to the old one (Rio et al., 2014a).

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\mathit{\%}\mathrm{VarU}={\displaystyle \frac{\sum _i{\left(u_{\mathrm{w}}-u_{\mathrm{w}}^{\mathrm{in}\mathrm{situ}}\right)}^{\mathrm{2}}}{\sum _i{\left(u_{\mathrm{w}}^{\mathrm{in}\text{-}\mathrm{situ}}\right)}^{\mathrm{2}}}}\text{(8)}& {\displaystyle }& {\displaystyle }\mathit{\%}\mathrm{VarV}={\displaystyle \frac{\sum _i{\left(v_{\mathrm{w}}-v_{\mathrm{w}}^{\mathrm{in}\mathrm{situ}}\right)}^{\mathrm{2}}}{\sum _i{\left(v_{\mathrm{w}}^{\mathrm{in}\mathrm{situ}}\right)}^{\mathrm{2}}}}\end{array}$$

The new surface model is fitted using the entire 2007–2016 period dataset, whereas the old surface model was based only on the 2007–2014 period, so that results might be biased toward the new solution. To remove this ambiguity, we have repeated the new model parameter calculation using data from 2007 to December 2014. We find (not shown) that the obtained parameter values are very similar to the values obtained using the entire period. We then calculate the percentage explained in independent data (2015–2016). It is much improved for both components of the surface velocity compared to the previous model from Rio et al. (2014a), especially at latitudes smaller than |5^∘| (Table 1). Similar results are obtained for the 15 m depth model (Table 2).

Table 2 Percentage of variance explained by the 15 m depth wind-driven model.

5.2 Wind slippage

The SVP-type drifters distributed by the SD-DAC consist of a surface float connected to a subsurface 7 m long holey-sock drogue centered at 15 m depth. Such drifters are designed to minimize the direct action of the wind on the buoy (wind slippage) so that the buoy is advected by the ocean currents at 15 m depth (Niiler et al., 1995). When the drogue is lost, the buoy trajectory is due to both the surface currents and the wind slippage.

In order to correct the undrogued drifting buoy velocities for the wind slippage, an updated version of the method described in Rio et al. (2011) is used. It consists of finding the optimal α coefficient so as to minimize the correlation between the wind stress and a residual velocity U_r as defined in Eq. (9):

$$\begin{array}{}\text{(9)}& U_{\mathrm{r}}={\left(U_{\mathrm{buoy}}-U_{\mathrm{ekman}}-U_{\mathrm{stokes}}-U_{\mathrm{geo}}\right)}_f-\mathit{\alpha }W,\end{array}$$

where U_buoy is the drifter velocity, U_geo is the geostrophic current derived from altimetry and W is the wind speed.

In Rio et al. (2011), the vectorial correlation is computed on a 100 d moving window. So only trajectories with more than 100 d can be processed, and no wind slippage values are computed for the first and last 50 d of the trajectory, resulting in the loss of a high number of velocity observations.

We therefore upgrade the method in order to calculate the wind slippage correction for the whole drifter's trajectories.

To evaluate the impact of using different correlation window lengths (100, 50, 30, 20 and 10 d), we use a selection of 21 trajectories longer than 300 d. Results show that below 30 d, the computed alpha coefficient is too noisy, so we choose 30 d as the lower bound of the correlation window T_L: 30 d ≤ T_L≤100 d.

In order to process the first and last $T_{\mathrm{L}}/\mathrm{2}$ d of the trajectory portions and also the trajectories shorter than 60 d, two methods have been investigated:

• the use of a mean α value, calculated from the α obtained along the drifter trajectory;

• the use of a climatological α value – a mean α is computed for each drifter's trajectories longer than 60 d, and values are then averaged into $\mathrm{4}{}^{\circ }\times \mathrm{4}{}^{\circ }$ spatial bins to yield climatological α values.

To validate the different wind slippage corrections, we compare the geostrophic altimetric velocity U_geo to ($U_{\mathrm{buoy}}-U_{\mathrm{ekman}}-U_{\mathrm{stokes}}{)}_f-{\mathit{\alpha }}_{xx}W$, where the subscript xx is set to

• CLIM for climatological α,

• MEAN for mean α value,

• 0 when no wind slippage correction is applied (α₀=0) and

• N for the nominal α (T_L=100 d).

We compute root mean square differences for both components of the velocity (RMSDU and RMSDV).

The first expected result is that applying a wind slippage correction always yields better results than not applying any correction: RMSDU (RMSDV) is reduced from 13 % to 7 % (from 10 % to 6.5 %). The second result is that using a mean α along the drifter trajectory yields better results than using the climatological α value. The RMSDU and RMSDV are reduced by 2.2 %.

Consequently, we finally apply the following procedure to remove the wind slippage from the drifting buoy velocities: the first and last $T_{\mathrm{L}}/\mathrm{2}$ d of the trajectories longer than 60 d are completed using α_MEAN. Once all drifter trajectories longer than 60 d have been processed, α_CLIM is used to correct the trajectories shorter than 60 d. Note that a higher error is associated with the velocity estimate at the beginning and end of the trajectory since the wind slippage correction is less accurate. It corresponds to 21 % of the total amount of drifters for the mean α and 1.5 % for the climatological α.

Figure 4 Mean zonal (a) and meridional (b) synthetic velocities.

5.3 Synthetic mean velocities calculation

The synthetic mean velocities are then calculated following Eq. (1): the 15 m drogued SVP velocities are corrected from the wind-driven component using the empirical model from Eq. (4) and the β(z=15 m), θ(z=15 m), c(z=15 m) parameters. Undrogued SVP-type drifters are first corrected from the wind slippage (Sect. 5.2) and then from the wind-driven component using the empirical model from Eq. (4) and the β(z=0 m), θ(z=0 m), c(z=0 m) parameters. In order to remove the tidal and inertial components from the drogued and undrogued SVP drifters, a low-pass filter is then applied along the drifter trajectory. In past work (Rio and Hernandez, 2004; Rio et al., 2011, 2014a) a unique 3 d filtering length was used. Here we refine the filtering length by considering for each drifter trajectory the lowest latitude sampled (absolute value) and filter the trajectory using a cutting period P_c=max(P_i, 24 h), where P_i is the inertial oscillation period at the lowest latitude, so as to be sure to filter the main tidal currents, the inertial oscillations and any residual high-frequency ageostrophic signal. Note that an upper bound is set at 6 d for the cutting period to avoid overly high values at low latitudes.

The Argo floats' surface velocities are corrected from the surface wind-driven model. Argo floats surface velocities being sampled only every 10 d, no further low-pass filtering can be applied on the Argo velocity dataset. Also, as discussed in Rio et al. (2014a) and mentioned earlier, no wind slippage correction is applied to Argo floats since, thanks to their design, they are less impacted by wind slippage than undrogued SVP.

Temporal variability is then removed by interpolating altimetric geostrophic velocity anomalies along the drifter trajectories and subtracting them from the in situ geostrophic velocities. The obtained mean geostrophic velocities are further averaged into $\mathrm{1}/\mathrm{8}$^∘ boxes (instead of $\mathrm{1}/\mathrm{4}$^∘ boxes in Rio et al., 2014a) and displayed in Fig. 4. Associated errors are obtained as described in Rio et al. (2011).

6 Dynamic heights data processing

The method used to calculate the synthetic MDT (MDT_synth) is the same as the one fully described in Rio et al. (2011) and also used in Rio et al. (2014a). We use the temperature (T) and salinity (S) profiles from the CORA database (Sect. 3) to calculate instantaneous dynamic heights h(t,r)_Pmax relative to a maximum profile depth P_max. Following the same philosophy as in the computation of synthetic velocities, the temporal variability is subtracted from $h(t,r{)}_{/P\max }$ using altimetric anomalies (SLA) referenced over the 1993–2012 time period to end up with a mean field $<h>(r{)}_{/P\max }$ also referenced over the 1993–2012 time period. The resulting mean dynamic height $<h>$ represents only the baroclinic processes from the surface down to P_max. We need to add the mean contribution of deep baroclinic and barotropic processes. The missing quantity is evaluated as the difference between a dynamic height climatology referenced at P_max ($h_{\mathrm{clim}/P\max }$) and the CNES-CLS18 MDT first guess (Eq. 10). To avoid any large-scale differences and time period discrepancy between $<h>(r{)}_{/P\max }$ and the climatology $h_{\mathrm{clim}/P\max }$, a specific climatology is computed using the same observations from CORA while a WOA climatology was used by Rio et al. (2011).

$$\begin{array}{}\text{(10)}& {\mathrm{MDT}}_{\mathrm{synth}}=<h>(r{)}_{/P\max }-h_{\mathrm{clim}/P\max }+\mathrm{first}\mathrm{guess}\end{array}$$

The calculation of the CNES-CLS18 MDT is based on a remove–restore technique where the MDT first guess is first removed and then added back to the optimally estimated field. Thus, the useful quantity in Eq. (10) is the difference between $<h{>}_{/P\max }$ and its climatology $h_{\mathrm{clim}/P\max }$, i.e., the small scales of the mean dynamic heights referenced at P_max. Indeed, the aim of the objective analysis is to optimally add scales smaller than those resolved by the MDT first guess. Most of the $T/S$ profiles go down to 2000 m; thus Fig. 5 illustrates the difference between $<h{>}_{/P\max }$ and its climatology $h_{\mathrm{clim}/P\max }$ for P_max=2000 m. The signal is high for instance in the Antarctic Circumpolar Current with anomalies higher than 10 cm and in the equatorial currents system area mainly in the western Pacific with anomalies around ±7 cm. We thus expect that these in situ data will bring useful information in these areas.

Figure 5 $<h>\left(r\right){/}_{\mathrm{2000}\mathrm{m}}-h_{\mathrm{clim}/\mathrm{2000}\mathrm{m}}$ (cm).

Finally, the synthetic mean heights are averaged in $\mathrm{1}/\mathrm{4}$^∘ boxes. Note that, unlike the drifters, we do not have enough observations of height to average in $\mathrm{1}/\mathrm{8}$^∘ boxes. The associated errors are computed as described in Rio et al. (2011). These are dominated by oceanic variability patterns, with errors up to 10 cm in high-variability areas (western boundary current, Antarctic Circumpolar Current – ACC) and less than 2 cm far from these areas.

Figure 6 (a) The CNES-CLS18 mean dynamic topography (cm) and (b) intensity of their associated mean geostrophic velocities (cm s⁻¹).

7 High-resolution MDT inversion

The synthetic mean geostrophic velocities and mean heights calculated in the previous sections are then used to improve the GOCE-based first guess using a multivariate objective analysis as in Rio and Hernandez (2004) and Rio et al. (2011, 2014a). The method relies on the prescription of the a priori MDT variance and the a priori zonal and meridional spatial correlation scales of the estimated field. We use the same a priori statistics as in Rio et al. (2014a). The obtained CNES-CLS18 MDT is displayed in Fig. 6a.

Estimates of the mean geostrophic velocities are obtained as output of the multivariate objective analysis. In the equatorial band, only the synthetic mean velocities are used for the inversion in order to circumvent the issue due to the failure of the geostrophic approximation. For the same reason, only the synthetic height is used to estimate MDT in the equatorial band. Figure 6b shows the norm of the obtained mean currents. All currents are significantly enhanced compared to the speed of the mean currents from the GOCE-based first guess (not shown).

A spectral analysis tool is used to investigate the spatial scales of the different MDT solutions: first guess, CNES-CLS13 MDT, new CNES-CLS18 MDT and Glorys12 numerical model MDT (a $\mathrm{1}/\mathrm{12}$^∘ numerical model from Mercator Ocean; Artana et al., 2019). An example is given in Fig. 7 for a 10^∘ box in the Agulhas Return Current area. We clearly see the increased energy at short scales contained in the new CNES-CLS18 MDT (in light blue) compared to either the GOCE-based first guess (in purple) or the previous CNES-CLS13 MDT (in red). The energy level is in good agreement with the GLORYS12 energy level (in green) until 50 km wavelength (25 km resolution). At shorter scales a flat response is obtained. This means that the shortest scales (between 12.5 and 25 km) of the CNES-CLS18 MDT might be contaminated by noise. Figure 7 shows also the power spectra density of the geodetic DTU19 MDT (Knudsen et al., 2019a, b) based on a combined geoid model. At scales shorter than 200 km, DTU19 is less energetic than CNES CLS18 MDT. A visual inspection of the geostrophic current associated with the MDTs confirms that the CNES-CLS18 is more energetic than the DTU19 (not shown). This underlines the value of using oceanic in situ data in the CNES-CLS18 MDT compared with a combined geodetic MDT.

Figure 7 Meridional spectral power density obtained for different MDT solutions in a 10^∘ box centered around (49^∘ S, 43^∘ E).

8 Validation

The validation of the obtained MDT-CNES-CLS18 solution uses three different approaches:

1. a qualitative validation in specific areas (Japan Sea, Gulf of Maine, Chilean coast, Agulhas Current) to highlight the enhancement of coastal currents and western boundary currents in those areas compared to the previous CNES-CLS13 solution;

2. a quantitative regional validation around Australia and in the Drake Passage;

3. a global, quantitative validation by comparison to an independent dataset of drifting buoy velocities – we use the 15 m drogued SVP drifter velocities available for the year 2017 (not used in the calculation of the synthetic velocities).

8.1 Qualitative validation

8.1.1 The Japan Sea

The Japan Sea is characterized by a number of thin coastal currents as depicted for instance in Fig. 1 from Lee et al. (2009). Most of these currents are poorly resolved (or not resolved at all) by the $\mathrm{1}/\mathrm{4}$^∘ resolution CNES-CLS13 MDT (Fig. 8c). In particular, the Liman cold current (LC), the Soya Current (SC), the nearshore branch (NB) of the Tsushima Warm Current and the East Korean Warm Current (EKWC) are nicely resolved in the new CNES-CLS18 MDT. None of these are resolved in the GOCE-based first guess, so this small-scale information is only due to the inclusion of the in situ information. The Liman cold current was not at all resolved in the previous CNES-CLS13 solution.

Figure 8 Mean geostrophic velocities in the Japan Sea in the CNES-CLS18 GOCE-based first guess (a), in the CNES-CLS18 MDT (b) and in the CNES-CLS13 MDT (c). Abbreviations for the main currents are as follows – LC: Liman Current; EKWC: East Korean Warm Current; NB: nearshore branch of the Tsushima Warm Current; and SC: Soya Current.

8.1.2 Gulf of Maine

The freshwater mass enters the eastern Gulf of Maine from south of the Nova Scotia shelf and follows a general counterclockwise circulation with a significant contribution to the Maine Coastal Current (MCC). MCC generally flows westward along the coast between 43.5 and 44.5^∘ N and creates signals in the MDT: the MDT increases toward the coast. As shown by Feng et al. (2018), this coastal current is not resolved by the MDT CNES-CLS13 while it is resolved by an independently developed MDT_RU from a regional ROMS numerical model by Rutgers University that assimilates various in situ measurements such as $T/S$ profiles, HF radar and drifters (Levin et al., 2018).

An update of this study is done here including MDT CNES-CLS18 to verify if the new MDT resolves this coastal current pattern (Figs. 9 and 10). Figure 9 shows a comparison, at six buoy sites, of the mean geostrophic currents inferred from MDT CNES-CLS13 and MDT CNES-CLS18 with ADCP mean upper ocean currents (averaged from the surface down to 50 m and over the whole buoy-operated time period: roughly July 2001 to July 2019). Specifically, at the offshore-most buoy N deployed in the eastern side of the northeast channel (NEC), both MDTs agree very well with the buoy mean in both magnitude and direction. Along the Maine coast at buoys I, E and B, the improvement of MDT CNES-CLS18 is obvious although disagreement still exists to some degree at buoys I and B at the western coast. At buoy L, MDT CNES-CLS18 looks reasonable in magnitude and direction, but not consistent with buoy measurements, most likely because at this location ADCP operates over a shorter time period (July 2001–August 2008) leading to larger uncertainty. At buoy M deployed at the deeper Jordan Basin (JB), MDTs and the buoy mean are not in agreement, maybe because at this location the ratio of “temporal variability” to “mean signal” is high, and thus the comparison with currents averaged over a different time period is not accurate. Even if the mean currents from ADCP represent a mean over a different time period than the MDTs referenced over 1993–2012, these comparisons show that the mean MCC estimated by MDT CNES-CLS18 is significantly improved compared with the MDT CNES-CLS13. Indeed, Fig. 10 shows that the MCC has a signature in the MDT CNES-CLS18 that increases toward the coast, while this is not the case in the MDT CNES-CLS13.

Figure 9 Map of Gulf of Maine with the mean upper ocean (50 m averaged) current vectors from in situ measurements at six buoy sites (N, L, M, I, E and B) in blue, and the mean geostrophic velocity vectors inferred from MDT_CNES_CLS13 in red and MDT_CNES_CLS18 in black. One Jason altimeter track 24 (T24) is also shown. (Abbreviations – JB: Jordan Basin; WB: Wilkinson Basin; NEC: northeast channel).

Figure 10 (a) MDT CNES-CLS18 and (b) MDT CNES-CLS 13 in the Gulf of Maine.

As highlighted in Figs. 11 and 12, this improvement is due to the better representation of this current both in the first guess and in the mean synthetic velocities used in the computation. First, thanks to the improvement of the processing, the first guess is improved close to the coast and the new first-guess MDT increases toward the coast, which is not the case for the MDT CNES-CLS13 first guess (Fig. 11). Second, the number of available drifters, the improved processing of this dataset and their averaging being performed in $\mathrm{1}/\mathrm{8}$^∘ boxes instead of $\mathrm{1}/\mathrm{4}$^∘ boxes (Fig. 12) are all key elements leading to the improvement of the CNES-CLS18 in this specific area.

Figure 11 First guess of the MDTs (a) CNES-CLS18 and (b) CNES-CLS13 in the Gulf of Maine.

Figure 12 Zonal component of the mean synthetic velocities from drifters used in the computation of the (a) MDT CNES-CLS18 and (b) MDT CNES-CLS13.

8.1.3 Chilean coast

Over the south Chilean coast, the eastward South Pacific Current (SPC) bifurcates to become the Humboldt Current (HC) that flows northward and the Cape Horn Current (CHC) that flows southward. The circulation in this area is fully described by Strub et al. (2019). The HC and CHC are much better resolved in the MDT CNES-CLS18 compared with the previous MDT CNES-CLS13. The circulation resolved by the MDT CNES-CLS18 is thus in better agreement with what is expected and with numerical models (Fig. 13). As in the Gulf of Maine, this improvement of the MDT CNES-CLS18 is due to the new drifters dataset and the improved processing of the first guess. The mean circulation along the Chilean coast is difficult to accurately resolve using geodetic MDT due to the influence of the Andes and the very rugged coastline. The improved processing of the first guess helps to overcome this issue.

Figure 13 MDT and associated mean geostrophic current along the south Chilean coast from the MDT CNES-CLS13, MDT CNES-CLS18 and the numerical model HYCOM interpolated on the MDT CNES-CLS18 grid ($\mathrm{1}/\mathrm{8}$^∘).

However, offshore of Chiloe Island, the MDT CNES-CLS18 shows an eddy that does not appear to be realistic. In this area, between the HC and the CHC, one or two observations that have an unrealistic, strong velocity may introduce some noise. We expect that this should be found in other areas and can explain the noise underlined by the spectral analysis at the smaller scales (see Sect. 7).

8.1.4 Agulhas Current

As first demonstrated by Chapron et al. (2005), highly valuable information on ocean surface currents can be retrieved by analyzing the Doppler shift measured between the signal emitted by a synthetic aperture radar (SAR) antenna and the signal backscattered by the sea surface.

Within the ESA GlobCurrent project, mean velocities were calculated in the Agulhas Current area from the ENVISAT SAR Doppler velocities. The ENVISAT data span the period from 2007 to 2012. The main limitation of this technique is that only the radial component of the velocity is retrieved. To deal with this limitation, the total current velocities are estimated by using the direction given by the altimeter-derived geostrophic currents, whose norm is corrected using a radial component measured by the SAR from either the ascending or descending passes. The corrected norm using the ascending ($V_{\mathrm{a}}^{*}$) or the descending ($V_{\mathrm{d}}^{*}$) passes are given by Eqs. (11) and (12), respectively:

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{a}}^{*}={\displaystyle \frac{V_{\mathrm{a}}^{\mathrm{SAR}}}{\cos \left({\mathit{\beta }}_{\mathrm{a}}\right)}},\text{(12)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{d}}^{*}={\displaystyle \frac{V_{\mathrm{d}}^{\mathrm{SAR}}}{\cos \left({\mathit{\beta }}_{\mathrm{d}}\right)}},\end{array}$$

where $V_{\mathrm{a}}^{\mathrm{SAR}}$ and $V_{\mathrm{d}}^{\mathrm{SAR}}$ are the SAR-derived radial velocities in ascending and descending passes, and β_a and β_d are the angles between the SAR range direction and the altimeter-derived current direction for ascending and descending passes.

The obtained mean velocities corresponding to the 1993–2012 time period are displayed in Fig. 14c. Velocities along the African coast are in good agreement with the mean velocities obtained from drifters (Fig. 14b), with maximum values reaching up to 2 m s⁻¹ in the core of the current. This is around twice the intensity obtained from the GOCE-derived MDT first guess (Fig. 14a).

Figure 14 Mean velocities in the Agulhas Current area from (a) the GOCE-only derived MDT first guess, (b) drifting buoy velocities and (c) SAR-derived mean velocities.

By using the drifter velocities to improve the GOCE-derived first guess, much higher velocity speeds are obtained (Fig. 15b). Maximum values reach 1.6 m s⁻¹, which is around 30 % more than the mean velocities obtained in the previous CNES-CLS13 MDT (maximum velocities are around 1.2 m s⁻¹). This is due to the higher resolution of the solution ($\mathrm{1}/\mathrm{8}$^∘ instead of $\mathrm{1}/\mathrm{4}$^∘), an improved processing of the drifters (less filtering is applied) and the removal of a threshold for a minimum number of data which was applied in the CNES-CLS13 MDT calculation, removing a number of very coastal in situ data for the inversion.

In the return current, a discontinuity is observed in the CNES-CLS18 MDT in the speed intensity east and west of around 24^∘ E. This is in qualitative agreement with both the drifter data and the SAR measurements (Fig. 14). Instead, due to a coarser resolution, the Agulhas return current intensity is very much continuous along its path both in the CNES-CLS13 MDT (Fig. 15a) and the GOCE MDT first guess (Fig. 14a).

Figure 15 Mean geostrophic velocities in the Agulhas Current area derived from the CNES-CLS13 MDT (a) and the CNES-CLS18 MDT (b).

Figure 16 Mean geostrophic velocities around Australia in the CNES-CLS18 GOCE-based first guess (a), in the CNES-CLS18 MDT (b) and in the CNES-CLS13 MDT (c).

8.2 Quantitative regional validation

8.2.1 Australia

A number of coastal currents flow around Australia, in particular the Leuwin Current (along the west coast), the East Australian Current and the Hiri Current (along the northern tip of Australia). None of these currents are resolved in the GOCE-based first guess (Fig. 16a) while they are nicely resolved in the CNES-CLS18 MDT (Fig. 16b). The Leuwin Current was hardly resolved in the previous CNES-CLS13 MDT (Fig. 16c).

Figure 17 Ocean current roses of the (a) 30 d low-pass-filtered ADCP observations, (b) Copernicus GlobCurrent based on MDT CNES-CLS13 and (c) Copernicus GlobCurrent based on MDT CNES-CLS18. The numbers in blue indicate the angle of the ocean current direction (0^∘ is eastward, 90^∘ is northward). The numbers in black and the concentric circles give the number of occurrences in the observations for each angle. Note that for the middle plot, the scale is different regarding the occurrences.

Figure 18 Mean surface velocity magnitude (in cm s⁻¹) and mean surface velocity vectors from MDT CNES-CLS18. Black lines are isolines of mean sea surface height corresponding to the main Antarctic Circumpolar Current (ACC) fronts and their branches: SAF-N (23 cm), SAF-M (0 cm), PF-N (43 cm), PF-M (−70 cm), PF-S (−79 cm), SACCF-N (−93 cm) and SACCF-S (−114 cm). White dots indicate the position of the DRAKE moorings.

Figure 19 Mean surface geostrophic velocities (MSGVs) derived from MDT CNES-CLS09 (blue), from MDT CNES-CLS13 (green), MDT CNES-CLS18 (cyan) and LUT (red) along Jason track no. 104. Location of the Drake moorings are indicated M1 to M10. Below bathymetry along the section.

Cancet et al. (2019) validate the Copernicus GlobCurrent product that is distributed by CMEMS (https://resources.marine.copernicus.eu/documents/QUID/CMEMS-MOB-QUID-015-004.pdf, last access: 31 May 2021). They compare the version of the Copernicus GlobCurrent currents that use the MDT CNES-CLS13 (Rio et al., 2014a) with three ADCP time series located within the East Australian Current close to the coast. Figure 17 is an update of the results from Cancet et al. (2019) at the station SYD100 (33.9438^∘ S, 151.3821^∘ E), including a comparison with a version of Copernicus GlobCurrent based on the MDT CNES-CLS18 (only the MDT field is changed). When the MDT CNES-CLS13 is used, the current does not seem to have a predominant direction. On the contrary, when the MDT CNES-CLS18 is used, the ocean current rose is in good agreement with the ADCP data, with the main current flowing with a 250^∘ angle.

8.2.2 Drake Passage

The persistent westerly winds over the Southern Ocean drive the Antarctic Circumpolar Current (ACC) flow around the Antarctic continent. The ACC flow is organized in three oceanic frontal systems, which correspond to water mass boundaries as well as deep-reaching jets of eastward flow: the Subantarctic Front (SAF), the Polar Front (PF) and the Southern ACC Front (SACCF) (Fig. 18). The Drake Passage, located between the South American and Antarctic continents, is the narrowest stretch of water (about 850 km wide), separating Antarctica from the northern continents; thus it forms the most practical location to monitor the ACC (Fig. 18). In the context of the DRAKE project (Provost et al., 2011) hydrographic data and current meter time series were collected during 3 years (2006–2009) in the upper 3000 m of the water column below the ground track 104 of the Jason-1 altimetry satellite (Fig. 18, black section and white dots). Combining the DRAKE in situ mooring and satellite altimetry velocities Koenig et al. (2014) build a look-up table (LUT) to compute the ACC transport. Koenig et al. (2014) adjusted the mean cross-track surface velocity at each mooring location using an iterative error correction scheme with the CNES-CLS13 mean surface geostrophic velocities as a first guess. The three mean surface geostrophic velocities from the CNES-CLS09 (Rio et al., 2011), CNES-CLS13 (Rio et al., 2014a) and CNES-CLS18 MDTs are compared to the mean surface geostrophic velocities from Koenig et al. (2014) in Fig. 19. The new CNES-CLS18 MDT provides mean surface geostrophic velocities consistent with those from Koenig et al. (2014). Both present a strong SAF (located at 56.2^∘ S) and SACCF-S (60.5^∘ S). They show two recirculation cells: negative velocities around 57^∘ S and at 59.7^∘ S as identified in Ferrari et al. (2012, 2013) with the DRAKE current-meter data. This is an improvement as the recirculation at 57^∘ S is not captured in the velocities derived from the CNES-CLS09 and CNES-CLS13 MDTs. The mean velocities derived from the CNES-CLS18 MDT further indicate a strong current close to the coast reaching 35 cm s⁻¹ at 55^∘ S. This feature is observed in all the L-ADCP sections (Renault et al., 2011) and in the GLORYS12 reanalysis (Artana et al., 2019). We propose an adjustment of the mean geostrophic surface velocities from Koenig et al. (2014) (LUT – updated, red dashed line) close to the South American continent to account for this coastal current associated with a northern branch of the SAF. This coastal current carries slightly less than 2 Sv, increasing the mean Koenig at al. (2014) ACC transport to 143 Sv. The MDT CNES-CLS18 is also shown to perform better than the MDT CNES-CLS13 in the SAF at 41^∘ S (Malvinas Current), downstream of Drake Passage, in Artana et al. (2019).

8.3 Quantitative validation through comparison to independent drifter velocities

The number of available independent 15 m drogued drifters for year 2017 is shown in Fig. 20. The world ocean is rather well sampled, apart from the Indonesian throughflow, the Arctic Ocean and the Southern Ocean (latitudes poleward of 60^∘ S).

Figure 20 Trajectories of the 15 m drogued SVP drifters available over the year 2017. Colors give the speed velocities deduced along the drifter trajectories.

Absolute dynamic topography values are calculated by adding gridded sea level anomalies from DUACS-2018 (Sect. 3) to the new CNES-CLS18 MDT. Corresponding geostrophic velocities were then derived and interpolated along the drifter trajectories. Root mean square differences (RMSDs) between the obtained geostrophic velocities and the drifter-derived geostrophic velocities for year 2017 are then calculated both for the zonal and meridional components. The same comparison is performed using the CNES-CLS18 first guess, the previous CNES-CLS13 MDT and its first guess.

Figure 21 Zonal (a) and meridional (b) RMS differences as a function of coastal distance between the independent drifting buoy velocities and the altimeter geostrophic velocities obtained using different MDT solutions – pink: the CNES-CLS13 MDT first guess; blue: the CNES-CLS18 MDT first guess; black: the CNES-CLS13 MDT; and red: the CNES-CLS18 MDT.

The RMSD of the independent velocity dataset are then calculated as a function of distance to the coast (Fig. 21). Everywhere, the RMSD is reduced with the MDT CNES-CLS18 and especially in coastal areas (distances from the coast less than 100 km). Part of this improvement is due to first-guess improvement, since the RMSD is also reduced when the CNES-CLS18 first guess is used rather than the CNES-CLS13 first guess. This was already pointed out in the Gulf of Maine and along the Chilean coast. As described in Sect. 8.1, this improvement is also due to in situ data, especially to processed drifters.

9 Conclusions

The new CNES-CLS18 MDT is presented, focussing on the innovative elements in the processing and datasets used compared to the previous CNES-CLS13 MDT solution, and on validation results at global, regional and coastal scales. Among the main improvements, a more recent GOCE geoid model is used and the calculation of the first guess is improved. More in situ hydrological profiles and drifting buoy data are used, and a new wind slippage correction and wind-driven model are computed to better extract the geostrophic current from the drifters. Grid resolution is enhanced from $\mathrm{1}/\mathrm{4}$^∘ for earlier solutions to $\mathrm{1}/\mathrm{8}$^∘. The effective resolution is also improved as shown by spectral analysis. Thus, this new MDT better retrieves the strongest ocean currents and shows huge improvement in coastal areas, as highlighted by a number of specific coastal and regional examples. Indeed, narrow and coastal currents are now resolved, which were not in the previous solutions.

The continuous improvement of MDT accuracy and resolution is necessary, all the more so in the context of the upcoming swath altimetry missions (e.g., SWOT, whose launch is planned in November 2022) that will be able to retrieve smaller spatial scales of the sea surface height. Moreover, as underlined by Hamon et al. (2019), an accurate MDT and its error are important to correctly assimilate altimetric observations in numerical models. Further processing improvements might be considered. First, the computation of the first guess could be improved by using a new method such as the one described by Siegismund (2020). In particular, this could help to better resolve MDT along the coast. Second, we should address the residual noise issue highlighted in some areas at scales shorter than 25 km. This might be done by refining the errors prescribed on the drifter-derived velocities. We need also to better understand and quantify the residual temporal variability and how the spatiotemporal coverage of the in situ data may impact the final mean estimate. Third, other kinds of observations could be used. In this study SAR-derived velocity are used for validation purposes, but this kind of measurement could be included in the MDT inversion. In order to improve coastal area results, HF radar measurements are valuable, as Caballero et al. (2020) show in the Bay of Biscay. They process HF radar velocities to extract only the mean geostrophic part and show that, when including this information in the MDT computation, the local mean circulation is better resolved. Finally, a realistic MDT error should be estimated. A formal error is given as an output of the objective analysis. It highly depends on the parameters prescribed in input and it is known to be underestimated. We plan to use new method to estimate this field, for instance an ensemble method.