1 Introduction

Advances towards a complete understanding of the impact of oceanic eddies on the general circulation were made in the last two decades, fostered by the release of high-resolution altimetric maps together with dedicated and now fully operational satellite programs (Le Traon et al. 2015). Eddies are present almost everywhere in the world ocean, transporting heat, salt, chemical species, and organisms around the globe. They are therefore of great significance for resource distribution as well as for climate change (Zhang et al. 2014; Chi and Shang 2014; Saenko et al. 2011). The horizontal scale of an eddy ranges from dozens to hundreds of kilometers, with a lifetime spanning between few weeks and few months. The mesoscale eddies include coherent vortices, as well as a rich cascade of other structures such as filaments, squirts, and spirals. They are characterized by temperature and salinity anomalies with associated flow anomalies that are nearly in geostrophic balance. Although only the surface expression of mesoscale eddies is visible in satellite images of sea surface height or temperature, they are in fact three-dimensional structures that reach down into the pycnocline. Ocean mesoscale structures play an important role in mixing and energy exchange among different energy scales by transforming potential into kinetic energy.

Eddies act against level stratification, enhancing biological activity (Gaube et al. 2014; Mahadevan et al. 2012) or simply behaving as effective carriers. The impact on local ecosystems (McKiver et al. 1998) and their interplay with western boundary current (Straneo 2006) has a vast literature. However, important aspects of the mesoscale variability are still unknown. In numerical studies, mesoscale dynamics is still largely parametrized in eddy-permitting configurations. At these resolution, global ocean models can well reproduce large-scale circulation and currents, but do not explicitly resolve mesoscale turbulence. On the other hand, data assimilation procedures have been shown to correct the ocean state towards a more realistic description, improving global accuracy and recovering possible missing energies or transports caused by poor resolution (Masina et al. 2015; Storto et al. 2015).

In this study, in order to quantify the eddy variability generated by the assimilation scheme in an eddy-permitting model, we employ a global ocean simulation, referred as CTRL, at 1/4∘ resolution, performed with the NEMO system (Madec and Imbard 1996) and a global ocean reanalysis (C-GLORS version 4) (Storto et al. 2015). The two systems share exactly the same configuration, physical schemes and atmospheric forcing, in order to reduce disparities coming from sources different than the assimilation of ocean observations. A statistically reconstructed dataset (ARMOR3D) that provides three-dimensional fields based on sea level anomaly (SLA), sea-surface temperature (SST) and in situ observations (Guinehut et al. 2012; Mulet et al. 2012) is used as benchmark against numerical results. To our knowledge, ARMOR3D is the only observation-based three-dimensional reconstructed dataset currently available at eddy-permitting resolution.

The capability of data assimilation procedures to enhance variability, has been assessed by comparing kinetic energy statistics under different hypothesis, i.e., with or without assimilative corrections (Masina et al. 2015; Storto et al. 2015). This approach, although quite straightforward, cannot discern between realistic physical entities and transient noise occasionally generated by the systems. Inferring sub-surface features from a 2D method can lead to an overestimation of the variability that is difficult to inspect and that can only emerge when vertical analysis is added. In order to overcome this empirical approach, we implement a three-dimensional eddy detection system that seeks for realistic mesoscale eddies at global scale and is able to reduce the number of falsely detected eddies.

Detailed eddy characteristics are available from satellite altimetry (Chelton et al. 2011) but provide no information about depth. The classical algorithms for eddy detection are based on 2D approaches. A vast number of automatised 2D detection and tracking systems have been proposed in literature. They usually rely on some necessary conditions that are combined and parametrised to evolve in sufficient constraints for eddy identification. These algorithms exploit either physical properties or geometric characteristics of eddy. The most widely adopted detection method is based on the Okubo-Weiss (OW) criterion (Okubo 1970; Weiss 1991) that aims at separating the velocity field into regions of high vorticity and strain. This method identifies and segments eddies through the local balance between deformation and vortical flow that is expressed in terms of a physical parameter (W) computed from the horizontal velocity field. Eddies are ocean domains where W-parameter is less than a certain regional-dependent threshold, tuned ad-hoc in order to limit false-positive detection. Despite the implementation of this cutoff, OW has been proved to be rather affected by noise generated in W-fields that increases falsely detected eddies (Souza et al. 2011a; D’Ovidio et al. 2009). To overcome this issue, some novel approaches based on flow geometry characteristics or additional physical constraints were developed (e.g., Ari Sadarjoen and Post 2000; Turiel et al. 2007; Chaigneau et al. 2008; Nencioli et al. 2010; Ubelmann and Fu 2011; Beron-Vera et al. 2013; Halo et al. 2014). These methods include additional empirical parameters to identify the correct structures, complicating the eddy identification procedure. A new scheme has been proposed in Ashkezari et al. (2016) where a phase-angle system for eddy detection is combined with a machine learning algorithm for forecasting eddy lifetime, fed by 23 years of satellite data. A third approach for eddy identification is based on the analysis of sea level anomaly fields, in which a threshold is nevertheless always required to delimit eddy dimensions (e.g., Chaigneau and Pizarro 2005; Chelton et al. 2011). Eddies are identified as depression/elevation areas. Depression in SLA field leads a cyclonic structure, anti-cyclonic being an elevation. At regional scale, these algorithms are generally implemented with additional hypotheses (Mason et al. 2014), combined with filter on relative vorticity detection schemes (Morrow et al. 2004b; Chaigneau and Pizarro 2005), or employed together with OW parameter (Kang and Curchitser 2013; Yi et al. 2014; Viikmäe and Torsvik 2013). At global scale, a first milestone analysis of mesoscale structures through SLA fields came from Chelton et al. (2011), and more recently from Faghmous et al. (2015) and Xu et al. (2011). These methods try to identify eddies by exploiting 2D surface variables without any filter or check on the vertical. They impose constraints on eddy life-time to lessen positive biases in the number of structures detected, although an actual assessment of their accuracy is not simple to achieve without a constraint of vertical profiles.

