A multi-year timeseries of observation-based 3D horizontal and vertical quasi-geostrophic global ocean currents

Estimates of 3D ocean circulation are needed to improve our understanding of ocean dynamics and to assess its impact on marine ecosystems and Earth climate. Here we present the OMEGA3D product, an observation-based timeseries of (quasi) global 3D ocean currents covering the 1993-2018 period, developed by the Italian Consiglio Nazionale delle Ricerche within the European Copernicus Marine Environment Monitoring Service (CMEMS). This dataset was obtained by applying 10 a diabatic quasi-geostrophic (QG) diagnostic model to CMEMS data-driven ARMOR3D weekly reconstruction of temperature and salinity and ERA-Interim fluxes. Outside the equatorial band, vertical velocities were retrieved in the upper 1500 m, at nominal 1⁄4° resolution, and successively used to compute the horizontal ageostrophic components. Root mean square differences between OMEGA3D total horizontal velocities and totally independent drifter observations at two different depths (15 m and 1000 m) decrease with respect to corresponding estimates obtained from zero-order geostrophic balance, meaning 15 that estimated vertical velocities can also be deemed reliable. OMEGA3D horizontal velocities are also closer to drifter observations than velocities provided by a set of re-analyses spanning a comparable time period, but based on data assimilation in ocean general circulation numerical models. The full OMEGA3D product (released on 31st of March 2020) is available upon free registration at 20 https://doi.org/10.25423/cmcc/multiobs_glo_phy_w_rep_015_007. The reduced subset used here for validation and review purposes is openly available at https://doi.org/10.5281/zenodo.3696885 (Buongiorno Nardelli, 2020).

levels between 2.5 m and 1482.5 m. This grid was specifically designed to get more accurate numerical solutions within the 95 ocean upper boundary layer (Kalnay de Rivas, 1972;Sundqvist and Veronis, 1970). Pre-processing of input data thus includes as a first step the vertical interpolation of ARMOR3D data on OMEGA3D vertical layers (using python class scipy.interpolate.interp1d (Virtanen et al., 2019) set to cubic spline interpolation) and the mapping of ERA-Interim data on OMEGA3D horizontal grid (using python class scipy.interpolate.griddata (Virtanen et al., 2019), set to fit data to a piecewise cubic, continuously differentiable, curvature-minimizing polynomial surface). 100 As ARMOR3D data may occasionally display density inversions along the water column, that are not compatible with QG Omega solution, vertical profiles of potential density are adjusted to impose static stability: moving from the surface to depth, whenever a density inversion is found, density is set to the upper level value plus a 0.0001 kg/m 3 increment.

Quasi-Geostrophic equations
A diabatic Q-vector formulation of the quasi-geostrophic Omega equation (Buongiorno Nardelli et al., 2018b;Giordani et al., 105 2006) is solved to get the OMEGA3D vertical velocity fields: In Eq. (1), w represents the vertical velocity (positive upwards), N 2 is the Brunt-Väisälä frequency, f is the Coriolis parameter, 110 h indicates the horizontal components, and the Q vector includes three components reflecting different processes (kinematic deformation, twg, turbulent buoyancy, th, and turbulent momentum, dm), as defined below:

565
(2) In the above definitions, r indicates the potential density, g is the gravitational acceleration, ) and ) represent the geostrophic velocities, while turbulent terms are defined following the KPP parameterizations, namely through non-local 120 effective gradients * and classical viscosity/diffusivity Kx. The function used here to estimate KPP vertical mixing coefficients (Smyth et al., 2002) (only modified to handle derivatives in non-staggered non-uniform vertical grids) is designed to account for Langmuir cells mixing by including an amplification of turbulent velocity scales, and includes a nonlocal momentum flux term and a parameterization of Stokes drift effects. Forcing terms are computed from ARMOR3D potential density and geostrophic velocity fields, and ERA-interim atmospheric re-analyses. 125 In one of the analytical steps to obtain Eq. (1), the details on which are given elsewhere (Buongiorno Nardelli et al., 2018b;Giordani et al., 2006), the following two equations are found: Once vertical velocities are retrieved through Omega solution, these two equations allow to estimate also the horizontal ageostrophic components ( + and + ) through a simple vertical integration of the right-hand side terms (assuming horizontal 135 ageostrophic velocities are negligible at the bottom boundary).

