1 Introduction

The South China Sea (SCS) is the largest marginal sea in the Northwest Pacific and has an average depth of more than 2000 m. Every year, a large number of oceanic eddies occur in the SCS. Because of their vital roles in ocean circulation and in the marine ecosystem, oceanic eddies in the SCS have drawn a large amount of attention from oceanographers and ecologists. Ocean eddies regulate the distribution of ocean properties by transporting oceanic energy, mass, and materials in both the horizontal and vertical directions (Chen et al. 2012; Zhang et al. 2014). The rich nutrients that are trapped by eddy-induced upwelling contribute extensively to marine productivity (Lee et al. 1991; Mahadevan 2014).

The upper layer ocean circulation in the SCS is mainly dominated by the monsoon and has a strong western boundary current and a large gyre. The boundary current, along with the gyre, reverses its direction during the summer monsoon and winter monsoon. In the northern SCS, one branch of the Kuroshio intrudes into the SCS through the Luzon Strait and results in exchanges of heat, salt, and moment (Wang et al. 2006; Nan et al. 2015). Wang et al. (2003) summarized several generation mechanisms of eddies in different regions of the SCS by analyzing 86 mesoscale eddies derived from altimeter data from 1993 to 2001. The frontal instability southwest of Taiwan, related to the Kuroshio, is thought to be a factor in shedding eddies. In addition, the strong eastward current jet during the southwesterly monsoon has been found to generate eddies offshore of Vietnam, particularly a stationary pair of an anticyclonic eddy (AE) and a cyclonic eddy (CE). The vorticity from the Kuroshio front, the wind stress curl, and the interaction of strong currents with the topography also influence the eddy generation mechanisms.

Many case studies of eddies in the SCS have demonstrated and revealed the details of these generation mechanisms. Most of the case studies focused on the entire lifetime of the eddies, including their formation, evolution, propagation, dissipation, and three-dimensional (3D) structure (Chu et al. 2014, 2017; Yuan et al. 2007; Wang et al. 2015; Zu et al. 2013; Geng et al. 2016, 2017; Li et al. 2015). Wang et al. (2008) found that the relaxation of Ekman transport anomalies may have contributed to the shedding of two AEs in the northeastern SCS. Zhang et al. (2016) captured the full-depth 3D structure of an AE and CE eddy pair near the Luzon Strait in 2013/2014 based on a multi-month “SCS Mesoscale Eddy Experiment” and suggested that the dominant dissipation mechanism of the eddy pair may have been the generation of submesoscale motions. In the summer, an eddy dipole, with one strong AE and one weak CE is always located near 11° N. The formation of this eddy pair is related to the vorticity transport from the western boundary current, which is driven by the wind stress curl (Chu et al. 2017).

In our paper, we examine the multidecadal climatology of travelling SCS eddies, specifically the annual cycle, interannual variability, interdecadal variability, and long-term trend in the number, intensity, and spatial size. For doing so, we identified and tracked eddies, and determine their intensities and sizes, in a multidecadal simulation found to properly reproduce the main features of the SCS circulation (Zhang and von Storch 2017).

Recently, several statistical analyses of the multiyear variability in SCS eddies have been performed (Wang et al. 2003; Sun et al. 2016; Chen et al. 2011, 2012; Xiu et al. 2010; Li et al. 2011; Nan et al. 2011; Feng et al. 2017). Xiu et al. (2010) investigated the eddy activity in water deeper than 1000 m from 1993 to 2007, using both 7-day interval AVISO altimeter data and the output from a Regional Ocean Model System (ROMS). The authors identified approximately 32.8 eddies each year from the AVISO data (of which 52% were cyclonic eddies) and found that the interannual variabilities in eddy number and the area occupied by eddies are not correlated with El Niño events. Chen et al. (2011) focused on mesoscale eddies and identified 827 mesoscale eddies (approximately 48.6 per year) in the SCS during 1993–2009 from 7-day interval AVISO satellite data with different detection parameters. Their study found that more AEs occurred than CEs. However, using the same data set, Nan et al. (2011) detected many more CEs (41) than AEs (27) in southwest of Taiwan. Feng et al. (2017) used daily AVISO data to identify mesoscale eddies in the SCS from 1993 to 2007 and found more CEs than AEs. The number of detected eddies is sensitive to the region and the identification parameters. In the analysis of Chen et al. (2011), the eddy intensities over 17 years exhibited a weak negative correlation with the sea surface temperature anomalies in the central Pacific (given by the Nino3 index), but no correlation was found between the eddy number and El Niño activity. However, Chu et al. (2017) found in a composite analysis that the El Nino–Southern Oscillation (ENSO) influences the stationary eddy pair off the Vietnam coast. This eddy pair is one of the most important phenomena off the Vietnam coast during the summer monsoon periods, but it was absent during the ENSO transition years. They suggested that the ENSO transition events led to changes in the southwesterly monsoon and then resulted in the eastward current jet turning northward and the eddy pair disappearing. Based on the output of an eddy-resolving ocean simulation for the Earth Simulator (OFES), Sun et al. (2016) investigated the interannual variability in the eddy kinetic energy (1980–2014) in the northeastern SCS and determined that the Luzon Strait transport exhibited a modulating effect, but their study did not consider the eddy activity.

An issue, which has not been completely resolved, is to what extent the statistics of eddy formation and lifecycle may be seen as being forced by large-scale conditions, say the seasonal mean barotropic or baroclinic state (or the atmospheric conditions, which, however, likely would not act directly but indirectly via changing currents). The experiment by Tang et al. (2019) demonstrated that at least a significant part, if not the dominant part of eddy activity, is reflecting internal dynamic processes, unprovoked by external causes.

The previous analyses of eddy activity, such as the eddy number distribution and the seasonal and interannual variabilities in the SCS, have suffered from several limitations. The limited accuracy of the satellite data used for the analyses has made it necessary to limit the analysis to mesoscale eddies with high intensity. In addition, only time scales up to a few years could be considered because of the relatively short time span of no more than 30 years, starting from 1993.

