A multi-year time series 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 their impact on marine ecosystems and Earth climate. Here we present the OMEGA3D product, an observation-based time series 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 a diabatic quasi-geostrophic (QG) diagnostic model to the data-driven CMEMS-ARMOR3D weekly reconstruction of temperature and salinity as well as ERA Interim fluxes. Outside the equatorial band, vertical velocities were retrieved in the upper 1500 m at 1∕4^∘ nominal 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 and 1000 m) decrease with respect to corresponding estimates obtained from zero-order geostrophic balance, meaning 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 reanalyses spanning a comparable time period but based on data assimilation in ocean general circulation numerical models.

The full OMEGA3D product (released on 31 March 2020) is available upon free registration at https://doi.org/10.25423/cmcc/multiobs_glo_phy_w_rep_015_007 (Buongiorno Nardelli, 2020a). The reduced subset used here for validation and review purposes is openly available at https://doi.org/10.5281/zenodo.3696885 (Buongiorno Nardelli, 2020b).

1 Introduction

The recognition of the key role played by the oceans in the Earth system led the United Nations to proclaim the Decade of Ocean Science for Sustainable Development (2021–2030). Major efforts will consequently be made in the next years to analyse state-of-the-art observations and models and provide the indispensable knowledge basis to preserve the marine environment through effective, science-informed policies. Providing accurate reconstructions of 3D ocean circulation time series is a fundamental part of this effort, aimed at better describing ocean dynamics and assessing their responses and feedbacks to natural and anthropogenic pressures. However, assessing the time-evolving lateral and vertical transport of energy, momentum, gases, nutrients, marine organisms and pollutants would require repeated synoptic observations of the 3D ocean state and surface forcings that cannot be presently achieved even with the most advanced technologies. Hence, a combination of measurements collected from in situ and remote-sensing platforms and proper modelling frameworks is needed to describe the ocean circulation both in the ocean interior and at the domain boundaries. Two main complementary approaches can be followed to this end: the assimilation of observations in global ocean circulation numerical models (Carrassi et al., 2018; Moore et al., 2019; Stammer et al., 2016) and the combination of diagnostic models and purely data-driven reconstructions. The latter is presently more widely used for surface circulation retrievals, but its extension to 3D ocean state reconstruction is generating growing interest as advanced statistical and machine learning tools are becoming more computationally efficient (Buongiorno Nardelli et al., 2012; Buongiorno Nardelli and Santoleri, 2005; Guinehut et al., 2012; Lopez-Radcenco et al., 2018; Mulet et al., 2012; Rio et al., 2016; Ubelmann et al., 2015; Yan et al., 2020). Both approaches present advantages and drawbacks: data assimilation in prognostic models may guarantee a description of the ocean state evolution that is fully consistent with the physics represented by the model, but the uncertainties in its initialization, the limits of the parameterizations of unresolved processes and the difficulties to properly represent model and observation errors and to account for their representativeness can significantly reduce models' ability to reproduce non-assimilated observations. Conversely, synergic use of satellite in situ observations and data-driven reconstruction methodologies, in combination with simpler dynamical models (often limited to zero-order balances, as in the retrieval of geostrophic currents from sea level data), can provide snapshots that better match independent observations (Mulet et al., 2012; Rio et al., 2016; Ubelmann et al., 2016).