Numerical solution
All equations used for the OMEGA3D retrieval are solved here numerically (Buongiorno Nardelli et al., 2018b). At each grid point in the interior domain (i.e. excluding the boundaries), the Omega equation is re-written substituting derivatives with central finite differences, considering a non-staggered grid. Vertical derivatives are computed considering a variable grid 140 spacing, increasing as the square of depth (Kalnay de Rivas, 1972), and adopting a second-order accuracy finite difference scheme (Sundqvist and Veronis, 1970).
At the surface and topographical boundaries Dirichelet conditions are imposed (namely zero vertical velocities) and Neumann conditions at the bottom and lateral boundaries (namely zero vertical velocity partial derivatives). These latter are imposed through forward/backward finite schemes and make the solution not suited to model current topography interactions along the 145 coasts. Grouping all factors multiplying w at each grid point one finally gets a set of equations that represent a closed linear system in w, which can thus be solved through a matrix inversion. Solving the linear system at once for the entire domain, however, would be computationally extremely demanding. As such, considering the elliptical nature of the Omega equation (which confines the impact of boundary conditions to a limited number of grid points), the original grid was split here in smaller horizontally overlapping sub-domains (each tile having a horizontal dimension of 75 grid points, overlapping by one 150 third). The inversion is carried out sequentially on these sub-domains, imposing the vertical velocity values that resulted from the previous step as lateral boundary conditions to the subsequent calculations. The algorithm used for the matrix inversion is Loose Generalized Minimum Residual, a.k.a., LGMRES (Baker et al., 2005) with incomplete lower-upper (LU) preconditioning (as implemented in python sparse linear algebra package scipy.sparse.linalg (Virtanen et al., 2019), imposing a tolerance for convergence of 10 -7 ). 155 Vertical velocities are finally used to integrate Eq. (3) by a simple trapezoidal rule to obtain ageostrophic horizontal velocities.

