Glider observations of interleaving layers beneath the Kuroshio primary velocity core east of Taiwan and analyses of underlying dynamics

Submesoscale interleaving layers are caused by lateral intrusions of dissimilar water masses in frontal zones, which are significant processes in shaping physical, biogeochemical, and ecological parameters in the ocean. Possible interleaving layers were sometimes observed by ship-based conductivity-temperature-depth (CTD) surveys with coarse spacing between adjacent stations in the Kuroshio region east of Taiwan but have never been examined dynamically. Here we show the characteristics of interleaving layers observed by a Seaglider with two repeated hydrographic surveys along a triangle track east of Taiwan from December 2016 to March 2017. Salinity profiles indicate that prominent interleaving layers appeared in the intermediate layer (approximately 500–800 m) with vertical and horizontal length scales of O(50) m and O(10–100) km, respectively, during our observations. A dipole eddy pair and a relatively large anticyclonic eddy impinged on the Kuroshio during the first and second surveys, respectively, which brought certain impacts on the interleaving motion as the eddy potentially altered the density slope across the Kuroshio. The associated instability analysis and the Turner angle suggest that the double diffusive instability is the primary driving mechanism for the development of interleaving layers.


Region Depth range (m) Thickness (m) Length (km)
Drake Passage -Arctic Polar Front 13 www.nature.com/scientificreports www.nature.com/scientificreports/ intrusions from one to the other water (see Supplementary Fig. S1). The following growth of the interleaving is attributed to either the double diffusion 21 or the advective driven intrusion through the cross-front velocity perturbations 22 . The associated driving mechanism of double diffusive intrusions is the salt fingering and diffusive convection generated subsequently with those initial disturbances. The alternative appearance of salt fingering and diffusive convection in the stratified layers generates convergence or divergence in the vertical buoyancy fluxes and, in turn, stretches the intrusions into distinct interleaving structures 1,14,15,21 (Supplementary Fig. S1). In contrast, the advective intrusion is dominated by the cross-front inertial velocity anomalies through the baroclinic instability 22 .
The importance of examining the interleaving in the Kuroshio is that it provides a mechanism for the transformation of water mass properties, which can alter the density slope across the Kuroshio and, in turn, the associated velocity structure and throughflow volume transport. However, the spacing of adjacent CTD stations 6 (~20 km) is too coarse to resolve the horizontal coherence of the interleaving layers and associated dynamics. The existence of this submesoscale structure indicates potential water mass exchange between South China Sea Water, Kuroshio Water or North Pacific Water in the Kuroshio. The characteristics of interleaving layers in the Kuroshio east of Taiwan and the formation of the interleaving have yet to be examined. To the best of our knowledge, this study serves as the first effort to quantify and examine the interleaving layers in the Kuroshio east of Taiwan using high-resolution hydrographic data collected by autonomous underwater vehicle (Kongsberg's Seaglider in this study). The Seaglider was navigated along the three sections of a triangle track east of Taiwan (Methods). Two repeated surveys, the first and second surveys, were completed from 15 December 2016 to 14 January 2017 and from 14 January 2017 to 5 March 2017, respectively (Methods).

Results
Mesoscale features of sea surface height during the observations. The 7-day composite sea surface height (SSH) obtained from the Archiving, Validation, and Interpretation of Satellite Data in Oceanography (AVISO) provides synoptic oceanic features during our glider observations (Fig. 1). After the launching of the glider near the southernmost of Taiwan on 8 December 2016, a dipole eddy pair (A1 and C1 in Fig. 1) encountered the offshore side of the Kuroshio (Fig. 1a). The hydrographic samplings along the transect N1 captured the influence of the anticyclonic eddy (A1 in Fig. 1) of the eddy pair. The dipole eddy pair migrated ~50 km westward during the sampling along the transect N2, where it was mostly influenced by the eastern half of the anticyclonic eddy A1 as well as the westward-flowing confluent current of the eddy pair at ~22°N (Fig. 1b). The confluent flow continuously affected the velocity and hydrography at the transect N3 (Fig. 1c). During the second repeated survey at N1, a small, short-lived cyclonic eddy (C2 in Fig. 1) with diameter of only ~100 km appeared on the eastern half of N1 (Fig. 1d). The rest of the second survey was primarily influenced by a larger anticyclone with a diameter of ~400 km (A2 in Fig. 1e,f).
The westward-propagating mesoscale eddy is a plausible albeit debated mechanism transporting North Pacific water from its origin to the eastern flank of the Kuroshio. The inner core of an eddy, defined by the contour of zero relative vorticity, may comprise only water mass from its origin, and the outer ring of the eddy may contain a mixture of ambient water from throughout the eddy's life 24 . Alternatively, the original water mass can be trapped and transported independently by propagating mesoscale eddies 25 . Regardless of this debate, the encountering of North Pacific water with Kuroshio water potentially creates fronts between the two water masses or alters the density slope across the Kuroshio, which provide basic conditions for the interleaving layers observed in this study. The T-S diagram suggests that intermediate waters primarily consisted of NPIW and KIW and their mixture in all transects during the first survey ( Fig. 2a-c). The zigzag T-S curve, which is a necessary anatomy of interleaving layers, is seen in many of the T-S profiles on the three transects, e.g., the gray curve in Fig. 2a. SCSIW was rarely observed during the first survey; instead, a mixed water mass of KIW and SCSIW appeared near the onshore end of N1 and N3, between σ θ = 26.0 and 26.8 kg m −3 (Fig. 2a,c).

Water masses in the
During the second survey, the water mass at N1 was mostly NPIW near the offshore end, tended to be KIW in the middle and was close to NPIW again near the onshore end of this section (Fig. 2d). It is plausible that the cyclonic eddy C2 in Fig. 1d played a certain physical role in affecting the spatial distribution of NPIW and KIW at N1. The salinity distribution in Fig. 3d (next subsection) shows that isohalines are elevated around 122.8°E because of the upwelling in the cyclonic eddy. Therefore, the intermediate water mass appeared to be KIW-like in the middle portion of N1 transect. It is clear that NPIW was the dominant water mass observed at the meridional transect N2. At N3, the dominant water mass was NPIW over the eastern half of the transect, and the water mass became the mixture of NPIW and KIW west of 122.5°E. The zigzag features remained in many T-S profiles, but the salinity anomaly value in a single zigzag was approximately 50% of that observed during the first survey. Since the surrounding outer ring of an anticyclonic eddy had reached the east and south sides of the triangle track www.nature.com/scientificreports www.nature.com/scientificreports/ during the second survey (Fig. 1e,f), the observation of very low S min there suggests that the eddy carried NPIW. This observation lends support to the viewpoint that an eddy carries water mass 25 . Interleaving layers in the salinity profiles. We only demonstrate glider observed salinity in the three transects in Fig. 3 because the interleaving feature is more distinct in the salinity than in the temperature profiles. Overall, the most noticeable interleaving appeared in 121.9-123°E and 123.3-123.9°E, within 500 to 800 m during the first survey at N1 (dashed red rectangle in Fig. 3a). The higher salinity layer (S > 34.32) intruded toward the east from ~122°E to 122.9°E (~100 km), and its thickness was ~180 m near 122°E and 20-50 m between 122.5 and 122.9°E, with smaller scale layered structures seen from the zigzag feature in salinity profiles. At the meridional transect N2, sporadic interleaving features with length scale smaller than those appeared at N1 are seen from Fig. 3b. At N3, no large scale interleaving layers such as those observed at N1 appeared in the salinity transect (Fig. 3c). Approximately one month later (Table 1), the glider reached the western corner of the triangle track again, and the repeated survey at N1 suggests that the salinity pattern in Fig. 3d differs from that in Fig. 3a, presumably due to the influence of the small cyclone (Fig. 1d). The isohalines within 200 and 450 m depths were heaved between 122.8° and 123.2°E. Comparing Fig. 3d with Fig. 3a, interleaving layers on the onshore half of N1 became blurred during the second survey. A very low salinity patch (S < 34.10) was observed at ~600 m near the easternmost part of N1, likely accompanied with shorter interleaving scale (Fig. 3d). The second survey at N2 and N3 observed that relatively low salinity waters occupied the intermediate layer from 550 to 750 m, which were presumably brought by the impinging anticyclone (Fig. 1e,f). Possible interleaving layers are seen in the intermediate layer between 22.2° and 22.5°N and between 23.3° and 23.5°N at N2 (Fig. 3e) and between 122° and 122.3°E and between 123.2° and 123.7°E at N3 (Fig. 3f). Overall, the salinity differences between these interleaving layers are from 0.05 to 0.10.
Note that the prevailing winds in our hydrographic survey region is the northeasterly monsoon in winter-like months from October to February and southwesterly monsoon in summer, with the former being stronger than the latter. During our glider observations, the mixing in the upper ocean was conceivably enhanced by the winter monsoon 26   www.nature.com/scientificreports www.nature.com/scientificreports/ interaction, which ranges from ~60 to 90 days 12 . Hence the upper bound for the uncertainty of salinity measurement in a zonal transect within 7 days is ~0.02, which is still less than the salinity differences between the observed interleaving layers 0.05-0.10.

Discussion
To objectively represent the locations of interleaving layers, we adopted the diapycnal spiciness curvature (DSC, τ σσ ) (Methods) as a measure of the lateral coherence and strength of interleaving layers. A high value of τ σσ represents a clear diapycnal interface between two different water masses. As τ σσ tends to be zero, it means that the diapycnal mixing diminishes the curvature of the vertical profile in the absence of horizontal advection 1 . Figure 4 shows τ σσ in longitude (or latitude) versus potential density (25.0-27.4 kg m −3 ) plots of each transect. It is easily seen that the corresponding locations of relatively high absolute values of τ σσ are primarily associated with the salinity transects in Fig. 3. The prominent interleaving mostly occurred along the σ θ = 27.0 kg m −3 isopycnal during the observations of the first survey at N1 (Fig. 4a) and N3 (Fig. 4c) and the second repeated survey at N3 (Fig. 4f).
In addition to the stratification and the differences between diffusivities of momentum, temperature, and salinity, the growth of interleaving layers is conceivably affected by the background shear and density gradient in a highly dynamic environment such as the Kuroshio off the east coast of Taiwan. The influence of baroclinic and shear effects on the growth of interleaving intrusions has been examined with the consideration of baroclinic and sheared barotropic geostrophic fronts containing constant horizontal and vertical shears and a horizontal density gradient 27 . The velocity shear could tilt the slope of the interleaving layer and, in turn, influence the growth of the interleaving. The horizontal density gradient in the front could enhance or suppress the growth of the www.nature.com/scientificreports www.nature.com/scientificreports/ interleaving, depending on the isopycnal slope relative to the interleaving layer slope. The growth of the interleaving is enhanced when the interleaving layer is along the isopycnal slope and is suppressed when the interleaving slope opposes the isopycnal slope 27 , which are schematically illustrated in Supplementary Fig. S3.
The effect of the baroclinicity on the growth of interleaving layers in the Kuroshio is examined using the inverse growth rates of the fastest growing perturbation in a low-shear and high-shear environment, λ −  (Fig. 1), which is close to transect N1, were used, and the shear was calculated using the barotropic (depth-averaged) velocity difference and the distance between two adjacent stations. The mean barotropic shear is thus 7.39 × 10 −7 s −1 with the local barotropic shear between two adjacent stations from −1.91 × 10 −6 to 2.64 × 10 −6 s −1 (Supplementary Fig. S4). Considering that diffusivities in the Kuroshio should be eddy-driven due to the dramatic influence of the frequent impinging eddies and the baroclinic shear, the turbulent Prandtl number is set to 1. With the horizontal and vertical properties of the observed interleaving layers in the intermediate layer at N1 in Table 2, the values of λ − l 1 and λ − h 1 are 2.1-3.7 and 42.0-145.5 d, respectively. Because the life time of the interleaving at each transect is hardly obtained by the glider observations, it is not easy to evaluate which of the two time scales is reasonable in this study. However, the two interleaving features observed during the two surveys at the same transect are dramatically different (e.g., Fig. 4a vs. 4d and 4c vs. 4f), and the times of the two surveys at the transect are ~30 days apart, which provides an upper bound on the resident time of an interleaving structure observed in a fixed transect. In addition to the time separation of the two surveys at a transect, the horizontal advection induced by the velocity in the intermediate layer (typically less than 0.1 m s −1 ) could shift interleaving features to the downstream for ~8.64 km per day, which suggests that the resident time of an interleaving structure in a transect may be even shorter. Hence, the time scale of 42.0-145.5 d for the growth of the interleaving layers at the N1 transect can be excluded, and it suggests that the high-shear limit is unlikely the condition in the intermediate layer of the Kuroshio. Moreover, the absolute value of mean barotropic shear v x (~7.39 × 10 −7 s −1 ) is much smaller than the fastest growth rate in the low-shear limit λ l (3.1-5.5 × 10 −6 s −1 ). This quantity of λ l is obtained using its inversion (λ − l 1 ) defined in Methods and the values quantified from our glider observations in Table 2. This relationship suggests that the development of interleaving layers (where k is cross-front wavenumber, and m is vertical wavenumber). The criterion is applied to the hydrography at N1 obtained from the first survey. The upper bound of the criterion ranges in 1.2-2.9 × 10 −3 , which was calculated using the quantities in Table 2. To estimate m, the wavenumber spectrum of the salinity anomaly was calculated. The salinity anomaly ′ S (z) was obtained by subtracting the background salinity S(z) from the observed salinity profile, i.e., ′ = − S (z) S(z) S(z). The background salinity profile was calculated by applying a fourth-order Butterworth filter with a cutoff wavenumber of 0.0125 m −1 (wavelength of 80 m) to the vertical salinity data 18 . After these calculations, the spectral peak of the salinity wavenumber spectrum is located at approximately 0.0167 m −1 (i.e., vertical wavenumber m). The horizontal wavenumber k estimated from the interleaving features in Fig. 1a is ~1.0 × 10 −5 m −1 . Therefore, the maximum interleaving slope k/m is ~6 × 10 −4 , which is within the range of 0 and the upper bound of 1.2-2.9 × 10 −3 . This result indicates that the driving mechanism of the observed interleaving layers in the N1 section (Fig. 3a) can be attributed to the double-diffusive instability.
The impingements of cyclonic and anticyclonic eddies on the Kuroshio during the observations altered the density slope across the Kuroshio 8,12 . The variation in the density slope may provide extra horizontal density gradients that are either positive or negative to enhance or suppress the growth of interleaving layers in the cross-front direction as depicted in Supplementary Fig. S3. During the first survey, the western half of the anticyclonic eddy A1 (Fig. 1a) potentially enhanced the density slope across the Kuroshio, which was a positive effect on the growth of the observed interleaving layers at N1 and N3 (Supplementary Fig. S3a). The eddy condition was complicated during the second repeated survey. The western half of the small cyclone C2 on the transect N1 (Fig. 1d), which tended to weaken the density slope of the Kuroshio, could possibly inhibit the development of interleaving layers at N1 (Figs 3d and 4d). The following anticyclonic eddy A2 conceivably increased density slopes to its surrounding waters, which further influence the growth of the interleaving at the three transects. The inference merits further comprehensive field observations to verify.
The Turner angle T u 28 (Methods) is helpful in analyzing stability and regimes of salt-fingering and diffusive convection in the observed interleaving layers. The values of T u at each transect during each survey are illustrated in Fig. 5. The geostrophic velocities calculated from the hydrographic data at each transect through the thermal wind relationship are overlaid on Fig. 5. The stratification is in the (strong) diffusive convection regime as T u is between −45° (−72°) and −90° and in the (strong) salt fingering regime between 45° (72°) and 90°. The stratification is stable as T u is between −45° and 45° and statistically unstable as T u is beyond ± 90°. Here, we focus on processes associated with the double-diffusive instability. The value of T u indicates that most of the regions in each transect are subject to salt fingering. In the primary interleaving layers described based on Fig. 3, the salt fingering regime and diffusive convection regime alternately appeared in the interfaces. Accompanying the hydrographic structure of each transect in Fig. 3, the Turner angle further helps locate where double-diffusive processes are effective and dominating the development of the interleaving layers. As previously described, how the double-diffusive instability causes interleaving layers is schematically illustrated in Supplementary Fig. S1. The isotachs of geostrophic velocities also suggest that the primary interleaving layers in the intermediate layer were in weak flow region (<0.1 m s −1 ) of the Kuroshio.
Note that the double diffusive instability, including salt fingering and diffusive convection processes, is a necessary but not sufficient condition for the generation of interleaving layers. The strong current with velocity greater than 0.2 m s −1 in the upper 500 m (Fig. 5) and the weak salinity variation in the layer between 200 and 500 m (around T~15 °C and S~34.6 in Fig. 2) restrict the development of interleaving layers in a fixed transect.
Finally, whether the observed interleaving layers were caused by the advection of the cross-front near-inertial velocity perturbations or the inertial instability similar to the observations in the western Equatorial Pacific 18 or in the Agulhas Current 22 is discussed. Since there is no direct measurement of velocity with our glider observations, we are not able to examine the relationship between the salinity perturbations and the cross-Kuroshio velocity anomalies as done for the data set collected in the Agulhas Current 22 . However, in our observations, the maximum interleaving slope (k/m) satisfies the criterion of the double-diffusive instability as the primary cause  Table 2. Horizontal and vertical properties quantified from the observations during the first survey at the intermediate layer of the N1 transect (500-800 m and 122.71°-122.51°E in Fig. 3a).
www.nature.com/scientificreports www.nature.com/scientificreports/ of the interleaving, which suggests that the cross-Kuroshio velocity anomaly may be not effective. To verify this inference, the criterion for the effectiveness of inertial instability on the growth of interleaving layers 18 , i.e., fQ < 0, 2 ) is the background potential vorticity, is further evaluated. The inertial frequency f is 5.82 × 10 −5 s −1 (listed in Table 2) at our glider survey region. To our estimate in Discussion, the order of magnitude of the barotropic shear v x is O(10 −7 -10 −6 ) s −1 , which is much larger than u y and fv N / z 2 2 in the components of potential vorticity Q. Therefore Q must be positive and fQ > 0 in the Kuroshio. In this light, the cross-front velocity perturbations and the inertial instability are excluded from the primary driving mechanism of the interleaving in the Kuroshio.
In conclusion, the interleaving layer between dissimilar water masses has been observed in many low energy frontal zones in the world ocean. Our glider observations off the east coast of Taiwan also observed prominent interleaving layers beneath the Kuroshio primary velocity core (poleward velocity > 0.2 m s −1 ). The results from the high-resolution hydrographic data quantified the characteristics of the interleaving beneath the primary velocity core of the Kuroshio east of Taiwan. Not surprisingly, the main feature of these interleaving layers primarily appeared in the cross-front direction of the Kuroshio. During the second survey at N2 and N3, the glider observed the western half of a westward propagating anticyclonic eddy that encountered the Kuroshio, and the characteristic of the water mass in the junction of the Kuroshio and the eddy more likely presented North Pacific Intermediate Water than that observed during the first survey on the eastern portion of the tracks. Although it is not the central focus of this study, the difference of water masses between the first and second surveys suggests that this anticyclonic eddy carried North Pacific water mass, which favors the concept of westward propagating mesoscale eddies that can trap and transport water masses 25 .
Since our glider survey was not a time-series observation, when and where the interleaving is initiated to perform is not resolved. It could be an evolving process during the observations or a result remaining after a previous www.nature.com/scientificreports www.nature.com/scientificreports/ event of eddy impingement on the Kuroshio. Moreover, the variations of the Kuroshio's velocity structure, baroclinicity, and water mass composition are relatively energetic due to the frequent impingement of mesoscale eddies east of Taiwan and the input of South China Sea water through the Luzon Strait to the onshore side of the Kuroshio, which creates different conditions for the initiation and growth of the interleaving. Comprehensive, multiple platform observations using ships, gliders, and moored instruments are needed to better resolve the evolution of the interleaving of water masses in the Kuroshio.