The OMEGA3D product, developed by the Consiglio Nazionale delle Ricerche within the European Copernicus Marine Environment Monitoring Service (CMEMS; http://marine.copernicus.eu/services-portfolio/access-to-products/, last access: 23 July 2020; product ID: MULTIOBS_GLO_PHY_W_REP_015_007), delivers the first observation-based retrievals of the global 3D vertical and horizontal ocean currents computed with a quasi-geostrophic (QG) model that explicitly considers the effect of both geostrophic advection and upper-layer turbulent mixing. The data are provided weekly at 0.25^∘ latitude–longitude resolution, with 75 non-uniformly spaced vertical levels between the surface and 1500 m depth, covering the period from January 1993 to December 2018 (with planned yearly extensions based on updated upstream datasets). Outside the equatorial band, vertical velocities are obtained by solving a diabatic Q vector formulation of the QG omega equation (Buongiorno Nardelli et al., 2018b; Giordani et al., 2006), with vertical mixing parameterized through the K-profile parameterization (KPP; Smyth et al., 2002). Only once vertical velocities are known can the horizontal ageostrophic components be retrieved.

The accuracy of the QG velocities depends on the input data and on the theoretical limits of the model and parameterization used. Omega forcings are estimated here from the multi-year CMEMS product ARMOR3D (Guinehut et al., 2012), providing a statistical reconstruction of 3D temperature and salinity fields from a combination of in situ profiles and satellite observations of sea surface temperature, salinity and topography. ERA Interim air–sea fluxes are used to evaluate the forcing terms due to vertical mixing (Dee et al., 2011). QG approximation implies that the omega equation cannot be solved at the Equator, and increased errors are expected in the low-latitude bands.

A direct validation of the vertical velocities is not possible due to the lack of direct reference observations. As such, OMEGA3D vertical velocity mean pattern and variability have been compared here with two global model reanalyses that include vertical velocity fields as disseminated output, namely Estimating the Circulation and Climate of the Ocean (ECCO; Forget et al., 2015) and Simple Ocean Data Assimilation (SODA; Carton et al., 2018).

Total horizontal and geostrophic components are instead compared with fully independent velocity estimates obtained from drifting-buoy and Argo float displacement. For reference, a similar comparison is carried out between two reanalyses, SODA and CMEMS GLORYS (Global Ocean Reanalysis and Simulation; Drévillon et al., 2018), and drifter data.

2 Methods

2.1 Input datasets

Two datasets are taken as input for OMEGA3D processing:

1. The global ARMOR3D reprocessed dataset (Guinehut et al., 2012). This is distributed by CMEMS within the product MULTIOBS_GLO_PHY_REP_015_002 (http://marine.copernicus.eu/services-portfolio/access-to-products/, last access: 23 July 2020, dataset ID: dataset-armor-3d-rep-weekly), which is one of the two data-driven ocean state reconstructions included in the Ocean Reanalyses Intercomparison Project (ORA-IP; Balmaseda et al., 2015). ARMOR3D is built in successive steps that include the retrieval of temperature (T) and salinity (S) profiles from gap-free surface temperature (Reynolds et al., 2007), surface salinity (Droghei et al., 2018) and sea level anomaly (AVISO+, 2015) fields, carried out through a multilinear regression of historical profiles (Cabanes et al., 2013), as well as the successive combination of 3D synthetic fields with in situ T and S profiles through an optimal interpolation algorithm. ARMOR3D provides weekly fields at 0.25^∘ nominal latitude–longitude resolution over 33 regularly spaced vertical levels with different spacing depending on depth between the surface and 1500 m.

2. The ERA Interim (Dee et al., 2011) surface air–sea fluxes. These are included in the global atmospheric reanalysis by the European Centre for Medium-Range Weather Forecasts (ECMWF; https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/, last access: 23 July 2020).

ERA Interim assimilates several observations of upper-air atmospheric variables (e.g. satellite radiances, temperature, wind vectors, specific humidity and ozone) through a four-dimensional variational (4D-VAR) system, running with a 12-hourly analysis cycle. OMEGA3D diabatic forcings take in input of the mean daily fields of the zonal and meridional components of the turbulent surface stress, the surface latent and heat flux, the surface net solar and thermal radiation, and total precipitation and evaporation (needed to estimate the equivalent surface salinity flux in KPP).

2.2 Input data preprocessing

The numerical tool used to retrieve the OMEGA3D product is designed to run on a non-uniform vertical grid that displays a refined mesh close to the surface. The vertical layer thickness increases with the square of depth, and the final grid includes 75 vertical levels between 2.5 and 1482.5 m. This grid was specifically designed to obtain more accurate numerical solutions within the ocean's upper boundary layer (Kalnay de Rivas, 1972; Sundqvist and Veronis, 1970). Preprocessing 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 set to cubic spline interpolation; Virtanen et al., 2019) and the mapping of ERA Interim data on OMEGA3D horizontal grid (using Python class scipy.interpolate.griddata set to fit data to a piecewise cubic, continuously differentiable, curvature-minimizing polynomial surface; Virtanen et al., 2019).

As ARMOR3D data may occasionally display density inversions along the water column that are not compatible with the QG omega solution, vertical profiles of potential density are adjusted to impose static stability: moving from the surface to depth, density is set to the upper level value plus a 0.0001 kg m⁻³ increment whenever a density inversion is found.

2.3 Quasi-geostrophic equations

A diabatic Q vector formulation of the quasi-geostrophic omega equation (Buongiorno Nardelli et al., 2018b; Giordani et al., 2006) is solved to get the OMEGA3D vertical velocity fields:

$$\begin{array}{}\text{(1)}& {\mathrm{\nabla }}_h^{\mathrm{2}}\left(N^{\mathrm{2}}w\right)+f^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}w}{\partial z^{\mathrm{2}}}}={\mathrm{\nabla }}_h\cdot \mathit{Q}.\end{array}$$