The data studied in our work have some advantages, such as no observational errors, and uninterrupted, homogeneous coverage across six decades. However, the data are the output of a model and may be affected by unknown shortcomings, related to the limited spatial and temporal resolutions, the incompleteness of the physical processes, and the absence of an interactive exchange of properties between the ocean and atmosphere. However, we have no indications that the STORM simulation is inconsistent with the limited observational data in the South China Sea; a systematic comparison of some large-scale features as derived from satellite data and the simulation revealed no such inconsistencies (Zhang and von Storch 2017).

Our paper is intended to improve our understanding of the multidecadal statistical characteristics and the variability of eddy activity by examining a 61-year time model hindcast. This simulation called STORM is part of the German consortium project STORM (von Storch et al. 2012; https://www.dkrz.de/redmine/projects/storm/wiki/STORM_list_of_experiments).

Our work investigates the statistical feature and the variability on different time scales of the travelling eddies in the SCS. It does not directly contribute to an improvement of our knowledge about their dynamics. However, the provision of the statistical characteristics and spatiotemporal variabilities of long-term oceanic eddy activities are needed as a prerequisite for assessing the roles and impacts of eddies, as well as projecting the future activities. In particular, such knowledge will allow building empirical downscaling models, which allow estimating future eddy statistics, given a change in the large scales.

This paper is organized as follows. The “Data sets” section describes the details of the eddy-permitting and multidecadal STORM simulation, as well as the AVISO altimeter data set for validating STORM data. To avoid problems from the calculation of differential and integral operators on discrete fields (Chelton et al. 2011), we developed an eddy detection and tracking method that was exclusively based on the discrete sea surface height anomaly (SSHA) fields. It is explained in the “Methods” section. In the “Simulation validation” section, the STORM simulation was assessed. The “Results and discussion” section focuses on the statistical characteristics and the variabilities on different time scales. In addition, the relationship with El Niño is briefly discussed. The “Summary and discussion” section provides a summary of our work and discusses some related aspects.

2 Data sets

2.1 Satellite observations

We use the “AVISO” data set (Archiving, Validation, and Interpretation of Satellite Data in Oceanography (AVISO) [AVISO 1996, 2015; https://www.aviso.altimetry.fr/en/home.html]) for validating our simulation with respect to the eddy activity. The “Delayed Time” and “All-Sat merged” AVISO has merged data from all altimeter missions available (up to four at a given time) (Faghmous et al. 2015).

We use the gridded AVISO data set, with a grid resolution of 0.25° and time resolution of 1 day for the time period 1993–2010. An issue is the accuracy of the data set. The information provided by the AVISO handbook (Taburet et al. 2018) is not really conclusive, but we may expect the root mean square error to be larger than 1 cm, possibly much larger. While the grid resolution is about 25 km, the effective resolution in describing phenomena is also expected to be larger, namely 65 km and possibly much more (Dufau et al. 2016). The study of Amores et al. (2018) suggested the gridded AVISO data set underestimated the number of eddies due to its limited spatial resolution.

2.2 STORM simulation

The quasi-realistic modelling of the South China Sea, resolving the specifics of the circulation in that regional ocean, was begun by Pohlmann (1987), who was the first to find an upwelling regime along the Vietnam coast. Today, doing so has become a standard routine with various global and regional models.

Our analyses are based on the daily global ocean model product from the German consortium Project STORM, covering a period of 1950–2010 (von Storch et al. 2012). This simulation employs a state-of-the-art ocean model—the Max Planck Institute Ocean Model (MPI-OM). For ensuring an isotropic horizontal resolution, the bipolar grid of the model has been replaced by a tripolar grid. In total, this simulation has 3600 × 2392 horizontal grid points and 80 uneven vertical levels. The grid resolution in our domain (Fig. 1) is approximately 0.1° so that eddies are resolved (Hallberg 2013). Tidal forcing was not activated in this simulation (von Storch et al. 2012). The model was forced by the 6-hourly National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis (Kalnay et al. 1996), after a 25-year spin-up phase using the German Ocean Model Intercomparison Project (OMIP) forcing. The kinetic energy reached a quasi-steady state in the deep ocean after the spin-up. For the following analyses, we refer to the simulation as STORM.

Fig. 1
figure 1

The domain in this work. The contour shows the 200-m isobath in the SCS

STORM performed well in the research on meso- or small-scale oceanic phenomena, including eddy-related heat and salt fluxes (von Storch et al. 2016), upwelling systems (Tim et al. 2015; Yi et al. 2017), and sea level along the Ghana coast (Evadzi 2017).

The skill of STORM in representing oceanic dynamics in the SCS has been validated by Zhang and von Storch (2017). STORM agrees well with AVISO and C-GLORS (Storto et al. 2016) in the seasonal means and interannual variabilities of SSHA in the SCS. The finer resolved STORM generates stronger currents and presents more details of the strong upwelling offshore the Vietnam coast during summertime. The assessment demonstrated the ability of STORM to capture the main features of the SCS oceanic hydrodynamics.

3 Methods

3.1 Eddy detection and tracking scheme

Our study focuses on travelling eddies irrespective of the size. Therefore, in this paper, an eddy is defined to be an extremum that moves in space.

The eddy-detection scheme begins with searching for extrema (minima and maxima) in the SSHA field. A center point in a box will be defined as an extremum if the corresponding SSHA is greater (or smaller) than all other points in the box. The difference between the extremum SSHA and the mean SSHA of all other points in the box defines the relative intensity (RI) of an eddy.

We employ boxes of 5 × 5 grid cells, as suggested by Faghmous et al. (2015), who showed that a box of this size is suitable for detecting extrema, while a box with 7 × 7 cells fails with small eddies, and a box with a 3 × 3 neighborhood leads to too many extrema.

The next choice is that of a minimum RI. Too strict thresholds may break an eddy track into several pieces. An automated tracking program cannot adequately address such a situation, as it will terminate a track in the middle and derive more than one track. To avoid this problem, a mild threshold is set at the moment. According to sensitivity tests (Zhang et al. 2017), only extrema with ∣RI∣ ≥ 3 mm will be considered. In addition, an “eddy core” with at least 5 pixels, with one extremum in the center and one pixel in each of the four directions (north, south, east and west), is requested. Our new method to determine an eddy core will be introduced in detail in the “Eddy size detection methodology” section.

The tracking procedure is used to connect the local extrema at consecutive time steps (days). For each extreme at time T, another extreme at time T + 1 will be considered as a member of the same track if it fulfills the following criteria:

  • It is the closest extreme to the extreme at time T when all extrema at time T + 1 are considered.

  • Taking eddy travel speed and the spatial resolution of our data into consideration, the distance between the extrema is less than 25 km.

  • The RI properties of the extrema differ by no more than a factor of 1.5.

As we focus on eddies with long track lengths and high intensities, we apply the following filters:

  • The accumulated track is required to be at least 100 km.

  • In case some eddies move back-and-forth, a restriction on the distance from the initial position to the final position is implemented (at least 50 km).

  • In the eddy detection part, all eddy centers with a |RI| ≥ 3 mm are kept. However, tracks are kept only if the strongest extreme RI surpasses the threshold RImax = 6 mm.

  • Many small disturbances occur near the shore (especially in shallow water). Since eddies need vertical space to form and develop, a depth criterion is implemented, namely, a track must travel for over 90% of its lifetime in water deeper than 200 m.

Then, for each eddy track that satisfies the criteria listed above, the track length, the lifetime, and the maximum strength along the track are derived.

3.2 Eddy size detection methodology

To properly describe eddies that form in the South China Sea, we need a measure of the sizes of such eddies. Most eddy-detection methods are based on a physical feature, a geometric feature, or a hybrid of them (Nencioli et al. 2010). Two typical methodologies are widely used, the W-based method and the winding angle (WA) method (Chen et al. 2011; Xiu et al. 2010; Zhan et al. 2014).

The W-based method is based on a physical parameter W, which describes the relative importance of rotation compared with the deformation in the flow (Chelton et al. 2011):

$$ W={\left({v}_x+{u}_y\right)}^2+{\left({u}_x-{v}_y\right)}^2-{\left({v}_x-{u}_y\right)}^2 $$

u and v represent the eastward and northward velocity components, respectively, and the subscripts x and y indicate partial differentiation. W is used to divide the SSHA field into two parts that are dominated by divergent flow (W > 0) or rotational flow (W < 0; suspected eddy area). Then, a threshold value of W (W0) is specified to identify an eddy core (Xiu et al. 2010). The eddy size is determined by the closed contour where W = W0. It is slightly problematic to obtain W through differentiation. In addition, if only SSHA data are available, as is the case for some satellite data, only geostrophic velocity can be derived from the SSHA field. Additionally, Chelton et al. (2011) found that the eddy area identified by the closed contour of W does not closely match the closed contour of the SSHA data.

Another popular method is the WA method, which is based on the geometric assumption that the streamline around the eddy core is close to a circle or a spiral. The accumulated angle of each consecutive segment of the streamline is computed. The streamlines with absolute accumulated angles larger than 2π are defined as closed streamlines (circular or spiral curve). An SSHA extremum in the center that is surrounded by a series of closed streamlines constitutes the eddy structure in the WA method (Zhan et al. 2014). The eddy margin is characterized by the outermost closed streamline. The WA method works better than the W-based method in the eastern South Pacific (Chaigneau et al. 2008). However, the WA method also has limitations. The need to approximate the geostrophic velocity field from the SSHA field results in a problem similar to that in the W-based method.

In this paper, we present a geometry-based method that exclusively relies on the discrete SSHA field to determine the size of an eddy. The outermost closed contour of the SSHA around the eddy core is defined as the eddy edge. According to the geostrophic balance theory, water tends to flow along the SSHA contour. Therefore, the flow along the outermost closed contour can also be considered closed. The SSHA distribution within an eddy is characterized by an extremum (minimum or maximum) in the center and a number of closed SSHA contours around the eddy center, whose SSHA displays a monotonic increase (or decrease) towards the center. Because of this monotonicity, we refer to this method as “M-based.”

An iterative detection procedure aims to identify the largest (if the center is a minimum; otherwise smallest) closed contour that fulfills the monotonicity condition. The eddy size estimation procedure extends the neighborhoods surrounding the eddy core iteratively until the monotonous condition is violated. The condition is fulfilled unless all neighbors around the eddy core are greater (or less) than one threshold. The initial threshold is set to the center value of one eddy. In the following steps, the threshold is always updated by the value of the most recently detected closed contour. As the detected eddy size grows in this way, the threshold increases (or decreases) step by step. Every time the monotonicity is extended, another larger (smaller) closed contour exists in these neighborhoods, and the contour corresponds to the minimum (or maximum) from these neighborhoods. Then, the eddy size will grow by including one or more points in this new contour, and the value of this contour will be used to update the threshold to check the monotonous condition in next step.

We define the eddy size as \( d=2\sqrt{A/\uppi} \), where A indicates the area of the eddy, and d is the diameter of the eddy if it is a circle.

This M-based method is similar to the method of Faghmous et al. (2015). However, a key criterion is different. In their method, the condition for an eddy size to stop growing is another extremum included in its interior, which differs from our monotonous condition. Their condition could result in an overestimation of the size. On the other hand, the increment of the threshold in their method is defined by the user in advance. Too coarse threshold step may fail to obtain an accurate eddy core, but a step that is too fine will require extensive computational time. The appropriate steps vary for different eddies, and it is difficult to select steps that are appropriate for all eddies.

Figure 2 shows the detected results from our M-based method on 2010-01-01 from STORM data. The sizes and shapes of the masked eddies coincide with the contour lines in the SSHA fields, which reveals the ability of the M-based method to capture eddy sizes. In some cases, the eddy size is underestimated due to the interruption of the eddy structure by islands or land. This algorithm also addresses the situation that monotonicity cannot be closed as a result of the islands. The M-based method can detect the smallest eddies with structures containing only five points. In this study, we keep all the small eddies.

Fig. 2
figure 2

a The SSHA distribution (contours; units: m) and the eddies detected by the M-based method (blue: CE; red: AE) on 2010-01-01 and b the single cyclonic eddies derived from a

4 Simulation validation

In the following, we compare our simulation results with data from the AVISO gridded data set, using the paper by Chen et al. (2011) as a reference. Therefore, we use the same criteria for detecting potential eddies as in Chen et al.’s work. In our main chapter on the long-term statistics of eddies in the SCS, we use a different set of parameters (see the “Eddy detection and tracking scheme” section), because our work is aimed at investigating the statistics of eddies with all sizes, and the absolute intensity limitation will lead the small eddies to be disregarded.

For assessing if STORM simulation could reproduce eddy activity in the SCS, we derived the minima and maxima in the SSHA fields with intensity ≥ 3 cm and diameter d ≤ 35 km (Chen et al. 2011) from STORM and AVISO during the joint period 1993–2010. These minima, and maxima, may be part of a migrating eddy, but also, in case of AVISO represent observational noise. For fair comparison, the 0.1-degree STORM in this section was interpolated onto the same grids with the 0.25-degree AVISO.

Figure 3 a and b show the occurrence frequency of such minima, or maxima, at all grid points. Both data sets show potential eddy centers occurring most frequently in the Luzon Strait. Starting from the Luzon Strait, the region with high occurrence frequency extends westward along the continental shelf, passing by the southeast of China and reaching Vietnam coast. Apart from these seas, some small areas with high occurrence frequency are scattered in the central and northern part of SCS.

Fig. 3
figure 3

Comparison of STORM with AVISO. The frequency of eddy centers in each grid box during 1993–2010 from AVISO (a) and coarsened STORM (b); c and d are the same with a and b but after eddy connecting and filtering. Note the difference with Fig. 4, where the full resolution of STORM is employed

It is notable that much more potential eddies show up in AVISO data than we have in the STORM simulation. However, errors of the magnitude of 1 and more cm prevail in AVISO data set (AVISO 2015), so that it is plausible that, given the minimum intensity of 3 cm and an RMS error of 1 and more cm, some of the extrema may be artifacts generated in the observational and gridding procedures. We assume that these artifacts are short-lived and do not travel consistently in space. Therefore, we focus the comparison sequences of such minima, or maxima, which form tracks.

It may be worthwhile to examine such cases of short-lived, not-travelling minima, or maxima, in the AVISO data set, but this effort would be beyond the framework of the present study.

Eddy propagation speeds are similar to the phase speed of long baroclinic Rossby waves (Faghmous et al. 2015) and are expected to be less than 30 cm/s. Thus, an eddy cannot travel longer than 30–40 km in 1 day. Given the coarse resolution of AVISO here, we use this maximum daily travel distance, but when deriving in the “Results and discussion” section the statistics of eddy migration for STORM alone, with its higher grid resolution, we use the 25 km. The maximum daily travel distance 40 km is different from 150 km in Chen et al.’s work based on the 7-day interval data (about 21.4 km per day) either. So, for the filtering step, there is no need to use the same filtering criteria as in Chen et al.’s paper any more. Thus, we use the filters in the “Eddy detection and tracking scheme” section. After this connecting and the filtering, many potential eddies disappear in AVISO. The patterns in AVISO and STORM match better and the difference are weaker (Fig. 3c, d).

STORM reproduces a too small frequency at 14° N, 119° E and 16° N, 110° E, and a too large one at 19° N, 115° E. These discrepancies are thought to be related to the simulated Kuroshio. In the SCS, ocean models often produce stronger Kuroshio intrusion near the Luzon Strait but a weaker effect from Kuroshio in the middle and southern SCS. We suggest that the too large frequency may be related to the simulated strong Kuroshio intrusion, whereas the small ones may be related to the simulated weak Kuroshio effect in the middle of SCS.

The general pattern (Fig. 3b) coincides with the eddy probability derived from 7-day interval AVISO shown in Chen et al. (2011) as well, but their region with high frequent occurrence is connected, larger, and concentrated. Their work got the eddy probability at each point by figuring out the time when the grid was covered by a vortex, different from our only using the time when grid point is occupied by a minimum or a maximum in SSHA. A vortex can cover many grid points, but an eddy center is only located in one grid point.

In addition, it is worth noting but ignored by the previous researchers that, for altimeter missions, the mesoscale resolution capability is limited (Dufau et al. 2016). Amores et al. (2018) analyzed the extent that the gridded altimeter products can characterize ocean eddies. Their results suggested that the gridded altimeter data set underestimated the eddy density and overestimated the amplitudes. The data set can capture less than 16% of the eddies, and the limited spatial resolution of the products is attributed to the underestimation. But no other better observations are available at this time for deriving eddies and assessing the realism of the STORM data.

An additional check of the realism of the detected eddies was made concerning the temperature distribution. Anticyclonic eddies (AEs) are supposed to have a warm core, while cyclonic eddies (CEs) a cold core. We examined the temperature gradient between the core and the outermost point of the eddies at the surface and at 100-m depth at the time of peaking (maximum sea surface height (SSH) deviation of the center from the surrounding). The result of the temperature gradient and the corresponding SSHA gradient is shown in Fig. 4 for AEs and CEs in 2010.

Fig. 4
figure 4

Scatter diagram of the SSHA gradients of the peak eddy points in 2010 and the corresponding seawater temperature gradients (a and b are SST gradients; c and d are the temperature gradients at a depth of 100 m). a and c are for anticyclonic eddies, and b and d are for cyclonic eddies

For the surface, the signal is weak (top diagrams of Fig. 4); for AEs, a relationship of sea surface temperature and SSHA is hardly visible, but obvious albeit not strong for CEs. At 100-m depth (lower diagrams in Fig. 4); however, the link between SSHA and temperature becomes clear, in particular for strong eddies (|SSHA| ≥ 10 cm). This result is consistent with earlier studies; for instance, Itoh and Yasuda (2010) found most anticyclonic eddies to be warm, but clearly not all (namely 85%).

5 Results and discussion

In this section, we consider parameters such as eddy intensities, diameters, track lengths, and lifespans, and examined climatological statistics, such as frequencies and temporal variability.

5.1 Climatological characteristics

During the period 1950–2010, a total of 62,317 AE points (maxima along a track) and 115,133 CE points (minima along a track) were detected from the STORM daily data, corresponding to 1709 AE tracks (AEs; 28.0 per year) and 3331 CE tracks (CEs; 54.6 per year), respectively.

The index “IN” is defined to measure the relative portion of CEs and AEs, given by:

$$ \mathrm{IN}=\frac{N_{\mathrm{CE}}-{N}_{\mathrm{AE}}}{N_{\mathrm{CE}}+{N}_{\mathrm{AE}}} $$

where NCE and NAE represent the genesis numbers of CEs and AEs. The IN index can vary between − 1 and 1, with 0 indicating that the same genesis numbers of AEs and CEs. IN = 1 (− 1) indicates that no AEs (CEs) were generated. The monthly IN has an average value of 0.29 and a skewness of − 0.67. Negative skewness indicates that the median is larger than the expectation of 0.29; the number of CEs is over 1.8 times greater than the number of AEs in most months.

The combined effect of the enhanced spatial and temporal resolutions of the STORM data, together with the relatively mild sizes and intensity criteria, lead to many more eddies in our results, compared with the work using 7-day interval and 0.25-degree AVISO (Chen et al. 2011; Xiu et al. 2010; see the “Simulation validation” section). The detected AEs are much less than CEs. It is related to the difference in the eddy size distribution between the two kinds of eddies, and will be discussed later, when the eddy diameter (ED) distribution is addressed.

Figure 5 presents the spatial distribution of the frequencies when an eddy center passes by each grid box during 1950–2010. Most eddy centers occur in the northern SCS. Note the difference to Fig. 3d, where the analysis was done with coarsened data and only the years 1993–2010 for allowing a fair comparison with the AVISO data.

Fig. 5
figure 5

The occurrence frequency of eddy centers in each grid box during 1950–2010. Note the difference to Fig. 3, where a coarse resolution of STORM is employed for allowing a fair comparison with AVISO

The regions with the highest frequencies are located near the Luzon Strait and extend southwestward to the Vietnam coast along the continental slope, which correspond to the strong western boundary currents in the northern SCS (e.g., Zhang and von Storch 2017). Frontal instability due to the Kuroshio intrusion and wind stress curl are the major factors to generate eddies near Luzon Strait and offshore Vietnam coast. On the basis of the westward propagation characteristics of eddies in the northern SCS, the high occurrence frequency along the continental slope might result from the eddy moving from Luzon Strait.

The maximum frequency along the Vietnam coast is not as high as that in Chen et al. (2011) from AVISO 7-day interval data. This difference may occur because our algorithm filters stationary eddies, so a stationary pair of an AE and a CE off the Vietnam coast associated with the coastal wind jet (Chu et al. 2017; Lin et al. 2017) is not counted in our analysis.

To describe the distributions of eddy track length and eddy lifespan, probability density functions of eddy lengths and lifespans have been calculated (Fig. 6). In our time period, AEs travelled no further than 1941 km with an average track length of 292 km, and CEs did not travel further than 1988 km with a mean length of 274 km. The eddy lifespans ranged from 6 to 240 days for AEs and to 293 days for CEs, with mean lifespans of 36.5 days and 34.6 days, respectively. The length and lifespan distributions do not differ much between AEs and CEs, but the ranges of the lengths and lifespans of CEs are wider than those of AEs.

Fig. 6
figure 6

The probability density function of eddy track length (a) and eddy lifespan distribution (b)

The distributions of the eddy intensities (EIs) and eddy diameters (EDs) of the peak points (those are the eddy points with highest EI along an eddy track) were investigated. The largest intensities range from 0.55 to 37.3 cm, with maxima of 3–4 cm for AEs and 2–3 cm for CEs (Fig. 7a). Our eddy detection and tracking algorithm removed the eddy tracks that have highest RImax < 6 mm. Therefore, many weak eddies with an EI <1 cm are deleted; thus, the frequency of eddies with intensity 0–1 cm is low. If all eddy points are considered (not shown), the EI values vary from 0.11 to 37.3 cm, with the mean EI of 6.2 cm for AEs and 4.9 cm for CEs. AE points tend to occupy higher percentages of eddies with EI >6 cm than CE points. More than 40% of the AE points and less than 30% of the CE points have EI >6 cm.

Fig. 7
figure 7

The probability density function of eddy intensity of peak eddy points with an interval of 1 cm (a; start from 0.0 cm) and of eddy diameter of peak eddy points with an interval of 10 km (b; start from 25 km)

The corresponding diameters at the strongest eddy points range from 43 to 631 km (Fig. 7b). The EDs most frequently occur in 95–105 km for the eddy peak points. No eddy peaks associated with a diameter less than 35 km, which is related to the size limitation of 5 pixels. Most of the peak CE points (over 41%) have EDs less than 135 km, which represent more than 29% of the AEs. It is worth noting that the situation reverses with peak AEs, which occupy over 13% more than the peak CEs when the ED increases over 205 km.

If all eddy points are considered (not shown), diameters vary from 30 to 572 km for AEs and to 633 km for CEs. The percentage of large AE eddies is larger than that of CEs: 41.3% of AEs have a diameter over 175 km, whereas the percentage of CEs is only 26.2%. Much more CE cases have small eddy sizes.

We attribute the difference in the detected eddy number of AEs and CEs to the presence of smaller size CEs. The previous work (Chen et al. 2011) was based on the AVISO altimeter data, with a coarse grid space of 0.25°, and detected more AEs than CEs. As suggested by Amores et al. (2018), AVISO often underestimates the eddy activity because its coarser spatial resolution for capturing small eddies. However, our work employs the high-resolution STORM simulation, which describes smaller eddies. On the other hand, our improved algorithm allows more smaller eddies to be kept. In those additional small eddies, CEs have a much higher percentage than AEs, finally leading to the number of CEs being higher than that of AEs.

The climatological analysis presents more active CEs, with larger range but a smaller mean value of EI and ED, which means more percentage of AEs have a larger size and higher intensity.

5.2 Dominant spatial patterns

Empirical orthogonal functions are a convenient and much-used tool to identify dominant spatial patterns of variability. We have used this approach to study the annual number of eddies (EN), the annual mean EIs, and annual mean EDs of the AEs and of the CEs. Thus, our data fields X consist of 61 time “points” and 27,271 spatial “sea-points” in the SCS. To calculate the eigenvalues of X’X is not possible because of the sheer size of the matrix, but the matrix XX’ has the same nonzero eigenvalues (von Storch and Hannoschöck 1984), and the eigenvectors are related through X. Since XX’ is much smaller than X’X, we solve the eigenproblem by the smaller matrix, without encountering numerical problems.

The result is surprising—when the eigenvalues, for instance for the spatial distribution of cyclonic eddy diameter, are plotted as a spectrum (the black curve in Fig. 8), the results is a distribution rather similar to what one would expect from a white noise analysis, using an estimate of the error in estimating eigenvalues following Lawley (1956). The same result is obtained for the other two variables, EI and EN (not shown).

Fig. 8
figure 8

Eigenvalue spectra of the covariance matrix of annual anomalies of the mean diameter of annual mean cyclonic eddies in the SCS on different spatial grids—the original 0.1-degree of STORM (black), and after coarsening to 1-degree (red) and 2-degree (blue). Since 61 temporal samples are available (61 years), all eigenvalues beyond number 61 are zero and not plotted

This is a rather rare case; in most cases, at least the estimates eigenvalue of the first EOF is well separated from the rest of the estimated eigenvalues, so that only the tail of the eigenvalues, sometimes beginning at the second eigenvalue but in most cases at higher-index eigenvalues becomes consistent with a white eigenvalue spectrum, and with arbitrary eigenvectors.

The first EOF is shown in Fig. 9a—the pattern is indeed quite fragmented, and a pattern can hardly be recognized. The next eigenvectors exhibit a similar texture of randomness.

Fig. 9
figure 9

The first EOFs of annual man diameters of cyclonic eddies in the South China Sea determined from spatial distribution on the original 0.1-degree grid (left; represents about 3% of overall variance), and on the coarsened 1-degree grid (right; about 8% of overall variance)

However, when the spatial field, of originally 0.1-degree resolution, is coarsened, by binning into 1-degree boxes, the situation improves. The spectrum of eigenvalues is displayed in Fig. 8, as red curve. This time, the difference between the first two eigenvalues is about the estimated error of the first eigenvalue, so that the first two eigenvalues are about separated, possibly also the second, whereas the tail, beginning with the third number is consistent with a white noise spectrum. Also, the pattern, shown in Fig. 9b, begins to show some structure. But, the percentages of variance represented by the first two eigenvalues are meager 8% and 6%. Thus, differently from most EOF analyses, the identified patterns do not describe a relevant part of the overall variability.

Indeed, even after coarsening to a 1-degree grid, the variability of the seasonal mean diameters of CEs in the SCS reflects little variations in large-scale conditions, such as currents or wind systems. Similar results are obtained for seasonal EIs and ENs, both CEs and AEs.

An even stronger coarsening was done to boxes of size 2°. The spectrum, shown in blue in Fig. 8, is getting a little steeper, and the first two EOFs represent about 25% of the overall variance. Maybe, even the third eigenvalue is getting a little more robust. The first eigenvectors for the ED, the EI, and the EN for CEs are shown in Fig. 10.

Fig. 10
figure 10

First EOFs of annual anomalies of mean diameter (ED; left, 14% of variance) and of mean intensity (EI, middle, 12% of variance), and number of CEs (EN; right; 22% of variance) on a coarsened 2-degree grid

Interestingly, the principal components, i.e., the time coefficients, of the EOFs derived on the three different grids of cyclonic eddy diameters, are positively correlated. For the very noisy EOFs on the 0.1-degree grid, the correlation to that on the 1-degree grid is only 0.44, whereas for the 2-degree grid, it is a meager 0.32. The correlation between the PCs derived in the 1-degree grid and the 2-degree grid is 0.78. Thus, even if the EOFs on the finest grid must be considered mostly useless because of noise contamination, it represents to a minor extent the signal found on the coarsened grids.

The first pattern for the diameter has its center of activity in the southern SCS, with a mostly uniform sign; the spatial variability of the intensity is given by a large area with more intense and a small area near Luzon strait of less intense (and vice versa) CEs; the number of CEs is dominated by a multipole with an overall spatial mean of zero.

Similar results hold for AEs (Fig. 11). The AE patterns are limited to the central and northern part of the SCS, reflecting the absence of AEs in the southern SCS. The first EOF of ED describes a general enlarging or shrinking of eddies, while the intensity and number are summarized by a dipole, favoring either the southern or the northern part of the Luzon strait throughflow.

Fig. 11
figure 11

First EOFs of annual anomalies of mean diameter (ED; left, 16% of variance) and of mean intensity (EI, middle, 18% of variance), and number of AEs (EN; right; 19% of variance) on a coarsened 2-degree grid

Thus, building, in the spirit of empirical downscaling, links between spatially disaggregated distributions of eddy parameters to large-scale patterns of, say, regional wind, barotropic stream function, current shear, or vertical stability, will be a challenge. However, this will be dealt with in a subsequent paper. We hypothesize that the internal variability, unprovoked by changing large-scale conditions—thus “noise” in the spirit of Tang et al. (2019)—is a significant factor for the variations of eddy activity. Of course, for verifying the hypothesis, more analysis is needed.

5.3 Seasonal variability

In this section, the seasonal features of accumulated EN, mean EI, and ED are discussed.

More CEs are generated than AEs in all the 12 months (Fig. 7a). The number of eddies generated in each month exhibits a significant seasonal variation. Most CEs are generated in the winter half year, with a peak of 364 eddies in March. More AEs are formed in the spring than in the autumn. The minimum number of eddies generated is only 105 in September.

For all 1709 AEs and 3331 CEs, the EIs and EDs of the peak eddy points are determined. The monthly means of EI and ED of peak eddy points (Fig. 12b, c) exhibit different variations. AEs always have stronger intensities and larger sizes than CEs. In addition, the monthly EI and ED data reveal an annual cycle, peaking in July with EI = 9.14 cm and ED = 205 km for AEs, while these values peak in October for CEs with EI = 7.93 cm and ED = 192 km. From July to October, the ED and EI values of peak AEs decline sharply, from the maxima hitting the minima. On the contrary, rapid growth occurs in both the ED and EI variability of CEs and increase to the maxima during these 4 months. In terms of the monthly means of all 62,317 AE points and 115,133 CE points, an annual cycle is also apparent (not shown), AEs peak in August, with EI = 7.01 cm and ED = 174 km. CEs reach peaks in November when EI = 5.58 cm and ED = 147 km.

Fig. 12
figure 12

The annual cycle of the multiyear sum number of eddy tracks (a) and the means of EI (b) and ED (c) of the peak eddy points

The presence of well-defined annual cycles points to the presence of a link between the eddy activity and large-scale conditions, in wind, currents of stability, which is undergoing annual variations (cf. Zhang et al. 2017), in contradiction to the above-formulated hypothesis that the variability, is mostly generated internally.

Therefore, we have examined the distributions of differences between months, when the diameter, the intensity EI, and the number EN have their multiyear maximum and multiyear minimum. This is done for all years, within which at least one eddy is found. The result is presented in Fig. 13 as Box-Whisker diagrams, with minimum, lower quartile, median (mean quartile), upper quartile, and maximum, and, sometimes outliers (as dots).

  • The annual cycle is not very stable. In up to 25% and more of all years, the normally positive difference between maximum and minimum is reversed.

  • The variability for anticyclonic eddies is stronger, in particular for diameter and intensity, so that referring to an annual cycle is of limited usage.

  • The annual cycle of cyclonic eddies is more stable, with less than 25% of years showing a reversed sign, in terms of number, diameter, and intensity.

Fig. 13
figure 13

Box-Whisker diagram of the annual differences of eddy parameters in the month of multiyear maximum and multiyear minimum, for describing the stability of the annual cycles. The gray box represents all cases between the upper and lower quartile, and the contraction the medium. The whiskers mark the minima and maxima, and the lone dot an outlier

On the other hand, the annual cycle, as being manifested in the large-scale monsoonal changes in winds and large-scale currents, is strong (not shown, but see also Zhang and von Storch 2017). Therefore, we conclude that also the annual cycle of the eddy properties is significantly affected by factors unrelated to the large-scale forcing which represent the annual march of the sun.

5.4 Interannual variability

To assess the long-term variability of the eddy properties, the annual time series of the number of eddies, EIs, and EDs have been plotted (Fig. 14). Interannual variability dominates the annual series of eddy number. There is some decadal variability with maxima in the 2000s and 1970s, and minima in the 1980s and 1950s.

Fig. 14
figure 14

The time series of the annual number of eddy tracks. a Number found in STORM using the standard set-up of RI = 3 mm and RImax = 6 mm. The axis on the left (right) is for AEs (CEs). b Series according to STORM with a standard set-up, AVISO satellite data (Xiu et al. 2010) and the simulation by Xiu et al. (2010). The axis for the STORM result is on the right

Xiu et al. (2010) published the annual number of generated eddies derived from the AVISO satellite data and their own modelled data during 1993–2007 by using the W-based eddy detection and tracking method. Figure 14 b combines our results from STORM and their results to further assess our results. Due to their strict criteria for eddy size, water depth, and so on, the number of eddy tracks identified in their study is less than ours. In terms of variability, STORM outperforms Xiu et al.’s simulation when compared with the AVISO satellite data. The correlation coefficient of the annual number of eddies during 1993–2007 between AVISO and STORM is 0.21, which is much higher than the correlation coefficient of 0.03 between AVISO and Xiu et al.’s simulation. This comparison suggests that STORM simulation describes more realistically the variability of the travelling eddies in the SCS, which should benefit from the realistic reproduction of the SCS dynamics (Zhang and von Storch 2017).

With regard to the annual mean EI and ED values of the peak eddy points (Fig. 15), the annual series also presents mostly interannual variability with very little decadal variability and no trends. AEs on average have higher EIs and larger EDs than CEs in most years, which is consistent with the climatological results. In the “Climatological characteristics” section, the skewness analysis and the eddy property distribution reveal much more CEs with small size and low intensity than such AEs in the SCS. This feature of the eddies in the SCS may be related to the different generation mechanisms of AEs and CEs.

Fig. 15
figure 15

The annual time series of EIs (a) and EDs (b) of peak eddy points

The AEs vary more strongly than the CEs. The annual mean EDs and EIs are strongly correlated, with correlation coefficients of 0.76 and 0.68 for AEs and CEs, respectively. Thus, strong eddies tend to be large, as is illustrated by the scatter diagrams of the two parameters (Fig. 16). The annual mean EDs and EIs of all eddy points (not shown) show similar features, with EI-ED correlation coefficients of 0.88 and 0.78 for AEs and CEs, respectively.

Fig. 16
figure 16

The scatter diagram between eddy diameter and eddy intensity of AEs (a) and CEs (b) derived for all points along all eddy tracks

6 Summary and discussion

In this work, the long-term variability and statistical features of the travelling eddies in SCS during the period 1950–2010 have been documented by means of the eddy-permitting and multidecadal output from the ocean global model simulation STORM, which was forced by a 6-hourly NCEP atmospheric reanalysis. STORM simulation has a horizontal resolution of 10 km on average, which enables it to resolve eddies in the SCS.

For a limited time, beginning in 1993, the satellite data (AVISO) of the SSHA in the SCS are available—we have used these data to find out if the STORM simulation data generate consistent numbers and characteristics of eddies in the SCS.

Sometimes the claim is made that “Since STORM is running in a realistic hindcast mode, in principle, it should be possible even to reproduce the eddy structures that are observed at a specific day.” This claim is wrong, since the model is run in the “climate mode,” i.e., driving only through upper and lateral boundary values with no deterministic influence by initial values. What happens in the interior of the South China Sea is in part related to the large-scale forcing of atmospheric states (in particular winds) and to lateral boundary conditions, but in part generated internally (“noise”; see Tang et al. 2019). Thus, we should be able to successfully hindcast the large-scale state of the dynamics in the South China Sea, but not the small scales (even if we have to add the caveat that it is unknown if periods of “intermittent divergence in phase space” take place in such an oceanic set-up; cf. Weisse et al. 2000). This is also the case in the South China Sea. To demonstrate this challenge, we have randomly selected a day, namely 31 December 2010, and have compared the total SSH anomalies as given by AVISO and by STORM, a spatial low-pass components (first 10 EOFs) and a spatial high-pass (difference between full and low-pass filtered fields). The first is quite similar on the large scale, but the difference shows small-scale patterns all over the South China Sea (Fig. 17).

Fig. 17
figure 17

a and b are the anomalous fields of the daily detrended SSHA on 2010-12-31 from AVISO and STORM. c and d are similar to a and b, but the anomalous fields are reconstructed from the first 10 EOFs. e and f is the difference between the anomalous fields and the reconstructed fields. (Note the same scales)

Unfortunately, the accuracy of local data in the AVISO gridded data set suffers from significant uncertainties, of the order of 1 cm and more (Taburet et al. 2018); also, the spatial resolution in AVISO is less than in STORM. A large underestimation of eddy activity in the AVISO due to its coarse spatial resolution is found by Amores et al. (2018). This contributes to differences between AVISO and STORM, but some of the differences may be due to the insufficiencies of the model simulations. Unfortunately, we cannot quantify the different contributions. But, given the large uncertainty of the AVISO data, we may diagnose a general consistency between AVISO and STORM. Unfortunately, in past studies, which have employed only AVISO data for describing cases and statistics of eddies, these inaccuracies have not always been recognized. Our analysis extends these earlier studies, providing more detail, such as the stability of the annual cycles, the strength of the link of warm/anticyclonic and cold /cyclonic eddies, and the organization of spatial co-variability (EOFs). By having a much longer time series (61 years), we were also able to examine multiyear variability.

For deriving the eddy activity from STORM SSHA data, the M-based eddy detection and tracking method, which makes use of the discrete SSHA fields, has been developed. This algorithm avoids the problems from the calculations of differential or integral operators. In addition, eddies are identified as the moving extremes in space and time. Stationary eddies are not considered in this work.

On average, 28.0 AEs and 54.6 CEs per year are found in the STORM data in the SCS according to the criteria and the parameters in the M-based detection and tracking procedure. Those tracks cover 62,317 points (for AEs) and 115,133 points (for CEs), which means that the average lifespan of AEs is 36.5 days and that of CEs is 34.6 days. More CEs are detected than AEs. The lifespans of eddies in the SCS range from 6 to 293 days and travel up to 1988 km. The spatial distribution results show that the most frequent regions for eddy tracks to travel are the Luzon Strait and towards the continental slope to the Vietnam coast. Those eddy points have an EI range of 0.11 to 37.3 cm and an ED range of 30 to 633 km, while the strongest eddy points along the eddy tracks have an EI range of 0.55 to 37.3 cm and an ED range of 43 to 631 km. For both kinds of eddies, the stronger eddies tend to be larger in size.

CEs are much more active than AEs, but the AEs with high intensities or large diameters are more frequent than similar CEs. Our algorithm and STORM simulation allow much more small eddies to be kept, compared with the previous work (Chen et al. 2011). In addition, AVISO often largely underestimates the eddy activity (Amores et al. 2018), because the spatial resolution is too coarse to capture the small eddies. Then many more small eddies appear in CEs instead of AEs in our work, which could contribute to the large difference in the detected eddy number of AEs and CEs.

When doing an EOF analysis, no significant spatial patterns are found, when the analysis is done on the original 0.1-degree grid. Since the numerical solution is obtained by the eigenanalysis of the smaller (61 × 61) XX’ covariance matrix, instead of the much larger X’X matrix, this strong contamination of the spatial co-variability by internal small-scale dynamics is not an unwanted artifact of solving large numerical problems but real, in consistence with the result of Tang et al. (2019). The absence of significant patterns seems inconsistent with the hypothesis that the spatial distribution of eddy activity would be noteworthily be constrained by some large-scale “drivers,” such as large-scale wind patterns or current patterns. Instead, the formation of eddies is more an issue of internal variability, largely unprovoked by external drivers. However, more analysis and numerical experimentation are needed to arrive at more robust results.

The STORM data suggest that the variability in the eddy properties in the SCS is dominated by intra-annual variability. Eddy number peaks in March and the fewest number of eddies occurs in September for AE and in August for CEs. For the EI and ED of the peak eddy points, the peaks occur in July for AEs and October for CEs. When comparing the annual cycle across different years, we find that they differ strongly from year to year, also pointing to the above-mentioned internal variability.

The interannual variability is strong, while the decadal variability is weak, and a long-term trend is not found.

One purpose of this analysis was to prepare a data set that was suitable for the statistical downscaling of eddy properties. Our results presented here indicate that a spatially resolving estimation of eddy properties as conditioned by the large scales in the South China Sea cannot be done straightforward. This aspect will be worked out in more detail in a forthcoming paper. However, it may be said already here, that an alternative may be to use spatially aggregated eddy properties, such as the total number of eddies.

As a first attempt to quantify links to large-scale conditions, we have tried to relate the eddy parameters to ENSO conditions. El Niño is one of the most important phenomena in the tropical ocean. For the SCS, El Niño is considered significant because of its impact on the wind stress curl over the SCS and the ocean circulation in the SCS (Fang et al. 2006; Wang et al. 2006). In addition, with regard to the annual time series of eddy number and EI when considering all eddy points, a simultaneous drop is observed during the strong El Niño year of 1998. The study of Tuo et al. (2018) found a changing influence of El Niño on the eddy activity in the SCS around 2004. Based on the altimeter data (1993–2013), they found significant and negative correlation between El Niño and the eddy activity showed up before 2004, but disappeared afterward.

However, we found that the variability of the eddy activity shows weak or no correlation with El Niño. When the El Niño3.4 index (https://www.esrl.noaa.gov/psd/data/climateindices/list/for info) is correlated with the number of generated eddies and the eddy parameters of all eddy points on a monthly basis, we find very small correlations, which seem to be irrelevant (see Table 1). To take into account the effect of serial correlations (Zwiers and von Storch 1995), we assume that a 12-month time difference is needed to have somewhat independent samples. Using this “efficient sample size,” we found that none of the links were significant.

Table 1 The correlation coefficients between El Niño3.4 index and the eddy properties

Tuo et al. (2018) investigated 10-year sliding correlations of eddy number and El Niño3.4 index and suggested significant negative correlations. However, in their work, after 2005, the negative correlation weakens and reach zero around 2008. After that, weak positive correlation emerged. In different time windows, their correlations are quite different and sometimes even opposite. Our work focuses on their relationship in the whole 61 years, which differs from Tuo et al.’s work (2018).

We conclude that a relevant link to the ENSO dynamics on the tropical Pacific does not exist for the formation of eddies in the SCS, again supporting our hypothesis that the main source of variability is internal.