Model re-analyses used for intercomparison
Three different ocean state reconstruction timeseries have been compared with OMEGA3D. All of them are based on ocean general circulation models assimilating both in situ and satellite observations, though significantly differing in term of 160 numerical schemes used, input data ingested and assimilation strategies.
The first dataset considered is the third release of version 4 of ECCO Fukumori et al., 2018), hereafter ECCOv4r3, covering the 1992-2015 period, and available at https://ecco.jpl.nasa.gov/products/all/. ECCOv4r3 is based on the MIT General Circulation Model (Adcroft et al., 2001) and applies a 4D-VAR assimilation scheme to a wide set of observations (including satellite altimetry, in situ T/S profiles, satellite sea surface salinity and temperature, ocean bottom pressure) 165 (Fukumori et al., 2017), minimizing the observation-analysis misfits in a least squares sense (Wunsch and Heimbach, 2013).
ECCOv4r3 model grid includes 50 vertical levels, with a zonal resolution of 1° and a variable meridional resolution, ranging from 1° to approximately 0.25° in the Equatorial band and near the poles. ECCOv4r3 is thus a relatively coarse resolution ocean circulation model and the effect of mesoscale dynamics on vertical velocities is parameterized by introducing a ''bolus'' vertical velocity (Danabasoglu et al., 1994;Gent and Mcwilliams, 1990;Liang et al., 2017), that needs to be added to the large-170 scale Eulerian vertical velocity diagnosed from volume continuity. ECCOv4r3 vertical velocity data are released only as monthly averages (https://ecco.jpl.nasa.gov/products/V4r3/user-guide/). The third product used for the comparison is the output of the first version of the 1/12° horizontal resolution GLORYS system (Drévillon et al., 2018), hereafter GLORYS12v1, covering the period 1993-2018. GLORYS12v1 is obtained by jointly assimilating along track altimeter data, satellite SST, sea ice concentration and in situ temperature and salinity vertical profiles into a global ocean eddy-resolving model with 50 vertical levels. GLORYS12v1 model component is NEMO (Madec, 2016) and data assimilation is carried out by means of a reduced-order Kalman filter, while a 3D-VAR scheme provides a correction 185 for the slowly-evolving large-scale biases in temperature and salinity.
GLORYS12v1 data can be freely downloaded at: http://marine.copernicus.eu/services-portfolio/access-to-products/ (id=GLOBAL_REANALYSIS_PHY_001_030). GLORYS12v1 distributed output does not include vertical velocities. For the sake of comparison, GLORYS12v1 daily files have been subsampled here over the same dates for which SODAv3.4.2 is available. 190 As for OMEGA3D, ECCOv4r3, SODAv3.4.2 and GLORYS12v1 surface forcings are all taken from ERA-Interim (Dee et al., 2011).
In order to minimize wind slippage, SVP drifters are drogued with a 7 m long holey-sock centered at 15 m depth and their 200 velocity estimates is considered representative of currents at 15 m depth (Lumpkin et al., 2017). Before carrying out the validation, individual SVP drifter 6-hourly records have been averaged over a running time window (inversely scaled with the Coriolis parameter) to remove the signal due to inertial oscillations (Buongiorno Nardelli et al., 2018b). YOMAHA velocities are estimated by measuring the displacement of profiling Argo floats during their submerged phase (Lebedev et al., 2007). Argo floats drift at a pre-defined parking pressure and emerge only for near-real-time data transmission 205 through ARGOS/IRIDIUM satellites. Most of these instruments follow a profiling cycle of approximately 10 days, and their parking level is set to 1000 m.

Vertical velocity mean patterns and resolved variability
Vertical velocities cannot be measured in the open ocean due to their relatively small magnitude (of the order of 1-100 m d -1 , depending on depth and processes involved). Consequently, OMEGA3D vertical velocity cannot be directly validated (no 210 reference datasets exist). However, the algorithm used to retrieve OMEGA3D horizontal velocities requires the vertical velocity in input, and improvements in quasi-geostrophic horizontal components with respect to standard geostrophic velocities would necessarily imply that vertical velocity is reliable.
OMEGA3D vertical velocities are thus compared here with the output of the only two ocean climate reanalysis systems that presently distribute vertical velocity timeseries of comparable length. Considering that vertical velocity fields are provided at  Given its 5-days sampling, SODAv3.4.2 could be expected to reveal a stronger variability than OMEGA3D (7-days sampling), and both are expected to display much higher values than ECCOv4r3 (providing monthly averaged fields). Conversely, though associated patterns display very similar features, OMEGA3D standard deviation maximum value exceeds SODAv3.4.2 by a 250 factor ~2. In both cases, intense maxima are associated with the main current systems.
For the sake of a more consistent comparison with ECCOv4r3, OMEGA3D and SODAv3.4.2 vertical velocity standard deviations have also been estimated also after low-pass filtering the latter two time series (by a 5-point and 7-point moving window, respectively) to keep only frequencies lower than monthly as those provided by ECCOv4r3 (figure 2 should thus be compared to figure 1f). Even in that case, the variability observed in ECCOv4r3 is sensibly lower than these retrieved from 255 higher spatial resolution products, likely revealing the limits of the mesoscale parameterization used in ECCOv4r3 in terms of vertical exchanges.

Horizontal velocity validation vs independent observations
OMEGA3D horizontal velocity accuracy has been assessed in terms of mean bias and root mean square differences (RMSD) with respect to space-time co-located in situ reference observations. Estimated metrics have then been compared to that 260 estimated for geostrophic velocities directly obtained from DUACS altimeter data (AVISO+, 2015) (when looking at SVP velocities) or from ARMOR3D (Mulet et al., 2012) (when looking at YOMAHA velocities), and successively also compared to similar metrics computed from SODAv3.4.2 and GLORYS12v1 output.
To build our matchup databases, OMEGA3D velocities have been interpolated at the same nominal depth of drifter measurements through a weighted average of the two closest levels. 265 The first assessment covered surface currents as measured by SVP drifters. As SVP may occasionally loose their drogue, failing to correctly represent 15 m depth currents, only drogued SVP drifters data, collected within ±12 hours from nominal reconstruction dates, have been included in our matchup databases ( Fig.3a and 3b). The same matchup procedure has been applied to DUACS geostrophic velocities, to SODAv3  fig.6b and 6d) present very minor differences, while GLORYS12v1 presents significantly higher differences (by a factor ~2) along all major currents ( fig.4f). Even stronger discrepancies affect SODAv3.4.2 estimates, displaying up to ~4 280 times higher RMSD values than OMEGA3D and DUACS along all major current systems ( fig.4h). It must be stressed that altimeter data are not assimilated in SODAv3.4.2, and modelled velocities are thus less constrained by observations. Directly comparing OMEGA3D and DUACS RMSD demonstrates that quasi-geostrophic velocities also improve with respect to geostrophic velocities (by a few cm s -1 ), mainly along the Antarctic Circumpolar Current and in the Western boundary currents ( fig.5a). 285  The second assessment focused on velocities provided by YOMAHA dataset, which are representative of currents at 1000 m depth. In that case, in order to increase the number of samples, a temporal window of ±1 day has been considered to build the matchup database ( Fig.3c and 3d).
A general overestimation of deep currents is revealed by looking at mean biases with respect to YOMAHA observations. The 300 mean bias attains around 5 cm s -1 in OMEGA3D, ARMOR3D and GLORYS12v1 ( fig.6a, 6c, 6e), while it gets up to >15 cm s -1 in SODAv3.4.2 ( fig.6g). GLORYS12v1 displays more spatially homogeneous values ( fig.6f), while purely observationbased values tend to overestimate mean western boundary currents but present slightly lower biases elsewhere (especially in the entire Atlantic Ocean). RMSD values computed from OMEGA3D, ARMOR3D and GLORYS12v1 at 1000 m depth ( fig.6b, 6d, 6f) show extremely similar patterns and values. As for surface values, much stronger discrepancies are found 305 between SODAv3.4.2 and YOMAHA estimates, reaching up to ~3 times higher RMSD values along the major current systems ( fig.6h). Specific comparison of OMEGA3D and ARMOR3D, RMSD shows also in this case that quasi-geostrophic velocities improve with respect to geostrophic velocities, even if by only <0.5 cm s -1 , along the Antarctic Circumpolar Current and in the Western boundary currents ( fig.5b).  OMEGA3D (a,b), ARMOR3D (c,d), GLORYS12v1 (e,f), SODAv3.4.2