In Eq. (1), w represents the vertical velocity (positive, upwards), N² is the Brunt–Väisälä frequency, f is the Coriolis parameter, 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 on the right:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathit{Q}& =\mathrm{2}{\mathit{Q}}_{\mathrm{twg}}+{\mathit{Q}}_{\mathrm{th}}+{\mathit{Q}}_{\mathrm{dm}}\\ {\mathit{Q}}_{\mathrm{twg}}& ={\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\displaystyle \frac{\partial u_{\text{g}}}{\partial x}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial x}}+{\displaystyle \frac{\partial v_{\text{g}}}{\partial x}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial y}},{\displaystyle \frac{\partial u_{\text{g}}}{\partial y}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial x}}+{\displaystyle \frac{\partial v_{\text{g}}}{\partial y}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial y}}\right)\\ {\mathit{Q}}_{\mathrm{dm}}& ={\displaystyle \frac{f}{{\mathit{\rho }}_{\mathrm{0}}}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\mathit{\rho }\left[K_{\text{m}}\left({\displaystyle \frac{\partial v}{\partial z}}-{\mathit{\gamma }}_v\right)\right]\right.,\\ & \left.-{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\left[\mathit{\rho }{K}_{\text{m}}\left({\displaystyle \frac{\partial u}{\partial z}}-{\mathit{\gamma }}_u\right)\right]\right)\\ {\mathit{Q}}_{\mathrm{th}}& =-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathrm{\nabla }}_h\left({\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathit{\rho }}\left({\displaystyle \frac{\partial \mathit{\rho }}{\partial z}}-{\mathit{\gamma }}_{\mathit{\rho }}\right)\right]\right)\\ & ={\mathrm{\nabla }}_h\left({\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathit{\rho }}\left(N^{\mathrm{2}}+{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{\gamma }}_{\mathit{\rho }}\right)\right]\right).\end{array}\end{array}$$

In the above definitions, ρ indicates the potential density, g is the gravitational acceleration, and u_g and v_g represent the geostrophic velocities, while turbulent terms are defined following the KPPs, namely through non-local effective gradients γ_x and classical viscosity and diffusivity K_x. 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 cell mixing by including an amplification of turbulent velocity scales and includes a non-local momentum flux term and a parameterization of Stokes drift effects. Forcing terms are computed from ARMOR3D potential density and geostrophic velocity fields as well as ERA Interim atmospheric reanalyses.

In one of the analytical steps to obtain Eq. (1), the details of which are given elsewhere (Buongiorno Nardelli et al., 2018b; Giordani et al., 2006), the following two equations are found:

$$\begin{array}{}\text{(3a)}& {\displaystyle }\begin{array}{rl}{f}^{\mathrm{2}}{\displaystyle \frac{\partial u_{\text{a}}}{\partial z}}& =-{\displaystyle \frac{\partial N^{\mathrm{2}}w}{\partial x}}-\mathrm{2}f\left({\displaystyle \frac{\partial u_{\text{g}}}{\partial z}}{\displaystyle \frac{\partial v_{\text{g}}}{\partial x}}+{\displaystyle \frac{\partial v_{\text{g}}}{\partial z}}{\displaystyle \frac{\partial v_{\text{g}}}{\partial y}}\right)\\ & +{\displaystyle \frac{f}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\left(\mathit{\rho }{K}_{\text{m}}\left({\displaystyle \frac{\partial u}{\partial z}}-{\mathit{\gamma }}_u\right)\right)\\ & -{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathit{\rho }}\left(N^{\mathrm{2}}+{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{\gamma }}_{\mathit{\rho }}\right)\right]\right)\end{array}\text{(3b)}& {\displaystyle }\begin{array}{rl}{f}^{\mathrm{2}}{\displaystyle \frac{\partial v_{\text{a}}}{\partial z}}& ={\displaystyle \frac{\partial N^{\mathrm{2}}w}{\partial y}}+\mathrm{2}f\left({\displaystyle \frac{\partial u_{\text{g}}}{\partial z}}{\displaystyle \frac{\partial u_{\text{g}}}{\partial x}}+{\displaystyle \frac{\partial v_{\text{g}}}{\partial z}}{\displaystyle \frac{\partial u_{\text{g}}}{\partial y}}\right)\\ & -{\displaystyle \frac{f}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\left(\mathit{\rho }{K}_{\text{m}}\left({\displaystyle \frac{\partial v}{\partial z}}-{\mathit{\gamma }}_u\right)\right)\\ & -{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathit{\rho }}\left(N^{\mathrm{2}}+{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{\gamma }}_{\mathit{\rho }}\right)\right]\right).\end{array}\end{array}$$

