**Deep Flow Variability O**ff**shore South-West Svalbard (Fram Strait)**

**Manuel Bensi 1,\*, Vedrana Kovaˇcevi´c 1, Leonardo Langone 2, Stefano Aliani 2, Laura Ursella 1, Ilona Goszczko 3, Thomas Soltwedel 4, Ragnheid Skogseth 5, Frank Nilsen 5, Davide Deponte 1, Paolo Mansutti 1, Roberto Laterza 1, Michele Rebesco 1, Leonardo Rui 1, Renata Giulia Lucchi 1, Anna Wåhlin 6, Angelo Viola 7, Agnieszka Beszczynska-Möller <sup>3</sup> and Angelo Rubino <sup>8</sup>**


Received: 22 February 2019; Accepted: 29 March 2019; Published: 2 April 2019

**Abstract:** Water mass generation and mixing in the eastern Fram Strait are strongly influenced by the interaction between Atlantic and Arctic waters and by the local atmospheric forcing, which produce dense water that substantially contributes to maintaining the global thermohaline circulation. The West Spitsbergen margin is an ideal area to study such processes. Hence, in order to investigate the deep flow variability on short-term, seasonal, and multiannual timescales, two moorings were deployed at ~1040 m depth on the southwest Spitsbergen continental slope. We present and discuss time series data collected between June 2014 and June 2016. They reveal thermohaline and current fluctuations that were largest from October to April, when the deep layer, typically occupied by Norwegian Sea Deep Water, was perturbed by sporadic intrusions of warmer, saltier, and less dense water. Surprisingly, the observed anomalies occurred quasi-simultaneously at both sites, despite their distance (~170 km). We argue that these anomalies may arise mainly by the effect of topographically trapped waves excited and modulated by atmospheric forcing. Propagation of internal waves causes a change in the vertical distribution of the Atlantic water, which can reach deep layers. During such events, strong currents typically precede thermohaline variations without significant changes in turbidity. However, turbidity increases during April–June in concomitance with enhanced downslope currents. Since prolonged injections of warm water within the deep layer could lead to a progressive reduction of the density of the abyssal water moving toward the Arctic Ocean, understanding the interplay between shelf, slope, and deep waters along the west Spitsbergen margin could be crucial for making projections on future changes in the global thermohaline circulation.

**Keywords:** Fram Strait; deep sea thermohaline variability; slope currents; wind-induced processes; shelf-slope dynamics

#### **1. Introduction**

Water masses in the eastern Fram Strait, strongly influenced by the interaction between Atlantic and Arctic waters and by local atmospheric forcing, substantially contribute to drive the global thermohaline circulation [1–4]. There is a remarkable variability in the system due to several forcing mechanisms (e.g., atmospheric, internal, tidal, shelf dynamics) that play an important role, especially in the upper layer [5–10]. On the contrary, it is not completely clear which processes are responsible for the inter-annual and seasonal deep flow variability in the western offshore Spitsbergen region [11,12]. Several studies, using both experimental and numerical modelling approaches, have addressed the role of interactions between Atlantic Water (AW) carried by the West Spitsbergen Current (WSC, [10,11,13–16]), and shelf and fjord waters [2,17–20] in the observed variability. Some of the processes that may be relevant for deep water circulation and variability in the western Spitsbergen region have been summarized in Figure 1.

**Figure 1.** (**a**) Map of the study region showing bathymetry and main currents in the Fram Strait and along the west Spitsbergen margin. Red dots indicate the location of moorings S1 and ID2. Blue dots indicate CTD (conductivity-temperature-depth) stations along transects S and P. (**b**) Schematic of the shelf-slope dynamics along the west Spitsbergen slope (figure modified from [6], © American Meteorological Society. Used with permission). (**c**) S1 and ID2 moorings configuration and specification of instruments. (AW = Atlantic Water; NSDW = Norwegian Sea Deep Water; WSC = West Spitsbergen Current; BSW = Brine-enriched Shelf Water; EGC = East Greenland Current; NwAC = Norwegian Atlantic Current).

With this study, we aim at exploring signals observed in time series of temperature, salinity, current velocity, and turbidity, acquired by two near-bottom moorings (S1, ID2; Figure 1) deployed to assess the deep flow variability on short-term, seasonal, and multi-annual time scales in the southwest region offshore the Spitsbergen margin [12]. We also discuss potential links between the local oceanographic and atmospheric forcing in regulating shelf/deep sea interactions, which in turn can trigger baroclinic (internal) waves that stimulate internal mixing in the ocean (Figure 1b). Finally, we put the observed variability in the context of larger-scale circulation changes.

The overall circulation in the eastern Fram Strait (Figure 1) includes the northward flow of the WSC, i.e., the northernmost extension of the Norwegian Atlantic Current [21,22]. Flowing at a steady pace of about 0.25 m s<sup>−</sup>1, following the 1000-m isobath [23], the WSC transports warm and saline AW into the central Arctic Ocean and, beneath it, colder and fresher Norwegian Sea Deep Water (NSDW, [15,24]) that occupies the deep local layer (Figure 1b). The NSDW is influenced by water contributions from the Greenland Sea (south of Fram Strait) and Eurasian Basin (north of Fram Strait, [15,24]).

Shoreward of the WSC, the Spitsbergen continental shelf is influenced by Arctic waters as well as by drifting sea ice [19]. The northward WSC is compensated by the southward East Greenland Current, which transports cold and fresh Polar Water across the western part of the Fram Strait. The WSC is topographically steered along the continental slope [2,25] with streamlines parallel to contours of constant f/H (Coriolis parameter/water column thickness). It has been traditionally considered as a barotropic flow [13,21]. However, more recent studies have outlined its baroclinic component, together with associated instabilities and eddy formations [14,26]. The AW occupies a large portion of the water column within the WSC (roughly between the surface and 600 m depth; [10]) and its properties undergo a strong interannual variability [27].

Relatively warm and cold periods of AW have been alternating in the last century: cold periods took place during 1900–1920 and 1955–1985, while a warm period occurred in the 1930s–1940s [27]. More recently, two warm AW anomalies passing through the Fram Strait occurred in 1999–2000 and 2005–2007 [22]. Temperature increase (0.06 ◦C year−1, [22]) was accompanied by salinity increase (0.003 year<sup>−</sup>1, [28]). However, based on monthly gridded fields of the meridional current component, no significant trend in the volume transport of AW was observed [22]. A warming of the Arctic Ocean, recorded since 1985 [27,29] became particularly evident in the WSC core after 2004, along with an intensification of the northward propagation of the AW warm signal [30]. Concomitantly with this evolution, the AW influence on temperature and salinity on the West Spitsbergen Shelf and in fjords of the Archipelago (e.g., Hornsund and Kongsfjorden, [20,31]) has become stronger in the last years, probably due to changes in the patterns of the dominating large atmospheric pressure fields [2]. It is not entirely clear to what extent the increased heat transport toward the Arctic is related to a strengthening of the Atlantic Meridional Overturning Circulation [32], to an increase in temperature or in volume of the AW [4,33], or, instead, to the variability of the AW transport along the two preferential pathways (Barents Sea and Fram Strait branches, [33]). Notably, the AW heat transport [29,34] can affect air temperature especially during winter [35], which in turn has direct effects on the dense water formation around the Spitsbergen margin.

In the west Spitsbergen region ocean-atmosphere interactions lead to multiple oceanographic processes, like shelf-slope dynamics, deep water variability through Polar and Atlantic waters interaction, as well as sea ice and dense water formation [1–3,7]. Dense water formation off the southwestern tip of Spitsbergen and in the Fram Strait depends on the rate of cooling (heat loss to the atmosphere) and homogenisation of the upper layers (i.e., water column stratification and mixing depth), sea ice growth, and brine rejection. Brine-enriched Shelf Waters (BSWs) are formed during ice formation in coastal polynyas [36–40]. Within fjords, and in particular in Storfjorden (the largest fjord in the Svalbard Archipelago, see Figure 1a) the accumulated water close to the freezing point has salinity ranging typically between 34.8 and 35.8 [38]. A 120 m deep sill separates this fjord from the shelf edge, and dense water overspills the fjord with strong inter-annual variability [38]. Eventually, it flows northwards along the shelf and the continental slope west of Spitsbergen, at depths where water of similar density is transported by the WSC [41,42]. Numerical models [17,42] simulated the evolution in time and space of the plume in this region, revealing that it can reach velocities as high as 0.6 m s−<sup>1</sup> over the continental slope.

Dense water spreading is also strongly influenced by tidally induced dispersion [8,9] and by bathymetric constraints [43]. Tides can drive barotropic and baroclinic water exchanges [44], and also generate shelf waves, which propagate as topographically trapped internal waves [13,45] stimulating mixing in the ocean interior. Inall et al. [44] reported about coastal trapped waves in Kongsfjorden, a western Spitsbergen fjord at latitudes (79◦ N) where the internal Rossby radius is small and the effects of trapped waves become more evident. Increased bottom shear velocities caused by the interaction of semidiurnal internal tides with the sloping bottom can cause sediment resuspension and prevent deposition, shaping thus the continental slope [46]. Topographically induced intrusions are also responsible for the generation of nepheloid layers (i.e., turbid layers with significant amounts of suspended sediment) transported along isopycnals [16,47]. Numerical models support geological observations showing that suspended sediment transported by bottom-arrested gravity plumes may have a role in the slope convection, by increasing the kinetic energy and the bulk density of the water (Figure 1b, [5,6]). Along their routes toward the abyss, offshore southwestern Svalbard, temperature of water plumes increases, while salinity decreases rather slowly due to entrainment of surrounding warm, but relatively saline AW [41,42,48,49]. The mixing and entrainment of these plumes with AW can reduce the maximum reached depth during their cascading [10,17,41,48]. According to Akimova et al. [48], BSW plumes leaving Storfjorden reached the bottom of the Fram Strait only three times, in 1986, 1988, and 2002, over a period from 1984 to 2003. Eddies can also break away from density-driven currents, as demonstrated in laboratory experiments and numerical models [42,50]. Intruding dense water can generate cyclones in the water above and anticyclones in the water below the depth of intrusion [43].

#### **2. Data and Methods**

#### *2.1. Oceanographic Moorings*

Various instruments attached to two oceanographic moorings (S1, ID2; Figure 1) collected data within a ~150 m thick near-bottom layer along the west Spitsbergen continental slope from June 2014 until June 2016. The payload of mooring S1 included a downward looking ADCP (Acoustic Doppler Current Profiler, Teledyne RD Instruments, Poway, CA, USA) RDI 150 kHz located ~140 m above the sea bottom, a CTD (Conductivity-Temperature-Depth) SBE16 (Sea-Bird Electronics, Bellevue, WA, USA) with Seapoint turbidity meter (Seapoint Sensors, Inc., Exeter, NH, USA) coupled with a sediment trap McLane (PARFLUX Mark 78H-21, McLane Res. Labs, East Falmouth, MA, USA) ~25 m above the sea bottom, and a single point Aanderaa current meter (Aanderaa Data Instruments, Bergen, Norway) RCM8 ~20 m above the sea bottom. Turbidity expressed as Formazin Turbidity Unit (FTU) was calibrated in the laboratory to obtain the corresponding values of suspended sediment concentration (mg L<sup>−</sup>1). Mooring ID2 was equipped with two Aanderaa current meters (RCM9 and RCM4, ~120 m and ~20 m above the sea bottom, respectively) and two CT SBE37-MicroCAT recorders below each current meter (substituted with T loggers SBE56 in June 2015). A complete scheme of each mooring line (position, bottom depth, instrument type, deployment depth, and sampling time interval) is presented in Figure 1c.