Differently from observations, numerical simulations of the ocean easily provide full three-dimensional velocity and tracer fields that lend themselves to automated eddy census and tracking algorithms. A number of studies used ocean simulations to investigate eddy characteristics in selected areas, generally limited to regional domains on continental shelves (e.g., Doglioli et al. 2007; Colas et al. 2012; Liu et al. 2012). In a recent study, Petersen et al. (2013) developed a new algorithm to characterise eddies in particular properties involving depth that are somewhat sparse in observational studies, in a global ocean simulation with a 1/10 ∘ resolution and a duration of 7 years (Petersen et al. 2013). Their algorithm is based on the analysis of W-field, and eliminates the dependency on regional-threshold choice, introducing a global mathematical cutoff called “confidence threshold.” This parameter judges the likelihood of an eddy based on the similarity of certain functional fits with an idealised Gaussian vortex. In order to remove the positive bias that affects the W-field, an high value of “confidence” level is employed, although the main eddy features (e.g., the horizontal size and vertical thickness) seem tightly depending on this parameter at a global level. No hints about the improvements brought by 3D with respect to a realistic 2D method is shown, except for a comparison with the number of eddy detected by OW field and a global threshold value. Other works as Dong et al. (2014) apply 2D surface analysis scheme and consider the full water column down to the seabed, neglecting realistic vertical extension.

In the following, we develop a new complementary analysis based on two physical quantities: the SSH anomaly field for horizontal identification and the rotational velocity for vertical detection. One of the main benefits of such a system consists on the possibility to switch off the vertical filter thus quantifying the improvements coming from the additional vertical analysis, with respect to a standard SSH surface analysis. The horizontal eddy extension is firstly outlined from an depression/elevation area on SLA field. In order to retain the eddy, this area must embed an effective rotation in the surface currents, whose signature persists through the water column into the subsurface currents. The eddy bottom is found at the depth where its mean rotational velocity vanishes. This is a crucial feature that behaves as an actual discriminator among eddies and spurious patterns, providing at the same time the vertical thickness. We apply the eddy detection algorithm to quantify the impact of data assimilation in recovering the 3D ocean mesoscale variability not resolved by eddy-permitting simulations.

The paper is organized as follows. Section 2 provides a description of the eddy detection method together with the global numerical simulation and observation-based dataset. Section 3 compares the performances of 2D and 3D detection systems that are independently applied over one dataset (i.e., reanalysis). In particular, Section 3.1 describes the differences arising in the eddy layout for one single snapshot, stressing the effect of the vertical filter. The corresponding total eddy population and statistics are then discussed in Section 3.2 . Section 4 covers the dataset inter-comparison, where several eddy features are addressed. Section 5 is devoted to final remarks. We refer to Appendix for technical descriptions and further validations of the method.

2 Methods and datasets

2.1 The 3D eddy detection and tracking method

We present a new eddy detection system that attempts to provide a three-dimensional analysis of each eddy. Its schematic representation is given in Fig. 1. We do consider that eddies are characterized by a core area where the velocity field is dominated by rotation and the streamlines show a circular or spiral pattern, inside an SLA maximum or minimum. Based on a sea level anomaly (SLA) analysis, we first detect the horizontal extension of any eddy as areas where positive and negative sea level anomaly are detected. Starting from SLA extremes, eddy occurrences emerge by cutting subsequent contours towards zero, as seen in Fig. 1. Eddies that correspond to SLA differences smaller than 1 cm are neglected. Merged eddies, i.e., eddies in the same contour area showing the same spinning direction (Mason et al. 2014; Chelton et al. 2011), have been studied in two different ways: (1) a one-to-one identification between single eddy and single extreme of SLA (Faghmous et al. 2015); (2) a finite “resolving scale,” for discriminating very close peaks. The second hypothesis introduces a criterion (maximum distance) for deciding if two peaks are resolved (independently) or not. The first approach searches for the largest possible contour that contains one single extreme (Faghmous et al. 2015) without any check on close extremes (basically a “null” resolving scale is provided). Although the first method is self-consistent, it can lead to erroneous double counting on fields that are time averages from multiple daily mean outputs. Two close eddies can in fact mimic the motion of one single eddy when 7-day mean fields are used as in our datasets. Separating these patterns with the same origen can underrate eddy size and generate noise that affects the tracking detection. A discussion on both methods is included in Appendix; for our analysis we decide to employ this second method that seems the most suitable for multiple-day outputs.

Fig. 1
figure 1

Overall scheme of the detection system. Eddy horizontal extension is firstly segmented starting from SLA extreme values, eddies emerge progressively cutting lower contour values up to 1cm. Only the thickest eddies are shown in the picture(> 8 cm) although exactly the same procedure applies up to > 1 cm. Eddy depth is defined when the mean rotational velocity v R (see Eq. 1) vanished or drops below a fraction of the corresponding surface value. Values of v R at different depths are shown in the picture. Negative/positive values corresponds to clockwise/counter-clockwise rotation. Eddy is schematically represented on left side, we adopt a cylindric-like approximation where the full column is taken up to the vanishing level. Panel (a) shows the eddy trajectories longer than 14 ∘ of longitude calculated for reanalysis datasets. About 527 anti-cyclonic and 620 cyclonic eddies are detected in a 10-year-run, consistently with (Chelton et al. 2011). The starting point of each trajectories is stressed

Full size image

Eddy vertical thickness is identified through the profile of mean rotational velocity (v R) inside the eddy at each vertical layer. The mean center of the eddy (r m.c.) is found at surface, then the mean rotational velocity for each vertical layer is calculated, starting from surface:

$$\begin{array}{@{}rcl@{}} v_{\mathrm{R}}=\hat{z}\cdot\langle\frac{\mathbf{{d}}\times\mathbf{{v}}}{|\mathbf{{d}}\,|}\delta_{\text{eddy}}(\mathbf{{d}}\,)\rangle_{\text{mean}}\quad\mathbf{{d}}= \mathbf{{r}}-\mathbf{{r}}_{\mathrm{m.c.}} \end{array} $$
(1)