Once vertical velocities are retrieved through the omega solution, these two equations allow the estimation of the horizontal ageostrophic components (u_a and v_a) through a simple vertical integration of the right-hand-side terms as well (assuming horizontal ageostrophic velocities are negligible at the bottom boundary).

2.4 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 rewritten substituting derivatives with central finite differences and considering a non-staggered grid. Vertical derivatives are computed considering a variable grid spacing, increasing with the square of depth (Kalnay de Rivas, 1972), and adopting a finite difference scheme of second-order accuracy (Sundqvist and Veronis, 1970).

At the surface and topographical boundaries Dirichlet conditions are imposed (namely vertical velocities of zero), and Neumann conditions are imposed at the bottom and lateral boundaries (namely partial derivatives of vertical velocity are set to zero). The latter are imposed through forward–backward finite schemes and make the solution unsuitable for modelling current topography interactions along the coasts. Grouping all factors and 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 for the entire domain at once, however, would be extremely computationally 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 into smaller horizontally overlapping subdomains (each tile having a horizontal dimension of 75 grid points, overlapping by one-third). The inversion is carried out sequentially on these subdomains, 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 (LGMRES; Baker et al., 2005) with incomplete lower–upper (LU) preconditioning (as implemented in the Python sparse linear algebra package scipy.sparse.linalg, imposing a tolerance for convergence of 10⁻⁷; Virtanen et al., 2019).

Vertical velocities are finally used to integrate Eq. (3) by a simple trapezoidal rule to obtain ageostrophic horizontal velocities.

3 Validation

3.1 Model reanalyses used for inter-comparison

Three different ocean state reconstruction time series 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 terms of numerical schemes used, input data ingested and assimilation strategies.

The first dataset considered is the third release of version 4 of ECCO (Forget et al., 2015; Fukumori et al., 2018), hereafter ECCOv4r3, covering the 1992–2015 period and available at https://ecco.jpl.nasa.gov/products/all/ (last access: 23 July 2020). ECCOv4r3 is based on the MIT General Circulation Model (Adcroft et al., 2004) and applies a 4D-VAR assimilation scheme to a wide set of observations (including satellite altimetry, in situ T and S profiles, satellite sea surface salinity and temperature, and ocean bottom pressure; Fukumori et al., 2017), minimizing the observation analysis misfits in a least squares sense (Wunsch and Heimbach, 2013). The 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-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/, last access: 23 July 2020).