The recorded data were post-processed after each recovery (June 2015, June 2016). Data with different sampling frequencies (15 min, 30 min, or 120 min) were averaged or interpolated to obtain homogenous hourly time series. Temperature, salinity, and turbidity data have been cleaned and despiked according to MyOcean in situ quality control standards and methodology (http://catalogue.myocean.eu.org/static/resources/user\_manual/myocean/QUID\_ INSITU\_TS\_OBSERVATIONS-v1.0.pdf). Moreover, temperature and salinity time series were checked by comparing them with data from CTD casts performed before and after each mooring deployment. Potential temperature (θ), salinity (S), and potential density anomaly (σθ) were calculated from in situ data. Note that only in situ temperature was recorded at ID2 between June 2015 and 2016, when thermistors (without conductivity sensors) were used, hence it was not possible to obtain θ, S, and σθ.

Each vertical cell of the ADCP was treated as an individual time series. Current vectors, both from ADCP and RCMs, were decomposed into u (eastward, zonal), v (northward, meridional), and w (vertical) components. Tidal constituents, obtained from harmonic analysis with a signal-to-noise ratio > 1, were subtracted from the original time series, to provide de-tided data [51]. For specific analyses, a low-pass filter with a cut-off period of 33 h was applied to the de-tided time series to remove inertial oscillations and obtain sub-inertial non-tidal flow (i.e., sub-tidal). The principal component analysis was applied to the time series to determine the direction of the major variance of the deep currents. As this direction almost coincided with the along-slope component, the coordinate system was rotated accordingly. The angles of rotation in the trigonometric system were between −34◦ and −52◦ at S1, and between −40◦ and −48◦ at ID2. By rotating the reference system two new components for the velocity, i.e., *ur* (along-slope) positive towards SE, and *vr* (cross-slope) positive towards NE, were obtained.

Correlation and cross-correlation coefficients calculated between different variables, and reported further in the text, are always within the significance level of 0.05, unless otherwise specified.

#### *2.2. Oceanographic Surveys*

CTD measurements during oceanographic cruises provided vertical profiles of temperature (T) and conductivity (C) approaching the seabed within ~ 7–10 m. θ, S, and σ<sup>θ</sup> were calculated from in situ data using MATLAB toolbox TEOS-10 (Gibbs SeaWater Oceanographic Toolbox) including thermodynamic equations of seawater (http://www.teos-10.org/software.htm). Oxygen concentrations were measured using a SBE43 sensor mounted on the CTD/Rosette Water Sampler while taking water samples at discrete depths for Winkler titrations [52]. Data were quality checked and averaged every 1 dbar with overall accuracies within ± 0.002 ◦C for T, ± 0.005 for S, and 2% of saturation for oxygen. Turbidity in the water column was detected by optical sensors mounted on the CTD probe.

CTD data used in this study come from the PREPARED cruise (r/v G.O. Sars, June 2014, PREsent and PAst flow REgime on contourite Drifts west of Spitsbergen), HH cruise (r/v Helmer Hanssen, June 2015), PS99.1 cruise (I/B Polastern, June 2016), and from oceanographic cruises carried out by the University Centre of Svalbard (UNIS) between 2014 and 2016 along section P (Figure 1). Long-term (1997–2017) yearly variability and linear trends were calculated from hydrographic data collected during AREX cruises carried out by the Institute of Oceanology Polish Academy of Sciences (IOPAN, [28]). In our study, we consider data collected along section S (~ 77◦30 N, see Figure 1) in the proximity of mooring ID2, from which we calculated the yearly values of θ, S, and σ<sup>θ</sup> averaged in the upper-intermediate (100–800 m depth) and in the deep layers (>800 m depth). Trends have been calculated using a linear regression model (Table 1). Standard error (SE) and Root Mean Squared Error (RMSE) for each statistically significant (p < 0.05) trend are also reported in Table 1. In order to focus on the variability of the AW in the upper layer, we only take into consideration CTD measurements that have temperatures compatible with those of the AW (i.e., θ > 2◦) in addition to the depth criterion [22]. Furthermore, we use the thickness of the AW layer as the weight while calculating the average θ, S, and σ<sup>θ</sup> values. Some figures presented in this paper were elaborated using the Ocean Data View software [53].


**Table 1.** Linear trends of θ, S and σ<sup>θ</sup> in the intermediate (100–800 m) and deep (>800 m) layers from CTD casts collected along section S (see Figure 1 and Figure 3). Data cover the period from 1997 until 2017. Significant linear trends (p < 0.05) are reported in red.

#### *2.3. Atmospheric Data*

To study atmosphere-ocean interactions between 2014 and 2016 we employed ERA-Interim dataset from the ECMWF (European Center for Medium-Range Weather Forecasts interim reanalysis) regularly spaced at 0.75◦ (interpolated at 0.25◦) of latitude and longitude and available at 6 h time intervals. We initially compared the ECMWF output (e.g., wind speed/direction, air temperature) with data collected at the Amundsen-Nobile Climate Change Tower (CCT, see http://www.isac.cnr.it/~{}radiclim/ CCTower/) in Ny-Ålesund. The 34 m high tower is equipped with a large set of instruments along the vertical to investigate physical properties of the atmospheric boundary layer over a long period and to study the turbulent exchange processes of energy and mass at the atmosphere-land interface. The comparison between land observations of wind speed and direction and ECMWF data revealed that both short-term and seasonal variabilities recorded at the CCT are reproduced coherently at S1 and ID2 locations (not shown). Hence, six meteorological parameters were selected from ECMWF dataset to characterize the air-sea interface: air pressure at sea-level, total cloud cover, 10 m zonal (*u*) and meridional (*v*) wind components, 2 m air and dew point temperatures, and sea surface skin temperature. The total heat flux (daily and monthly mean values) at the air-sea interface was computed taking into account four heat flux components also downloaded from ECMWF dataset: solar radiation, net longwave radiation, latent and sensible heat. The heat loss and the relative sea water temperature variation occurring during winter months were estimated using the methodological approach described in [54]. To analyse the periodicity and recurrence of intensive events, a Morlet wavelet analysis [55] was applied on the horizontal principal components both of the ECMWF wind data and of the current velocity at moorings S1 and ID2. The obtained results enabled studying non-stationary signals and investigating the dynamic response of the marine currents to the atmospheric forcing.

#### **3. Results**

#### *3.1. Thermohaline Patterns on the West Spitsbergen Margin during Oceanographic Cruises between 2014 and 2016*

Vertical distribution of hydrographic data along the zonal section P (Figure 1) on the West Spitsbergen Shelf north of Hornsund is shown in Figure 2. It refers to data collected in April 2014, June 2014, May 2015, and August 2016. Thermohaline properties on the shelf, usually occupied by Arctic Water [2,56], show a large seasonal and interannual variability. They depend on the variability of waters transported by the coastal current itself, on the seasonal variability in the contribution of fresh water from the main fjords, and on the variability of the AW inflow from offshore, whose properties and extension vary with large-scale circulation patterns in the area.

**Figure 2.** θ ( ◦C), S, σ<sup>θ</sup> (kg m<sup>−</sup>3), DO (mL L−1), along the W–E section P (see Figure 1) in April 2014 (**a**), June 2014 (**b**), May 2015 (**c**), and August 2016 (**d**) on the West Spitsbergen Shelf area. The data are from the PREPARED cruise and from several UNIS cruises.

In April 2014 (Figure 2a), the westernmost part of the section was occupied by warm and saline AW (3◦C < θ < 4◦C and S > 34.9), while its easternmost part was filled by cool and fresh water (θ < 1 ◦C and S < 34.9). AW could also be detected by its lower oxygen content with respect to Arctic Water. In June 2014 (Figure 2b), during the PREPARED cruise, a cold and fresh waterfront extended offshore. At that time, the westernmost part of the shelf was occupied by AW. The latter was also present on the continental slope, roughly down to 800 m depth, above the NSDW in the deep layer (>800 m, not shown). At depths between 100 m and 800 m, where AW was present, values of θ ~ 0–4.5◦C, S ~ 34.95–35.1, σ<sup>θ</sup> ~ 27.85–28.05 kg m<sup>−</sup>3, and dissolved oxygen (DO) ~ 6.4–6.6 mL L−<sup>1</sup> were found, while at depths larger than 800 m θ was below 0 ◦C, S ~ 34.91, σ<sup>θ</sup> ~ 28.05–28.08 kg m−3, and DO ~ 6.2–6.4 mL L−1. Light transmission (not shown) decreased when approaching the seabed, revealing high suspended matter concentrations in the bottom layer of the shelf and continental slope. In May 2015, AW retreated offshore (Figure 2c) like in August 2016 (Figure 2d) when, however, θ was everywhere larger than 3 ◦C and S values were relatively low (<34.85) in the shallow part of the section. In August 2016, the AW was confined in the western part of the section and it was characterized by the highest θ and S values (θ > 6 ◦C and S > 35.15). Overall, a large variability characterizes the temporal and spatial evolution of the observed parameters along the section [2]. In some cases sharp thermal and saline fronts tend to compensate, and no density structure was found (Figure 2b,c).

#### *3.2. Multiannual Variability of the Thermohaline Properties along the West Spitsbergen Margin*

CTD casts carried out during annual IOPAN summer cruises aboard the r/v Oceania between 1997 and 2017 help to understand the interannual variability of the thermohaline properties of the along-slope flow (i.e., WSC). Data collected yearly along section S (Figure 1) at latitude ~ 77◦30 N, close to the ID2 mooring, were spatially averaged within two layers: the upper-intermediate layer, between 100 and 800 m, where the AW is detected by the additional θ > 2 ◦C criterion, and the deep layer (>800 m), mainly occupied by NSDW (Figure 3). Linear trends were calculated for each of the three variables, θ, S and σ<sup>θ</sup> (Table 1). They revealed long-term positive trends of θ and S in the upper layer, while no significant trend was found in the σ<sup>θ</sup> (Figure 3). However, the upper layer went through relatively warm and saline periods (i.e., 2005–06 and 2011–12) as well as cold and fresher periods (i.e., 2003–04, 2008, and 2013). Noteworthy, starting from 2013 an overall θ increase was not accompanied by S increase. In the deep layer, slightly positive long-term trends in θ and S were also found, and they became more evident between 2009 and 2017, when also a small but significant trend in σ<sup>θ</sup> decrease was detected. Years characterized by higher standard deviations for S in the deep layer (i.e., 2000, 2002, 2006, 2008) may indicate the occurrence of advection/intrusion phenomena of different water masses into the layer. No significant correlations (p < 0.05) were obtained by comparing detrended time series between the upper and deep layers for each parameter. Observed trends are in agreement with results published by, e.g., Chatterjee et al. [57], where authors point out how changes in the AW properties depend on the strength of the subpolar gyre in the North Atlantic, which in turn is strongly influenced by the wind stress curl. The analysis of the multi-annual variability of the average thermohaline properties gives a general framework in which we undertake a more detailed deep flow variability analysis (including high-frequency signals) presented in the following section. The trends exposed here also provide an indication that the abyssal layer west of Svalbard is slowly becoming warmer and saltier on long temporal scale.

**Figure 3.** Weighted average of thermohaline properties (θ, S, σθ) calculated from CTD stations along the section S at ~ 77◦30 N (between 8 and 10◦30 E, see Figure 1): upper layer (**a**–**c**) (100–800 m with θ > 2 ◦C) and deep layer (**d**–**f**) (>800 m). Weights are defined by the layer thickness at each station where bottom depth varied between 1000 and 2300 m. Data are from the hydrographic measurements carried out during summer AREX cruises aboard r/v OCEANIA (IOPAN).

#### *3.3. Temporal Variability in Oceanographic Parameters at Moorings S1 and ID2*

Time series data recorded at moorings S1 and ID2, at depths >900 m, highlighted the presence of water with properties (<sup>θ</sup> ~ <sup>−</sup>0.90 ◦C, S ~ 34.91, and <sup>σ</sup><sup>θ</sup> ~ 28.07 kg m−3) typical of the NSDW [46]. However, perturbations caused by occasional "thermohaline intrusions" of water, which is warmer (up to ~2 ◦C), saltier (up to ~35), and less dense (down to 27.98 kg m−3), are evident in Figures 4 and 5. Noteworthy, temporal fluctuations of θ, S, and σ<sup>θ</sup> were particularly apparent between October and April and occurred almost simultaneously at both moorings, ~170 km apart from each other (see Figure 1a). A larger variability was recorded during winter 2014–15 at both sites. The duration of the observed intrusions of relatively warm and salty water varied from a few hours to several days (up to 10–15 days), during which thermohaline variations oscillated mainly at a diurnal (24 h) frequency.

**Figure 4.** Time series (hourly data) of θ (**a**), S (**b**), and σ<sup>θ</sup> (**c**) at S1 (1017 m depth). Green squares show data extracted from CTD casts (θ and S) at depths close to the moored instrument. Panel (**d**) shows normalized values of current magnitude (1022 m) and turbidity (1017 m). Data span from June 2014 to June 2016.

The normalized time series of current magnitude and turbidity at S1 are shown in Figure 4d (normalization was obtained by subtracting the average of each deployment phase and dividing by the corresponding standard deviation). Turbidity peaks were observed generally from December to July. During the first year of measurements, high turbidity values were found around mid-May 2015, while in 2016 several episodes characterized by high turbidity were observed almost each month (Figure 4d).

Overall, data recorded at mooring ID2 show a larger variance with respect to S1 data, except for θ at 1025 m depth (Figure 5). During the first year of measurements, events associated with temperature increase occurred mainly between October and April (Figure 5), while during the second year they were limited to the period October–February (Figure 6).

**Figure 5.** Time series (hourly data) of θ (**a**), S (**b**), and σ<sup>θ</sup> (**c**) at ID2 between June 2014 and June 2015 at 922 m and 1025 m depth. Green and blue squares indicate data extracted from CTD casts (θ and S) at depths close to the moored instruments. Note that after June 2015, only temperature sensors were deployed at ID2 (see in situ temperature in Figure 6).

**Figure 6.** In situ temperature time series at S1 (red) and ID2 (blue and green). Data cover the period from June 2014 to June 2016 (for the position of the moorings refer to Figure 1).

The correlation between θID2 at two depths (1025 m and 922 m) is rather high (0.75), while the correlation between SID2 at the two depths is smaller, 0.28. Short-term oscillations at 922 m were larger than at 1025 m (Figure 5). Moreover, the difference between S at 1025 m and S at 922 m diminished from ~0.01 to almost zero after one year. Accordingly, the difference in σ<sup>θ</sup> progressively diminished as well (Figure 5c). This would be consistent with the effect of mixing in the ~100 m thick near-bottom layer.

A choice to bring together and compare the T (in situ temperature) time series at S1 and ID2 (Figure 6) rises from the fact that only temperature sensors were deployed at ID2 after June 2015. The lack of conductivity sensors prevented calculation of θ, S, and hence, of σ<sup>θ</sup> as well. As said, T peaks appeared almost concomitant at the two stations. A significant cross-correlation, ~0.43–0.45, between T recorded at S1 and ID2 was found with the maximum value corresponding to a time lag of

5 h. Similarly, but not shown in Figure 6, S peaks were almost simultaneous at the two stations, despite small cross-correlations (0.06 between S1 and ID2 at 1025 m; 0.17 between S1 and ID2 at 922 m).

In general, large current velocities were observed during the winter period. To study more in detail the possible linkage between the ocean currents and the thermohaline variability in the deep layers at S1 and ID2, we focus on the most energetic period, 1 October 2014–30 April 2015 (Figure 7). A similar variability (but less energetic) was observed also in the period October 2015–April 2016.

**Figure 7.** Time series at S1 and ID2. (**a**) Stick diagram of deep currents; (**b**) current magnitude; (**c**) current component along the direction of maximum variance at S1 (1017 m depth, *ur* positive towards SE); the colorbar refers to vertical velocity (*w*) from the ADCP at the same depth. The angle of maximum variance is referred to the trigonometric system; (**d**,**e**) θ and S at 1017 m (S1) and 1025 m (ID2) depth. Temperature (T) from the ADCP at 909 m and turbidity are also shown (green lines in panels d and e). Data refer to the period 1 October 2014–30 April 2015.

The principal components determined from current data revealed that the main direction of variability was along isobaths. In this frame, the appropriate velocity components are the along-slope (*ur*, positive towards SE) and the cross-slope one (*vr*, positive towards NE, i.e., toward the coast). Major events, mainly associated with anticlockwise (cyclonic) rotations, occurred at the beginning of January, beginning of March, and at the end of March 2015. They could be attributed to the passage of eddies at the mooring location. Sometimes these eddies were detected both at S1 and ID2, but with a certain time delay. This fact was particularly evident at the beginning of January 2015 (Figure 7a), when a cyclonic eddy lasted more than 10 days at S1, and it appeared at ID2 after approximately 30 hours. Strong (>40 cm s<sup>−</sup>1) positive along-slope currents, and low (<10 cm s<sup>−</sup>1) positive cross-slope currents (Figure 7c) characterized the first phase of this eddy passage, together with enhanced positive vertical velocities (upward). The second phase of the eddy passage was instead characterized by negative along-slope currents and slightly negative vertical velocities. In general, the eddy was not associated with significant θ and S variations. On the contrary, main episodic events with enhanced θ and S values (Figure 7d,e) were associated with negative vertical velocity (i.e., downward flow, correlation was 0.13), but with relatively small horizontal velocities (Figure 7b–d). Such episodic increases in θ and S seem to occur coherently throughout the 100 m thick bottom layer. The correlation between *ur* (*vr*) and θ was small and around −0.14 (−0.14). However, the cross-correlation between θ and *vr* reached a maximum value (−0.27) with a time lag of 14 h (θ lags *vr*), probably suggesting that θ increases could be associated, to some extent and delayed in time, with enhanced negative cross-slope currents (i.e.,

directed offshore). It is noteworthy that the correlation between turbidity and *vr* (*ur*) at S1 was not statistically significant considering the period June 2014–June 2015, but it resulted significant, –0.41 (–0.13), by restricting the analysis to the period April–June 2015. These facts suggest that the seasonal increase of turbidity (Figure 4c) was mainly related to a major availability of sediment transported by cross-slope currents from the shelf to the deep sea during spring.

Progressive vector diagrams (Figure 8) show that sub-tidal currents at S1 were directed prevalently northwards, following, as expected, the bathymetric constraints. However, the prevalent current direction slightly changed with increasing depth: it tended to rotate clockwise approaching the seabed (Figure 8a). Periodical current reversals (with the current direction changing from NW to SE) were also evident with a periodicity of about 15 days, usually accompanied by stronger currents, especially during the winter season. Currents, in particular those at 1022 m depth (Figure 8b), were prevalently NE oriented between June 2014 and December 2015 (mean u = 0.5 cm s<sup>−</sup>1, mean v = 4 cm s<sup>−</sup>1), while they were prevalently NW oriented between January 2016 and June 2016 (mean <sup>u</sup> <sup>=</sup> <sup>−</sup>0.3 cm s<sup>−</sup>1, mean v = 3.6 cm s<sup>−</sup>1). In this last period, the prevalent direction of the flow resulted coherent throughout the deep layer at different depths (not shown) with only a very slight clockwise rotation with increasing depth. Moreover, less current reversals occurred, especially at the deepest measured level. This fact suggests the temporary presence of a more homogeneous deep flow during the last phase of measurements.

**Figure 8.** Progressive vector diagram of sub-tidal flow at S1. Panel (**a**): data from ADCP at different depths (909 m, 944 m, and 1014 m) and from current meter RCM8 at 1022 m within the period 1 October 2014–30 April 2015. The colorbar refers to θ at 1017 m depth. Panel (**b**): data from RCM8 at 1022 m within the whole 2-year period (June 2014–June 2016). Black dots correspond to the beginning of each month.

At ID2, the main direction of the sub-tidal currents followed the isobaths (not shown), like those observed at S1. It was prevalently northwestwards (mean u = <sup>−</sup>6.3 cm s−1, mean v = 7.6 cm s−<sup>1</sup> at 921 m, and mean u <sup>=</sup> <sup>−</sup>9.6 cm s<sup>−</sup>1, mean v <sup>=</sup> 7 cm s−<sup>1</sup> at 1024 m). The direction of the current changed slightly with increasing depth, but in the opposite way than at S1: indeed, while the main flow at 921 m depth was NW, it progressively became W–NW approaching the seabed, hence rotating anti-clockwise. Moreover, from April 2016 to June 2016, at the deeper level (1024 m) the orientation of the main flow changed from NW to N–NE.

#### *3.4. Local Wind Variability and Dynamic Response of the Ocean*

To investigate the relation between atmospheric forcing and dynamic response of the deep layer, we compared the ECMWF time series with oceanographic data recorded at S1 and ID2. In general, in both years of measurements air temperature values decreased below 0 ◦C starting from September–October (Figure 9), with negative heat fluxes at the air-sea interface (heat loss from the ocean) until April. Maximum daily heat losses around <sup>−</sup>640 W m−<sup>2</sup> occurred between December and March during intense storms characterized by strong (>10–15 m s<sup>−</sup>1) northeasterly winds (not shown). Wintertime (December–March) monthly mean heat losses were between <sup>−</sup>200 and <sup>−</sup>260 W m−2. Starting from May, most of the area southwest of Svalbard gained heat from the atmosphere (peak daily net heat flux values up to 200 W m<sup>−</sup>2). Heat losses were larger during winter 2014–15 when a larger thermohaline and current variability at S1 and ID2 was observed as well, compared to the winter 2015–16. The correlation between ECMWF wind at S1 and ID2, using data smoothed with a daily moving average, is large (>0.84). It confirms the synopticity of meteorological events in the study region. Moreover, the correlation increases up to 0.91 using data smoothed with a seven-days moving average, which still reproduces accurately the features of the strongest events. Overall, by comparing ECMWF wind and currents at S1 and ID2 a coherent variability emerges (Figure 9). Indeed, we observed intensification of the deep currents during winter months when the wind was stronger. Cross-correlations revealed that current peaks lag those of wind by about one day at S1 (0.13), and two days at ID2 (0.08).

**Figure 9.** ECMWF wind speed (m s<sup>−</sup>1) and deep current (cm s−1) time series at S1 (**a**). Data are smoothed with a daily moving average. Air temperature (◦C) and net heat fluxes (W m<sup>−</sup>2) at the air-sea interface obtained from ECMWF ERA-Interim dataset (**b**).

In order to explore the temporal evolution of the energetic events at different time scales, the wavelet analysis was applied on the ECMWF wind and horizontal current decomposed into along-slope and cross-slope principal components. As far as the wind is concerned, there were no prominent differences between outcomes at the two mooring sites, hence, for simplicity, we present only results for S1. We found a pronounced intra-annual variation with the most energetic events occurring during winter months and characterized by periodicities around 500 h (~20 days). In general, the wind in the summer period (i.e., June–August) was more quiescent (Figure 9) and the energy concentrated at a short time scale was much lower with respect to the energy concentrated at longer time scales (Figure 10e,f).

**Figure 10.** Wavelet power spectrum of the along-slope (**a**,**c**,**e**) and cross-slope (**b**,**d**,**f**) components (obtained from principal component analysis) of deep currents at ID2 and S1, and ECMWF wind data at S1 (not filtered). Periodicities of 12 and 24 h are indicated by black dashed lines.

The wavelet power spectrum of the horizontal currents at S1 and ID2 (Figure 10a–d) shows that the along-slope component was characterized by larger energy than the cross-slope one. However, the time evolution of the energetic events associated with the cross-slope current component (Figure 10b,d) reveals that major events at S1 occurred more or less regularly all over the study period, while at ID2 higher energy was found during the second year, with a peak between April and June 2016.

Referring to both components, most of the variance (i.e., kinetic energy) was concentrated in the low-frequency (long period) domain with typical periodicity centred around 500 h (~20 days) and ranging roughly between 200 and 1000 h (i.e., 1–6 weeks), similarly to what observed for the wind. Part of the energy was also distributed in the time domain between 50 and 200 h (i.e., 2–8 days), as well as in the high-frequency domain (~12, 24 h). As for the energy concentrated in the semidiurnal (12 h) and diurnal (24 h) time scales, which contain tidal signals, we found the former negligible compared to the latter. Noteworthy, the energy concentrated at the diurnal time scale was comparable for both components, although it was larger at S1 than at ID2, where it seems shifted to about 35 h. Moreover,

we found that the energy at the diurnal time scale in the cross-slope current component, in addition to being larger at S1 than at ID2, decreased with increasing depth (not shown). Another peculiarity of this signal is that it was not distributed regularly over time, contrary to what we may expect for tidally induced oscillations. In fact, there were periods of the year when the diurnal signal appeared more energetic, mostly during winter, while there were other periods (e.g., August–September 2014) when it was less pronounced or almost absent.

Summarizing, we found a correspondence between wind and current energetic events, especially at longer time scales. The energy induced by the wind and transferred to the water column can trigger fluctuations in the current field at different time scales, even at high frequencies. However, the fact that high energy over broad time scales appear also in periods when the wind is more quiescent, like during summer, suggests that the deep flow is influenced by various concomitant phenomena.

#### **4. Discussion and Conclusions**

Our measurements revealed the presence of vigorous deep current activity on the continental slope southwest of the Svalbard region. Several times, at depths of ~1000 m, we observed along-slope jets exceeding 40 cm s<sup>−</sup>1, with highest values reaching ~60–70 cm s−<sup>1</sup> (Figures 7 and 9). The prevalent flow was directed northward (from NE to NW), but periodical inversions of the flow were observed (Figure 8). In particular, the period December–March was the most energetic, as it was characterized by current magnitudes almost constantly larger than 30 cm s<sup>−</sup>1. Thermohaline variations, with large peaks of temperature and salinity, were also observed mainly between October and April. They appeared almost concomitant at the two stations, which are ~170 km apart along the mean pathway of the WSC. In fact, the cross-correlation between temperature at S1 and ID2 revealed a time lag of about 5 h only. However, water parcels travelling with an average velocity of 30 cm s−<sup>1</sup> would cover the distance between the two moorings in little more than one week. From our analyses, a strong eddy activity at both moorings was also evident. In particular, an energetic eddy that lasted more than 10 days at S1 in January 2015, was observed after about 30 h also at ID2 (Figure 7). Moreover, fluctuations of the thermohaline properties, as well as of the currents, included frequently diurnal time scale. According to these findings, we can exclude that such fluctuations observed at the two sites are part of localized phenomena transported along the west Spitsbergen slope by the WSC. In addition, we can exclude the presence of one large single eddy at both mooring locations since the typical diameter of mesoscale eddies in this polar region is around 10–20 km [58]. On the other hand, eddies may be a manifestation of the mesoscale variability that propagates along the shelf break and along the continental slope, and they are able to induce vertical water displacement. Similarly, the passage of trains of topographically trapped waves, with a predominantly diurnal periodicity, travelling along the continental slope may result in internal oscillations, highlighted by the observed thermohaline variability. This interpretation is in agreement with Nilsen et al. [45], whose findings revealed the presence of both mesoscale variability and topographically trapped waves at shallower depths (shelf break) in a wider zone including our study area, resulting from rotating wind field during passing storms. In particular, beside the mesoscale variability they found internal waves with wavelength λ = 32 km for periods of 24 h and 35 h. The passage of trains of internal waves would partly explain the thermohaline oscillations observed at the mooring sites since these waves can induce changes in the vertical distribution of the AW along the entire west Svalbard margin. A similar phenomenon was also observed in numerical simulations along the Norwegian continental shelf [59]. Lien et al. [59] found that general intensifications of North Atlantic Oscillation and consequent atmospheric events, are able to amplify the AW depth variability during winter, due to the Ekman cross-shore transport induced by the along-slope (i.e., along-shore) wind component. This fact seems in agreement with the enhanced winter energy found from the wavelet analysis applied to our data. Indeed, we found an overall consistency between the variability of the deep sea currents and the local wind signal (see Figure 9), especially on seasonal and low-frequency scales. In particular, we observed the occurrence of the most energetic events during winter both for the deep current flow and for the wind speed with typical

periodicities of ~20 days. Consequently, our interpretation is that the passage of high-low atmospheric pressure systems can explain the origins and amplification of internal oscillations observed along the west Spitsbergen margin at 1000 m depth. Meteorological forcing would hence modulate the current flow, extending its influence from the surface down to deep layers. Von Appen et al. [15] also found variability of deep ocean flows with periodicities of 1–3 weeks in the northern Fram Strait, at depths exceeding 2000 m. They explained such variability by the passage of basin-scale topographic Rossby waves excited by synoptic atmospheric forcing.

Sanchez-Vidal et al. [11] observed similar sporadic intrusions of relatively warm and salty water at their Stations A (close to S1) and B (1500 m depth, slightly deeper on the continental slope) in 2010–11. They suggested that such modifications of the ambient water at large depths were likely related to the influence of shallower AW. Moreover, they also found a 100 m thick bottom layer characterized by relatively high turbidity values, likely associable to a nepheloid layer. Increased turbidity along the west Spitsbergen continental slope can result from resuspension of sediment due to enhanced bottom currents (i.e., passage of trains of internal waves), but also from particles transported during cascading events, for example. According to this, we suppose that observed water intrusions at depths larger than 1000 m have two possible causes: (i) they originate from internal oscillations triggered by the passage of trains of internal waves or eddies, as discussed above; (ii) they occasionally originate from slope currents as a result of dense water formation and vertical mixing events occurring on the shelf, the latter being strongly influenced by the progressive intensification of the AW signal on the shelf and within fjords [20,31].

In support of the second hypothesis, we discuss here a combined analysis of meteorological and oceanographic data. Air-sea interaction and strong vorticity in the wind field lead to Ekman pumping and vertical convection over the West Spitsbergen Shelf [1,23,60]. Nilsen et al. [2] demonstrated that in the period autumn–winter (between September and May) low-pressure atmospheric systems influence the West Spitsbergen Shelf area. They also found a negative correlation with a zero time lag between along-shore wind stress and ocean temperature at their mooring I1, at the mouth of Isfjorden on the shelf, at 50 m and 190 m depths. Subsequently, northerly wind events (and associated positive wind stress curl field) could cause upwelling of AW that can be cooled after reaching the surface and become denser. The wintertime heat fluxes in our study had daily peak heat loss from the ocean to the atmosphere reaching <sup>−</sup>640 W m−2, while average monthly in the period December–March had values around <sup>−</sup>200 W m−2. Such heat losses could cool down by 2–3 ◦C the intruding warm and saline AW all the way to the bottom on the ~200 m deep continental shelf, in agreement with findings by Hakkinen and Cavalieri [61]. Cooling of this magnitude would increase the density of the AW (θ > 2 ◦C, S > 34.92, 27.70 <σ<sup>θ</sup> < 27.97 kg m<sup>−</sup>3) on the shelf by reaching values of ~28.04 kg m−3, triggering the sinking of shelf water downslope to depths > 900 m (Figure 1b). Consequently, thermohaline and current variability observed in the deep layer at S1 and ID2 could also be sporadically caused by the arrival of gravity currents driven by dense water plumes formed over the shelf during intense meteorological events (i.e., those producing shelf convection) or formed after sea ice formation within fjords (i.e., BSW). Dense plumes descending as bottom-arrested currents follow preferential routes constrained by bathymetry, undergo a strong entrainment of AW that occupies the intermediate and upper layers of the WSC [41]. These dense plumes can collect sediments (Figure 1b) that increase their kinetic energy and bulk density [6,62,63]. In this regard, Fohrmann et al. [6] pointed out how a volumetric concentration of 1000 mg L−<sup>1</sup> of suspended quartz particles increases the bulk water density of 0.6 kg m−3. They also demonstrated how turbidity plumes in the Svalbard region are 10 times faster and can reach deeper layers with respect to mere temperature-salinity plumes. However, at S1, we found turbidity values up to 6 FTU in late winter season, corresponding to a concentration of suspended sediment of about 16–18 mg L−1, which would not be sufficient to compensate the low density deriving solely from temperature and salinity. Hence, an input of sediments from the shelf/slope that would stimulate water cascading is plausible, but not so explicit from our data. In any case, at S1 we found a larger correlation between the cross-slope velocity component directed

offshore and turbidity during spring. This fact may indicate cross-slope currents transporting more sediments toward the deep layers. Since slope currents can affect ocean stratification through baroclinic instabilities, they can themselves generate oscillatory signals in the deep layer [64] and induce current reversals, which in turn, cause vertical displacements of the water as well as sediment resuspension, interacting with the continental slope. Hence, the two phenomena, i.e., internal oscillations and gravity currents, can coexist although hardly discernible from our data.

Finally, our study shows that the West Spitsbergen Shelf experiences high seasonal and interannual thermohaline variability associated with the inflows of AW and outflows of cold and fresh water from fjords. A large thermohaline variability can strongly influence the generation of dense water plumes on the shelf, as well as their properties. At the same time, a warming trend of AW emerges in agreement with previous studies [22,28,65]. Similarly, a slightly positive temperature trend appear also in the deep layer, which is consistent with the observed continuous warming of the deep waters in all sub-regions of the Nordic Seas [65].

In conclusion, our observations suggest that shelf-slope dynamics modulated by synoptic atmospheric forcing can increase the mixing rate between upper and deep layers along the west Spitsbergen continental slope, contributing to the slow modification of the deep layer (>800 m depth) along the west Svalbard margin, which is experiencing a slight tendency to become warmer and saltier. Further analyses are required to understand if prolonged injections of relatively warm water within the deep layer along the west Spitsbergen margin could potentially be responsible for future modifications of the abyssal waters in the Arctic Ocean.

**Author Contributions:** M.B. set the leading hypotheses of the work, performed data analysis, and wrote the main manuscript text. V.K., L.L., S.A. co-designed the research, participated in data analyses and interpretation and contributed to writing the text. A.R., L.U., I.G., T.S., R.S., F.N., contributed to data analyses and writing the manuscript. D.D., P.M., R.L., M.R., L.R., R.G.L., A.V., A.W., and A.B.-M. participated in data acquisition and interpretation. All authors reviewed the manuscript.

**Funding:** This research was funded by the FP7-EU/Eurofleets2 initiative through the PREPARED (PREsent and PAst flow REgime on contourite Drifts west of Spitsbergen) project and by the Italian Ministry of University and Research through the National Programme of Research in Antarctic (PNRA), project DEFROST (DEep Flow Regime Off SpiTsbergen). This study was partially supported by the Svalbard Integrated Arctic Earth Observing System (SIOS) through the first SESS report pilot call (SOA project, contract n◦2017\_0007).

**Acknowledgments:** We acknowledge Matthias Forwick (UiT), as well as Captains and Crews of the r/v Helmer Hanssen, r/v G.O. Sars, I/B Polastern, and r/v Alliance (Italian Hydrographic Institute, "High-North" programme) for moorings maintenance and assistance in the field works. We thank Daniela Accettella (OGS) for preparing the bathymetry map shown in Figure 1a. We also acknowledge the anonymous reviewers for their valuable comments.

**Conflicts of Interest:** The authors declare no competing interests.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Census of the 1993–2016 Complex Mesoscale Eddy Processes in the South China Sea**

**Huimeng Wang 1,2, Yunyan Du 1,2,\*, Fuyuan Liang 3, Yong Sun <sup>4</sup> and Jiawei Yi 1,2**


Received: 16 March 2019; Accepted: 5 June 2019; Published: 10 June 2019

**Abstract:** Mesoscale eddy process with at least one splitting and/or merging event can be defined as either a complex process or a simple process. Investigation of the difference between these two categories could provide new insights into how different factors, such as the seabed topography, Kuroshio intrusion, and winds, affect the origin, migration, and decay of the mesoscale eddies. This study compared the characteristics of the complex against the simple eddy processes in the South China Sea (SCS) from 1993 to 2016. We comprehensively analyzed the eddy processes with regards to their characteristic points, trajectories, and networks. The simple and complex processes share many similarities but do show significantly different behaviors. Both the simple and complex processes mainly start from the eastern SCS. However, the complex processes mainly vanish in the western SCS whereas the simple processes disappear almost everywhere across the SCS. The complex processes last longer and migrate more than the simple processes. Lastly, the complex processes mainly move westward within the community. The complex processes can be further categorized into complex anticyclonic and cyclonic eddy processes. Spatially, the splitting and merging events mainly occur in the southwest of Taiwan, northwest of the Luzon Island, and the southeast of Vietnam. Temporally, the merging and splitting events mainly occur in the fall. The interaction among the communities reveals the different migration patterns of the complex anticyclonic and cyclonic eddy processes in the SCS.

**Keywords:** complex processes; eddies; mobility indicators; splitting and merging; community division; South China Sea

#### **1. Introduction**

Mesoscale eddies are a common dynamic phenomenon in the ocean, which is constantly moving around and changing with respect to their geometric and thematic characteristics. Previous studies have shown that mesoscale eddies migration facilitates transport and exchange of energy and substance in the oceans, and thus plays a significant role in the marine ecosystem, atmospheric environment, and surface ocean circulation [1–5]. Mesoscale eddies even show their effects on the migration and biogeochemical processes of the deep-ocean biological communities [3–5]. The heat exchange between the eddies and the atmosphere also affects the local wind field, clouds, and precipitation [6].

Many studies have been done on mesoscale eddies in different sea areas (such as the Pacific Ocean, the Indian Ocean, the Mediterranean Sea, Peru, etc.) to develop methods for identifying and tracking eddies [7–9], to examine eddies statistical characteristics [10–15] and the driving mechanisms behind

them [16–19], and to study the influences of eddies on current circulation [2,20,21], oceanic ecology, and biogeochemical processes [5,22]. As in other waters, mesoscale eddies are also very common and active in the South China Sea (SCS), which is a marginal sea of the Pacific Ocean and also the largest semi-enclosed sea in the tropics. Different methods and multiple-source data sets have been used to examine the eddies in the SCS. For example, Yang et al. [23] and Yuan et al. [24] discovered a seasonal cyclonic (CE) and an anticyclonic eddy (AE) in the northern SCS from the climatological Levitus data and satellite altimeter data. From the Argo float data, Chow et al. [25] found perennial cold eddies in the south of the Dongsha Islands, mainly in winter or spring. Observational studies in the east of Vietnam [26] show that during the summer monsoon season, more anticyclonic eddies are found in the southern SCS whereas more cyclonic eddies are found in the northern SCS. Chen et al. [27] used in situ hydrographic data, 375 Argo profiles, and sea level anomaly (SLA) data to find relatively large current heat exchange both along the parallel and the meridian directions. Polar eddy current heat transfer was found in summer in the east of Vietnam. In the west of Luzon Island, polar heat exchange was found in winter. In the western Luzon Strait, a large equatorial heat transfer was found in winter.

Statistics has been widely used to depict the overall characteristics and movement patterns of the mesoscale eddies in different regions, from which efforts have been made to infer the driving mechanism and provide support for numerical simulation of the eddies [11–14,26]. Statistical analysis of the eddies in the SCS has been conducted from different perspectives by using the satellite altimeter data and observational data. Many studies focused on statistically analyzing the total number of eddies over a certain period and their lifespans, propagation directions, movement speeds, emerging and vanishing locales, spatiotemporal distribution patterns of their attributes, and even the correlation with the El Nino or La Nina events. For example, Wang et al. [11] examined the 1993 to 2000 sea surface height anomaly data and delineated the SCS into four regions based on the statistical analysis results of emerging locales of the eddies. By using the 1993 to 2001 TOPEX/Poseidon merged ERS-1/2 altimeter data, Lin et al. [26] showed that about 80% of eddies moved westward from the southwest of Taiwan to the east of southern Vietnam. Xiu et al. [12] conducted an eddy census in terms of the number, size, and lifespan from 1993 to 2007 and showed that the radius of the eddies ranged from 46.5 km to 223.5 km, with an average of 87.4 km. More than 70% of the eddies have a radius smaller than 100 km. They also found that the eddy activities in the SCS are not directly related to the El Niño–Southern Oscillation events. Chen et al. [13] examined the 1992 to 2009 SLA data and statistically summarized the seasonal and inter-annual variability of the eddies. They found that the eddies in the SCS mainly develop in an elongated northeast–southwest zone in the southwest of the Luzon Strait, and argued that eddy activity is sensitive to the wind stress curl in the northern SCS. Du et al. [14] investigated the general characteristics of the eddies in the SCS, particularly the spatial distribution patterns of eddy disappearing, reappearing, splitting, and merging activities.

Nowadays, mesoscale eddies are observed splitting and merging frequently thanks to the availability of high-quality spatial data, including satellite observation data, in situ observation data and the emergence of various model data [28–32]. In this study, we define a mesoscale eddy evolution process with at least one or without any splitting and merging event as a complex or simple process, respectively. It is valuable to study complex eddy processes as the splitting and merging events tend to impact seawater interaction and marine ecosystems more than the simple eddy processes [29,30]. For example, eddy splitting can cause eddy deformation, substance dispersion, and subsequent energy attenuation, affecting the temperature, salinity, water mass and microorganisms in the marine ecosystems [31]. Eddy merging may strengthen seawater interaction and the energy and entropy transfer across multiple scales in turbulent flows [29,32]. In fact, some studies have examined the complex processes and already recognized the great value of comparing simple against complex processes [28–32]. For example, Li et al. [31] argued that an eddy merging event may trigger a tropical storm, whereas a splitting event can initiate and enhance energy and material exchange among different regions. Not much statistical analysis has been conducted on the complex eddy processes in the SCS,

and it remains unclear how the spatiotemporal characteristics are different between the simple and complex eddy processes.

Various factors, such as Kuroshio intrusion, wind stress, and coastal jet have been frequently cited to explain the origin, migration, and extinction of the eddies in the SCS. Unfortunately, the driving mechanism is still pretty vague. For example, the northern SCS is where eddies mainly start. It has been attributed to the shedding of the Kuroshio intrusion [24,33], wind stress [11,34] or influence of other eddies that invaded from the Pacific Ocean through the Luzon Strait [35]. In the east of Vietnam, the unstable strong currents along the coast may trigger more births of eddies [36]. A comprehensive statistical census of the complex eddies processes, in comparison with the simple processes, is thus needed and the results may shed new light on the mechanisms that drive the splitting, merging, and even the migration patterns of the eddies in the SCS.

This paper reports a comprehensive census of the complex eddy processes in the SCS from 1993 to 2016. The aim of this study is to examine the spatiotemporal characteristics of the complex eddy processes in the SCS to better understand their spatial distribution and mobility patterns in comparison with those of the simple processes. The movement patterns of eddy processes not only reflect their evolution characteristics but also imply the spatial interaction between the different geographic locations where the eddies pass. In this study, we used a variety of Geographic Information System (GIS)-based methods and network data mining method, including kernel density, hot spot analysis, and community detection, to examine the complex eddy processes.

The kernel density analysis calculates the density of point or linear features within a neighborhood around those features. The hot spot analysis identifies where features with either high or low values cluster spatially. Both kernel density and hot spot analysis are widely used in crime, utility, disease, rainfall, disasters, and wildlife mapping [37–39]. We applied these two methods to the SCS eddies to identify where the starting, demising, splitting, and merging locales are most clustered. Identifying the clustered locales would help us understand why certain locales are clustered at certain places. Community detection is an important method of network data mining, which is widely used to analyze the regional and mobile characteristics of objects such as mobile phone users and web graph in space [40,41]. The community detection method is used to delineate the SCS into different communities. We then examine the activities of the eddies within a community and between two adjacent communities to identify interaction among different communities (regions) in the SCS.

The paper is organized as follows. Section 2 introduces the data sets and methods used in this study. Section 3 focuses on the analysis results and we also compare and discuss the characteristics of the complex and simple eddy processes in the SCS. Section 4 focuses on the influence of topography on the splitting and merging of the complex eddy processes. Finally, we summarize this study and outline our future work in Section 4.2.

#### **2. Data and Methods**

#### *2.1. Data*

This study used the SLA data from January 1993 to December 2016. The SLA dataset is widely used in studying eddies [38] and has a temporal resolution of 1 day and a spatial resolution of 1⁄4 degree [42]. We developed and improved the eddy identification and tracking methods and built a spatiotemporal database, from which every eddy in the SCS could be queried and studied for different purposes. In our dataset, there were 5773 eddy processes (3452 simple and 2321 complex processes). We provided a brief description of the data processing and data accuracy assessment results below. More details are available in Yi et al. [9,43–45].

The eddies in our database were identified using a hybrid detection (HD) method [43], which takes the advantages of the two widely-used eddy identification methods: the Okubo–Weiss (OW) [46,47] and the Sea Surface Height (SSH) [8] methods. An eddy truly exists only if it meets the OW criterion (W < −0.2σw, where σw is the spatial standard deviation of W, and W is computed from the horizontal velocity field asW = Sn<sup>2</sup> +Ss<sup>2</sup> <sup>−</sup> <sup>ω</sup>, where Sn, Ss, and <sup>ω</sup> are the straining deformation rate, the shearing deformation rate, and the vorticity, respectively) and includes either a local maximum or minimum SLA value. The identification results derived from the HD, OW, and SSH methods were compared, showing success rates of 96.6%, 96.6%, 100% and failure rates of 14.2%, 70.3%, 36.8%, respectively. Overall, the HD method outperformed the other two methods in correctly identifying eddies in the SCS [43].

The eddies identified were then concatenated using the method proposed by Yi et al. [44]. It is very common that an eddy is concatenated to a wrong predecessor or successor, particularly when an eddy splits, merges, disappears, and reappears. A global nearest neighbor filter (GNNF) approach [44], which combines the Kalman filter and optimal data association technologies to recursively recover the predecessor and successor of a specific eddy. The GNNF tracking results of the eddies in the SCS were compared with those derived from the distance-based search [48] and the overlap-based search [49] methods. The mismatching rates were 0.19%, 0.30% and 1.17% for these three methods respectively. The three methods were also used to track synthetic eddy trajectories with mismatching rates 0.2%, 0.4% and 0.5%, respectively. The comparison shows that the GNNF approach had its advantages in concatenating eddies.

#### *2.2. Methods*

In this study, an eddy process is represented as a trajectory with multiple points in chronological order (p1, p2, p3 ... ... pn). Every point along the trajectory was defined with its longitude, latitude, and a timestamp (pi = (xi, yi, ti)). The trajectory could bear with or without branches. A structure containing no branches would represent a simple eddy process and otherwise a complex process (Figure 1). Any branch along the trajectory would represent either a splitting or a merging event.

The complexity of a trajectory was first quantitatively measured by a complexity index η (0 < η ≤ 1), which is the ratio between the number of the points along the branch(es) and all points on the trajectory. Figure 1a illustrates a complex trajectory. At a specific time-stamp, the eddy represented by point p3 splits into two eddies p4 and p5, which then merge into eddy p7. At another time-stamp, eddies p8 and p9 merge into eddy p10, which then splits into eddies p11 and p12. The total number of points along this trajectory is 13 and nine of them are on the branches. Therefore, the η is 0.69. A simple eddy trajectory (Figure 1b) is a linear structure without any branches thus has a complexity index 0 (η = 0).

**Figure 1.** A complex trajectory (**a**) and a simple trajectory (**b**).

In this study, conceptually the eddy trajectories could be represented as different objects: feature points, lines, and then networks. A variety of methods were used to statistically examine these objects (Figure 2). We first conducted kernel density and hot spot analysis on the feature points. Then we used a series of activity parameters to describe the movement characteristics of the feature points. The networks, formed by the points and lines, were examined using the fast folder community detection method [41]. The sections to follow provide more details about these methods.

**Figure 2.** A diagram showing the structure of an eddy trajectory, which is represented as different objects (points, lines, and networks) (**a**), which were then examined using different methods (**b**).

#### 2.2.1. Kernel Density and Hot Spot Analysis

We first conducted spatial kernel density and hot spot analysis on the feature points, which include the starting, ending, splitting, and merging points along the trajectories of cyclonic and anticyclonic eddies. Kernel density shows the concentration of points over an interpolated and smooth surface. In the places where points are clustered, the surface would have a higher density, which then gradually diminishes with increasing distance from the point clusters. The kernel density of the starting and ending points showed the places where eddies mainly originated and disappeared. By contrast, the kernel density of the splitting and merging points revealed the locales where the eddies tended to split and merge, probably due to the influence of different ocean environmental and climatic factors.

The Gi\* method [50] was used to identify the hot and cold spots of the feature points by examining the z-score and *p*-value. A hot spot is a place where high-value features are surrounded by other high-value features. By contrast, a cold spot is a place where low-value features are clustered together. First, a 1◦ × 1◦ square fishnet was generated to cover the whole study area. We chose the 1◦ × 1◦ grid (about 108 km × 108 km in the SCS) as the radii of most eddies in the SCS are between 100 km and 200 km, and the average moving distance is 74 km [12,14]. Secondly, we counted the total number of every type of feature points in each grid. The counts and the Euclidean distance among the points were then used to calculate the Gi\* values. Lastly, statistical significance was evaluated by looking at the z-score, *p*-value, and confidence level. A hot spot was identified when a z-score was greater than 1.65 and a *p*-value less than 0.10. In contrast, a low negative z-score (<−1.65) and a small *p*-value (<0.10) indicated a cold spot clustering, i.e., few eddies nearby.

#### 2.2.2. Eddy Mobility Index

Activity space in geography is mainly used to examine human activities within a space [51,52]. We employed the same idea and computed four indexes, namely the complexity, activity duration, activity scope, and the number of activity points, to study how an eddy migrates within an activity space. The complexity, as discussed in Section 2.2, showed how complex an eddy trajectory is. The activity duration referred to the lifespan of an eddy, indicating how long an eddy process lasted. The activity scope was defined as the maximum distance between any two points along a specific track, indicating the maximum span within which an eddy was moving around. The number of activity points was the number of different positions that an eddy may occupy over its life span. More activity points suggested a more active eddy process.

#### 2.2.3. Community Detection

A geographic area could be delineated into a set of communities, each of which has similar attributes and/or functions in a network graph [53]. Many partitioning methods have been developed for community detection. This study used the fast unfolding method [54]. Eddy movement essentially is inhomogeneous as eddies are more active in certain regions across the SCS [13]. From a GIS perspective, any movement from one position to another could be represented as a connection (i.e., an edge). Movements of a large number of eddies could be translated and represented as a complex network consisting of nodes and edges. The network could then be delineated into different communities (Figure 3). Within a specific community, the eddies were similar in terms of their characteristics. The eddies' characteristics showed homogeneity within a community but heterogeneity among communities.

**Figure 3.** Community partition workflow.

The fast unfolding method includes two major steps in a sequence (Figure 3). The first step is to construct the movement network. We generated tessellated square grids of 1◦ × 1◦ to cover the whole study area and the center of each grid was represented as a node. We then counted the total number of eddies within each grid (node), respectively. The connections between any two adjacent grid nodes in the network were weighted using the total number of eddies that crossed their boundary. A network thus was constructed with the nodes representing the grids and edges the total number of eddies that crossed the boundary of two adjacent grids.

The second step was to partition the communities. Each node was regarded as a community and would be grouped with its neighboring nodes (communities) one by one. A modularity gain ΔQ was calculated as below [41].

$$
\Delta \mathbf{Q} = \left[ \frac{\sum\_{\text{in}} + \mathbf{k}\_{\text{i,in}}}{2\mathbf{m}} - \left[ \frac{\sum\_{\text{tot}} + \mathbf{k}\_{\text{i}}}{2\mathbf{m}} \right]^2 \right] - \left[ \frac{\sum\_{\text{in}}}{2\mathbf{m}} - \left[ \frac{\sum\_{\text{tot}}}{2\mathbf{m}} \right]^2 - \left[ \frac{\mathbf{k}\_{\text{i}}}{2\mathbf{m}} \right]^2 \right] \tag{1}
$$

where m represents the total connection weights across the entire network, in is the total connection weights within a specific community, tot is the sum of all connection weights associated with nodes within the specific community, ki is the sum of the weights of all the connections to node i, ki, in denotes the sum of the connection weights between nodes i and all the other nodes in the community. If the maximum ΔQ was greater than 0, the node would be grouped with its neighboring node. Otherwise, it stayed as it was. If a node was grouped into another community, a new network was constructed and the afore-mentioned procedure was repeated until no more nodes shifted to another community.

#### **3. Results and Analysis**

Over the study period from 1993 to 2016, 5773 tracks were identified and stored in our database. Among them, 2321 and 3452 were complex and simple eddy processes, respectively. The 2321 complex processes included 1120 and 1201 anticyclonic eddies (CPAEs) and cyclonic eddies (CPCEs), respectively. By contrast, 1682 and 1770 out of the 3452 simple processes were anticyclonic eddies (SPAEs) and cyclonic eddies (SPCEs), respectively.

#### *3.1. Feature Points Analysis*

#### 3.1.1. Kernel Density Analysis of the Simple and Complex Eddy Processes

The kernel density difference of the starting and ending points of the SPAEs (Figure 4a) and SPCEs (Figure 4b) showed that the eastern SCS witnesses more births than deaths of SPAEs and SPCEs, respectively. The death places of SPAEs and SPCEs were scattered across the whole SCS.

The kernel density difference of the starting and ending points of the CPAEs (Figure 4c) and CPCEs (Figure 4d) showed that more CPAEs and CPCEs were born in the eastern SCS and disappeared mainly in the western SCS. Such a finding is consistent with previous studies conducted by Wang et al. [11], Xiu et al. [12], and Chen et al. [13], in which the simple and complex eddy processes were not differentiated and studied as a whole. In this study, we examined the simple and complex eddy processes separately and found that the complex processes mainly disappeared in the western SCS whereas the simple processes disappeared almost everywhere across the SCS.

**Figure 4.** The difference of the kernel density of the starting and ending points of simple processes of anticyclonic eddies (SPAEs) (**a**), simple processes of cyclonic eddies (SPCEs) (**b**), complex processes of anticyclonic eddies (CPAEs) (**c**), and complex processes of cyclonic eddies CPCEs (**d**).

#### 3.1.2. Hot Spot Analysis of Simple and Complex Eddy Processes

Figure 5 shows the hot and cold spots of where the simple and complex eddy processes originated, respectively. The hot spots were mainly located in the eastern SCS whereas the cold spots were mainly located in the western SCS. Such a pattern indicates that the eastern SCS was where most eddies start from. However, the hot spots of the complex eddy processes were more extensive and continuous in the eastern SCS, covering the area from the southwest of Taiwan, west of Luzon Strait and Luzon, Huangyan Island, Nansha Islands, northwest of the Kalimantan. By contrast, the hot spots of simple eddy processes were more isolated in the eastern SCS and could only be found in the southwest and northwest of Luzon and Palawan, southwest of Nansha Islands (7◦ N–9◦ N, 110◦ E–113◦ E). The same pattern was also found for the cold spots. Extensive cold spots were found for the complex eddy processes in the northwestern SCS. The cold spots of simple eddy processes were isolated and found in only a couple of grids in the east of Hainan Island and the southwestern SCS.

The difference between the hot and cold spots of the simple and complex eddy processes indicates that the complex processes originated from an extensive area in the eastern SCS. Very few complex eddy processes started from a quite extensive area in the northwestern SCS. By contrast, the simple processes mainly originated from three concentrated regions in the eastern SCS as shown in Figure 5a. There are only a couple of places in the western SCS that few simple processes originated from.

Both the kernel density and hot spot analysis indicated that the southwest of Taiwan Island, the west of Luzon Strait (Luzon strait between Taiwan Island and the Luzon Island) and the northwest of Luzon Island were where significant numbers of simple and complex processes started. Previous studies have shown that the local wind stress curl and Kuroshio intrusion contributed to the generation of both cyclonic and anticyclonic eddies in this area [33,55]. This study shows that origination of complex processes in this region could also be attributed to the local wind stress curl and Kuroshio intrusion. In Section 4.1.1, we discussed and analyzed how wind stress curl and Kuroshio intrusion affect the eddy evolution processes.

**Figure 5.** The hot and cold spots of the simple (**a**) and complex eddy processes (**b**).

3.1.3. Splitting and Merging of Complex Eddy Processes

Unlike the simple eddy processes, the complex eddy processes split and merge. On average, the complex eddy processes split and merged 75 and 45 times every year in the SCS (Figure 6a). The CPAEs and CPCEs showed no significant difference between the numbers of merging and splitting events (Figure 6a). However, the numbers of splitting and merging events varied significantly in different seasons, with 908 (542), 898 (539), 1013 (612), and 809 (495) splitting (merging) events occurring in the spring, summer, fall, and winter, respectively (Figure 6b). Fall was the most favorable season that the complex eddy processes tended to split and merge.

**Figure 6.** The statistics of splitting and merging events of CPAE and CPCE by year (**a**) and season (**b**).

The kernel density difference between the splitting and merging feature points shows that CPAEs merged more in the southwest of the Taiwan Island and the east of Vietnam, whereas they split more in the southeastern SCS (Figure 7a). By contrast, the CPCEs split more in the east of Luzon Island and the southwest of Taiwan Island (Figure 7b). In the southern SCS, the CPCEs split more except southeast of Vietnam. These regions were also where the eddies were concentrated, the differences in splitting/merging events between different areas will be discussed in Section 4.1.2.

**Figure 7.** The kernel density difference between the splitting and merging feature points of CPAEs (**a**) and CPCEs (**b**).

#### *3.2. Line Analysis*

Figure 8 shows the trajectory density of all (a), complex (b), and simple eddy processes (c) in the SCS, respectively. There was no significant difference between the trajectory distribution of all and complex processes, with higher density in the southwest of Taiwan and northwest of Luzon, west of Mindoro, and east and south of Vietnam. The trajectories of simple processes also concentrated in the southwest of Taiwan and northwest of Luzon. However, more trajectories of simple processes were found across an extensive region in the southwestern SCS.

**Figure 8.** The grid-based trajectory density for all (**a**), complex (**b**), and simple (**c**) eddy processes.

The simple and complex eddy processes were also different in terms of the mobility indexes, including the activity duration, activity scope, and the number of activity points (Table 1). The complex eddy processes had a much longer activity duration (33.5 days) than that of the simple eddy processes

(8 days). Both the activity scope and the number of activity points were over three times more extensive than those of the simple eddy processes. About 90% of the simple eddy processes have less than 10 active points and the travel range was less than 150 km in a lifespan less than 28 days (Figure 9). By contrast, about 90% of complex eddy processes had less than 43 active points and the travel range was less than 503 km in a lifespan less than 123 days. More than 80% of the complex eddy processes had a complexity index between 0.3 and 0.8 and the complexity index of the simple eddy processes, of course, was 0. All the differences between the means of all three mobility indexes were statistically significant (Table 1). However, there was no significant difference between the mobility indexes of the SPAEs and SPCEs. There were significantly more CPCEs than CPAEs with a complex index of 0.5 (Figure 9).

**Table 1.** Comparison of the means between the simple and complex eddy processes. The last column shows the significance test results of the Mann–Whitney U test for the mobility indexes.

**Figure 9.** Statistics of the mobility indexes of CPAEs and CPCEs.

#### *3.3. Network Analysis*

The SCS beyond the 200-m bathymetry contour line could be divided into eight communities based on the networks built from all, simple, and complex eddy processes, respectively (Figure 10). The communities were roughly similar. Particularly, the community delineation results based on the complex- and all-eddy networks were very similar (Figure 10a,c). Roughly, the northern, central, and southern SCS could be divided into three, two, and three communities, respectively. However, difference does exist. For example, only one community was partitioned from the all-eddy networks in the southern part of the southern SCS (C8 in Figure 10a), whereas this region was divided into two simple- and complex-eddy communities, separately. A community was identified in the east of Hainan Island and Vietnam based on the simple-eddy networks (C3 in Figure 10c). This region was divided into two communities based on the complex- and all-eddy networks.

**Figure 10.** The community delineation results based on the all- (**a**), simple- (**b**), complex-eddy processes (**c**).

We then mainly focused on discussing how the SCS could be divided into different communities based on the networks constructed with the CPAEs (complex processes of anticyclonic eddies) and CPCEs (complex processes of cyclonic eddies) separately. The delineation results were very similar for the northern and southern SCS (Figure 11) but not for the central SCS. For example, the communities could be delineated for the central SCS based on the complex anticyclonic eddy networks (C4–C6 in the Figure 11a). However, no such communities were identified in the central SCS based on the networks of complex cyclonic eddy processes (Figure 11b). The complex cyclonic eddy processes-derived communities (C6–C7 in Figure 11b) in the central SCS tended to extend into either the northern or the southern SCS.

Eddy processes frequently migrated between adjacent communities. The interaction patterns of the complex anticyclonic and cyclonic eddy processes shared many similarities yet also with a significant difference. The CPAEs and CPCEs moved more frequently in the southwest of Taiwan Island (CPAEs-C1, CPCEs-C1), northwest of Luzon Island (CPAEs-C3, CPCEs-C3) and southeast of Vietnam (CPAEs-C7, C8, CPCEs-C6, C8) and much less in the southeastern SCS (CPAEs-C9, CPCEs-C7).

In the northern SCS, the CPAEs mainly migrated from C3 to C2 along the 18◦ N parallel from the northwest of the Luzon Island to the east of the Hainan Island. The CPAEs processes also frequently flowed from communities C1 to C2 along the shelf slope of the northern SCS from the southwest of Taiwan to the southeast of Hainan [56,57].

In the central SCS, the CPAEs mainly migrated westward from communities C6 to C4 between the 12◦ N and 16◦ N parallels (from the southwest of the Luzon Island to the northeastern Vietnam). In the southern SCS, there were strong interactions of CPAEs between communities C7 and C8, and from C9 to C7. In general, CPAEs showed strong interactions along the parallels between communities in the northern (C3 and C2), central (C6, C5 and C4), and southern (C8 and C7) SCS. Stronger regional interaction in the southeastern Vietnam Sea and the dominant westward migration across the central SCS are consistent with what was reported in previous studies [12–14].

The CPCEs in the northern SCS (Figure 11b) showed a strong interaction between the communities in the west of the Luzon Strait (C1) and the northwest of Luzon (C3). In the central SCS (C4, C5, C6, C7), interaction was more notable between C5 and C7 in the southwest of Luzon. In the southern SCS, interaction was more prominent in southeastern Vietnam between communities C6 and C8, C9 and C8.

Compared to the CPAEs, the CPCEs also showed significant interactions among communities along the meridian direction in the eastern SCS (C1, C3, C5, C7), which accounted for about 36% of the total CPCEs interaction in the SCS. Modeling results about how an island/seamount splits eddies showed that, after splitting, CPAEs tended to propagate southwestward while CPCEs did so toward the northwest [32]. The same patterns were not only observed in the splitting but also the merging processes in this study.

**Figure 11.** The community division results of CPAEs (**a**) and CPCEs (**b**). The numbers show the percentages of the eddies that migrate along with the dominate flow directions, which are shown by the blue arrows.

Within the community (Figure 12), both the CPAEs (a) and CPCEs (b) showed a dominant westward propagation direction, mainly due to the beta effect [55]. In the northwest of Luzon, quite large portions of the CPAEs and CPCEs (CPAEs-C3, CPCEs-C3) moved northward. In the southeast of Vietnam offshore, there were slightly more north–south than east–west eddy migration, probably because of the uneven Coriolis force [36].

**Figure 12.** The inter-community movement directions of CPAEs (**a**) and CPCEs (**b**).

#### **4. Discussion and Summary**

#### *4.1. Discussion*

4.1.1. The Generation and Extinction of Simple and Complex Eddy Processes

Our results show that both complex and simple eddy processes are mainly born in the eastern SCS, including the southwest of Taiwan Island, the west of Luzon Strait, the northwest of Luzon Island. This is mainly due to the strong wind stress curl variations and the intrusion of the Kuroshio in these areas [33–55].

The local wind stress curl is mainly controlled by the Asian monsoon, which can be divided into the winter (November to April) and summer (mid-May to mid-September) monsoon. Due to the topographic and narrow tube effects, the winds in the southwest of Taiwan, west of Luzon Strait and northwest of the Luzon Island are stronger than those in other regions in the SCS, especially during the winter monsoon period [34]. Under such a constraint, the surface water tends to move as a fast jet at the place with stronger wind shear potential. The jet flow would be blocked once it reaches the edge of the water and the water would reverse and move along the regions with weaker wind shear potential. Such an unstable water flow field make this region an active place of eddy formation [11,13,26]. The local wind stress curl may cause Ekman effect in this region [34]. The Ekman downwelling prevents the deep cold water from advecting upward and is favorable for the formation of anticyclonic eddies. By contrast, the Ekman upwelling is favorable for the formation of cyclonic eddies.

In addition to the local wind stress curl, the Kuroshio intrusion [33,58] is another significant contributor to eddy generation in the Luzon Strait. Previous studies have shown that Kuroshio intrusion can cause baroclinic instability and subsequently a horizontal density gradient. A large enough horizontal density gradient would provide enough energy and mass for the generation of mesoscale eddies in the SCS [33]. In the Luzon Strait, the Kuroshio may form a "current jacket", which could shear mesoscale eddies [58,59].

For the extinction of eddy processes, both simple and complex eddy processes mainly migrate toward the west because of the beta effect [55]. The complex eddy processes, with a longer life span than the simple eddy processes tend to travel longer toward the west (Table 1). As a result, they mainly disappear in the western SCS. However, the simple eddy processes travel a shorter distance due to their shorter life spans and smaller ranges of activity (Table 1) and consequently show no concentrated extinction zones in the SCS.

#### 4.1.2. Variations in Splitting and Merging Events in Different Regions in the SCS

Interactions among eddies through splitting and merging could be another factor that leads to new eddy formation. Our study shows that the regions in the SCS with significant splitting and merging events also witness the significant formation of eddies. The regions with significant eddy splitting and merging tend to be packed with more eddies and therefore more frequent or even stronger interactions among the eddies. Such interactions may form more new-born eddies.

Spatial distribution in eddy splitting and merging events probably indicate how the factors locally impact the eddy evolution processes. More splitting in the western Luzon Island may be attributed to the seabed topography, which is more complex in the west of Luzon Island [33]. Results from this study are consistent with previous studies [31], showing that topography is the major factor the triggers eddy splitting and merging events. This study also shows the west of Luzon Island probably is an ideal place for scientists to model the relationship between the eddy splitting activities and the topography, which, however, obviously is not the only factor leading to eddy splitting and merging.

We also found that, in the southwest of Taiwan, CPAEs tend to merge whereas CPCEs tend to split more. It is likely that this phenomenon could also be attributed to wind stress and the Kuroshio intrusion, which significantly affect the oceanic environment of this area [56]. Our knowledge about the SCS is not enough to elaborate on why CPAEs and CPCEs are prone to merging and splitting in this area, respectively. In fact, dense observation stations have been deployed in this area, providing enough data for physical oceanographers to investigate the mechanism behind this phenomenon. A more thorough discussion about the effects of topography and wind stress curl on eddy splitting and merging is available in the discussion Section 4.1.3.

In the southeast of Vietnam, both CPAEs and CPCEs tend to merge more, though they both also split. The strong eastward jet shooting away from the coast in this region pushes anticyclones to move southward and cyclones move northward [27]. It also significantly disturbs the water masses in this region, which forces CPAEs and CPCEs to either merge or split.

#### 4.1.3. Topography and Wind Stress Curl Effects

Previous studies have shown that, in other oceans, seabed topography may cause an eddy to split or merge and thus significantly affects the evolution and structure of eddies [28,29]. To understand how the seabed topography in the SCS affects the eddy evolution, we used GIS spatial analysis tools to calculate the slope of the seabed topography in the SCS. Here the slope refers to the maximum rate of change in elevation from a cell to its immediate neighbors. In this study, the cell was defined as a grid of 1◦ × 1◦ (about 108 km × 108 km in the SCS). We chose this resolution as the radii of most eddies in the SCS are between 100 km and 200 km, and the average moving distance is 74 km [12,14]. The slope was calculated for each cell from the seabed topographic data that was obtained from National Centers for Environmental Information. We then calculated the number of splitting and merging events within each grid. The calculation results, including the total number of grids and the number of grids that splitting and merging events occurred at least four times, were summarized by slope ranges (Table 2).


**Table 2.** The percentage of splitting and merging events in different slope ranges.

The eddies were more prone to splitting or merging in areas with a steeper slope. The percentages of the CPAEs and CPCEs that merged or split increased four times and more significantly increased with the increased seabed slopes, indicating the rugged seabed tends to disturb the eddies more significantly (Table 2). All (100%) of the CPAEs and CPCEs passing through the grids with a slope over 50 degrees either split or merge at least four times over our study time period. By contrast, roughly half (41%–47%) of the CPAEs and CPCEs merge or split four times within the cells with a seabed slope less than 10 degrees. Eddies in the southwest of Taiwan and the west of Luzon Island were also prone to splitting and merging as these two areas have an average slope of 37 and 39 degrees, respectively.

In addition to the topography, other factors such as wind stress, background currents, and Kuroshio invasion may also influence eddy splitting and merging. Wang et al. [11] pointed out that the wind stress curl formed by the interaction between wind and topography may be one of the important mechanisms for the formation of mesoscale eddies in the SCS. The southwest of Taiwan and the northwest of Luzon are the active areas of eddy formation, splitting, and merging. This study also showed that the eddy splitting and merging events in the SCS more frequently occurred in autumn. This is the reason when summer monsoon starts to weaken and winter monsoon starts to strengthen. Under such an unstable wind field, multi-eddy structures frequently appear in the upper circulation. The multi-eddy structure generally is a transitional state before eddy splitting, and/or merging. Once the winter monsoon dominates the northern SCS, wind stress significantly increases in the southeast of Taiwan and the west of Luzon Island due to the narrow tube effect [34]. The interaction between the wind stress and the topography probably further enhances eddy splitting and merging in these areas. However, further studies are needed to clearly illustrate how the topography and wind combine to drive the eddies to split or merge.

#### 4.1.4. The Interaction of Complex Eddy Processes among Communities in the SCS

In the northern and central SCS, both the CPAEs and CPCEs processes mainly propagate westward and the CPCEs also northward and southward in the west of Luzon. The remarkable westward migration trend is mainly due to the beta effect [55]. In the southern SCS, CPAEs and CPCEs show no uniform propagation direction but the north–south interactions are more frequent in the regions next to southeastern Vietnam. The strong eastward jet shooting away from the Vietnam coast pushes anticyclones to move southward and cyclones move northward in the southeastern of Vietnam [27].

Previous studies have also shown that significant poleward heat transfer occurs in the east of Vietnam in summer and the west of Luzon Island in winter, whereas large equatorward heat transfer occurs in Western Luzon Strait in winter [60]. However, it is not clear whether the same patterns of heat transfer could be found in other areas in the SCS due to the limitation of time scale and space range of the observation data [27,60]. Our results show that the interactions among communities show very similar heat transfer patterns across the whole SCS, as reflected by the eddy migration patterns (Figures 11 and 12). Further studies are needed to use the observation data, when they become available, to confirm the heat transfer patterns in the other SCS regions except for the east of Vietnam, west of Luzon Island, and Western Luzon Strait. Once the heat transfer patterns are identified for the whole SCS, we may better unravel the underlying physical mechanisms that drive the eddy-induced heat transfer in the SCS.

#### *4.2. Summary*

This study examined the mesoscale eddies in the SCS from a GIS perspective. As a real-world phenomenon, eddies could be represented either as a point, a line, or a polygon. The point and line then form a network, showing the connection and interactions among the eddies. In this study, we translated the eddies in the SCS into points, lines, and networks. Different GIS analysis methods were used to examine these three types of objects and unravel the kernel density and hotspots of the starting, ending, splitting and merging points, eddy migration patterns, eddy exchange, and interactions among communities in the SCS. We studied the difference of all the afore-mentioned aspects between the simple and the complex eddy processes, as well as between the anticyclonic and cyclonic eddies in the SCS.

Our results show that both the complex and simple eddy processes are mainly born in the eastern SCS, including the west of Luzon Strait, the northwest and southwest of Luzon Island. This is mainly due to the strong wind stress curl variations and the intrusion of the Kuroshio in these areas [33,55]. Both simple and eddy processes mainly migrate toward the west due to the beta effect [55]. The complex processes mainly disappear in the western SCS whereas the simple processes disappear almost everywhere across the SCS because of the different life span and scope of activity.

The CPAEs show more merging than splitting events in the southwest of Taiwan Island and the east of Vietnam. By contrast, the CPCEs mainly split in the southwest of Taiwan Island and the west of Luzon. Merging also frequently occurs in the southeastern SCS for both CPAEs and CPCEs. More rugged topography tends to trigger more splitting and merging events, particularly when the slope is more than 50 degrees (Table 2). However, it is not clear why the eddies merge or split more in some specific regions and less in other regions. Temporally, most splitting and merging events occur in the fall, when the wind direction starts to reverse in the SCS.

The SCS could be divided into the different number of communities based on the networks built from all-, simple-, and complex-eddy processes, separately. The community delineation results, as well as the eddy interaction and exchange among the communities share similarities but difference does exist. The most active inter-community interactions are mainly found in the southwest of Taiwan Island, the west of Luzon Island and the southeast of Vietnam. These are the regions where most eddies were born and then start to travel westward. The interactions between adjacent communities reveal the dominant migration pattern of both CPAEs and CPCEs from a totally different perspective. In the northern and central SCS, both the CPAEs and CPCEs processes mainly propagate westward. In the southern SCS, CPAEs and CPCEs show no uniform propagation direction but the north–south interactions are more frequent in the regions next to southeastern Vietnam due to the strong eastward jet shooting away from the Vietnam coast [27].

In short, multiple factors including the Kuroshio intrusion, local wind stress, topography, and the beta effect could all affect the behaviors of complex eddy processes in the SCS. Contributions of these factors tend to vary in different regions in the SCS. More comprehensive studies are needed to better reveal their different contributions in different regions.

**Author Contributions:** Conceptualization, H.W. and Y.D.; methodology, H.W.; software, Y.S.; validation, H.W., Y.D. and F.L.; formal analysis, H.W.; investigation, H.W.; resources, J.Y.; data curation, J.Y.; writing—original draft preparation, H.W.; writing—review and editing, F.L.; visualization, Y.S.; supervision, Y.D.; project administration, Y.D.; funding acquisition, Y.D.

**Funding:** This work was supported in part by a grant from the National Science Foundation of China (41671445), National Key R&D Program of China (2017YFB0503605).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Four Patterns of the East Branch of the Kuroshio Bifurcation in the Luzon Strait**

**Ruili Sun 1,\*, Fangguo Zhai <sup>2</sup> and Yanzhen Gu <sup>2</sup>**


Received: 17 November 2018; Accepted: 6 December 2018; Published: 10 December 2018

**Abstract:** Based on the self-organizing map (SOM) method, a suite of satellite measurement data, and Hybrid Coordinate Ocean Model (HYCOM) reanalysis data, the east branch of the Kuroshio bifurcation is found to have four coherent patterns associated with mesoscale eddies in the Pacific Ocean: anomalous southward, anomalous eastward, anomalous northward, and anomalous westward. The robust clockwise cycle of the four patterns causes significant intraseasonal variation of 62.2 days for the east branch. Furthermore, the study shows that the four patterns of the east branch of the Kuroshio bifurcation can influence the horizontal and vertical distribution of local sea temperature.

**Keywords:** Kuroshio bifurcation; Luzon Strait; SOM; temporal and spatial variation; sea surface temperature; mesoscale eddy

#### **1. Introduction**

The Kuroshio, which is the strongest western boundary current in the Northwest Pacific (NWP), originates from the North Equatorial Current (NEC). Encountering Luzon Island, the NEC bifurcates into the northward Kuroshio Current and the southward Mindanao Current. Then the Kuroshio passes to the east of the Luzon Strait (LS), flows northward along the Taiwan Island into the East China Sea, and returns to the Pacific Ocean through the Tokara Strait. Because the Kuroshio transports a lot of heat and matter from low latitudes to mid-latitudes, it plays an important role in regional and global climate, and the marine ecosystem [1–5]. Meanwhile, the Kuroshio carries enormous amounts of momentum and kinetic energy, so it strongly interacts with mesoscale eddies and influences regional circulation such as the one in the LS [6–11].

The LS, located between Taiwan Island and Luzon Island, is an important gap in the NWP (as shown in Figure 1). Affected by the Kuroshio, monsoon, mesoscale eddies, and topography, the circulation in the LS is very complicated, and has been studied extensively in previous literatures [8,9,12]. These studies have specifically discussed the interaction between the Kuroshio and mesoscale eddies: (1) the Kuroshio can generate mesoscale eddies by its own variation [11,13]; (2) the Kuroshio can affect mesoscale eddies in the NWP into the South China Sea [7,9]; (3) mesoscale eddies can alter the Kuroshio current field [8,10,14]; (4) interaction between the Kuroshio and mesoscale eddies can influence the local marine environment [14,15]. Although there are many related studies, some questions involving small and moderate ocean phenomena in the LS remain unclear.

**Figure 1.** Climatology of geostrophic current and sea surface temperature (SST) around the LS. The black line represents 27.2 ◦C, in order to emphasize the spatial patterns of two warm tongues. The white box covers 20.3◦ N–22◦ N, 121.875 ◦E–123.125 ◦E. Geostrophic current (m/s, in vector); SST (◦C, in color); TWI: Taiwan Island; BIS: Batanes Islands; LZI: Luzon Island; SCS: South China Sea. The extent of the main map is shown as a black box in the inset.

The accumulation of high-resolution satellite remote sensing data enables us to better study small and moderate ocean phenomena. Based on high-resolution satellite remote sensing data and hydrological data along with numerical model simulations, the Kuroshio bifurcation phenomena in the LS has been put forward: The Kuroshio in the LS bifurcates into two branches located on the western and eastern sides of the Batanes Islands, and produces two warm tongues east of the LS [14]. The western branch is the main part of the Kuroshio, and its temporal and spatial variability has been widely studied [16–19]. However, the temporal and spatial variability of the eastern branch has hardly been studied and remains unclear, although its intensity is related to mesoscale eddies in the NWP [14]. The Kuroshio transports warm water from low latitudes to middle latitudes, so it can cause temporal and spatial variability of local sea temperature in most cases [14,16]. The east branch of the Kuroshio bifurcation, as a part of the Kuroshio, can produce a warm tongue [1] (Figure 1). However, it is unclear whether and how local sea temperature patterns are related to the patterns of the east branch. The purpose of the present study is to assess the variability of the east branch of the Kuroshio bifurcation and its effect on the local sea temperature. The rest of the paper is organized as follows: Section 2 briefly introduces the data and methods. Section 3 shows the temporal and spatial variability of the four patterns of the east branch and their effect on local sea temperature. Section 4 is a summary of this paper.

#### **2. Data and Method**

#### *2.1. Data*

Absolute dynamic topography (ADT) data and geostrophic current anomalies data were processed by Ssalto multimission ground segment (SSALTO)/Data unification and Altimeter combination system (DUACS) (AVISO, Ramonville Saint-Agne, France) and distributed by Archiving Validation and Interpretation of Satellite Data in Oceanography (AVISO) with support from CNES (Centre National d'Etudes Spatiales). The product is merged from Jason-1/2, Topex/Poseidon, Envisat, Geosat, and other altimeter missions. The raw satellite measurements have applied atmospheric corrections and a sea-land transition. The product is computed on a Cartesian 1/4◦ × 1/4◦ spatial resolution and provides daily data from 1 January 1993 to the present [20,21].

Sea surface temperature (SST) is provided by Remote Sensing Systems (RSS, Santa Rosa, CA, USA). The product combines the through-cloud capabilities of the microwave SST data with the high spatial resolution and near-coastal capability of the infrared SST data. The microwave SST data is merged from Advanced Microwave Scanning Radiometer for Earth Observing System, WindSat, Tropical Rainfall Measuring Mission Microwave Images and Giant magnetoimpedance. The infrared SST data is merged from Moderate Resolution Imaging Spectroradiometer and Visible Infrared Imaging Radiometer Suite. The product has applied atmospheric corrections. It provides daily data from 1 July 2002 to present and its spatial resolution is 9 km × 9 km.

Hybrid Coordinate Ocean Model (HYCOM) reanalysis data is provided by HYCOM organization, which provides access to near-real-time global data-based ocean prediction system output. The HYCOM organization provides four HYCOM data-assimilation product, and we choose the product with the longest time span, which is from 2 October 1992 to 31 December 2012. The product has uniform 0.08 degree latitude/longitude grid between 80.48 ◦S and 80.48 ◦N and is interpolated to 40 standard z-levels. It provides sea surface height (SSH), temperature, salinity, meridional flow, and zonal flow.

#### *2.2. Method*

SOM can be classified as a type of clustering technique [22], and is an effective method for feature extraction and pattern classification [22,23]. A best-matching unit (SST), which records the category of each pattern and can reflect the evolution of each pattern, is obtained based on the minimum Euclidean distance from the input data. Previous studies have demonstrated that the SOM is more useful than conventional EOF (Empirical Orthogonal Function) to extract nonlinear information [24–27], and showed that the SOM has been applied in oceanography and atmosphere science [15,26,28–31]. It needs to be pointed out that the SOM can give the occurrence time of the patterns, which can be used to calculate the period of the patterns [15,32].

Tunable parameters need to be specified in the process of SOM training. The parameters such as training method, lattice, initialized weights, and map shape are chosen according to Liu et al. [27], who gave a practical method for choosing parameters. The parameter of map size defines the number of patterns and should be specified subjectively according to the variability of the shift and strength of the east branch in the LS. The shift of the east branch can generally be classified into two types: westerly or easterly, and the strength of the east branch can also be classified into two types: strong or weak. Thus, a map size of 2 × 2 is chosen. The method we use to choose the map size is similar to the one described in Jin et al. [22]. We run a series of tests with different map sizes such as 2 × 3, 3 × 3, 3 × 4 and 4 × 4. The east branch patterns with these larger map sizes are very similar to the ones with the map size of 2 × 2 although the patterns in the larger SOM provide more details about the east branch variability.

#### **3. Results**

#### *3.1. The Four Patterns of the East Branch of the Kuroshio Bifurcation*

Applying the SOM method to the corresponding AVISO ADT data in the white box of Figure 1, we can obtain the dominant variation patterns of the east branch. Figure 2 shows that the east branch has four coherent patterns: anomalous southward (P1), anomalous eastward (P2), anomalous northward (P3), and anomalous westward (P4), which correspond to a cyclonic eddy (CE), CE-anticyclonic eddy (AE), an AE, and AE-CE to the east of the LS, respectively. Briefly, the four patterns of the east branch are dominated by single and dipole eddy structure, which is very similar to the SSH variability east of the Taiwan Island [33,34]. To the authors' knowledge, the four patterns of the east branch in the LS have not been previously reported. P1 shows significantly southward velocity anomaly, which is caused by a CE to the east of LS. P2 shows significantly eastward velocity anomaly, which is accompanied by an eddy pair consisting of a CE on the north side of 21.2 ◦N and an AE on the south side of 21.2 ◦N. P3 is the opposite of P1, and is caused by an AE to the east of the LS. P4 is the opposite of P2, and is accompanied by an AE on the north side of 21.2 ◦N and a CE on the south side of 21.2 ◦N. The four patterns of the east branch also are identified in the HYCOM SSH data. Figure 3 shows that, P1 (P3) corresponds to a CE (an AE), and P2 (P4) corresponds to a CE (an AE) on the north side of 21.2 ◦N and an AE (a CE) on the south side of 21.2 ◦N, which corresponds well to the ones in Figure 2, although the occurrence rate of corresponding pattern of Figures 2 and 3 is not exactly the same, which may be related to HYCOM model bias.

**Figure 2.** The 2 × 2 SOM patterns derived from the AVISO ADT fields for the east branch. Geostrophic current anomaly (m/s, in vector). The number on each panel represents the occurrence rate of corresponding pattern.

**Figure 3.** The 2 × 2 SOM patterns derived from the HYCOM SSH fields for the east branch, and corresponding sea temperature anomaly in 300 meters deep. Sea Surface Height Anomaly (SSHA) (m, in contour); Sea temperature (◦C, in color). The number on each panel represents the occurrence rate of corresponding pattern. The red dotted line represents meridional section of sea temperature, and its longitude is 123.2 ◦E.

The temporal evolution of the four patterns can be seen from the time series of the BMU (Figure 4). The cycle characterized as a clockwise trajectory in the SOM space is very noteworthy, which has been marked by red lines in Figure 4. We take the cycle of P1, P2, P3, P4 and P1 as an example. First, a CE approaches from the east and weakens the east branch (P1). Second, as the CE flows northward and gradually weakens, an AE begins to form in the south (P2). Third, with the AE eddy strengthening, it approaches and enhances the east branch (P3). Finally, as the AE moves northward and gradually weakens, a CE in the south begins to form (P4). Note P1 as the initial pattern in the sequence is arbitrary. We define the occurrence of the preferred cycle by counting occurrences of clockwise trajectories that pass through all four patterns at least once. The percentage of time that follows the preferred cycle (red lines) to the total time analyzed is 39.2%.

Power spectrum analysis shows that the most significant cycle of the BMU is 62.2 days (Figure 5), which is consistent with the significant cycle of upper-layer and intermediate-layer transports based on observation data in the LS [35], and Zhang et al. pointed out that the ~60 days of variability might be attributed to the impinging mesoscale eddies from the Pacific [35]. The average consecutive days are 19.2, 6.0, 18.8, and 6.5 for P1, P2, P3, and P4, respectively, causing the total time length of a cycle of P1, P2, P3, and P4 to be between 56.5 days and 69.7 days, which is very close to 62.2 days. Therefore, we think that the cycle of the 62.2 days of the BMU is mainly caused by the cycle of P1, P2, P3, and P4.

**Figure 4.** Time series of the BMU for the four patterns, as shown in Figure 2. On the y axis, 1, 2, 3 and 4 correspond to P1, P2, P3, and P4, respectively. The subgraphs (**a**), (**b**), (**c**), and (**d**) correspond to the period from 1993 to 1998, from 1999 to 2004, from 2005 to 2010, and from 2011 to 2015, respectively.

**Figure 5.** Power spectrum analysis of time series of the BMU for the four patterns. The red dotted line represents the 95% significant level.

#### *3.2. The Effect of the Four Patterns on the Local Sea Temperature*

SST around the LS has obvious seasonal variability. Figure 6 shows that there is a good correspondence between the Kuroshio and warm tongues in most months, especially in winter, However, the Kuroshio fails to correspond well with warm tongues in some months such as in June and August. The high temperature in June and August is not conducive to the formation of warm tongue around the LS [32].

**Figure 6.** Seasonal variability of SST and geostrophic current around the LS. SST (◦C, in color); Geostrophic current (m/s, in vector).

We apply a 4◦ moving average for RSS SST field for the purpose of removing the background state to highlight the effect of mesoscale process. Figure 7 shows that the four patterns of the east branch have a significant effect on SST around the LS. Sun et al. used a mixed-layer model based on temperature equation to demonstrate that the Kuroshio affects SST mainly through geostrophic thermal advection in the LS [14]. Through thermal advection, the east branch of the anomalous southward causes SST reduction in the location of the east branch (P1); the east branch of the anomalous eastward causes an eastward warm tongue (P2); P3 is the opposite of the P1, and SST increases in the location of the east branch (P3); P4 shows that the east branch of the anomalous westward causes a westward cold tongue. Figure 7 embodies that mesoscale eddies can influence sea temperature distribution by modulating the east branch of the Kuroshio of the western boundary current.

**Figure 7.** The spatial pattern of high-pass-filtered RSS SST corresponding to 2 × 2 SOM patterns derived from the AVISO ADT fields for the east branch. Geostrophic current anomaly (m/s, in vector); SST (◦C, in color). The spatial filtering is done by subtracting a 4◦ moving average from the original RSS SST data to remove the large-scale background state. The black line represents a zero contour line.

Besides affecting SST, the four patterns of the east branch can also affect deep sea temperature. Figure 3 shows that there is a good correspondence between sea temperature in 300 meters depth and the four patterns of the east branch: the CE (AE) corresponds to a cold (warm) center of sea temperature anomaly because of the following reasons. When the CE (AE) exists, SSHA is negative (positive) and upwelling (downwelling) occurs, and then the temperature in the deep ocean correspondingly decreases (increases) [36], which indicates that the four patterns of the east branch can affect the vertical distribution of sea temperature. Figure 8 gives vertical profile of sea temperature anomaly along the 123.2◦E, which is marked as a red dotted line in Figure 3, and shows negative (positive) temperature anomaly corresponding to a CE (an AE) can reach 5000 meters deep. Based on Conductivity Temperature Depth (CTD) data, Zhang et al. found that mesoscale eddies can affect sea temperature in 2530 meters depth [36], and the conclusion in this paper is consistent with it in magnitude. However, based on HYCOM data in this paper, we found that the effect of the four patterns of the east branch on sea temperature can reach 5000 meters in depth, or even deeper [37].

**Figure 8.** The profile of HYCOM temperature anomaly corresponding to the four patterns of the east branch. Subgraphs (**a**), (**b**), (**c**), and (**d**) correspond to P1, P2, P3 and P4 of the east branch patterns, respectively. The black line represents the zero contour line. It is noted that the zonal mean has been removed. Sea temperature (◦C, in color).

#### **4. Summary**

Spatial patterns of the east branch of the Kuroshio bifurcation in the LS are extracted from the AVISO data and HYCOM data using the method of SOM. Four coherent patterns are used to summarize the variation of the east branch: anomalous southward, anomalous eastward, anomalous northward, and anomalous westward, which corresponds to a CE, CE-AE, an AE, and AE-CE to the east of the LS, respectively. The robust clockwise cycle of the four patterns causes a 62.2 days cycle of the east branch.

The four patterns of the east branch modulated by mesoscale eddies can influence horizontal and vertical distribution of local sea temperature. East branch of anomalous southward (northward) causes SST reduction (increase) in the location of east branch. East branch of anomalous eastward causes an eastward warm tongue while east branch of anomalous westward causes a westward cold tongue. The effect of the four patterns of the east branch on sea temperature can reach 5000 meters deep.

The patterns of the east branch of the Kuroshio bifurcation and their temporal evolution, as revealed by the SOM, provide new insights into the variability of the Kuroshio and their effect on the local sea temperature. Detailed dynamic mechanism studies with numerical models and more in situ observations are warranted as part of future work.

**Author Contributions:** Conceptualization, R.S.; Methodology, R.S.; Software, R.S.; Validation, R.S., and F.Z.; Formal Analysis, R.S.; Investigation, R.S.; Resources, R.S.; Data Curation, R.S.; Writing-Original Draft Preparation, R.S.; Writing-Review & Editing, R.S., F.Z. and Y.Z.; Visualization, R.S.; Supervision, R.S., F.Z. and Y.Z.; Project Administration, R.S.; Funding Acquisition, R.S.

**Funding:** This research was funded by National Natural Science Foundation of China grant number (41806019), Natural Science Foundation of Zhejiang Province grant number (LY18D060004), Foundation of Department of Education of Zhejiang Province grant number (1260KZ0417982), National Natural Science Foundation of China grant number (41506008), Open Fund of the Key Laboratory of Ocean Circulation and Waves, Chinese Academy of Sciences grant number (KLOCW1503), Talent Start Foundation of Zhejiang Gongshang University grant number (1260XJ2317015 and 1260XJ2117015), and The APC was funded by National Natural Science Foundation of China grant number (41806019).

**Acknowledgments:** The authors thank AVISO (Archiving Validation and Interpretation of Satellite Data in Oceanography, available for download at: ftp://ftp.aviso.oceanobs.com/), RSS (Remote Sensing System, available for download at: ftp://ftp.remss.com), and HYCOM (Hybrid Coordinate Ocean Model, available for download at: https://hycom.org/dataserver/) for providing the data used in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