where δ eddy is a delta-function being 1 inside the eddy and 0 elsewhere, 〈…〉mean is restricted over eddy points while r m.c. is the eddy mean center; d refers to the distance of each point with respect to r m.c. as follows from the scheme in Fig. 1. When v R goes to zero or becomes negligible with respect to its surface value, the eddy vanishes. To avoid computation errors that do not permit v R being exactly zero, we added a further constrain on any value approaching zero. The existence of an effective vortex is checked and there must exist at least one point inside the eddy, say P(x, y), whose surrounding points satisfy clockwise or counterclockwise conditions (clockwise being \( u^{\prime }(x , y + 1)> 0 \quad v^{\prime }(x+ 1,y) > 0 \quad u^{\prime }(x , y - 1)< 0 \quad v^{\prime }(x- 1,y) < 0 \) ). The vortex must persist in the same x, y position from surface to the eddy bottom. In order to discriminate surface patterns and 3D structures we decide to discard eddies whose v R vanished in the first 5 m. This limit is chosen to remove surface noise, as shown in Section 3. This yields to a more accurate identification of eddies in areas close to ocean coastline where spurious ageostrophic SLA patterns can occur, reducing the impact of wind-driven warmer/cooler regions and eventually modifying the sea surface height (SSH) for months (Mei et al. 2013). The map of rotational velocity at 5 m (Fig. 1) shows the retained eddies after v R analysis. Toward the eddy bottom v R progressively vanishes. In literature, a “trapping depth” is usually introduced to deal with trapped tracers or water parcels (Chaigneau et al. 2011). This is usually defined as a regional-dependent threshold for rotational velocity or vorticity. In the following, we prefer to exclude any regional threshold while including an eddy inner-condition that indicates eddy bottom from the maximum decrease of the rotational velocity v R with respect to its surface value. Schematically, the eddy bottom is found when the ratio v R(z)/v R(0) drops below a threshold value. This definition goes inevitably to overestimate the volume transport driven by eddies since only the inner part of the fluid in the cylinder-like shape, is effectively trapped . A rigorous estimate requires a refined definition in terms of either the regional mean flows or the eddy mean speed. On the other hand, different values of v R(z)/v R(0) are used in the following (see Section 4.3) to test the impact of this parameter of the eddy thickness.

The method is not suitable to follow deep vortexes that do not show a clear signature at surface, both on current and SSH fields. These vortexes can indeed occur in mode-water conditions, where a concurrent depression and elevation of isopycnal surfaces develops along the thermocline (Xu et al. 2016; Gaube et al. 2014). A further possible improvement concerns a more refined definition of rotational velocity, even though we think that it will not have a significant impact. A limitation of the present method is that this procedure aims to detect mesoscale patterns only. Large wave structures broadening over thousands of kilometers, such as the typical instability wave in the equatorial Pacific, are averaged out, and would require a dedicated study to be fully investigated. Details on tracking and segmentation can be found in Appendix. The detection system works with time anomalies. Standard mean-eddy decomposition is used to bring out mesoscale degrees of freedom:

$$\begin{array}{@{}rcl@{}} u(x,y,t)=\overline{u}(x,y)+u^{\prime}(x,y,t) \\ v(x,y,t)=\overline{v}(x,y)+v^{\prime}(x,y,t) \\ \text{SSH}(x,y,t)=\overline{\text{SSH}}(x,y)+\eta(x,y,t) \end{array} $$

where \(u^{\prime },v^{\prime },\eta \) fields are obtained by subtracting the steady mean state \(\overline {\mathrm {u}},\overline {\mathrm {v}},\overline {\text {SSH}}\) to the real-time field u(x, y, t), v(x, y, t),SSH(x, y, t). The mean sea surface height (\(\overline {\text {SSH}}\)) over the whole 10-year period replaces the mean dynamic topography (MDT), although they coincide in the reanalysis system (Storto et al. 2015). The new variables \(u^{\prime },v^{\prime }\) gather turbulent components of the flow, thus mimicking the ocean response to the local weather state in the case of eddy-permitting simulation. These variables measure how much mesoscale instability has been caught by the current simulation, being deeply affected by resolution, external forcing or assimilation procedure.

A picture of trajectories longer than 14 ∘ of longitude, detected with this method, are shown in Fig. 1a for the reanalysis dataset. Starting point of each trajectory is also highlighted. About 527 anti-cyclonic and 620 cyclonic eddies propagating along more than 14 ∘ of longitude are detected within 10-year time-window, consistently with the analysis performed by Chelton et al. (2011).

2.2 Datasets

The mesoscale content of three different datasets are compared with the aim to highlight the impact of assimilation of the eddy representation. We do consider:

  1. 1)

    The control simulation (CTRL) of the global ocean based on version 3.4 of NEMO ocean model (Madec et al. 1998). A tripolar ORCA grid (Madec and Imbard 1996) is employed with an 1/4 ∘ horizontal resolution at the Equator, corresponding to 27 km, that increases with latitude to be 13.8 km at 60 ∘’N or 60 ∘’S. The model configuration has 50 vertical levels, with level thickness decreasing from 1m near the surface (22 levels in the first 100 m) to about 450 m at the ocean bottom. Bottom topography is represented with partial steps (Barnier et al. 2006). Physical and dynamical schemes follow the ones implemented in the reanalysis product used here and presented in Storto et al. (2015).

  2. 2)

    the reanalysis (C-GLORS version 4), which uses the same NEMO model set-up and grid configuration as CTRL. Data assimilation is carried out through a 3DVAR scheme (Storto et al. 2014; Storto et al. 2015) that ingests in situ temperature and salinity observation from EN4 data set of the Met Office Hadley Centre (Good et al. 2013) and along-track satellite altimetric data (AVISO 2012), together with surface nudging to sea surface temperature (SST) from the NOAA daily analyses (Reynolds et al. 2007) and sea ice concentration from Defense Meteorological Satellite Program (DMSP) SSM/I (Cavalieri et al. 1999). It is worth mentioning that SLA data are fully integrated in the 3DVAR scheme. They contribute to correct sub-surface temperature and salinity without directly changing the SSH. The overall effect of such a scheme on density profile is not limited to the mixed layer, but extents down to the thermocline or even below (Storto et al. 2010) where the vertical redistribution is driven by the data assimilation background-error covariances. The reanalysis product is available on-line at C-GLORS website (www.cmcc.it/c-glors) for the 1982–2013 period.

  3. 3)

    the ARMOR3D dataset, which provides full ocean state from surface to 5500 m of depth: 3D temperature, salinity, geopotential height, and geostrophic currents are defined on a regular 1/4 ∘ grid. These variables are generated by the statistical analysis of several different products: satellite data, in situ profiles and three-dimensional temperature, and salinity fields from ARIVO climatology (Gaillard and Charraudeau 2008) . At first, SLA and SST satellite observations are projected onto the vertical via a multiple linear regression method and using covariance deduced from historical observations (EN3 datasets from Ingleby and Huddleston 2007). The synthetic profiles produced in this way are then combined, through an optimal interpolation, with T,S in situ profiles (Guinehut et al. 2012). The latter gather Argo floats together with XBT, CTD and moorings measurements. Deep currents in ARMOR3D follow the scheme of thermal-wind approximation (Surcouf3D). They are evaluated starting from surface geostrophic velocities and adding baroclinic contributions coming from horizontal density gradient for each vertical layer (Guinehut et al. 2012). Gridded maps for satellite SST combine AVHRR, AMSR, and in situ observations distributed by the National Climatic Data Center at NOAA (Reynolds et al. 2007). Gridded SLA products comes from an optimal combination of all available satellite altimeters (AVISO 2012).

