Observed 3D Structure, Generation, and Dissipation of Oceanic Mesoscale Eddies in the South China Sea

Oceanic mesoscale eddies with horizontal scales of 50–300 km are the most energetic form of flows in the ocean. They are the oceanic analogues of atmospheric storms and are effective transporters of heat, nutrients, dissolved carbon, and other biochemical materials in the ocean. Although oceanic eddies have been ubiquitously observed in the world oceans since 1960s, our understanding of their three-dimensional (3D) structure, generation, and dissipation remains fragmentary due to lack of systematic full water-depth measurements. To bridge this knowledge gap, we designed and conducted a multi-months field campaign, called the South China Sea Mesoscale Eddy Experiment (S-MEE), in the northern South China Sea in 2013/2014. The S-MEE for the first time captured full-depth 3D structures of an anticyclonic and cyclonic eddy pair, which are characterized by a distinct vertical tilt of their axes. By observing the eddy evolution at an upstream versus downstream location and conducting an eddy energy budget analysis, the authors further proposed that generation of submesoscale motions most likely constitutes the dominant dissipation mechanism for the observed eddies.


Results
Observations. To investigate the 3D structure and generation and dissipation processes of oceanic eddies in the SCS, we designed and conducted the S-MEE in the northern SCS. As part of S-MEE, two bottom-anchored subsurface mooring arrays were deployed along the historical pathway of mesoscale eddies in this region (Fig. 1a,  Supplementary Fig. 1). The first array, cross-shaped and consisting of 10 moorings, was deployed west of the Luzon Strait in the generation region of eddies in late October 2013, and was successfully recovered in early June 2014. Before this experiment, mooring M1 at the western end of this array had been deployed and recovered two times, and it provided continuous measurements for a total length of more than 2 years. The second array, consisting of 7 moorings along a line normal to the local isobaths, was deployed in the dissipation region of eddies in mid-January 2013, and successfully recovered in mid-June 2014. All the moorings were equipped with acoustic Doppler current profilers (ADCPs), recording current meters (RCMs), conductivity-temperature-depths (CTDs) and temperature chains to monitor horizontal current velocity, pressure, and temperature/salinity (T/S) in the whole water column (see Methods and Supplementary Table S1). The moored velocity and T/S data are used to analyze the 3D structure and generation and dissipation processes of oceanic eddies.
In order to identify and track the oceanic eddies, AVISO merged satellite altimeter sea level anomaly (SLA) and absolute surface geostrophic velocity data during the experiment period (October 2013-June 2014) are used 34 . The temporal and spatial resolution of the AVISO data is daily and 1/4°, respectively. From the 9-month SLA and surface geostrophic velocity maps, a total of 5 distinct oceanic eddies were observed to cross the mooring arrays. Among these eddies, one anticyclonic and cyclonic eddy pair was fully captured by the first mooring array from mid-November to early January (Fig. 1b-e). After crossing the first array, the anticyclonic eddy (AE) propagated southwestward and was well measured by the second array from mid-January to early February (Fig. 1f,g). The AE during its two stages (at the first and second arrays) and the cyclonic eddy will be referred to as AE0, AE1 and CE, respectively, below. The analysis of this study will primarily focus on this eddy pair. On the basis of real-time SLA information, we also conducted two transects across the center of AE in early December 2013 and mid-January 2014, respectively ( Supplementary Fig. 2). Along these two transects, high-resolution hydrographic and turbulent mixing measurements were made. The obtained turbulent kinetic energy dissipation rate ( Supplementary Fig. 3a,b), CTD and shipboard ADCP velocity data are used to analyze the dissipation of AE.  Fig. 1b-g). Before AE was generated, the Kuroshio, the western boundary current in the North Pacific subtropical gyre, appeared to intrude into the northeastern SCS in a form of loop current. After the generation of AE, the Kuroshio returned to its normal path, leaping across the Luzon Strait. This result suggests that the origin of AE is from the shedding of the Kuroshio current. For the CE, it was generated right behind the AE and seems to originate from the western part of the Taiwan Strait. The upper-layer (< 350 m) water within the AE0 shows properties close to the Kuroshio water that is warmer and saltier than the SCS water (Fig. 2). The upper-layer water within the CE, on the other hand, is colder and fresher than the SCS water (Fig. 2). In addition, the high-resolution satellite sea surface temperature (SST) reveals a warm-and a cold-core within the AE and CE, respectively ( Supplementary Fig. 4), and the historical hydrographic observations reveal a low temperature and salinity tongue in winter in the western Taiwan Strait 35 . All these results verified the AVISO observations that the AE was shed from the Kuroshio, while the CE originated from the western part of the Taiwan Strait which trapped cold and fresh near-shore water into the northern SCS. Based on 22 years (1993-2014) of the AVISO data, a total of 19 events of anticyclonic eddy shedding associated with the Kuroshio intrusion were detected. For more than 80% of these cases, a CE is immediately generated behind the AE to its northeastern side, suggesting the common occurrence of such eddy pairs ( Supplementary Fig. 5). Formation of the trailing CE is likely due to the vortex tube stretching induced by the offshore transport on the northeast side of AE combined with the cyclonic velocity shear associated with the intruding Kuroshio south of Taiwan. Further study, however, is needed to verify this hypothesized formation mechanism for CE. 3D eddy structures. By combining the mooring-array and AVISO data, we constructed 3D structures of the eddies over the full water column (see Methods). Figure 3a-c show the 3D structures of horizontal velocity for AE0, CE and AE1, respectively. Vertically, the eddies are surface-intensified, with the velocity sharply decreasing with depth. At the same time, the eddies are found to be deep-penetrating, extending from the surface to near bottom (with depth over 2000 m). The maximum velocity magnitude at the 2000 m depth and deeper exceeds 3.5 cm/s and 5.0 cm/s for the AE0 and CE, respectively, which is much stronger than the time-mean current detected in the deep SCS circulation 36,37 . The deep signals of the eddy pair may be due to quick geostrophic adjustment in response to pressure anomalies caused by the upper-layer mass loading 38 . Horizontally, the velocity structure is axis-asymmetric and the velocity is enhanced on the northeastern side of AE0 and southwestern side of CE due to the lateral alignment of the eddy pair. The large current anomalies associated with the eddy pair are expected to play an important role in deep-ocean material transport as proposed by previous studies 32, 39,40 . Another notable feature for the 3D structures is that the axis of the eddies (purple dotted lines in Fig. 3) strongly tilts southwestward from surface to bottom. From the surface to 2000 m, the tilting distance reaches up to ~100 km for both AE0 and CE. For AE1, the tilting distance extends to ~150 km from the surface to 1500 m. These tilting structures of eddies can also be confirmed from the depth-time plots of velocity at each mooring ( Supplementary Fig. 6e,f). Because the tilting direction is generally coincident with the propagation direction of the eddies (yellow and magenta lines in Fig. 1a), the phase of velocity in the deep layer is detected to lead that in the upper layer ( Supplementary Fig. 6e,f). Quantitatively, the velocity phase at 2000 m leads that at 100 m by about 15 days, which, if multiplied by the eddies' propagation speed of 0.08-0.11 m/s, agrees well with the aforementioned tilting distance. Although the reconstructed 3D thermohaline structures are confined to the upper 1000 m due to the limit of temperature chain data (Supplementary Table S1), tilting features similar to the velocity structures can also be clearly seen in them ( Supplementary Fig. 7), i.e., eddy center tilts southwestward with increasing depth.
Dissipation of the eddy. After its generation in early November, the AE quickly grew, and its amplitude (defined as the maximum SLA, see Methods) reached the maximum in early December ( Fig. 1b-g). Then the AE propagated southwestward along the continental slope and its amplitude gradually decayed. The AE has a lifespan and travelling distance as long as 4.5 months and 900 km, respectively, and it finally disappeared north of the Xisha Islands (approximately at 112.5°E, 17.7°N) in mid March 2014. Compared to the AE, the CE is very short-lived and only has a lifespan of one and half months. After reaching its maximum amplitude in early   Table 1. Terms contributing to the dissipation of AE. All energy values are volume-integrated over the AE. Contribution is the ratio of each energy term to dEE. Resi denotes the residue of the eddy energy budget. See detailed definitions and calculation of each term in Methods. Values to the right of the ± symbol denote the uncertainty of each term. The uncertainty is primarily associated with the fact that there were no in-situ observations between the first and second mooring arrays and we assumed a linear evolution for AE between them (from AE0 to AE1). Here, it is estimated as the half value of the difference between that we integrated the observed terms at AE0 and AE1, respectively, during the whole period.
Scientific RepoRts | 6:24349 | DOI: 10.1038/srep24349 January, the CE was quickly dissipated. Finally, it was observed to disappear just southwest of the first mooring array in mid January. Because the AE was well captured by both of the mooring arrays, the eddy dissipation analysis hereafter is primarily focused on the AE. From AE0 to AE1, the AE's amplitude, kinetic energy (EKE) and available potential energy (EPE) had a prominent decrease by 40%, 47% and 35%, respectively (Figs 1b-g and 4a,b). In order to investigate the dissipation mechanism of AE, we carried out an eddy energy budget analysis based on available observational data. The detailed quantification of each eddy energy budget terms is given in the Methods section. As summarized in Table 1 and schematically depicted in Fig. 5, summation of all the terms gives rise to a small residue (11% of the observed dissipated eddy energy (dEE) from AE0 to AE1), suggesting that the eddy energy budget is largely closed when considering the estimation uncertainties. During the decay process, the eddy-mean flow interaction still fed the AE by providing it eddy energy at the expense of the energy of mean circulation and accounted for 10% (± 5%) of dEE. Both the surface wind stress and bottom drag did negative work on the AE (i.e., weakening the AE), but they accounted for only 3% (± 2%) and 18% (± 8%) of the dEE, respectively. Vertical eddy diffusion, with an observed vertical eddy diffusivity on order of 10 −4 m 2 s −1 (Supplementary Fig. 3c,d; cf. background level is 10 −5 m 2 s −1 in the ocean interior), explained 20% (± 8%) of the AE's dissipation. The enhanced vertical mixing is primarily caused by breaking of strong internal tides radiated from the Luzon Strait 41,42 and it may also be partially due to the breaking of near-inertial waves radiating from the eddy itself 43 .
In the presence of AE, submesoscale motions (see definition in Methods) are prominently enhanced (Fig. 4c,  Supplementary Fig. 8a), with the submesoscale eddy energy within the AE more than double that outside of the AE (Fig. 4c). The close correlation between the submesoscale EKE (KE sm ) and mesoscale EKE (KE ms ) is confirmed by examining the 2-year moored observations ( Supplementary Fig. 9), implying that the mesoscale eddies can feed the generation and growth of submesoscale motions. If we assume that the increased submesoscale energy within the AE all derives from AE, this part of energy would account for 58% (± 12%) of the dEE over the period from AE0 to AE1 (see the detailed calculation in Methods). This implies that the generation of submesoscale motions, or downscale energy transfer, likely plays a dominant role in the dissipation of AE.
Taking into account that the submesoscale motions act on horizontal velocity shear of AE, this amount of dissipated energy (D h ·dt in Table 1) yields an "effective horizontal eddy viscosity" (i.e., A h ) of 75-110 m 2 s −1 for submesoscale motions (see Methods), which is within the range of previous estimates based on observations and numerical simulations 28,44,45 . The reasonable A h value estimated here in turn supports the notion that the downscale energy cascade/transfer is likely the primary dissipation mechanism for AE.