The second dataset is version 3.4.2 of SODA (Carton et al., 2018), hereafter SODAv3.4.2, which covers the 1991–2017 period and can be downloaded from https://www.atmos.umd.edu/~ocean/ (last access: 23 July 2020). SODAv3.4.2 reanalysis is based on the ocean component of the NOAA–Geophysical Fluid Dynamics Laboratory CM2.5 coupled model (Delworth et al., 2012), namely version 5.1 of the Modular Ocean Model, and was designed with a 0.25^∘ horizontal-resolution and 50-level vertical-resolution grid; thus it is improved with respect to previous versions of the same system (Carton et al., 2018) and provides the same nominal resolution of OMEGA3D. SODAv3.4.2 assimilates basic hydrographic data from the World Ocean Database (Boyer et al., 2013) and level 3 (night-time) sea surface temperature from different sources (Carton et al., 2018) through a linear deterministic sequential filter with a 10 d cycle. SODAv3.4.2 reanalysis output is provided on its native grid at 5 d resolution.

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 sea surface temperature data, sea ice concentration, and in situ temperature and salinity vertical profiles into a global ocean eddy-resolving model with 50 vertical levels. The GLORYS12v1 model component is Nucleus for European Modelling of the Ocean (NEMO; Madec and the NEMO team, 2016), and data assimilation is carried out by means of a reduced-order Kalman filter, while a three-dimensional variational (3D-VAR) scheme provides a correction 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/ (last access: 23 July 2020) (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.

As for OMEGA3D, ECCOv4r3, SODAv3.4.2 and GLORYS12v1, surface forcings are all taken from ERA Interim (Dee et al., 2011).

3.2 In situ validation data

Two fully independent in situ datasets have been considered for the validation of OMEGA3D horizontal velocities: Surface Velocity Program (SVP) data (Lumpkin et al., 2013) from the NOAA Global Drifter Program (covering the period 1993–2018 and freely available at https://www.aoml.noaa.gov/phod/gdp/, last access: 23 July 2020) and the YoMaHa'07 (Yoshinari Maximenko Hacker, hereafter YOMAHA) database (covering the period 1997–2018 and freely available at http://apdrc.soest.hawaii.edu/projects/yomaha/, last access: 23 July 2020; Lebedev et al., 2007). Both datasets provide velocity estimates obtained from the displacements of drifting platforms along a Lagrangian trajectory.

In order to minimize wind slippage, SVP drifters are drogued with a 7 m long holey sock centred at 15 m depth, and their velocity estimates are considered representative of currents at 15 m depth (Lumpkin et al., 2017). Before carrying out the validation, individual 6-hourly SVP drifters were 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 predefined parking pressure and emerge only for near-real-time data transmission through ARGOS–IRIDIUM satellites. Most of these instruments follow a profiling cycle of approximately 10 d, and their parking level is set to 1000 m.

3.3 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⁻¹, depending on depth and processes involved). Consequently, OMEGA3D vertical velocity cannot be directly validated (no reference datasets exist). However, the algorithm used to retrieve OMEGA3D horizontal velocities requires the vertical velocity as 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 time series of comparable length. Considering that vertical velocity fields are provided at different space–time resolutions, this comparison only describes the mean patterns and the amount of variability captured by each product.

Figure 1 Mean vertical velocity at 100 m (a, c, e) and related standard deviation (b, d, f) computed from OMEGA3D (a, b), SODAv3.4.2 (c, d) and ECCOv4r3 (e, f) data on the OMEGA3D domain over their 23-year overlapping period. Areas shallower than 100 m or contaminated at least once by sea ice within the averaging period are masked in light grey.

Specifically, mean vertical velocity patterns at 100 m depth and associated variability (standard deviation) are computed here from OMEGA3D, SODAv3.4.2 and ECCOv4r3 over their 23-year overlapping period (1993–2015), focusing on the domain covered by the OMEGA3D product and thus excluding the 5^∘ N–5^∘ S band and coastal areas. The large-scale patterns and range of values found in the averaged velocities are quite similar among the three reconstructions (Fig. 1). Maximum absolute mean values reach around 2 m d⁻¹, and areas dominated by large-scale, wind-driven upwelling at high latitudes and by downwelling at mid-latitudes are consistently identified in the three products, with values rarely exceeding 0.5 m d⁻¹. In the intertropical band, OMEGA3D vertical velocities display slightly higher values than the other two reconstructions (especially in the central Pacific). More substantial differences are found along the major western boundary current systems and along the Antarctic Circumpolar Current, where the different nominal resolution and the different dynamical representation of mesoscale features in the three systems are reflected in terms of averaged vertical transport. Specifically, even if OMEGA3D and SODAv3.4.2 display very similar patterns, the former displays stronger values in the Agulhas Return Current and along the Gulf Stream meanders as well as around the northern branch of the anticyclonic gyre around the Zapiola Rise, while the latter presents intensified exchanges in the Pacific–Antarctic Ridge and South Indian Ridge areas. The differences between OMEGA3D and SODAv3.4.2 in the Zapiola Anticyclone, an eddy-driven flow controlled by bottom friction, could be associated with the difficulties in accurately representing that circulation in many global reanalyses. As expected, ECCOv4r3 patterns do not resolve any of the alternated upwelling and downwelling patterns found along the main currents' meanders in the two 0.25^∘ resolution products, and its representation of mean vertical velocities in all western boundary currents basically consists of uniform upwelling and downwelling associated with the parameterization of baroclinic instabilities along steep isopycnal slopes. OMEGA3D mean vertical velocity patterns look very similar to SODAv3.4.2 in the northern part of the Atlantic Ocean as well, while ECCOv4r3 presents quite different large-scale patterns.

Given its 5 d sampling, SODAv3.4.2 could be expected to reveal a stronger variability than OMEGA3D (7 d sampling), and both are expected to display much higher values than ECCOv4r3 (providing monthly averaged fields). Conversely, though associated patterns display very similar features, the maximum OMEGA3D standard deviation value exceeds SODAv3.4.2 by a factor of ∼ 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 after low-pass filtering the latter two time series (by a five-point and seven-point moving window, respectively) to only keep frequencies lower than monthly, like those provided by ECCOv4r3 (Fig. 2 should thus be compared to Fig. 1f). Even in that case, the variability observed in ECCOv4r3 is sensibly lower than those retrieved from higher-spatial-resolution products, likely revealing the limits of the mesoscale parameterization used in ECCOv4r3 in terms of vertical exchanges.

Figure 2 Mean monthly vertical velocity patterns and standard deviations computed from OMEGA3D (a) and SODAv3.4.2 (b) after low-pass filtering the time series to remove signals above monthly frequency.

3.4 Horizontal velocity validation vs. independent observations

OMEGA3D horizontal velocity accuracy has been assessed in terms of mean bias and root mean square differences (RMSDs) with respect to space–time-co-located in situ reference observations. Estimated metrics have then been compared to those estimated for geostrophic velocities directly obtained from the Data Unification and Altimeter Combination System (DUACS) altimeter data (when looking at SVP velocities; AVISO+, 2015) or from ARMOR3D (when looking at YOMAHA velocities; Mulet et al., 2012) and also successively 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.

The first assessment covered surface currents as measured by SVP drifters. As SVP drifters may occasionally lose their drogue, thus failing to correctly represent 15 m depth currents, only drogued SVP drifter data collected within ±12 h from nominal reconstruction dates have been included in our matchup databases (Fig. 3a and b). The same matchup procedure has been applied to DUACS geostrophic velocities, SODAv3.4.2 and GLORYS12v1 (all of which share the same nominal horizontal resolution).

Figure 3 Number of matchups within $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ bins between OMEGA3D (a) as well as GLORYS12v1 and SODAv3.4.2 (b) velocities (at 15 m depth) and SVP observations and between OMEGA3D (c) as well as GLORYS12v1 and SODAv3.4.2 (d) velocities (at 1000 m) depth and YOMAHA velocities.

Figure 4 Mean and root mean square differences between OMEGA3D (a, b), DUACS (c, d), GLORYS12v1 (e, f) and SODAv3.4.2 (g, h) velocities at 15 m depth and co-located SVP observations.

Table 1 Main characteristics of CMEMS OMEGA3D product.

Mean biases between SVP and OMEGA3D velocities (Fig. 4a) display similar values and patterns to what obtained from altimeter-derived geostrophic velocities (Fig. 4c), with a slight underestimation of the current intensities. OMEGA3D actually appears more biased than DUACS geostrophic velocities close to the tropical band, likely due to the fact that the omega equation is derived from the f plane, and the forcing cannot be correctly estimated there by definition (as accurate horizontal velocities are needed to compute all its terms; see also Sect. 2.3). Overall, their mean biases do not exceed 10 cm s⁻¹. GLORYS12v1 slightly overestimates mean western boundary currents, though it displays very small biases elsewhere (Fig. 4e). Conversely, mean differences between SODAv3.4.2 at 15 m depth and SVP velocities (Fig. 4g) reveal a more significant underestimation of surface currents, with biases reaching up to 20 cm s⁻¹over wide portions of the domain. Similarly, OMEGA3D and DUACS RMSD values (Fig. 6b and d) present very minor differences, while GLORYS12v1 presents significantly higher differences (by a factor of ∼ 2) along all major currents (Fig. 4f). Even stronger discrepancies affect SODAv3.4.2 estimates, displaying RMSD values up to ∼ 4 times higher 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 centimetres per second), mainly along the Antarctic Circumpolar Current and in the western boundary currents (Fig. 5a).