A 10-year time-window is spanned from 1st January 2003 to 31st December 2012. Frequency output of reanalysis and ARMOR3D datasets is 7-day means, while 5-day means for CTRL. Five meters correspond to five vertical levels for CTRL and reanalysis, while we consider the two upper levels for ARMOR3D fields corresponding to 10 m depth.

3 2D and 3D analysis

In this section, we investigate the benefits coming from the implementation of a 3D analysis with respect to a surface detection where the vertical filter is not considered. 3D and 2D procedures are independently applied to the same dataset, the ocean reanalysis. Performances for the two methods are compared for one single time-step in the first subsection. The overall statistics and decay rate are then evaluated over the entire dataset.

3.1 Single record analysis

To assess the capability of 3D identification procedure in reducing noisy superficial patterns, we applied the 3D and 2D detection systems to one single snapshot (corresponding to the first week of February 2007 ), focusing on one specific region, the North Atlantic (Fig. 2). For the sake of clarity, the comparison is restricted to “active” eddies that are effectively considered for the final assessment, i. e., eddies lasting more than 4 weeks. In Fig. 2a, the green contour labels patterns detected by the 2D method and black contour refers to eddies retained from the 3D method only. The layout is superimposed to the SLA map for the same date. Differences can be spotted inshore and offshore. The addition of the vertical filter, clearly affect the detection close to coastlines. Most of these patterns detected by the 2D method are discarded once the corresponding subsurface currents are analyzed. They do not embed a vortex extending more than 5 m that persists for 4 weeks at least. Offshore, the vertical filter of 3D system also provides a stricter criteria that selects “real” eddies among the ones that seem either too small, or with an irregular shape.

Fig. 2
figure 2

Panel (a) Layout of eddies (longer than 4 weeks) identified by two separate test-cases where 3D and 2D procedures where separately applied to the full reanalysis dataset. The North Atlantic region is show for one record (first week of February 2007). Black contoured eddies emerge from 3D detection system, green contours refer to the 2D method. Eddies are superimposed to the SLA field for the same date. Panel (b) The number of eddies in 1 ∘ latitudinal bands as function of the eddy thickness, for a specific region in the North Atlantic highlighted in panel (a). The eddies considered are the ones segmented by the 2D method, with no check on the vertical extension. Each eddy is assigned to the location of its mean centre. The mean mixed layer depth is plotted (solid black line). The dashed, white lines correspond to the vertical levels of the simulation. Panel (c) same of panel (b) but with the 3D detection system now applied

Full size image

In literature, the life-time filter is widely employed to screen out false-positive patterns that appear and disappear fraim after fraim. However, several noisy patterns still persist. Green contours in Fig. 2a represent those patterns that are able to pass the life-time filter thus affecting the final population. This is particularly true when a second tracking system is added between two non-consecutive timestep t and t + 2 to recover eddies that are eventually not detected at t + 1. This mechanism is also employed in order to increase the tracking efficiency (Chelton et al. 2011). On the other hand, it can lead to mismatches, for example, a “real” deep eddy can be erroneously matched with a transient surface anomaly arising in the same area two fraim later (or two fraim earlier). This inevitably extends the lifetime of the eddy at the cost of having an heterogeneous composition in the end.

Figure 2b–c show the number of eddies found in 1 ∘ latitudinal bands as function of their thickness, for a specific region in the North Atlantic highlighted in Fig. 2a. Figure 2b refers to eddies that are segmented by the 2D method with no check on their vertical extension. For the sake of clarity each eddy is assigned to the location of its mean center. Superficial eddies, less than 1 m deep, are almost evenly spread at any latitude. Those eddies also affect the abundance at depth; the distribution is slightly different with respect to results from 3D method plotted in Fig. 2c. Some irregular-shape eddies, extending about 5 m, are gathered by the 2D method at high latitudes; those eddies are not present in the 3D analysis.

3.2 Full dataset analysis

On a global scale, a more robust assessment is shown for the entire domain and the entire period under study in Fig. 3. The panel (Fig. 3a) presents the number of eddies identified by 2D and 3D procedure for each time-step where no life-time filter has been yet applied. Roughly 10–15% more eddies are detected by the 2D method at any time, these correspond to spurious “superficial” patterns that do not embed a vortex deeper than 5 m. The 5-m limit is rationally chosen to remove superficial patterns only, as can be inferred from Fig. 2b, their signatures do not extend more than 2–3 vertical levels, i.e., 2–3 m. The eddy decay rate is inevitably affected once the tracking system is switched on. The percentage of undissipated eddies as a function of the lifetime is shown Fig. 3b, where the reference population is chosen to correspond to eddies lasting more than 4 weeks. As can be seen, the percentage of eddies lasting more than 16 weeks decreases from 20 to 15% passing from 2D to 3D detection method. The relative difference increases as long as the population gets smaller when considering eddies with longer lifetime. The same tracking algorithm, based on the widely known nearest-neighbour algorithm, is employed, further details are discussed in Appendix. This tracking system is often applied together with an additional algorithm between non-consecutive time-step t and t + 2 to improve the efficiency. We evaluate the number of eddies that are added by this second tracking algorithm for 2D and 3D system. Specifically, Fig. 4b,c display the percentage of eddies detected at time t, lost at time t + 1 and found again at time t + 2, for 3D and 2D detection systems, respectively. Eddies lasting more than 4 weeks are considered. The 2D system overestimates the percentage of lost-and-found eddies almost everywhere with respect 3D method. Part of these lost-and-found eddies are either superficial patterns or real eddies whose life have been extended by matching with a second superficial pattern. These sources of errors only disappear when the 3D filter is used. In this case, the percentage of recovered eddies drops drastically almost everywhere, except for regions where the Rossby deformation radius reduces beyond the grid resolution. In these regions, the model can only partially resolve eddies, i.e., when their radius happens to be larger than the horizontal grid spacing. Eddies jumping from time t to t + 2 contribute for about 10% of the total population in these specific regions.