Discussion
Based on observations from individual moorings in the northern SCS, previous studies found that, during the eddy pair events the deep current (> 2000 m) was intensified and tended to flow in opposite directions against the surface current 32,40 . The deep-penetrating and vertical-tilting eddy structure found in this study provides a good explanation for this previously observed phenomenon. Because the eddy propagation in the northern SCS follows the regional sloping bottom topography, the topographic β effect, which exerts more influence to the water column near bottom in the stratified ocean 46 , is likely the cause for the observed vertically tilting structures. Indeed, previous modelling studies showed that for an isolated equivalent barotropic eddy overriding a sloping topography, the topographic β effect works to shift the eddy core in the lower layer more to the right in the down-slope direction than in the upper layer 47 . It is important to emphasize that the mesoscale eddies reconstructed in the open ocean exhibited no such strongly tilted 3D structures 4,8-10 and, as such, the observed tilting structure is likely a unique characteristic of eddies in the northern SCS or other marginal seas with broad-scale continental slopes.
In defining submesoscale motions in the present study, we have adopted the separation method based on time scales (see Methods). It is worth pointing out that the submesoscale motions thus derived show characteristics consistent with the classical submesoscale dynamics. Specifically, different from the quasi-geostrophic mesoscale eddies whose EPE ms is larger than EKE ms (Fig. 4a,b), EKE sm is comparable with EPE sm for the submesoscale motions (Fig. 4c). This results in an order one Rossby number = Ro EKE EPE 38 , highlighting the ageostrophic feature of the submesoscale motions. Compared to the mesoscale eddies, submesoscale energy decreases more rapidly with depth and is primarily concentrated within the surface mixed layer (Fig. 4, Supplementary Fig. 10). All these characteristics indicate that the submesoscale motions identified in this study fall in the classical submesoscale regime 26,48,49 .
In addition to the perspective based on the time scale separation, submesoscale motions can also be seen from a spatial scale perspective. High-resolution (1 km) SST images reveal fine submesoscale structures (of horizontal scales under 50 km) along the peripheries of the cores of AE and CE ( Supplementary Fig. 4). From the high-resolution (1 km) shipboard velocity measurements ( Supplementary Fig. 2c,d) we also found that the submesoscale motions with horizontal scales smaller than 50 km show enhanced EKE sm on the peripheries of the AE's core. In addition, the high-resolution velocity data produced an energy spectrum with a spectral slope of k −2 (k is the horizontal wavenumber), which characterizes the submesoscale regime 26,50 . Based on these results, the timescale-based submesoscale motions defined in this study show no substantial difference with those defined based on spatial scales. Although the eddy energy budget analysis was applied to only the AE case, enhanced submesoscale energy was also found in other mesoscale eddy cases based on the 2-year moored observations at M1 (Supplementary Fig. 9). Based on the above observational evidence, we conclude that the downscale energy transfer to submesoscale motions likely plays an important role for the dissipation of oceanic eddies in the northern SCS. Finally, we note that the downscale energy cascade derived in this study is for the total eddy energy (EKE + EPE) signals. An inverse energy cascade is likely to occur if the EKE only is considered.