The second assessment focused on velocities provided by the 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 d).

Figure 5 RMSD of OMEGA3D quasi-geostrophic and geostrophic horizontal velocities vs. drifters in $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ bins at 15 m (a) and 1000 m (b) depth. Negative values indicate an improvement with respect to geostrophy.

A general overestimation of deep currents is revealed by looking at mean biases with respect to YOMAHA observations. The mean bias attains around 5 cm s⁻¹ in OMEGA3D, ARMOR3D and GLORYS12v1 (Fig. 6a, c, e), while it reaches up to > 15 cm s⁻¹ in SODAv3.4.2 (Fig. 6g). GLORYS12v1 displays more spatially homogeneous values (Fig. 6f), while purely observation-based 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, d, f) show extremely similar patterns and values. As for surface values, much stronger discrepancies are found between SODAv3.4.2 and YOMAHA estimates, reaching RMSD values up to ∼ 3 times higher along the major current systems (Fig. 6h). Specific comparison of OMEGA3D and ARMOR3D RMSD values also shows that quasi-geostrophic velocities improve with respect to geostrophic velocities in this case, even if by only < 0.5 cm s⁻¹, along the Antarctic Circumpolar Current and in the western boundary currents (Fig. 5b).

Figure 6 Mean and root mean square differences between OMEGA3D (a, b), ARMOR3D (c, d), GLORYS12v1 (e, f) and SODAv3.4.2 (g, h) velocities at 1000 m depth and co-located YOMAHA observations.