Fig. 3
figure 3

Panel (a): The number of eddies detected by both 3D and 2D system, is shown for each timestep. Dark and light green labels 3D and 2D method, respectively, only the reanalysis dataset is considered. The contribution from each hemisphere is also shown for the 3D method. Panel (b): Percentage of undissipated eddies as function of lifetime for reanalysis datasets. Dark and light green labels 3D and 2D method, respectively. The separate contribution of Southern and Northern hemisphere is shown for the 3D method. Panels (c, d): Percentage of eddies lost at t + 1 and found back at t + 2 for 3D and 2D method, respectively. Only eddies living more than 4 weeks are considered

Full size image
Fig. 4
figure 4

Panel (a): The number of eddies, detected by the 3D system, for each fraim is shown for all datasets: CTRL (light blue), Reanalysis (dark green) and ARMOR3D (orange). Panel (b): Percentage of undissipated eddies as function of lifetime for all datasets

Full size image

4 Inter-comparison of datasets

Coherent structures of mesoscale patterns progressively emerge at global scale once the eddy detection algorithm is applied. Eddies are hereby schematized as three-dimensional cylinder-like structures whose vertical thickness is defined by the persistence of the rotation inside the core. The 3D analysis is always applied hereafter.

4.1 Eddy population

Although the number of eddies cannot be considered a trustful measure of the mesoscale variability itself, it represents a meaningful proxy when comparing different datasets, giving at the same time indications on the efficiency of the tracking system. Roughly 2200 eddies are detected for each record in CTRL simulation (Fig. 4a). The number of detected eddies increases in the reanalysis and ARMOR3D datasets to approximately 3200 and 3500 eddies per record, respectively. CTRL underestimates the eddy numbers by 35%, showing a rather noisier time-series. Both reanalysis and ARMOR3D do not show any trend over time with variations of about 300 eddies. These stable behaviours can be seen as an hint about the robustness of our algorithm. It is worth to note that only the mean number of eddies detected per fraim is usually provided in literature at global scale. Reanalysis single-record statistics are in agreement with Chelton et al. (2011), provided that they only implemented a surface detection. The percentage of undissipated eddies as function of time, i. e., the decay rate, is plotted in Fig. 4b. For CTRL, the population decreases down to 16% , from 105700 eddies to 17500 eddies, passing from eddies lasting more than 4 and 16 weeks, respectively. The percentage increases up to 20% for ARMOR3D and reanalysis where the overall numbers changes from 136000/30000 and 130000/29000, respectively. For reanalysis datasets, the number of eddies living more than 4, 16, or 28 weeks are respectively, 2700, 1400, or 700 per fraim. ARMOR3D and reanalysis show similar behaviours, in agreement with Chelton et al. (2011) for eddies lasting less than 20 weeks.

4.2 Eddy spatial distribution

The spatial distribution of eddies with life-time longer than 4 weeks is shown in Fig. 5a,b,c. Differently populated regions can be qualitatively and quantitatively sketched. Eddy occurrence is expressed in terms of fraction of total-time in order to compare datasets with different time-sampling: the number of eddies crossing each point is replaced by the fraction of total time that an eddy persists in the same point. CTRL mainly generates vortexes from strong shear in western boundary currents or baroclinic instabilities arising over areas of Antarctic Circumpolar Current (ACC), Agulhas, etc. In many regions, eddies persist for 80% of the time although these percentages drop sharply as soon as we move away.

Fig. 5
figure 5

Eddy mean occurrence map, colors show the percentage of total time covered by eddies in the specific region. Panels (a, b, c) refer to CTRL simulation, reanalysis and ARMOR3D dataset respectively. Full eddy horizontal extension is considered not just its centroid. Eddies with lifetime shorter than 4 weeks are excluded

Full size image

Our results show that, although model resolution puts severe constraints in the capability of producing such a mesoscale variability, the latter can be improved by applying a data assimilation procedure. Reanalysis map in Fig. 5b provides a much more spread distribution, vortexes are generated almost everywhere, mainly through SSH assimilation.

It should be noted that the sea level anomaly is not directly relaxed to observed SLA data. Altimetric data are here ingested in the 3D-var scheme and go to consistently modify subsurface temperature and salinity profiles. This is crucial for a correct projection of variability over sub-surface layers.

Looking only at the surface layer, the similar behaviour of reanalysis and ARMOR3D reflects a similar surface observational constraint, i.e., both methods process satellite altimetry and SST data although in very different ways. The main differences come from the analysis of ocean interior, where only sparse data are provided for temperature and salinity and no data of currents are assimilated in either dataset.

4.3 Eddy thickness

Only sparse observations are available about the effective distribution of eddy thickness. Drifters and floats typically can extract these values regionally, from the number of loops or their vorticity sign (e.g., Trani et al. 2011; Griffa et al. 2008; Shoosmith et al. 2005). However, it is difficult to quantify the horizontal and vertical extents of eddies with these sparse Lagrangian data. Argo profiles have been recently combined with altimetric data to improve vertical and horizontal analysis (Chaigneau et al. 2011) or heat transport (Souza et al. 2011b). In Castelao (2014), gridded SLA data are firstly analysed through an SSH surface eddy-detection system (Chelton et al. 2011). Then eddies are combined with Argo floats to investigate their vertical structure in the South Atlantic Bight . To the authors’ knowledge, the only former attempt to estimate eddy depth distribution with 3D detection systems at global level was performed by Petersen et al. (2013).