Methods
In-situ moored data. To acquire the full-depth current velocity and T/S data, self-contained instruments including ADCPs, RCMs, CTDs and temperature chains were heavily mounted on the moorings (Supplementary Fig. 11 and Table S1). The instrument configuration on each mooring is similar but not identical due to the different water depth at the mooring locations. Nominally, the moorings were equipped with an upward-looking and a downward-looking 75 kHz ADCPs (sometimes only one upward-looking ADCP) at ~500 m to measure the velocity in the upper 1000 m (or 500 m) water column. Below that, discrete RCMs provided point measurements of velocity every 400-600 meters, with the upmost one at the nominal depth of 1500 m (or 550 m) and the lowermost one at depth 50-200 meters above the ocean bottom. Immediately beneath each RCM (for most cases), a CTD was mounted on the mooring to measure the T/S and pressure. The moorings were also equipped with temperature chains extending from tens meters beneath the sea surface to ~1000 m (or ~500 m), consisting of 2-6 CTDs and dozens of thermometers. Detailed mooring configurations can be found in Table S1 of the supplementary information.
The high-resolution raw velocity data were firstly hourly averaged. For each mooring, the hourly velocity data from both ADCPs and RCMs were then linearly interpolated onto the uniform 10 m levels between near surface (40-60 m) and near bottom. In order to remove the major tidal and inertial signals (the inertial period is 31-38 hours in the study region), all the hourly velocity time series were 2-day low-pass filtered with a fourth-order Butterworth filter. The low-pass filtered velocity data were further averaged on a daily basis. The processing procedure as described above was similarly applied to the temperature chain data, resulting in the daily-averaged temperature data at 10 m depths between the near surface and ~1000 m (or ~500 m). The aforementioned daily velocity and temperature data are used for analysis in this study.
Mesoscale and submesoscale motions. The observed temperature and velocity time series show strong intraseasonal variations with periods between 2 and 100 days (Supplementary Fig. 8). The intraseasonal variations are associated not only with the mesoscale eddies, but also other high-frequency motions embedded in the eddies (Supplementary Fig. 8a). The periods of these high-frequency motions are generally shorter than 10 days, which can be clearly seen from the power spectra ( Supplementary Fig. 8). To be distinguished from the oceanic eddies, the high-frequency motions with periods shorter than 10 days are regarded as submesoscale motions in this study. The temperature anomaly (T′ ) on the intraseasonal timescale is similarly separated into two terms, filter. The same separation was also applied to the velocity data. Throughout the paper, the subscripts 'ms' and 'sm' are used to denote variables associated with the mesoscale and submesoscale motions, respectively.
Construct the 3D structures. The moored and altimeter data are combined to construct the 3D structures of the oceanic eddies (AE0, CE and AE1). Before the construction of eddy's ′ v ms and ′ T ms field, the Okubo-Weiss parameter was first used to identify the eddies from the SLA maps 51,52 . The criterion of W < − 0.2σ w is used to define the eddy element within the eddy core, where W is the Okubo-Weiss parameter and σ w is the spatial standard deviation of W in the study region in each day. W is defined as where ′ u g and ′ v g are geostrophic velocity anomalies calculated from the SLA through geostrophy. We use the mean longitude and latitude of the eddy elements to represent the eddy center (x 0 , y 0 ), and the radius of an equal-area circle to represent the eddy radius R 0 . We define the amplitude of the eddy (H 0 ) as the maximum SLA around the eddy center. With the above altimeter-derived surface information, we tracked the eddies based on a previously used eddy tracking algorithm 30 . Then, we derived the trajectory (Fig. 1a) and temporal evolution of each oceanic eddy.
For eddy AE0, its shape was relatively steady from 14 November to 27 December 2013. Based on the moored data during this period, the 3D structure of AE0 was constructed as follows. Taking the AE0 center as the origin, we constructed a new cartesian coordinate in each day. Then, all the ′ v ms and ′ T ms profiles located within 200 km of the eddy center were projected onto this new coordinate system. Considering the variation of the eddy size during the given period, we normalized the coordinates (x, y) of the profiles using , where (x n , y n ) denote the normalized coordinates, and R 0 and R 0 are the instantaneous and mean eddy radius, respectively. ′ v ms and ′ T ms were normalized by H 0 of the eddy with the similar method. At each depth, the obtained ′ v ms and ′ T ms were then mapped onto a 10 km by 10 km grid using an objective interpolation scheme. The weight function for the mapping has an inverse distance form, where L = 50 km is the interpolation radius, and d is the distance between the data point and grid point. For each grid point, only the data points with d < L are used. Finally, the 3D structure of ′ v ms in the full water column and ′ T ms in the upper 1000 m were obtained for AE0. Using the data from 14 December 2013 to 8 January 2014 (from 14 January to 10 February 2014), the 3D structure of CE (AE1) was also constructed based on the same method. The similar method to construct eddy's 3D structure was also used by the previous studies 4, [8][9][10] . The difference is that the previous studies reconstructed the mean structure of many eddies because there were very few observations for individual eddies. To get the mean structure, those studies had to assume that all eddies have similar 3D structures.
Eddy energy budget. From the primitive equations for Boussinesq fluid we derived the eddy kinetic energy (EKE) and eddy available potential energy (EPE) equations 53,54 . Adding the EKE and EPE equations and integrating over the whole water column, we obtained Here, primes denote the mesoscale anomalies; overbars denote the time average; ∇ h is the horizontal gradient operator; τ v ( ) top w is the horizontal velocity (wind stress) at the sea surface; τ v ( ) b b is the horizontal velocity (drag stress) near the sea bottom; ρ = − 1030kgm 0 3 is the mean sea water density; N is the buoyancy frequency; A h , A v and K v denote the horizontal eddy viscosity, vertical eddy viscosity and vertical eddy diffusivity, respectively; the other symbols and notations are standard. The above equation can be written symbolically as The rhs terms, respectively, are rate of the wind stress work (W w ) and the bottom drag work (W b ), baroclinic (BC) and barotropic (BT) conversion terms of the eddy-mean flow interaction (if positive, the mean flow transfers energy to the eddy), divergence of perturbative pressure work (P d ), and terms of energy dissipation caused by horizontal eddy viscosity (D h ), vertical eddy viscosity (D Av ) and vertical eddy diffusivity (D Kv ). The last two terms (D Av and D Kv ) are termed D v for simplicity.
AE was observed two times by the mooring arrays, one during its mature stage (i.e. AE0) and the other during its decay stage (i.e. AE1). The energy budget of AE, therefore, can be analyzed in a Lagrangian framework. Using the formulas in equation (1) ρ′ ms , respectively. The ρ′ ms here was calculated based on ′ T ms assuming that contribution of salinity anomaly to ρ′ ms is negligible 32 . The eddy energy (EE, i.e. EKE + EPE) was then integrated over the eddy: vertically, from 500 m to the surface; horizontally, from the center to the edge of the eddy. Here, the eddy center at each depth is defined as the point with the minimum EKE or maximum EPE; the eddy edge here is defined as the points where EE is 5% of its maximum value. The 500 m was chosen as the lower limit of integration because the temperature observations below this depth are very scarce for AE1 (Supplementary Table S1). For AE0, the integrated EE over the upper 500 m accounts for more than 95% of the total EE, and as such, we believe that this choice does not affect the conclusions. In the Lagrangian framework, EE of AE1 minus that of AE0 equals to time integral of the rhs terms. Given that there were no observations between the two mooring arrays, we assumed that the rhs terms had a linear evolution from AE0 to AE1. As a result, time integral of the rhs terms equals to the averaged values of AE0 and AE1 multiplied by the time for AE to move from AE0 to AE1 (~60 days).
To quantify the rhs terms of AE0 and AE1, we utilized their 3D structures constructed using the available data, and integrated them over the eddy entity. In the calculation a 3 is the density of air at sea surface, C ds is the drag coefficient as a function of wind speed 56 , v w is the 10-m wind velocity of the ASCAT wind product 57 , and ′ v top is the mesoscale sea surface velocity derived from the altimeter data. The bottom drag stress τ b is formulated as , and ′ v b is the mesoscale velocity measured by the deepest RCMs (50-200 m above the sea floor). The horizontal density gradient in the BC term is approximated by the thermal wind relation where v is the mean velocity. In the BT term, there exists horizontal gradient of mean velocity, which cannot be determined by the available data except at station M6 (using centered difference). The calculated BT at M6 was one order of magnitude smaller than the BC term and it therefore was assumed to be negligible. Because the streamlines of geostrophic oceanic eddies are generally along the P′ (pressure anomaly) contours, the P d term is commonly small and will also be neglected. In the quantification of D v term, the mean profiles of A v and K v calculated based on the observed turbulent kinetic energy dissipation rate ( Supplementary Fig. 3) were used 58 . Because of the unknown parameter A h , the D h term cannot be directly calculated. Here, we estimated this term in an indirect way as below. Physically, D h is the energy dissipation caused by submesoscale motions that act on horizontal velocity shear of mesoscale eddies, and as a result, it represents the energy mesoscale eddies transfer downscale to the submesoscale motions. Therefore, this term can be estimated if we know how much submesoscale energy was generated by the AE over the period from AE0 to AE1. As seen from Fig. 4c (also Supplementary  Figs 8a and 9), in the presence of AE, submesoscale motions are prominently enhanced and the submesoscale eddy energy within the AE is more than double that outside of the AE. Based on the instability theory of the Eady model 50,59 , the observed mixed-layer mean stratification (N = 3.4 × 10 −3 s −1 ) and geostrophic shear (Λ = 1.2 × 10 −3 s −1 ) within and around the AE give rise to a mean growth rate (0.3fΛ /N, f is the local inertial frequency) of 1/2.2 day −1 for the most unstable mode of submesoscale motions in the mixed layer. If we assume that the increased submesoscale energy within the AE all derives from AE and that the submesoscale motions have a mean growth rate of 1/2.2 day −1 , it would result in a total of 11.5 × 10 14 J eddy energy transferring downscale from the AE to the submesoscale motions over the period from AE0 to AE1 (~60 days). Here, the calculation formula is  Fig. 4c, and integrated over the depth) and dt is 60 days. Based on this D h value, the horizontal eddy viscosity can be estimated with