4 Data availability

The OMEGA3D product is distributed as part of the CMEMS catalogue (https://resources.marine.copernicus.eu/?option=com_csw&task=results?option=com_csw&view=details&product_id=MULTIOBS_GLO_PHY_W_3D_REP_015_007, last access: 23 July 2020, https://doi.org/10.25423/cmcc/multiobs_glo_phy_w_rep_015_007; Buongiorno Nardelli, 2020a). The reduced subset used for validation and review purposes is openly available at https://doi.org/10.5281/zenodo.3696885 (Buongiorno Nardelli, 2020b). Access to the full product is granted after free registration as a user of CMEMS at https://resources.marine.copernicus.eu/?option=com_sla (last access: 23 July 2020). Once registered, users can download the product through a number of different tools and services, including the web portal Subsetter, DirectGetFile (DGF) and FTP. More information can be found at http://marine.copernicus.eu/services-portfolio/technical-faq/ (last access: 23 July 2020).

The basic characteristics of the OMEGA3D product are summarized in Table 1.

5 Conclusions

The 1993–2018 OMEGA3D time series provides weekly observation-based estimates of the 3D vertical and horizontal ocean currents in the upper 1500 m of the global oceans. The product is obtained by applying a quasi-geostrophic diagnostic model (based on the omega equation) that includes the effect of both geostrophic advection and upper-layer turbulent mixing and delivers, for the first time, estimates of the vertical velocities based on a combination of satellite and in situ observations. The OMEGA3D time series is provided over a 1∕4^∘ horizontal grid, and inter-comparison with model reanalyses of comparable length and different space–time resolutions indicates that OMEGA3D resolves the highest amount of vertical velocity variance. Both OMEGA3D and model reanalyses present slightly negative biases with respect to SVP and positive biases with respect to YOMAHA, likely reflecting data representativeness differences and missing processes (or only partially parameterized processes) in models (e.g. sub-mesoscale processes) as well as potentially inaccurate or biased fluxes in the input. OMEGA3D displays the smallest horizontal velocity root mean square differences with respect to independent estimates obtained from drifting-buoy or float displacements (when compared with reanalyses), mostly improving surface current estimates at the mid–high latitudes along the Antarctic Circumpolar Current and in the western boundary currents. Due to the approach followed to retrieve the horizontal ageostrophic components (once the omega equation is solved), the observed improvements (also with respect to simpler geostrophic estimates) mean that estimated vertical velocities should also be deemed reliable.

The model cannot be applied in the equatorial band (where QG approximation fails), and it does not include any parameterization of bottom-boundary-layer mixing. Dirichlet conditions are thus set at the bottom. As a consequence, considering that the domain is also limited to the upper 1500 m, OMEGA3D is not suited for studies of bottom-boundary dynamics or equatorial dynamics. Dirichlet conditions are also applied at coastal boundaries. Even if the effect of lateral-boundary conditions only propagates a few grid points due to the elliptical nature of the omega equation (Buongiorno Nardelli et al., 2001, 2012, 2018a), this also makes the OMEGA3D product unsuitable for studies of coastal dynamics (the product is considered reliable approximately 100 km away from masked coastal areas). As such, OMEGA3D is mostly suited to describe the role and long-term variability of open-ocean, large mesoscale dynamics and air–sea interactions (here parameterized through KPP), for example regarding the vertical exchanges and water mass transformation outside the equatorial band.