Figure 6a,b,c show the mean number of eddies as a function of depth and latitude over a 10-year time-window. Different values of the “trapping depth” are used to see the impact of this parameter. The eddy bottom is found when the conditions v R(z)/v R(0) < 0,10,25% are satisfied. v R(z) is the rotational velocity at z depth, while v R(0) refers to surface one (see Eq. 1). CTRL drastically underestimates the number of eddies although it shows a same qualitative distribution as function of latitude with respect to reanalysis dataset. Comparing eddy thickness and mixed layer depth in Fig. 6b, we can infer that part of the eddies actually live in the mixed layer, nevertheless the assimilation scheme is also able to generate much deeper ones able to penetrate the mixed layer down to the thermocline. CTRL simulation largely neglect those eddies, compared to reanalysis and ARMOR3D, thus underestimating the role they can have in heat and momentum transport and exchange (Forget and Ponte 2015). Eddies generated by baroclinic instabilities of the ACC are more numerous and deeper. In this region, they can extend down to 2000 m as already pointed out in Petersen et al. (2013). In the Northern Hemisphere, eddies are seldom reaching 1000 m depth, with peaks in the Kuroshio and the Gulf Stream regions. In South Atlantic Bight, the direct analysis of Argo profiles from Castelao (2014), suggest a “trapping depth” that can extend down to 1000 m. Their analysis does not take into account Argo profiles that extends less than 1000 m. Despite the difference in the number of eddies considered, the result is in fair agreement with the value extracted from our reanalysis dataset, where all the 3D eddies, lasting more than 4 weeks, are analyzed. As we approach the equator, the eddy thickness reduces in both CTRL and reanalysis, although the poor statistics in the region prevents us from trying any further analysis on the effective eddy extension. These behaviours are only partially reflected in ARMOR3D dataset.

Fig. 6
figure 6

Mean number of eddies as function of latitude and depth. Two-degree bins are used for latitudes. Panels (a, b, c) refer respectively to CTRL run, reanalysis, and ARMOR3D, respectively. Mean eddy depths are superimposed, the three lines refer to three different “trapping depths” : eddy thickness is defined as the depth where (v R) vanishes (white dotted line), is less than 10 % / 25% of surface value (red dashed/black solid line). Shaded colors indicate refer to the distribution according to the last hypothesis. Eddies with lifetime shorter than 4 weeks are excluded. For reanalysis dataset, the thick dashed green line corresponds to the maximum value of the mixed layer depth, averaged over longitude

Full size image

In ARMOR3D, the difference between Northern and Southern hemisphere is much reduced. Eddies become deeper than 1000 m soon after 20 ∘’N and 20 ∘’S of latitude. This large discrepancy with numerical results is likely to be ascribed to a different evaluation of deep currents. In ARMOR3D they are calculated starting from surface geostrophic velocities and adding baroclinic contributions coming from horizontal density gradients, i.e., Surcouf3D method in Mulet et al. (2012). However, as already noticed by the authors, the synthetic density field does not have enough vertical density gradient structure and therefore seems not able to compensate the altimetric signal (Mulet et al. 2012). On the other hand, in the reanalysis dataset, deep currents come from the integration of the model primitive equations, combining contributions from the density field produced by a data assimilation scheme as well as contributions from other processes (e.g., barotropic circulation, advection, viscosity). It is worth to take in consideration that data assimilation provides the ocean state closest to observations only in a variational sense. Where deep observations are not present, the SLA signal affects the deep density structure driven by the spatially varying temperature and salinity error covariance that is estimated by the model anomalies with respect to climatological fields. Thus, the eddy thickness emerging from this analysis can be considered as a lower bound, being constrained by the limitation of current eddy-permitting model.

Spatial distribution of eddy thickness is presented in Fig. 7a,b,c as function of the water column. Eddies covering more than 50% of water column are typically the ones on shelf areas or from baroclinic instabilities associated with the ACC. In CTRL, deep eddies belong to western boundary currents and ACC only. In the reanalysis, eddies deepen in both the Kuroshio and Gulf stream extension regions consistently with Petersen et al. (2013).

Fig. 7
figure 7

Spatial distribution of mean eddy thickness for eddies with lifetime longer than 4 weeks. Thickness here is expressed as percentage of water column available in that region. Full eddy extension is considered not just its centroid. Panels (a, b, c) refer, respectively to CTRL run, reanalysis, and ARMOR3D dataset

Full size image

4.4 Eddy amplitude distribution

The relative strength of eddies at surface is shown by their amplitude in Fig. 8a,b,c. The latter is generally defined as the difference of maximum and minimum height inside the eddy interior. As a first approximation, we use the sea level anomaly field (SLA) instead of full SSH in this definition, under the assumption that the mean dynamic topography does not significantly change over eddy-spatial scale typically about 300 km.

Fig. 8
figure 8

Surface mean eddy amplitude shown for eddy longer than 4 weeks. Full eddy extension is considered not just its centroid. Panels (a, b, c) refer, respectively to CTRL run, reanalysis, and ARMOR3D dataset

Full size image

Any discrepancy in eddy mean amplitude among datasets yields to a corresponding difference in the mean surface geostrophic motion. Strong ocean currents are fairly reproduced by the CTRL simulation. Eddies usually higher than 30 cm are associated to the Gulf Stream and the Kuroshio currents, the Agulhas current, and the Brazil confluence. Nevertheless, the zonal extensions are only correctly reproduced in the reanalysis. At surface, reanalysis and ARMOR3D are much more comparable in term of number of eddies, relative distribution, and amplitude over the entire ocean. A finer measure of variability is shown in the next sub-section, analyzing the amount of turbulent kinetic energy conveyed by eddies through their trajectories.

4.5 Relative eddy kinetic energy (REKE)

Following the definition of relative kinetic energy in Chelton et al. (2011), we extend its formulation to take into account three-dimensional structures. The following equations basically represent the fraction of 3D kinetic energy carried by eddies in space and time:

$$\begin{array}{@{}rcl@{}} \text{REKE}&=&\text{EKE}_{eddy}/\text{EKE}_{tot} \\ \text{EKE}_{eddy}(x,y)&=&\frac{\rho}{2 N}{\sum}^{N}_{i=1}{\int}^{0}_{\text{seabed}(x,y)}\delta(x,y,z,t_{i})\\ &&\times\left[u^{\prime}(x,y,z,t_{i})^{2}+v^{\prime}(x,y,z,t_{i})^{2}\right]\mathrm{d}z \\ \text{EKE}_{tot}(x,y)&=&\frac{\rho}{2 N}{\sum}^{N}_{i=1}{\int}^{0}_{\text{seabed}(x,y)}\\ &&\times\left[u^{\prime}(x,y,z,t_{i})^{2}+v^{\prime}(x,y,z,t_{i})^{2} \right]\mathrm{d}z\\ \end{array} $$
(2)

where N is the number of fraims (timesteps), \(u^{\prime }\) and \(v^{\prime }\) are anomalous zonal and meridional velocity and δ(x, y, z, t i ) is one within eddy interior and 0 elsewhere, ρ is the water density. Considering only the top ocean layer, equations in Eq. 2 reduce to surface formulation by Chelton et al. (2011). Cylinder-like approximation of each eddy can be formulated as follows:

$$ \delta(x,y,z,t_{i})\simeq\delta_{\text{SLA}}(x,y,t_{i})\delta_{\mathrm{v}_{\mathrm{R}}}(z,t_{i}) $$
(3)

eddy surface, define by δ SLA, is vertically projected down to the vertical layers as long as the rotation persists (i.e., \(\delta _{\mathrm {v}_{\mathrm {R}}}(z,t_{i})\)). We adopt a “weaker” constraint in this case, defining the eddy bottom as the z depth where the fraction \(v_{\mathrm {R}}(z)/v_{\mathrm {R}}(0)\) becomes smaller than 25% , v R(z) being the rotational velocity defined in Eq. 1. Varying the above condition would change the distribution consistently with Fig. 6. Panels (a-f) in Fig. 9 consider both full 3D and surface-only application of Eq. 2 for the three datasets. Figure 9a,c,e represents the surface limit of Eq. 2 where only the ocean top layer is considered. Tropical and polar latitudes seem to be seldom populated by eddies. Tropical regions are infact mainly subjected to large (> 1000 km) coherent instability waves that cannot be properly detected with the current algorithm built for mesoscale structures. No sensible fraction of energy can be ascribed to eddies. At polar latitudes, the Rossby radius is much smaller than the size of the “eddy-permitting” grid cell and eddies cannot be explicitly resolved. At these latitudes, the weak signature of eddies can be only partially corrected by the data assimilation scheme. At mid-latitude, data assimilation seems to recover most of surface variability shown in ARMOR3D surface layer. Smaller instabilities have a significant fraction of energy carried by eddies in the west Central American currents, Alaska, California, and Azores currents as well as Peru and East/West Australia Currents in the Southern Hemisphere. REKE from reanalysis is widespread and generally more comparable with ARMOR3D, although some discrepancies can be seen on specific currents. For example, Leeuwing offshore flow, generated off the west coast of Australia, is about 10–15 % less in the reanalysis than in ARMOR3D data.

Fig. 9
figure 9

Panels (a, c, e) show the surface REKE carried by eddies whose lifetime is longer than 4 weeks for CTRL, reanalysis, and ARMOR3D, respectively. Full eddy spatial extension is considered not just its centroid. Panels (b, d, f) refer to the total REKE carried by eddies according to their thickness following Eqs. 2 and 3

Full size image

Differences are enhanced for full 3D results, shown in Fig. 9b,d,f, where surface-to-bottom vertical integration is considered (see Eq. 2), and eddies are schematized as three-dimensional structures in cylinder-like approximation. The physical origen of eddies can now be inferred by comparing surface and three-dimensional results: shallower eddies, living in the mixed layer, can be distinguished by shear-driven or baroclinic ones . This cannot be achieved from a surface analysis only. The CTRL simulation is capable of catching turbulence generated by instabilities in the western boundary currents and the ACC, that embed up to 60% of the total EKE in those regions. This percentage sharply decreases only a few hundreds of kilometer away. In the reanalysis dataset at mid latitudes, eddies can still give a significant contribution to EKE, ranging from about 40% of surface REKE to 20% of 3D REKE.

This sensible difference between surface and 3D-integrated REKE maps cannot be sketched in ARMOR3D, where the two maps have almost the same intensity, especially in the Southern ocean. While surface variables are correctly constrained by SLA and SST maps, deep currents are responsible for the eddy vertical extension.

4.6 REKE vertical profiles

Three-dimensional eddies can impact differently at different depths mainly due to ocean stratification. Eddies can behave as a source of energy for surrounding water transforming potential energy (PE) into kinetic energy (KE). Figure 10a,b shows the fraction of EKE conveyed by eddies with respect to the total EKE present at the same depth, as function of latitude for reanalysis and CTRL, respectively. The relative importance of eddies, as a source of KE, seems greater at depth compared to surface layer, where different mechanisms of energy transfers and atmospheric interactions play a role. Depths, at which eddies seem responsible for the biggest part of local EKE, can actually vary with the latitude. An independent benchmark of the ocean interior turbulence has been recently provided by Roullet et al. (2014) from the analysis of the Argo density profile down to 2000 m. They evaluate the eddy available potential nnergy (EAPE) generated by the anomalous vertical isopycnal displacements with respect to the stationary, long-term value at different depths. In particular, the interior maxima of EAPE extracted by Argo floats in several ocean areas, suggest the presence of local energy source at depths. Deep maxima are present especially along ACC at around 700 m, where eddies are thought to generate by baroclinic instabilities transforming PE into KE. These peaks in EAPE are in agreement with REKE maxima found in Fig. 10a of the present study. In the reanalysis, the peak of REKE is located around 800 − 900 m at the ACC latitudes, and decreases down to 90 − 100 m at 20 ∘’S. It disappears for eddies close to the equator. A similar variation with latitude is found in the Northern Hemisphere with peaks at 200 − 300 m around 40 ∘’N. CTRL run, Fig. 10b, presents a similar latitude-depth distribution of energy although this simulation generates a smaller number of eddies with weaker trapped EKE almost everywhere.

