This paper proposes a new algorithm for parallel identification of mesoscale eddies from global satellite altimetry data. By simplifying the recognition process and the sea level anomaly (SLA) contours’ search range, the method improves identification efficiency compared with the previous SSH-based method even in the single-threaded process. The global SLA map is divided into several regions. These regions are identified simultaneously with a new SSH-based method. All the eddy identification results of these regions are merged seamlessly into a global eddy map. A β-plane approximation is used to calculate the geostrophic speed in the equatorial band. Compared with the computation complexity of the previous SSH-based method, which is , the computation complexity of the new method is , where K is the number of threads and L is the number of regional SLA maps. When applying the new method to the global SLA map, the computation is ~100 times faster than the previous SSH-based method on an average computer. The new method characterizes an eddy structure by radius, amplitude, eddy core, closed SLA contour, and closed SLA contour with maximum average geostrophic speed. In situ data and another global eddy dataset are applied to validate the reliability of eddies detected by the new algorithm. Global eddy mean properties, variability, and the geographical distribution of both datasets are analyzed to demonstrate the performance of this new method and to help understand eddy activities on a global scale.
Mesoscale eddies are circular currents of water that widely exist in the ocean, and they spanning tens to hundreds of kilometers over tens to hundreds of days. Eddies play a significant role in the transport of momentum, mass, heat, and nutrients, as well as salt and other seawater chemical elements, effectively impacting the ocean’s circulation, large-scale water distribution, and biology (Chelton et al. 2011; Chaigneau et al. 2011; Chen et al. 2011, 2012; Faghmous et al. 2012). There are thousands of eddies in the global ocean that can be generally classified as either cyclonic if they rotate counterclockwise (in the Northern Hemisphere) or anticyclonic otherwise (Faghmous et al. 2013).
There has been much research on the statistical analysis of regional eddy characteristics. Chaigneau et al. (2008) investigated the average characteristics of mesoscale eddies in the eastern South Pacific Ocean; Souza et al. (2011) made an analysis of eddy properties in the South Atlantic Ocean; Chen et al. (2011, 2012) used 17 years of satellite altimetry data to analyze eddy mean properties, spatiotemporal variability, and heat and salt transports in the South China Sea; and Zheng et al. (2015) inferred eddy characteristics in the south Indian Ocean from surface drifters. However, only a few works (Chelton et al. 2007, 2011; Faghmous et al. 2015) made long-term sequence estimates of global eddy characteristics.
Identification of mesoscale eddies on a global scale can ensure the continuity of tracking, and reveal the laws of ocean mass and energy exchange between oceanic basins. As demonstrated by Zhang et al. (2014), who used the Okubo–Weiss parameter method (Isern-Fontanet et al. 2003) to get the surface localization of eddies, eddy-induced zonal mass transport is comparable in magnitude to that of the large-scale wind- and thermohaline-driven circulation. Therefore, eddy research of long time series on a global scale has great significance on analyzing and understanding ocean energy and mass transport, as well as the spatiotemporal variability of global eddies.
A large number of automatic eddy identification algorithms have been developed to study eddy activity. These algorithms can be divided into three categories: 1) the physical-parameter-based method, including the Okubo–Weiss parameter method (Chelton et al. 2007; Henson and Thomas 2008), the winding angle method (Chaigneau et al. 2008; Sadarjoen and Post 2000), and the 2D wavelet method (Doglioli et al. 2007); 2) the flow-direction-based method (Nencioli et al. 2010; Williams et al. 2011); and 3) the SSH-based method (Chelton et al. 2011; Mason et al. 2014; Faghmous et al. 2015). Another modern method that is based on the instantaneous Lagrangian flow geometry (Haller 2005; Beron-Vera et al. 2013; Haller and Beron-Vera 2013, 2014; Wang et al. 2015; Onu et al. 2015; Beron-Vera et al. 2015; Haller 2015) is proposed to identify eddies in turbulent flows.
Despite the major circular structures that can be detected by existing eddy detection algorithms, there is still more work to be done. Some algorithms lack the computational efficiency of contour iterations (Wu 2014) or have complex calculation processes (Onu et al. 2015). In addition, because of low Coriolis effect strength when calculating the geostrophic speed of eddies, eddy detection becomes unreliable in the equatorial region (Chelton et al. 2007, 2011; Wu 2014).
Compared with the other two kinds of eddy identification methods, the SSH-based method presented the best performance. Chelton et al. (2011) has evidenced that the SSH-based method is superior to the Okubo–Weiss parameter method due to its ability to avoid extra noise and excess eddy detections. Moreover, the flow-direction-based algorithm (Nencioli et al. 2010; Williams et al. 2011) is primarily for simulated data with higher resolution. It cannot be applied to satellite altimeter data (e.g., AVISO 2014), because the algorithm requires higher-resolution data to get accurate flow rotation values for the classification of eddies (Wu 2014).
Previous SSH-based methods (Chelton et al. 2011; Mason et al. 2014) work well at basin scale but slow down significantly at the global scale, because of a high-order computation complexity on contour iterations. With the improvement of data resolution, existing hardware does not meet the requirement of global eddy recognition due to the increase of the number of SLA contours. Thus, a more efficient global eddy identification algorithm is required. In this paper, our method offers several improvements over existing eddy detection methods:
Both anticyclonic and cyclonic eddies are identified simultaneously when SLA contours are iterated from minimum SLA value of all cyclone seeds to maximum SLA value of all anticyclone seeds.
Space segmentation of global SLA map reduces the computation complexity of eddy identification by fewer SLA contour iterations. Regional eddy results are merged seamlessly into a global eddy map to avoid redundancy and loss of eddies.
A robust parallel computing method is adopted to improve eddy recognition efficiency on a global scale.
A β-plane approximation (Lagerloef et al. 1999) is applied to compute geostrophic velocity in the equatorial band.
This paper is organized as follows: In the next section, we describe the satellite dataset and introduce methods used in this paper. The global eddy dataset generated by our method is verified by in situ data and compared with the Faghmous et al. (2015, hereafter FA15) dataset in section 3. In section 4, we perform a global eddy statistical analysis of the two datasets on a global scale. In section 5, we conclude the paper with a discussion of the study’s contributions and future work.
2. Data and method
The altimeter data used in this study are Archiving, Validation, and Interpretation of Satellite Oceanographic Data (AVISO 2014), which were generated from a combination of data from TOPEX/Poseidon, Jason-1 and Jason-2, and Environmental Satellite (Envisat). AVISO delayed time altimeter data are global daily mean sea level anomalies at a spatial resolution of ¼° × ¼° and provide satellite images of SLA data spanning roughly 23 years, from 1993 to 2015.
a. Eddy identification
1) Complexity analysis of algorithm
The computation complexity of previous SSH-based methods (Mason et al. 2014; Chelton et al. 2011) is linear in the number of SLA contours. When generating SLA contours in the process of global eddy recognition, the SLA value is divided into a number (M) of levels for each level has several number (N) of SLA contours. Thus, we need to iterate the SLA contours times, which would result in a computation complexity of . To improve the iteration speed of the SLA contours for the eddy identification method, decreasing the SLA contour iterations is the most effective way. Space segmentation and the parallel computing method are used in this paper to achieve this purpose. First, the global SLA map is divided into the proper number (L) of regions. For each regional SLA map, there are levels of SLA contours and SLA contours at every level. Next, the regional SLA maps are processed simultaneously by parallel computing. Let K be the number of threads, , , and the SLA contour iterations of our method are calculated using
Thus, the computation complexity of our method is .
2) Eddy identification process
Figure 1 shows a flowchart detailing the parallel eddy identification algorithm. It first divides the global SLA map into several regions (1. Segmentation). Second, each regional SLA map is spatially high-pass filtered (2. High-pass filtered). Third, it scans for local extreme SLA pixels and designates these extreme pixels as eddy seeds (3. Eddy seed identification). Fourth, SLA contours are computed and searched to extract the eddy boundary (4. SLA contours computation and iteration). Finally, all the eddy identification results of these regional SLA maps are merged seamlessly into a global eddy map (5. Regional eddy results merged).
Previous SSH-based methods (Chelton et al. 2011; Mason et al. 2014) work well on basin scale, but they have low computational efficiency on global scale. So, a global SLA map is partitioned into appropriate regions on basin scale. Each region is divided into two blocks: one is the major part of the region, referred to as the “inner block” (boxes A and B in Fig. 2); the other block is the overlap of two adjacent regions, referred to as the “outer block” (box C in Fig. 2); this extended area is introduced here to avoid missing eddies in the boundary of the inner block due to segmentation. For example, in Fig. 2, eddy 2 and eddy 4 are two eddies in the same region (A + B + C), but they encapsulate to eddy 1 and eddy 3, respectively, because of segmentation.
(ii) High-pass filtered
To facilitate eddy detection, a spatially high-pass filter with half-power filter cutoffs of 20° longitude × 10° latitude is applied to remove steric heating and cooling effects, as well as other large-scale variability from the AVISO SLA map (Chelton et al. 2011).
(iii) Eddy seed identification
Eddy seed identification is an essential procedure that can effectively avoid multicore eddies. Potential cyclonic (or anticyclonic) eddy cores are detected by finding the local minima (or maxima) in the SLA map. An extremum is defined as a pixel in the SLA map whose value is higher or lower than all of its eight neighboring pixels in a 3 × 3 neighborhood. Figure 3 shows the result of eddy seed identification of the South China Sea on 8 August 2004.
(iv) SLA contour computation and iteration
Closed contours of SLA approximately correspond to streamlines surrounding an eddy (Chelton et al. 2011). SLA contours are computed at 0.25-cm intervals for levels of minimum SLA value of all cyclone seeds to maximum SLA value of all anticyclone seeds. It is different from previous SSH-based methods (Chelton et al. 2011; Mason et al. 2014), which use a range of SLA thresholds from −100 to +100 cm. Although 1-cm increments had been confirmed good enough to resolve the eddies with compact interiors (Chelton et al. 2011), it was found that the smaller increment is preferred, since it provides a better approximation of eddy dimensions (Yi et al. 2014; FA15). While the increment does affect the algorithm’s computational runtime, our algorithm allows for use of fine threshold steps with reasonable compute time. Furthermore, SLA contours are sequentially identified and analyzed to detect anticyclonic and cyclonic eddies simultaneously. Thus, our method improves identification efficiency by simplifying the recognition process and the SLA contours’ search range. Figure 4 shows eddy seeds and SLA contour plots of the South China Sea on 8 August 2004. To classify the closed contours into eddies and noneddies, closed contours of SLA must satisfy the following criteria:
Contain no more than one seed (local maxima/minima):
This is different from Chelton et al. (2011), who permit multiple local maxima/minima.
Eddy structures with a radius smaller than 0.4° are weakened, caused by the preprocessing of the AVISO Reference Series (AVISO 2014). Since eight pixels is about 0.4° radius on the AVISO SLA map, it is adopted to be the minimum size threshold of eddies. The eddy radius is far smaller than 4.5° (radius of ~1000 pixels), so we choose 1000 pixels as the maximum area for defining an eddy (Chelton et al. 2011).
Eddy shape test with an error of 55%:
Here, , where is the area of the closed SLA contour deviations from its fitted circle (circular shape) and is the area of that fitted circle.
Amplitude (A): 1 cm 150 cm:
Here, , where is the SLA value at the eddy seed, and is the average SLA value around the outermost closed SLA contour that defines the eddy perimeter.
When a closed SLA contour meets the abovementioned tests, it is identified as an eddy, and the seed inside it is identified as the eddy center. The outermost closed SLA contour is referred to as the effective perimeter of the eddy (). Then, a closed SLA contour with the maximum average geostrophic speed () is estimated by iterating from inward over all closed SLA contours whose pixel count ≥ 8 (Fig. 5).
The geostrophic speed (U) of the eddy corresponds to the geostrophic velocity amplitude (Chaigneau et al. 2008): . Here, geostrophic velocity components () can be computed from the SLA gradients:
where g is the acceleration due to gravity; and are the eastward and northward distances, respectively; is the Coriolis parameter; and is the rotation rate of the earth, with φ being the latitude.
However, it is important to note that f tends to 0 near the equator and that the geostrophic approximation will not be valid between ~5°N and ~5°S (Chaigneau et al. 2008). The identification of eddies therefore is compromised near the equator (Chelton et al. 2011; Wu 2014). Although FA15 makes when , it leads to inaccuracy in the geostrophic speed.
It has been shown that a β-plane geostrophic approximation [, where R is the earth’s radius] provides excellent agreement with observed velocities in the equatorial undercurrent (Lukas and Firing 1984; Picaut et al. 1989). Geostrophy varies smoothly from β-plane formulation at the equator to an f-plane formulation in the midlatitudes. The transition between the two approximations is done using a Gaussian weight function having a meridional decay scale that is found to be approximately the Rossby radius (~2.2° latitude) (Lagerloef et al. 1999). Thus, when we are searching for by iterating from inward, the geostrophic velocity components () of every closed SLA contour are calculated followed the stencil-width methodology (Arbic et al. 2012) using f-plane formulation, β-plane formulation, and ,
This does not affect significantly the statistical results of global eddies in this study. However, it is important to study the kinetic energy variation of mesoscale eddies in the equatorial region for further work.
(v) Regional eddy results merged
Not until all regions of a global SLA map finish eddy detection will we have judgment between adjacent regions: if any two eddies in the outer block have the same location of the eddy core, then they are supposed to be the same eddy, and the smaller one of the two will be removed (Fig. 2). After the above-mentioned processing, the eddy identification results of all regions are merged seamlessly to get a global eddy map (Fig. 6).
3) Eddy identification performance
According to Eq. (1), two things should be considered when the number (L) of regional SLA maps is set for partitioning the global SLA map: 1) it should not be too small. Because the size of the inner block and the outer block should be larger than the size of a normal eddy, to ensure that eddies were not shrunk or missed because of segmentation; 2) it should not be too big that the advantage of the parallel computing method will not be fully exploited and that the computation complexity of our method will not be fully reduced. In this study, the size of the outer block for each regional SLA map is set to 10°, which is roughly equal to the diameter of 1000 pixels. The size of the inner block is different, as the number of regions dividing the global SLA map changed, and it is not smaller than 10° to ensure the eddy will not be separated because of segmentation.
Table 1 shows the performance of global eddy detection by the previous SSH-based method and our method with different number (K) of threads while global SLA map is divided into different number (L) of regional SLA maps. The two algorithms have a similar amount of eddy numbers (~4800) but a different runtime. When , the previous SSH-based method spent ~2.5 times as much time as our method spent on detecting global eddies. There are two reasons that improved the identification efficiency of our method compared to the previous SSH method: 1) SLA contours are generated from the minimum SLA value of all cyclone seeds to the maximum SLA value of all anticyclone seeds, which caused fewer SLA contour iterations than the range from −100 to +100 cm; and 2) SLA contours are sequentially identified and analyzed to detect anticyclonic and cyclonic eddies simultaneously.
As L and K increase, the runtime of our method decreases. When and , our method spent only 8.4 min on detecting global eddies, which is about 145 times faster than the previous SSH-based method. Because of the cost of merging eddies of the outer block, practical increased efficiency compared with the previous SSH-based method is lower than the theoretical one [Eq. (1)]. All measurements were conducted with an Intel Core i7 with a 4790 (3.60 GHz) CPU and four-processor core, eight threads, running under Windows 7 64-bit Ultimate. The efficiency of our method for global eddy detection will be improved if the machine performance is better.
b. Eddy tracking
The eddy-tracking algorithm used in this study is adapted from Chelton et al. (2011). For each eddy identified at time step , the eddies identified at the next time step are searched to find the closest eddy lying within a restricted search region. In order for an eddy to be associated with another in the subsequent time frame, both the amplitude and area at must be within 0.25–2.75 times that of . Considering eddies may disappear between consecutive maps, we allow eddies to be unassociated for one day before terminating the track. Thus, searching distance D is set to 1.75 times the distance that a long baroclinic Rossby wave would propagate in one 2-day time step based on the long Rossby wave phase speed, calculated from the Rossby radius of deformation as presented in Chelton et al. (1998).
3. Verification with in situ data
As Argo floats with a short repeating cycle can track eddies for a long time, they have great prospect in monitoring eddies in the open ocean (Zhang et al. 2015). The satellite-tracked drifters have their drogues centered at 15-m depth to measure surface currents, and a semirigid material throughout the drogue can provide support for the drifter to maintain its shape in high-shear flows (Lumpkin and Pazos 2007). Thus, the trajectory data of Argo and drifters can be used to test the accuracy of eddy identification results. To demonstrate the performance of our method, verification figures prepared from FA15 are included for comparison.
a. Validated by Argo float trajectory data
Three Argo floats (2901380, 2901381, and 2901382, named A1–A3, respectively) observed the activities and evolution of an anticyclonic eddy (AE). The AE was generated on 7 October 2010 and then moved from 17.75°N, 120°E to 19.25°N, 114.25°E, and ended on 19 January 2011. We also get a similar eddy trajectory from the FA15 (Fig. 7) dataset, which has a period from 10 November 2010 to 22 January 2011. Float A2 was released on 23 October 2010, and floats A1 and A3 were launched on the following day in the eddy near Luzon Island. Float A1 was thrown out from the AE on 3 December. And floats A2 and A3 spun and moved with the AE until it disappeared. The trajectories of these three floats not only verified the existence of the AE but also successfully tracked its movement (Fig. 7), which matches well with the eddy-tracking result of our method.
The Physical Oceanography Laboratory of Ocean University of China launched 17 Argo floats (2901550–2901566) in March 2014 to observe an AE in the western Pacific area. These floats are different from other regular floats, having a parking depth at 500 m and a cycle of once a day. We select three Argo floats (2901551, 2901559, and 2901556, named A4–A6, respectively) to verify eddy identification and the tracking result of our method. Floats A4 and A5 were released on 28 March 2014, and float A6 was launched on the following day. Floats A4–A6 continuously tracked an AE in the western Pacific area from 28 March 2014 to 30 September 2014. Figure 8 shows three Argo float trajectories around an AE. The red line is shorter than the blue line, because the FA15 dataset was generated only until 1 May 2015. About a half-a-year eddy observation reflects the high accuracy of our algorithm.
b. Validated by surface drifter trajectory data
We also chose some satellite-tracked drifters in different oceans to validate the eddy detection results of both datasets. In the eastern North Pacific, a satellite-tracked drifter (43677, named D1) that was deployed on 21 January 2011 moved around a cyclonic eddy (CE) from 20 March 2011 to 31 August 2011 before leaving this area (Fig. 9a). In the Indian Ocean, a CE was generated on 3 September 2001, while a drifter (9525824, named D2) moved around it at the same time. D2 followed the CE for 108 days and moved northward after being thrown out of it (Fig. 9b). As shown in Fig. 9c, an AE was generated on 3 May 2005 and died on 7 October 2006, tracked by our method in the Atlantic Ocean. A drifter (36929, named D3) started moving around the AE on 30 May 2005 and followed it for 453 days until 25 August 2006.
4. Global ocean eddy statistics
In this part we compare the global eddy results from our method and FA15. The eddy identification conditions are different between the two methods. Area restriction of eddies in the FA15 dataset is ≥ 4 pixels, and there is no limit to the amplitude of eddies in FA15. In this case, we can provide an explanation for the differences between the two datasets on eddy-tracking results in section 3: having a stricter eddy identification threshold (size, amplitude), our method is more likely to terminate a track than FA15. To get more intuitive comparison results of the two datasets, we select eddies with sizes of at least eight pixels and 1-cm amplitude from the FA15 dataset.
a. Eddy mean properties
Histograms of global eddy mean properties (eddy radius, amplitude, and rotational speed) are shown separately for cyclones and anticyclones of the two datasets in Fig. 10. Cyclones and anticyclones of both datasets have similar presentation on global eddy mean properties.
As presented in Fig. 10a, both datasets show a radius with a peak at 50–60 km, and more than half of eddies have a radius between 50 and 100 km. However, eddy radii of the FA15 dataset are larger than our dataset, because the unconstrained approach to determining the eddy contour of FA15 method may sometimes lead to overestimating the eddy’s size (FA15).
Figure 10b also shows some differences of eddy amplitudes between the two datasets. Although the eddy amplitudes of both datasets peaked at 1–2 cm, they have different probability: about 32% in the FA15 dataset and 24% in our dataset. The factor that influences the amplitude bias may be the 0.05-cm increment in the FA15 method.
There are significant differences between the two datasets on the rotational speed U of eddies as presented in Fig. 10c. Globally, about 70% of the observed eddies in the FA15 dataset have 5 cm s−1, while it is ~10% in our dataset. Most observed eddies in our dataset have U between 5 and 20 cm s−1. The rotational speed of an eddy in our method is characterized by the maximum of the average geostrophic speed around all of the closed contours of SLA inside the eddy. However, it is defined as the average geostrophic speed of all pixels inside the eddy in FA15.
b. Variability of eddies
There is no significant seasonal variation in terms of global eddy number. About 2400 cyclones and 2400 anticyclones formed per day as detected by our method for each global SLA map in different seasons. This is close to the result of FA15, which identifies approximately 2300 cyclones and 2300 anticyclones for each daily SLA snapshot in different seasons.
There is good general agreement between the two datasets in interannual variations of mean number of global eddies generated per day as presented in Fig. 11. It shows some sharp decreases of mean number of eddies generated per day during the periods of December 1993–March 1994, 1997–98, 2002–03, 2009–10, and 2012–14. Within the period 1993–2015, a strong El Niño occurred in 1997/98, followed by two moderate events in 2002/03 and 2009/10. An interannual variation of eddies may be related to the El Niño event; however, further research is needed to verify a connection.
c. Geographical distribution
Figures 12a and 12b are similar in the geographical distribution of global eddies. There are ~260 000 eddy trajectories with lifetimes ≥ 4 weeks in our dataset from January 1999 to December 2015, and ~310 000 eddy trajectories with lifetimes ≥ 4 weeks in the FA15 dataset from January 1993 to May 2014. Figure 12 clearly shows several areas, including most of the western continental boundary, the Antarctic Circumpolar Current, the Agulhas Current, the Eastern Australian Current, the North Brazil Current, the Gulf Stream, and the Kuroshio Extension, where eddies are passed through more often than any other area of the ocean.
This paper proposes a new algorithm for parallel identification of mesoscale eddies from global satellite altimetry data. By using space segmentation and the parallel computing method, our algorithm is ~100 times faster than the previous SSH-based method on global eddy detection. Moreover, our algorithm improves recognition efficiency compared with the previous SSH-based method even in the single-threaded process. The β-plane geostrophic approximation is introduced to compute geostrophic velocity components accurately in the equatorial band, which is ignored or miscalculated in other methods. Besides, eddies detected by our algorithm and FA15 are verified by in situ data. Global eddy statistics of both methods are compared in terms of eddy mean properties, variability, and geographical distribution. Discrepancies between the two methods relating to eddy radius, amplitude, and rotational speed are identified and discussed. There are broadly consistent results on variability and geographical distribution of eddies on the global scale.
With the improvement of the observation technique and computational hardware, higher-resolution observation data or model data are generated, and our method is capable of supporting oceanic eddy identification of such data. In addition, our method not only improves the efficiency of eddy identification but also can be applied to accelerate other object recognition of ocean phenomenon globally.
The eddy identification algorithm applied in this paper reveals only a two-dimensional structure of eddies on the sea surface; the lack of knowledge about the vertical structure of eddies prevents people from figuring out the estimation of their role in the oceanic mass and heat transport (Souza et al. 2011). Global in situ observations, such as the Argo profiling floats and satellite-tracked drifters, can be used to study the vertical structure of eddies. The combination of satellite altimeter data and vertical profiles should provide new estimations for oceanic eddies in future work. Moreover, considering the potential benefits of the Lagrangian flow geometry method of eddy identification on transport assessments, we will combine the method with the parallel computing method to estimate global eddy transport more accurately and efficiently in future work.
This research was jointly supported by the Natural Science Foundation of China under Grants 41331172, U1406404, and 61361136001, and the Fundamental Research Funds for the Central Universities under Grant 201562024. The Argo data were obtained freely from the Coriolis Argo Global Data Assembly Centre (ftp://ftp.ifremer.fr/ifremer/argo/). The satellite-tracked drifter data were provided by the GDP Drifter Data Assembly Center (DAC) at NOAA’s Atlantic Oceanographic and Meteorological Laboratory (http://www.aoml.noaa.gov/envids/gld/dirkrig/parttrk_spatial_temporal.php). The authors thank the anonymous reviewers for their thorough reviews and constructive suggestions on a previous version of this manuscript.