Methods
Seaglider. The glider used in these observations is a Kongsberg Seaglider, which is equipped with conductivity, temperature, pressure, dissolved oxygen, fluorescence (chlorophyll), and backscatter optical sensors. Only the data collected by the first three sensors (i.e., CTD) were analyzed in this study. Seaglider is a buoyancy-driven autonomous underwater vehicle that is designed to glide from the ocean surface to 1000 m depth in a sawtooth pattern 29 . It steers through the water by controlling its pitch and roll and can navigate between waypoints to execute survey sections. Mission durations depend on payload, stratification, and profile depth but typically range from two to six months. Commanded remotely, gliders are able to acquire GPS position fixes and report their measurements via Iridium satellite telephone while at the sea surface (at the end of each dive cycle). Comparing distance travelled through the water with surface GPS position fixes allows estimation of depth-averaged current velocity. According to the report of Seaglider data quality control process (available at https://gliderfs2.coas.oregonstate.edu/sgliderweb/Seaglider_Quality_Control_Manual.html), the temperature measurement ranges from −2.5 to 43 °C and salinity ranges from 19 to 45. After quality control 30 , the accuracy is 0.001 °C for temperature and 0.01 (0.03 in strong thermocline regions) for salinity. The resolution of the sensor is 0.0001 °C for temperature and 0.001 for salinity. The salinity differences in the interleaving layers (Fig. 3) 0.05-0.10 are thus significant. estimate of glider's underwater location. A Seaglider also comprises built-in motion sensors, which provide pitch, roll, and heading angles while the glider is navigating in the water 29 . During the operation of the glider, a built-in microcomputer will calculate expected horizontal speeds and headings from its current position to the next target point. With the expected quantities of these parameters and the time interval between two adjacent CTD samplings, we can easily calculate the expected north-southward and east-westward displacements of the glider in one dive. When there is a background current, the glider can be drifted by the current and may therefore cause a drifting distance between the designate target point and the realistic position after one dive. The drifting speed, which is equivalent to a depth-averaged current, is obtained using the drifting distance and the time span of this dive. Furthermore, the drifting distance during two adjacent CTD samplings computed from the depth-averaged speed is added to the expected distance to obtain the coordinate of each CTD sampling. With this procedure, the integrated final coordinate should match the last surface position provided by the Global Positioning System (GPS). The coordinate of each salinity sample is used to plot the salinity transects in Fig. 1.

Data.
A Seaglider was launched off southernmost Taiwan on 8 December 2016 and was navigated to (121.51°E, 23.08°N) for hydrographic sampling along a triangular survey track covering the Kuroshio off the east coast of Taiwan (Fig. 1). The glider observations and associated study were sponsored by the Ministry of Science and Technology (MOST) of Taiwan under the Study of the Kuroshio -II (SK-II) project. After launching, the gilder was navigated to the western corner of the triangle track, and from there, the glider began to collect hydrographic data along the three sections (N1, N2, and N3) of the triangle successively (Fig. 1). Two repeated surveys for each transect were completed from 15 December 2016 to 5 March 2017, and 392 dives of hydrographic sampling (1 dive = 1 continuous dive and climb of the glider) were obtained during that period. The sampling period and number of dives at each transect are summarized in Table 3.
The LADCP used in the ship measurement at the KTV1 line (Fig. 1a) comprised two 300 kHz ADCPs with a downward-looking master and an upward-looking slaver, which was deployed with the CTD from the sea surface down to near the ocean bottom and then to the surface. The raw data were processed to be 8 m interval in the vertical for the whole water column (see Jan et al., 2015 for details). The processed data were depth-and time-averaged to obtain the time mean of barotropic velocity components in the zonal and meridional directions, and were further used to calculate horizontal shear of barotropic velocity (v x ) in the Kuroshio region east of Taiwan.
Near-real-time satellite altimeter sea surface height (SSH) and associated absolute geostrophic current data, which are produced by Ssalto/Duacs and distributed by Archiving, Validation, and Interpretation of Satellite Data in Oceanography (AVISO) at http://www.aviso.oceanobs.com/en/data/ products.html, were collected to supplement the interpretation of the glider data.
Diapycnal spiciness curvature. In the water mass analysis, the spiciness is used as a state variable with the unit of kg m −3 , which is a robust and sensitive indicator of how warm (spicy) and salty the water of a known density is, and a diagnosis of double diffusive instability and associated interleaving activity 1,31 . The spiciness is thus widely used to characterize water masses and their diffusive stability. The diapycnal spiciness curvature τ σσ (DSC) is a second derivative of spiciness (τ) with respect to potential density along a profile 1 , which is calculated using where τ is spiciness as defined by the relationship dτ = ρ(αdθ + βdS), σ is potential density, R ρ is the vertical density ratio calculated from d dz dS dz α β θ , α is the thermal expansion coefficient (1.5 × 10 −4 °C −1 ), β is the saline con- Table 3. Periods and sampling cycles of the Seaglider observations along the three sides of the triangular survey track.