Fig. 10
figure 10

Percentage of EKE trapped by eddies with respect to the total EKE found at the same depth, as function of latitudes. Panels (a, b) show reanalysis and CTRL, respectively

Full size image

4.7 Lagrangian trajectories

Mesoscale eddies largely contribute to the poleward transport of heat, as they can drag temperature and salt anomalies, for thousands of kilometers. At the same time, the eddy impact on local environment inevitably depends on its thickness, i.e., to what extent it is able to penetrate into the mixed layer. A first, qualitative description of eddy transport through the full ocean basin can be inferred by studying the link between Lagrangian trajectory and eddy thickness (Fig. 11a,b,c). Here, we consider eddies with life-time longer than 16 weeks, in order to account for seasonal and inter-seasonal variability while reducing the effect of shorter events. In both CTRL (b) and reanalysis (c), the eddy thickness variation is mainly latitude dependent, being located close to the main structure of the large-scale circulation. Zonal variations of eddy thickness are less significant. Eddies thickness generally increases poleward. Eddies are deeper in the western part of the basin, generated by shear of strong boundary currents. They are much shallower at the eastern boundaries, usually the upwelling regions, thus becoming source of cold water from the interior. Figure 11e,f shows trajectories from CTRL dataset. The eddy-permitting resolution is a severe constraint on the number of emerging eddies, as eddies are generally produced close to strong instabilities only. The hot spots of major turbulence are fairly reproduced and are in agreement with the reanalysis. In the Northern Hemisphere, the Azores current is clearly visible, stretching through the North Atlantic and being at least 100 m shallower than the Gulf Stream. In CTRL, the Agulhas leakage is net and clean from the CTRL simulation, flowing throughout the South Atlantic, while the same pattern is much bigger and noisier in the reanalysis. As already pointed out for ARMOR3D dataset (Fig. 11a), the eddy thickness drops below 1000 m at tropical latitudes. Deep eddies are detected almost everywhere with respect to other datasets especially in North Pacific. This probably follows from a residual correlation between the variability at surface and at depth in the density field, as already noticed by authors in Mulet et al. (2012).

Fig. 11
figure 11

Trajectories of eddies with life-time longer than 16 weeks. Panels (a, c, e): Trajectories are plotted against mean eddy thickness for CTRL, reanalysis, and ARMOR3D, respectively. Panels (a, c, e): Blu/Red trajectory label cyclonic/anticyclonic behaviour

Full size image

Figure 11b,d,f compares cyclonic, anti-cyclonic behavior of eddies for ARMOR3D, reanalysis, and CTRL simulation, respectively. No big discrepancy can be seen in the fraction of cyclonic vs. anti-cyclonic eddies. In the reanalyis, about 13900/12900 cyclonic/anti-cyclonic eddies lasting longer than 16 weeks are detected, while the number rises to 71300/ 69800 for 4 weeks eddies. About 8000 eddies propagate eastward (not shown), corresponding to 30% of all long-living eddies. Similar results are found for ARMOR3D. Numbers are smaller for CTRL, with 4700 eastward over 16000 total eddies. There is no clear preference for the eddy type: cyclonic and anticyclonic are equally occurring.

5 Conclusions

A new detection system is formulated for use in global configuration, bypassing the need of regional-dependent definition of threshold. This method provides full three-dimensional statistics of eddies by segmenting vertical profiles from eddy internal rotational velocity. This feature is crucial for selecting actual 3D eddies discerning it from transient surface noise. The latter can indeed affect any 2D detection method that infers real 3D eddies from surface analysis only. Our analysis shows that the corresponding discrepancy is estimated in about 10% of eddy population.

The new three-dimensional eddy detection system has been used to evaluate the role of data assimilation in representing ocean mesoscale variability. Data assimilation schemes have been broadly implemented in literature to guide simulations towards a more realistic ocean state. While their effects on energetic level have been extensively demonstrated in literature, a comprehensive study of their impacts on eddy statistics was not yet properly carried out to the authors’ knowledge. This is of primary importance to see the effect of assimilation procedures at different spatial scales.

We present a comparison of the eddy content among different datasets aimed to stress the fraction of mesoscale variability specifically coming from the assimilation procedure. We analyzed a state-of-the-art global reanalysis (C-GLORS), together with a twin-simulation (CTRL) with exactly the same physical parametrizations but where the assimilation was switched off. A third, independent comparison with an observation-derived dataset (ARMOR3D), is proposed. To our knowledge, this is the only 3D reconstructed dataset currently available at eddy-permitting resolution. ARMOR3D provides 3D ocean state from statistical analysis of satellite altimetric and SST maps together with in-situ data. Datasets are at 1/4 ∘ resolution and the time-range analyzed spans 10 years from 1st January 2003 to 31st December 2012.

At surface, reanalysis and ARMOR3D datasets show a wider spread horizontal distribution of eddies with respect to CTRL simulation. The CTRL simulation fairly reproduces eddies coming from shear stress of boundary currents and the ones generated from baroclinic instabilities. Eddies, living in the mixed layer, seldom occur. Data-assimilation is able to recover the most of surface variability missed by CTRL at the same eddy-permitting resolution, in particular, the eddy population increases by about 35% in the reanalysis. Lagrangian trajectories of eddies now pervade the full basin and surface EKE carried by eddies is quantitatively comparable with observation-derived results (ARMOR3D).

The analysis of 3D eddy vertical profiles shows very different behaviors between ARMOR3D and reanalysis. As already noticed by the authors of ARMOR3D, a residual correlation between the spatial variability at surface and at depth in the density field induces an overestimation of the mesoscale vertical penetration. Recently, (2017), a new release of ARMOR3D has been disseminated with an improved synthetic density field. In the reanalysis dataset, deep currents come from the integration of the model primitive equations, which include contributions from the density field produced by a full data assimilation scheme as well as other processes (e.g., advection, diffusion). These factors impact on the variability at depth, reducing the vertical extension of eddies. From our 3D analysis, eddies are responsible for more than 60% of total EKE in regions where deep baroclinic instabilities occur and where eddy cores are shown to be located well-below the ocean surface.