Variability of Kuroshio Surface Axis Northeast of Taiwan Island Derived from Satellite Altimeter Data

: The spatial and temporal variability of the Kuroshio surface axis northeast of Taiwan Island is investigated using 24 years of surface geostrophic currents derived from satellite altimeter data from 1993 to 2016. The Kuroshio surface axis is derived by an extraction method with three selected parameters, including the length of the subsidiary line, the intervals between two adjacent points, and the distance between the two adjacent subsidiary lines. The empirical mode decomposition analysis on the 24-year Kuroshio axes reveals that the mean periods of intra-seasonal and inter-annual variability, which are the two dominant components, are about 3.2 months and 1.3 years, respectively. The self-organizing map analysis reveals that the variation of Kuroshio axis northeast of Taiwan Island has four best matching unit (BMU) patterns: straight-path (BMU S ), meandering-path (BMU M ) and two transition stages (BMU T1 and BMU T2 ). The straight-path pattern shows strong seasonality: more likely occurring in summer. The meandering-path pattern is less frequent than straight-path pattern. During a typical period from November 26, 2012 to January 27, 2013, which is chosen as an independent example, the analysis on the satellite altimeter and sea surface temperature data shows that the patterns of the Kuroshio axis change successively in order of BMU T1 → BMU M → BMU T2 → BMU S , i.e., the Kuroshio axis migrates from the meandering-path to the straight-path pattern. During the typical period the warm water intrusion and a mesoscale eddy occur at the second stage corresponding to BMU M and migrate northwestward gradually at the last two stages corresponding to BMU T2 and BMU S . The transient order appears only during this typical period but it is not common for the whole study period. The monthly mean relatively vorticity is calculated and analyzed to evaluate the impact of the eddies on the Kuroshio surface axis variability, the results show that the anticyclonic (cyclonic) eddies can promote the Kuroshio surface axis to present the meandering-path (straight-path) pattern because of the potential vorticity conservation. The impacts of the anticyclonic eddies and the cyclonic eddies on the variability of the Kuroshio surface axis are opposite. The long-term day-to-day detection contributes to improving understanding the variability of Kuroshio surface axis northeast of Taiwan Island.

The KC transports relatively warm water, therefore, variations in the KC position can affect sea surface temperature distribution that moderates the air-sea heat and moisture fluxes [1]. Previous investigations were performed to study the variability of the Kuroshio path northeast of Taiwan Island [5][6][7][8][9]. The Kuroshio surface path always has an anticyclonic meander pattern because of monsoon forcing [5,6]. In spring, the KC branch enters into the ECS continental shelf through the ETC [7,8] and mostly remerges to the Kuroshio main stream [6]. The spatio-temporal variability of the entire Kuroshio path in the ECS from 7-day averaged absolute geostrophic velocity (AGV) data was derived in [9], and they showed that the Kuroshio surface axis moves closer to the 200 m isobath, mostly in winter, and away from it in summer. The intra-seasonal variability of the Kuroshio subsurface water (KSSW) intruding onto the ECS shelf was analyzed based on hydrographic observations and dissolved inorganic iodine species [10], the results showed that the KSSW intrusion northeast of Taiwan Island is weaker in autumn than in spring.
The shoreward intrusion of the KC northeast of Taiwan Island greatly influences the circulation over the ECS shelf and is one of the key mechanism that forms the Taiwan warm current [11][12][13].
Owing to the baroclinic instability, abundant mesoscale eddies are generated and propagate to the western boundary and affect the Kuroshio path and volume transport [14,15]. Delman et al. [16] demonstrated that the vertical stretching of relative vorticity in eddies indicated a deceleration (acceleration) of the Kuroshio Extension coincident with northward (southward) quasi-permanent meanders. On the northeast region of the Taiwan Island, the onshore intrusion of KC is associated with the Kuroshio volume transport and the interaction between the westward propagating openocean eddies and the Kuroshio east of Taiwan Island [17][18][19][20]. Gawarkiewicz et al. [18] and Velez-Belchi et al. [19] reported that the cyclonic eddies can contribute to the intrusion of KC onto the ECS shelf, in contrast, less onshore intrusion occurred during the period that the anticyclonic eddies The KC transports relatively warm water, therefore, variations in the KC position can affect sea surface temperature distribution that moderates the air-sea heat and moisture fluxes [1]. Previous investigations were performed to study the variability of the Kuroshio path northeast of Taiwan Island [5][6][7][8][9]. The Kuroshio surface path always has an anticyclonic meander pattern because of monsoon forcing [5,6]. In spring, the KC branch enters into the ECS continental shelf through the ETC [7,8] and mostly remerges to the Kuroshio main stream [6]. The spatio-temporal variability of the entire Kuroshio path in the ECS from 7-day averaged absolute geostrophic velocity (AGV) data was derived in [9], and they showed that the Kuroshio surface axis moves closer to the 200 m isobath, mostly in winter, and away from it in summer. The intra-seasonal variability of the Kuroshio subsurface water (KSSW) intruding onto the ECS shelf was analyzed based on hydrographic observations and dissolved inorganic iodine species [10], the results showed that the KSSW intrusion northeast of Taiwan Island is weaker in autumn than in spring.
The shoreward intrusion of the KC northeast of Taiwan Island greatly influences the circulation over the ECS shelf and is one of the key mechanism that forms the Taiwan warm current [11][12][13]. Owing to the baroclinic instability, abundant mesoscale eddies are generated and propagate to the western boundary and affect the Kuroshio path and volume transport [14,15]. Delman et al. [16] demonstrated that the vertical stretching of relative vorticity in eddies indicated a deceleration (acceleration) of the Kuroshio Extension coincident with northward (southward) quasi-permanent meanders. On the northeast region of the Taiwan Island, the onshore intrusion of KC is associated with the Kuroshio volume transport and the interaction between the westward propagating open-ocean eddies and the Remote Sens. 2020, 12, 1059 3 of 22 Kuroshio east of Taiwan Island [17][18][19][20]. Gawarkiewicz et al. [18] and Velez-Belchi et al. [19] reported that the cyclonic eddies can contribute to the intrusion of KC onto the ECS shelf, in contrast, less onshore intrusion occurred during the period that the anticyclonic eddies present on the east of Taiwan Island. The mesoscale eddies, which propagate westward to the east region of Taiwan Island, bring potential vorticity flux and modulate local vertical stratification, therefore, they are able to favor or inhibit the Kuroshio intrusion onto the ECS shelf [20].
While these studies have provided important information on the variability of the Kuroshio surface axis and the impact of the eddies, the analysis of the daily detection results of the Kuroshio surface axis using satellite AGV data and a vorticity analysis of eddy activities affecting the Kuroshio axis variability over a long period still need further research. The objectives of this study include: the variability of the Kuroshio axis northeast of Taiwan Island from 1993 to 2016 is analyzed using empirical mode decomposition (EMD) and self-organizing map (SOM) based on long-term day-to-day axis detection results; the impact of the eddies on the Kuroshio surface axis variability is investigated based on the vorticity analysis. The next section describes the data and methods used in this study. In Section 3, the variability of the surface Kuroshio axis northeast of Taiwan Island is quantified from satellite data using extraction method, EMD and SOM; the monthly mean relative vorticity is calculated and analyzed to evaluate the impact of the eddies on the Kuroshio surface axis variability. Section 4 presents the conclusions.

Data
24-year time series of merged AGV data from 1993 to 2016, which are now distributed by the Copernicus Marine Environment Monitoring Service (CMEMS; http://marine.copernicus.eu/), are used in this study. The data are derived from the Ssalto multi-mission ground segment/Data Unification and Altimeter Combination System (Ssalto/Duacs) products which integrates the measurements from altimeter missions including: HY-2A, Saral/AltiKa, Cryosat-2, OSTM/Jason-2, Jason-1, Topex/Poseidon, Envisat, GFO, ERS-1&2. The AGV was calculated from the absolute dynamic topography (ADT), which is the sum of the sea level anomaly (SLA) and the mean dynamic topography (MDT), using the 9-point stencil width method [21]. More detailed information about the altimeter data can be found in [22]. In this study, the delayed-time gridded "allsat" product is chosen. The product is the time series of the merged maps computed with all the satellites available and provide the maps with the best mapping as possible but not homogeneous in time. All of the data can be downloaded from ftp://my.cmems-du.eu/. All users need an account on the website after a registration to obtain an access to a given list of altimetry data. The daily series of the satellite AGV data with a horizontal resolution of 0.25 • × 0.25 • are used to detect the Kuroshio axis; the ADT, AGV, SLA, and the geostrophic velocity anomaly data during a typical period from November 26, 2012 to January 27, 2013 are used for detailed SOM analysis to analyze the transient evolution of the Kuroshio surface axis. In addition, the sea surface temperature (SST) analysis products with a spatial grid resolution of 0.25 • × 0.25 • and temporal resolution of 1 day are also used to study the variation of the Kuroshio surface axis. The products are provided by National Centers for Environmental Information (NCEI; https://www.nodc.noaa.gov/) of National Ocean and Atmospheric Administration (NOAA; https://www.noaa.gov/) and use advanced very high resolution radiometer (AVHRR) satellite data from the Pathfinder AVHRR SST dataset available from September 1981 through December 2005, and the operational Navy AVHRR Multi-Channel SST data from 2006 to the present. More detailed information about the SST products can be found at https://data.nodc.noaa.gov/cgi-bin/iso?id=gov. noaa.nodc:AVHRR_Pathfinder-NCEI-L3C-v5.3#. The AVHRR SST data can be downloaded from ftp://ftp.nodc.noaa.gov/pub/data.nodc/pathfinder/Version5.2/. The structures of the Kuroshio axis and the intrusion are recognized based on the AVHRR SST with higher gradients.

Extraction Method for Kuroshio Surface Axis
The positions of the Kuroshio axis are estimated by tracking the locally strongest AGV in the subsidiary lines. Along each of the subsidiary lines, the local maximum velocity (V 0 ) corresponds to the Kuroshio axis, and the velocities often decrease gradually from the axis to two sides. The method is established according to [23], which has been widely used [20,24,25]. Figure 2 shows the schematic diagram of the extraction method for the Kuroshio axis. At a certain point within the Kuroshio, a subsidiary line is provided crossing the Kuroshio almost perpendicularly. The AGV data are interpolated on the subsidiary line, and the position of the axis is determined as the point corresponding to the local maximum velocity V 0 . A new subsidiary line is then drawn at the new position moved from the upstream axis point by ∆r along the direction of V 0 . The direction of the subsidiary line is further adjusted to be perpendicular to the direction of the new V 0 . The subsidiary line is not always centered on V 0 , but the distance between the midpoint of the subsidiary line and the V 0 is often small, because the relative location of V 0 for the new subsidiary do not change significantly within ∆r. In this study, selection of several parameters is particular; the length of the subsidiary line l is chosen as 140 km to avoid that the subsidiary lines gradually deviate from the Kuroshio axis (locations of V 0 ); the velocities are interpolated to the points of the subsidiary lines, and the intervals between two adjacent points on the subsidiary line ∆l are chosen as 7 km; the distances between two adjacent subsidiary lines ∆r are chosen as 15 km. The extraction method is applied on the 24-year daily gridded AGV data from 1993 to 2016 to identify the locations of the Kuroshio surface axes. the Kuroshio axis, and the velocities often decrease gradually from the axis to two sides. The method is established according to [23], which has been widely used [20,24,25]. Figure 2 shows the schematic diagram of the extraction method for the Kuroshio axis. At a certain point within the Kuroshio, a subsidiary line is provided crossing the Kuroshio almost perpendicularly. The AGV data are interpolated on the subsidiary line, and the position of the axis is determined as the point corresponding to the local maximum velocity V0. A new subsidiary line is then drawn at the new position moved from the upstream axis point by Δr along the direction of V0. The direction of the subsidiary line is further adjusted to be perpendicular to the direction of the new V0. The subsidiary line is not always centered on V0, but the distance between the midpoint of the subsidiary line and the V0 is often small, because the relative location of V0 for the new subsidiary do not change significantly within Δr. In this study, selection of several parameters is particular; the length of the subsidiary line l is chosen as 140 km to avoid that the subsidiary lines gradually deviate from the Kuroshio axis (locations of V0); the velocities are interpolated to the points of the subsidiary lines, and the intervals between two adjacent points on the subsidiary line Δl are chosen as 7 km; the distances between two adjacent subsidiary lines Δr are chosen as 15 km. The extraction method is applied on the 24-year daily gridded AGV data from 1993 to 2016 to identify the locations of the Kuroshio surface axes.

Figure 2.
Schematic diagram showing the method how to detect the Kuroshio axis by tracking the maximum speed of the surface velocity. V0 denote the local maximum velocity, Δr is the distance between two adjacent subsidiary line (before the new subsidiary line being adjusted), l is the length of the subsidiary line, Δl is the interval between two adjacent points.

Empirical Mode Decomposition
The time series of the locations of the Kuroshio surface axis along a transect is analyzed by the EMD method [26]. The EMD was adopted to extract the inter-annual variability of the surface layer heat variance of the Kuroshio intrusion region in the South China Sea and the impact index of the Kuroshio water in southeastern Taiwan Strait [27,28], and used to extract a sea-level index for the Kuroshio Intrusion [29].
A brief introduction of EMD is given below for completeness. The EMD is used to decompose a long-term continuous time series dataset into a finite number of components, the so-called "intrinsic mode functions" (IMFs) and residual (mean trend), in the form of where cj(t) is the jth IMF and r(t) is the residual. Each IMF satisfies two conditions: (1) The number of extrema and the number of zero crossings must be equal or differ by one; (2) the mean value of the upper and lower envelopes connecting all local maxima and minima, respectively, is zero at any data point. The extraction of the IMFs can be simply presented as a sifting process: (1) Locate all local

Empirical Mode Decomposition
The time series of the locations of the Kuroshio surface axis along a transect is analyzed by the EMD method [26]. The EMD was adopted to extract the inter-annual variability of the surface layer heat variance of the Kuroshio intrusion region in the South China Sea and the impact index of the Kuroshio water in southeastern Taiwan Strait [27,28], and used to extract a sea-level index for the Kuroshio Intrusion [29].
A brief introduction of EMD is given below for completeness. The EMD is used to decompose a long-term continuous time series dataset into a finite number of components, the so-called "intrinsic mode functions" (IMFs) and residual (mean trend), in the form of Remote Sens. 2020, 12, 1059 5 of 22 where c j (t) is the jth IMF and r(t) is the residual. Each IMF satisfies two conditions: (1) The number of extrema and the number of zero crossings must be equal or differ by one; (2) the mean value of the upper and lower envelopes connecting all local maxima and minima, respectively, is zero at any data point. The extraction of the IMFs can be simply presented as a sifting process: (1) Locate all local maxima and minima of a data series and connect all maxima and all minima with two cubic spline lines separately to obtain the upper and lower envelops of the data series; (2) calculate the difference between the time series and the mean of the upper and lower envelops to yield a new series; (3) repeating steps (1) and (2) until the upper and the lower envelopes are symmetric with respect to zero mean under the certain criteria, then an IMF, c j (t), is derived; (4) subtract c j (t) from original time series to yield a residual r(t) and treat r(t) as the original time series and repeat steps (1)-(3) until r(t) becomes a monotonic function or a function with only one extremum. The EMD is adaptive and consists of independent IMFs, each of which has an orthogonal amplitude and frequency. For interpretation of temporally changing phenomena, the IMFs are more physically meaningful than the empirical orthogonal function (EOF) modes [26,29]. After the IMFs are obtained using the EMD method, we apply the Hilbert Transform (HT) to each component to calculate the instantaneous frequency by where X(t) is an IMF vector, ⊗ denotes convolution operator, P represents the Cauchy principal value. With this expression, the corresponding analytic signal, Z(t), can be expressed as in which a(t) is the amplitude and can be written as θ(t) = tan −1 HT(t) X(t) is the phase. The HT provides a unique way for defining the imaginary part, so that the result is an analytic function [26]. The instantaneous frequency for the IMF can be defined as In this study, the EMD method is performed on the 24-year daily series of the distance, which is between the locations of the daily axes and the location of the mean axis along Transect A, to generate variability components with different time scales.

Self-Organizing Map Analysis
With the derived Kuroshio surface axis locations, the SOM is used to cluster the Kuroshio axis. On the basis of unsupervised neural network, the SOM is an effective method for feature extraction and classification, and maps high-dimensional input data onto the elements of a regular, low-dimensional array to extract data features [30]. It is able to effectively preserve input data information because of the topology-preserving technique [31].
The SOM has been demonstrated to be powerful for feature extractions, especially when the signal is highly nonlinear [32,33]. It has been applied in the study of the interaction between the East China Sea Kuroshio and the Ryukyu Current [34], the variability of the Kuroshio intrusion through Luzon Strait [35], and the seasonal variability of Kuroshio intrusion northeast of Taiwan Island [31]. These applications suggest that the SOM may be a better way to cluster the Kuroshio axis information.
Prior to the neural network training process, the tunable SOM parameters need to be specified. Based on a brief analysis of the Kuroshio axis information, the parameters such as map shape, initial neighborhood size, layer topology function neuron distance function are specified. The SOM analysis is performed using selforgmap neural network toolbox in Matlab. Different SOM array sizes (1 × 2, 1 × 3, 2 × 2, 2 × 3, 3 × 3, 4 × 4, 5 × 5, and 6 × 6) are tested to evaluate the stability of the quantization error in the SOM analysis. A most optimized SOM array size can be determined based on the comparisons. The number of training steps and the initial neighborhood size are chosen as 100 and 3, respectively. The neurons in the layer of a SOM are arranged originally in positions according to a topology function. There are three functions, including grid, hexagonal, and random topology, in the Matlab selforgmap toolbox [36]. The hexagonal layer topology function is adopted in the classification process. The distance between neurons are computed from their positions with four distance functions, i.e., Euclidean, Manhattan, link, and box functions. The link function computes the number of links as distance, and the box function computes the layer neuron position vector based on distance. In previous studies on the SOM clustering technique application, the distance from a particular neuron to its neighbors has often been calculated using Euclidean [34,37,38] and link functions [39,40]. Comparisons of SOM results using the four distance functions [38,41,42] show that there is no "golden rule" to select a best distance function because each of the three functions, i.e., Euclidean, Manhattan, and link distance, is able to work better than others in different experiments. According to [36,43], the Euclidean and link distance functions are two of the most common for SOM analysis. In this study, the link function is adopted.
In this study, the input data for the SOM model is a matrix D M×N = {d 1 , d 2 , · · · , d N }, where d i denotes the vector containing the latitude information of the daily Kuroshio surface axis. The latitude of each axis is interpolated to the same meridional intervals from 122 • E to 128 • E. The spacing length is 0.1 • . In this study, the length of the vector d i (M) is 61, the number of the vectors (N) is 8737 because 41 failed extraction results are removed. In the learning stage of the SOM process, the data is introduced to SOM lattice with k neurons and the winning neurons (C 1 . . . C k ) are updated to match the input data [44]. For this reason the winning neurons are often referred to as best matching units (BMU). According to [38,44], the winning element error (WEE i ) can be expressed as where .× stands for the Hadamard (elementwise) product of vectors, S i is the neuron which is associated with the data vector d i with the shortest Euclidian distance and expressed as S i = arg min The SOM method is used to cluster the derived 24-year daily series of the Kuroshio surface axis location information.

General Features of the Kuroshio Mean Axis
The extraction method is applied on the 24-year AGV data to identify the location of the Kuroshio surface axis. Figure 3 shows mean ADT (color-filled contours), the detected Kuroshio mean axis (solid black curve), standard deviation (STD, dashed blue curves), envelopes (dashed cyan curves), 90 percentile (dashed green curve), and 10 percentile (dashed purple curve) from 1993 to 2016. The envelopes represent the farthest (closest) location that the KC reached the open ocean (shelf) over the study duration. The STD is a measure of how scattering the data are and can be calculated as the square root of the variance, which is the averaged squared deviations from the mean values.
The extraction method is applied on the 24-year AGV data to identify the location of the Kuroshio surface axis. Figure 3 shows mean ADT (color-filled contours), the detected Kuroshio mean axis (solid black curve), standard deviation (STD, dashed blue curves), envelopes (dashed cyan curves), 90 percentile (dashed green curve), and 10 percentile (dashed purple curve) from 1993 to 2016. The envelopes represent the farthest (closest) location that the KC reached the open ocean (shelf) over the study duration. The STD is a measure of how scattering the data are and can be calculated as the square root of the variance, which is the averaged squared deviations from the mean values.  From east of Taiwan Island to the Kuroshio middle segment in the ECS at about 128 • E, the Kuroshio mean axis generally moves between 1.0 m and 1.2 m contours of the mean ADT. The STD is obviously smaller than the envelopes, implying that relatively small shifts of the Kuroshio axis occur frequently. For the Kuroshio middle segment in the ECS between about 125 • E and 127.5 • E, the Kuroshio axis appears nearly as a straight line, because the KC flows more stably and more straightly compared with other segments in the ECS [6,45]. The comparisons among STD, envelope, 90 and 10 percentile indicate that most of the Kuroshio axis shift with relatively small amplitude because 90 and 10 percentile is significantly closer to the STD than envelope. This implies that 80% of the daily detection results of the Kuroshio axis (between 90 and 10 percentile curves) are far away from the envelopes, while only 10% of those (between 90 percentile and upper envelope curves) correspond to the large anticyclonic turning, which indicate the Kuroshio intrusion into the ECS [7,31,46]. The reason why the 90 and 10 percentile curves are far away from the envelope curves is that relatively small axis fluctuations appear in most of the study period from 1993 to 2016. Furthermore, the envelopes represent extreme events and are part of 10% of the data in this study, so they do not represent most of the Kuroshio axis variation. Another reason is that the extraction method in this study detects the Kuroshio axis based on the local maximum velocity along a subsidiary line, it is hard to show the whole Kuroshio surface current field. Therefore, the detected Kuroshio axes will more likely be stable, which makes the fluctuation amplitudes of 80% of the Kuroshio axis (between 90 and 10 percentile curves) significantly smaller than the envelopes.

Variation of the Kuroshio Surface Axis
To further quantify the variation of the Kuroshio axis, the daily distances away from the location of the mean axis along Transect A are calculated. The positive values denote shoreward distances and negative values denote seaward distances. The distances of zero correspond to the intersection between the mean Kuroshio axis and Transect A. The top panel of Figure 4 shows the Kuroshio axis location relative to the mean axis along Transect A. The EMD method is performed on the 24-year time series of the distance to generate 12 IMFs and a residual (Figure 4), which represent the variability components of the surface Kuroshio axis with different time scales. We apply the discrete Hilbert transform (using Matlab HILBERT toolbox, i.e., "Hilbert" function) to calculate the instantaneous frequency for all of the IMFs by Equation (5). The mean period T m can be expressed as where ω(t) is the instantaneous frequency at time t, ω m represents the mean frequency. The energy E that represents the intensity of the signal can be calculated as where a(t) is the instantaneous amplitude calculated by Equation (4).  The error bars in Figure 5a represent the mean value and STD of the frequency and amplitude for each of the IMFs. The mean periods and energy of the IMFs are shown in Figure 5b. From Figure  5 one can see that the mean amplitude and energy of IMF4-IMF9 are larger than other IMFs, which means that the variability of the Kuroshio axis corresponding to IMF4-IMF9 is dominant during the study period from 1993 to 2016. IMF1-IMF4 have mean periods of 5.4, 6.6, 10.1, and 20.5 days, respectively, i.e., short-period fluctuations below the monthly timescale. As the ratio of large fluctuations for IMF4 is relatively higher compared with IMF1-IMF3, this means that the intensity of The error bars in Figure 5a represent the mean value and STD of the frequency and amplitude for each of the IMFs. The mean periods and energy of the IMFs are shown in Figure 5b. From Figure 5 one can see that the mean amplitude and energy of IMF4-IMF9 are larger than other IMFs, which means that the variability of the Kuroshio axis corresponding to IMF4-IMF9 is dominant during the study period from 1993 to 2016. IMF1-IMF4 have mean periods of 5.4, 6.6, 10.1, and 20.5 days, respectively, i.e., short-period fluctuations below the monthly timescale. As the ratio of large fluctuations for IMF4 is relatively higher compared with IMF1-IMF3, this means that the intensity of the Kuroshio axis variation with a mean period of 20.5 days is stronger than the other shorter-period variations. IMF5-IMF6 correspond to the monthly variation with mean periods of 31.6 and 51.8 days, respectively. The mean periods of IMF7-IMF8 are 96.4 and 191.3 days, which indicate the intra-seasonal and inter-annual variations. IMF9-IMF11 have mean periods of 1.3, 2.1, and 6.1 year, respectively, i.e., the inter-annual variations. The mean amplitude and energy of IMF9 are much larger than IMF10 and IMF11 ( Figure 5), implying that the 1.3-year variation is the relatively important component among the three inter-annual variations. The longest mean period as long as 11.6 year occurs in IMF12 corresponding to decadal variation. The residual defines the variation trend of the Kuroshio surface axis, one can see that the Kuroshio axis northeast of Taiwan Island started with a seaward pattern in 1993, and gradually reversed toward the ECS continental shelf. The shoreward distance reached a maximum during the period from 2006 to 2008. Then the Kuroshio axis began to move back to the Pacific again. As the mean amplitude and energy of IMF7 and IMF9 are the largest of all the IMFs, implying that the intra-seasonal and inter-annual variations with mean periods of about 3.  To examine the intra-seasonal variation of the Kuroshio surface axis, 24-year daily mean series of IMF4-IMF9, which are the dominant components of the variation, are shown in Figures 6. From Figures 6 one can see that remarkable seaward-shoreward oscillations of IMF4-IMF9 occur throughout the year. For the variation corresponding to IMF4, the fluctuation is the strongest in winter, next in spring and autumn, and weakest in summer for the Northern Hemisphere. For the monthly variations, the amplitude of IMF6 is generally larger than IMF5. IMF6 shows that the Kuroshio axis presents shoreward state distinctly in first half of June, second half of November and early December, while seaward state in most of March, late June, most of July, September and second half of December. Moreover, there are three crosses, at which the oscillation direction reverses, i.e., from shoreward to seaward in June and December, and from seaward to shoreward in March. For the intra-seasonal variation corresponding to IMF7, the fluctuation occurs frequently in winter, and the seasonal characteristics are weaker than that of IMF8, implying an extremely distinct intraseasonal variation with a mean period of almost half a year. The intra-seasonal variability of IMF8 and IMF9 occurs quasi-periodically, implying the periodical variability of Kuroshio intrusion into the ECS. From IMF8, it can be seen that the Kuroshio axis shifts seaward in summer and shoreward in other seasons, because the KC, on the one hand, flows eastward directly with weak intrusion onto the shelf in summer, on the other hand, the Kuroshio intrusion is strong in other seasons [6]. The intra-seasonal variability of IMF9 is similar to IMF8. The intra-seasonal variation of the Kuroshio axis is consistent with that of the Kuroshio intrusion, which is related to the disposition of the current fields induced by monsoon [6,9]. To examine the intra-seasonal variation of the Kuroshio surface axis, 24-year daily mean series of IMF4-IMF9, which are the dominant components of the variation, are shown in Figure 6. From Figure 6 one can see that remarkable seaward-shoreward oscillations of IMF4-IMF9 occur throughout the year. For the variation corresponding to IMF4, the fluctuation is the strongest in winter, next in spring and autumn, and weakest in summer for the Northern Hemisphere. For the monthly variations, the amplitude of IMF6 is generally larger than IMF5. IMF6 shows that the Kuroshio axis presents shoreward state distinctly in first half of June, second half of November and early December, while seaward state in most of March, late June, most of July, September and second half of December. Moreover, there are three crosses, at which the oscillation direction reverses, i.e., from shoreward to seaward in June and December, and from seaward to shoreward in March. For the intra-seasonal variation corresponding to IMF7, the fluctuation occurs frequently in winter, and the seasonal characteristics are weaker than that of IMF8, implying an extremely distinct intra-seasonal variation with a mean period of almost half a year. The intra-seasonal variability of IMF8 and IMF9 occurs quasi-periodically, implying the periodical variability of Kuroshio intrusion into the ECS. From IMF8, it can be seen that the Kuroshio axis shifts seaward in summer and shoreward in other seasons, because the KC, on the one hand, flows eastward directly with weak intrusion onto the shelf in summer, on the other hand, the Kuroshio intrusion is strong in other seasons [6]. The intra-seasonal variability of IMF9 is similar to IMF8. The intra-seasonal variation of the Kuroshio axis is consistent with that of the Kuroshio intrusion, which is related to the disposition of the current fields induced by monsoon [6,9]. It is worth noting that the EMD is applied only for the time series of the Kuroshio axis along Transect A. The variability along Transect A cannot represent the whole Kuroshio axis. However, Transect A is defined in Figure 3 for one of the main objectives of this study, which is studying the relation between the Kuroshio surface axis variation northeast of Taiwan Island and the onshore intrusion. Therefore, the EMD analysis along Transect A should be relatively representative in this study.
Two transects (Transects B and C shown in Figure 3) are used to analyze the relationship between the relatively upstream Kuroshio axis east of Taiwan Island (along Transect B) and the Kuroshio axis northeast of Taiwan Island corresponding to Kuroshio intrusion (along Transect C). The two transects are independent with Transect A. Each of the daily Kuroshio axis intersects Transect B and C at two points, the longitudes along Transect B versus the latitudes along Transect C are shown in Figure 7. The color scale represents the number of intersections, namely, the data density. One can see that there are two high-value centers of about 50 and 60 in the south and southeast of the intersection between the mean Kuroshio axis and the two transects (black pentagram in Figure 7). The results suggest that the Kuroshio flows directly eastward in northeast of Taiwan Island with weak intrusion for most the time. The intersections in the red ellipse shown in Figure 7 should correspond to the strong Kuroshio intrusion, because the intersections appear in the quite far north regions. Furthermore, it can be seen that when the Kuroshio intrudes into the ECS, the longitudes of the intersections along Transect C are generally in an interval from 122.4° E to 122.6° E, in other words, the KC flows closely to Taiwan Island in the upstream region east of Taiwan Island when the strong Kuroshio intrusion occurs. It is worth noting that the EMD is applied only for the time series of the Kuroshio axis along Transect A. The variability along Transect A cannot represent the whole Kuroshio axis. However, Transect A is defined in Figure 3 for one of the main objectives of this study, which is studying the relation between the Kuroshio surface axis variation northeast of Taiwan Island and the onshore intrusion. Therefore, the EMD analysis along Transect A should be relatively representative in this study.
Two transects (Transects B and C shown in Figure 3) are used to analyze the relationship between the relatively upstream Kuroshio axis east of Taiwan Island (along Transect B) and the Kuroshio axis northeast of Taiwan Island corresponding to Kuroshio intrusion (along Transect C). The two transects are independent with Transect A. Each of the daily Kuroshio axis intersects Transect B and C at two points, the longitudes along Transect B versus the latitudes along Transect C are shown in Figure 7. The color scale represents the number of intersections, namely, the data density. One can see that there are two high-value centers of about 50 and 60 in the south and southeast of the intersection between the mean Kuroshio axis and the two transects (black pentagram in Figure 7). The results suggest that the Kuroshio flows directly eastward in northeast of Taiwan Island with weak intrusion for most the time. The intersections in the red ellipse shown in Figure 7 should correspond to the strong Kuroshio intrusion, because the intersections appear in the quite far north regions. Furthermore, it can be seen that when the Kuroshio intrudes into the ECS, the longitudes of the intersections along Transect C are generally in an interval from 122.4 • E to 122.6 • E, in other words, the KC flows closely to Taiwan Island in the upstream region east of Taiwan Island when the strong Kuroshio intrusion occurs.

Self-Organizing Map Analysis
With the derived Kuroshio axis location information, the SOM method is used to cluster the Kuroshio axis over the last 24 years. First, the adequacy of the SOM array size is evaluated using the relative error, i.e., the coefficient of variation, of the intra-class sum of squares (CVSSIntra) [44,47]. Evolution of CVSSIntra shows one minimum (Figure 8) for 2×2 SOM array, implying that the most optimized SOM array size is 2×2. Therefore, the SOM analysis identifies four BMU patterns.  Figure 9 shows the mean axes for the four BMU patterns. BMUM and BMUS correspond to meandering-path and straight-path patterns. BMUT1 and BMUT2 are regarded as transition stages between BMUM and BMUS. The mean axis of BMUM shows a northward flow across the isobaths from 122.1° E to 123.0° E and entering the ECS, then deflects southeastward at 123.5° E. BMUT2 is similar to BMUM at the beginning, but the mean axis deflects southeastward as soon as it flows across the 200 m isobaths. Subsequently, the mean axis of BMUT2 remerges to that of BMUS corresponding to straight-path state. BMUT1 is closer to the Kuroshio mean axis than other three from 122.1° E to 122.8° E, continuously flows northwestward and finally remerges to the axis of BMUM. For BMUS, the mean  Figure 7. Scatter diagram of the longitudes of the intersections along Transect B versus latitudes of the intersections along Transect C for the long-term day-to-day detection of the Kuroshio axis from 1993 to 2016, Transects B and C are shown in Figure 3, the color codes represent the numbers of intersections, i.e., the data density. Black pentagram corresponds to the mean axis at (25.43 • N, 122.59 • E). Two black dashed lines are zonal and meridional extension lines. Latitudes along Transect C within the red ellipse represent the Kuroshio intrusion events.

Self-Organizing Map Analysis
With the derived Kuroshio axis location information, the SOM method is used to cluster the Kuroshio axis over the last 24 years. First, the adequacy of the SOM array size is evaluated using the relative error, i.e., the coefficient of variation, of the intra-class sum of squares (CV SSIntra ) [44,47]. Evolution of CVSSIntra shows one minimum ( Figure 8) for 2 × 2 SOM array, implying that the most optimized SOM array size is 2 × 2. Therefore, the SOM analysis identifies four BMU patterns.

Self-Organizing Map Analysis
With the derived Kuroshio axis location information, the SOM method is used to cluster the Kuroshio axis over the last 24 years. First, the adequacy of the SOM array size is evaluated using the relative error, i.e., the coefficient of variation, of the intra-class sum of squares (CVSSIntra) [44,47]. Evolution of CVSSIntra shows one minimum ( Figure 8) for 2×2 SOM array, implying that the most optimized SOM array size is 2×2. Therefore, the SOM analysis identifies four BMU patterns.  Figure 9 shows the mean axes for the four BMU patterns. BMUM and BMUS correspond to meandering-path and straight-path patterns. BMUT1 and BMUT2 are regarded as transition stages between BMUM and BMUS. The mean axis of BMUM shows a northward flow across the isobaths from 122.1° E to 123.0° E and entering the ECS, then deflects southeastward at 123.5° E. BMUT2 is similar to BMUM at the beginning, but the mean axis deflects southeastward as soon as it flows across the 200 m isobaths. Subsequently, the mean axis of BMUT2 remerges to that of BMUS corresponding to straight-path state. BMUT1 is closer to the Kuroshio mean axis than other three from 122.1° E to 122.8° E, continuously flows northwestward and finally remerges to the axis of BMUM. For BMUS, the mean  Figure 8. Coefficient variation of intra-class sum of the squared errors (CV SSIntra ) with different self-organizing map (SOM) array sizes. Figure 9 shows the mean axes for the four BMU patterns. BMU M and BMU S correspond to meandering-path and straight-path patterns. BMU T1 and BMU T2 are regarded as transition stages between BMU M and BMU S . The mean axis of BMU M shows a northward flow across the isobaths from 122.1 • E to 123.0 • E and entering the ECS, then deflects southeastward at 123.5 • E. BMU T2 is similar to BMU M at the beginning, but the mean axis deflects southeastward as soon as it flows across the 200 m isobaths. Subsequently, the mean axis of BMU T2 remerges to that of BMU S corresponding to straight-path state. BMU T1 is closer to the Kuroshio mean axis than other three from 122.1 • E to 122.8 • E, continuously flows northwestward and finally remerges to the axis of BMU M . For BMU S , the mean axis shows a northeastward flow at the beginning and deflects eastward with weak intrusion onto the shelf. For each of the BMU patterns, the temporal evolution of the WEE i is shown in Figure 10, where the subscript i denotes the day number from 1993 to 2016. From Figure 10, one can see that the errors are significantly larger between 122 • E and 123.5 • E than other intervals for all of the four BMU patterns, indicating the high-amplitude variation of the Kuroshio surface axis northeast of Taiwan Island. At the end of 2001 and the beginning of 2005, there are relatively high errors for BMU S between 122 • E and 123 • E because the daily Kuroshio axes are located distantly to the south of the mean axis and close to the lower envelope curve (Figure 9). It indicates that the KC deflects directly northeastward at about 23 • N during the particular time period (the end of 2001 and the beginning of 2005). However, the special situation occurs in only less than 20 days during the whole 24-year study time period and is located on the east of the Taiwan Island which is not mainly focused in this study, more detailed analysis for the situation and mechanism will be carried out in future. . However, the special situation occurs in only less than 20 days during the whole 24-year study time period and is located on the east of the Taiwan Island which is not mainly focused in this study, more detailed analysis for the situation and mechanism will be carried out in future.    (Figure 9). It indicates that the KC deflects directly northeastward at about 23° N during the particular time period (the end of 2001 and the beginning of 2005). However, the special situation occurs in only less than 20 days during the whole 24-year study time period and is located on the east of the Taiwan Island which is not mainly focused in this study, more detailed analysis for the situation and mechanism will be carried out in future.  To quantify the occurrence percentage of each pattern, an index of the occurrence frequency (IOF) is calculated by summing the occurrence number of each pattern divided by the total record length [33]. Over the 24 years, the IOFs of BMU S , BMU M , BMU T1 , and BMU T2 are 61.7%, 10.0%, 17.6%, and 10.7%, respectively ( Figure 9). One can see that the IOF of meandering-path pattern BMU M occupies only about 1/6 of that of straight-path patterns BMU S . In other words, the strong deflection of the Kuroshio axis corresponding to BMU M does not occur frequently. This conclusion is consistent with the analysis results from Figure 7.
To examine the most probable case of BMU S (large IOF of 61.7%), the principal component analysis (PCA) [48,49] is used for the detected results of the Kuroshio surface axis. The detailed analysis results based on the PCA are described in Appendix A, which is able to corroborate the SOM results that the straight-path pattern more likely occurs than other three patterns.
The monthly IOFs for the four patterns are shown in Figure 11a. The monthly IOFs are calculated by summing the occurrence number of each BMU pattern per climatological months divided by the total record length. For each of the BMU patterns, the sum of the monthly IOFs over 12 months is equal to the IOF values, as shown in Figure 8. The monthly IOFs of BMU S are higher than those of other patterns in all twelve months, and show strong seasonality: higher values in summer, and lower values in spring and autumn for the Northern Hemisphere. The IOFs of BMU M and BMU T2 are highest in winter and spring, while lowest in summer. Unlike BMU T2 , the other transition pattern, BMU T1 , has larger IOFs from April to November. BMU T1 shows relatively weak seasonality, because there is noticeable increase from January to May and from July to October and decrease from May to July and from October to December. To quantify the occurrence percentage of each pattern, an index of the occurrence frequency (IOF) is calculated by summing the occurrence number of each pattern divided by the total record length [33]. Over the 24 years, the IOFs of BMUS, BMUM, BMUT1, and BMUT2 are 61.7%, 10.0%, 17.6%, and 10.7%, respectively (Figure 9). One can see that the IOF of meandering-path pattern BMUM occupies only about 1/6 of that of straight-path patterns BMUS. In other words, the strong deflection of the Kuroshio axis corresponding to BMUM does not occur frequently. This conclusion is consistent with the analysis results from Figure 7.
To examine the most probable case of BMUS (large IOF of 61.7%), the principal component analysis (PCA) [48,49] is used for the detected results of the Kuroshio surface axis. The detailed analysis results based on the PCA are described in Appendix A, which is able to corroborate the SOM results that the straight-path pattern more likely occurs than other three patterns.
The monthly IOFs for the four patterns are shown in Figure 11a. The monthly IOFs are calculated by summing the occurrence number of each BMU pattern per climatological months divided by the total record length. For each of the BMU patterns, the sum of the monthly IOFs over 12 months is equal to the IOF values, as shown in Figure 8. The monthly IOFs of BMUS are higher than those of other patterns in all twelve months, and show strong seasonality: higher values in summer, and lower values in spring and autumn for the Northern Hemisphere. The IOFs of BMUM and BMUT2 are highest in winter and spring, while lowest in summer. Unlike BMUT2, the other transition pattern, BMUT1, has larger IOFs from April to November. BMUT1 shows relatively weak seasonality, because there is noticeable increase from January to May and from July to October and decrease from May to July and from October to December. To explore a possible relationship between the upstream surface current transport along Transect B and the Kuroshio intrusion, the daily surface current transport time series are calculated from 24-year AGV data and expressed as where TC is the surface transport along Transect B, U is the AGV and actually denote the meridional component because of the zonal Transect B, ds is the distance between two neighboring points across the width of the Kuroshio, p1 and p2 represent two endpoints of Transect B. At the monthly scale, the IOF of BMUS (straight-path pattern) is positively correlated with the monthly mean value of TC (R=0.78; Figure 11b). The statistical results indicate that TC is larger in summer months, and the To explore a possible relationship between the upstream surface current transport along Transect B and the Kuroshio intrusion, the daily surface current transport time series are calculated from 24-year AGV data and expressed as where T C is the surface transport along Transect B, U is the AGV and actually denote the meridional component because of the zonal Transect B, ds is the distance between two neighboring points across the width of the Kuroshio, p 1 and p 2 represent two endpoints of Transect B. At the monthly scale, the IOF of BMU S (straight-path pattern) is positively correlated with the monthly mean value of T C (R = 0.78; Figure 11b). The statistical results indicate that T C is larger in summer months, and the corresponding IOF of BMU S is higher, i.e., the straight-path pattern more likely occurs when the upstream surface transport increases especially in summer.
To examine the transient evolution of the four Kuroshio axis patterns, we choose a typical period from November 26, 2012 to January 27, 2013. An important reason for choosing this period is that the four Kuroshio axis BMU patterns occur alternately during a short time (two months). Figure 12 shows the time series of the BMU patterns during the typical period. The period is divided into four stages, in which the patterns of the Kuroshio axis change successively in order of BMU T1 (from November 26 to December 11, 2012) → BMU M (from December 12, 2012 to January 1, 2013) → BMU T2 (from January 2 to January 20, 2013) → BMU S (from January 21 to January 27, 2013). It is worth noting that the evolution order appears only for the typical period and is not common during the whole study period. One of the limitation of the SOM analysis in this study is that it is too hard to make generalization about the transient evolution based on only one single cycle. Therefore, more detailed analysis for the transient evolution needs to be performed in future. The corresponding ADT and geostrophic currents, SLA and geostrophic velocity anomaly, and AVHRR SST are shown in Figures 13-15 to analyze the continuous evolution of the Kuroshio intrusion and the variability of the surface axis. The geostrophic currents show circulation structures in good agreement with the SST images ( Figure 14). The results suggest that sometimes there are significant anticyclonic intrusions (Figures 13b and 14b)  corresponding IOF of BMUS is higher, i.e., the straight-path pattern more likely occurs when the upstream surface transport increases especially in summer.
To examine the transient evolution of the four Kuroshio axis patterns, we choose a typical period from November 26, 2012 to January 27, 2013. An important reason for choosing this period is that the four Kuroshio axis BMU patterns occur alternately during a short time (two months). Figure 12 shows the time series of the BMU patterns during the typical period. The period is divided into four stages, in which the patterns of the Kuroshio axis change successively in order of BMUT1 (from November 26 to December 11, 2012) → BMUM (from December 12, 2012 to January 1, 2013) → BMUT2 (from January 2 to January 20, 2013) → BMUS (from January 21 to January 27, 2013). It is worth noting that the evolution order appears only for the typical period and is not common during the whole study period. One of the limitation of the SOM analysis in this study is that it is too hard to make generalization about the transient evolution based on only one single cycle. Therefore, more detailed analysis for the transient evolution needs to be performed in future. The corresponding ADT and geostrophic currents, SLA and geostrophic velocity anomaly, and AVHRR SST are shown in Figures  13-15   corresponding IOF of BMUS is higher, i.e., the straight-path pattern more likely occurs when the upstream surface transport increases especially in summer.
To examine the transient evolution of the four Kuroshio axis patterns, we choose a typical period from November 26, 2012 to January 27, 2013. An important reason for choosing this period is that the four Kuroshio axis BMU patterns occur alternately during a short time (two months). Figure 12 shows the time series of the BMU patterns during the typical period. The period is divided into four stages, in which the patterns of the Kuroshio axis change successively in order of BMUT1 (from November 26 to December 11, 2012) → BMUM (from December 12, 2012 to January 1, 2013) → BMUT2 (from January 2 to January 20, 2013) → BMUS (from January 21 to January 27, 2013). It is worth noting that the evolution order appears only for the typical period and is not common during the whole study period. One of the limitation of the SOM analysis in this study is that it is too hard to make generalization about the transient evolution based on only one single cycle. Therefore, more detailed analysis for the transient evolution needs to be performed in future. The corresponding ADT and geostrophic currents, SLA and geostrophic velocity anomaly, and AVHRR SST are shown in Figures  13-15 Figure 13, but for AVHRR multi-day mean SST. Figure 14. Same as Figure 13, but for AVHRR multi-day mean SST.   For the first stage of the time interval corresponding to the transient pattern BMU T1 , small meander occurs at about 123 • E ( Figure 13a); however, the Kuroshio axis is generally close to straight-path pattern BMU S , and the intrusion of the surface warm water into the ECS between 25 • N and 26 • N and between 121.5 • E and 122.5 • E is not distinct (Figure 14a). From Figure 15a one can see that there is a cyclonic eddy between 26.5 • N and 27.5 • N and between 121.5 • E and 122.5 • E, of which the southmost edge is adjacent to the KC, the northeastward current anomalies northeast of Taiwan along the Kuroshio axis indicate that the meander of the Kuroshio axis, or the Kuroshio intrusion, may be strengthened.
For the second stage of the time interval corresponding to meandering-path pattern BMU M , the large meandering state and the corresponding surface geostrophic currents (Figure 13b) indicate that a great volume of Kuroshio warm water enters into the ECS shelf (Figure 14b), meanwhile, an anticyclonic eddy occurs between 25.5 • N and 27 • N and between 122 • E and 123.5 • E (Figure 15b).
Next, another transient pattern BMU T2 presents at the third stage of the period. The Kuroshio axis between about 123 • E and 125 • E starts to migrate southward compared with BMU M (Figure 13c) and becomes close to the axis of straight-path pattern BMU S , implying that the Kuroshio intrusion becomes weaker. The transition of the Kuroshio axis and the geostrophic currents from BMU M to BMU T2 is consistent with the SST images, which indicates that the warm water intrusion (Figure 14b,c) migrates westward. Similarly, the center of the anticyclonic eddy of BMU T2 (Figure 15c) indicates a westward migration, implying the state of the Kuroshio axis change from meandering-path (BMU M ) to transient (BMU T2 ).
The final stage of the time interval corresponds to straight-path pattern BMU S . The Kuroshio axis migrates southward generally compared with other patterns. At this stage, the Kuroshio surface water still intrudes onto the ECS shelf which is shown as a warm water tongue in Figure 14d. Meanwhile, the eddy appearing at the third stage gradually flows northward and away from the main stream of the Kuroshio (Figure 13c,d). It means that the Kuroshio intrusion still occurs in boreal winter when the surface axis corresponds to straight-path pattern. The reason should be that there is the Kuroshio branch current (KBC) on the northeast of Taiwan Island especially in winter. The KBC is referred to as a quasi-steady current [31]. On intra-seasonal time scale, the intrusion events of the Kuroshio are more frequent in winter than in summer [50,51]. Therefore, the straight-path pattern in winter should not entirely correspond to the KC without intrusion or even weak intrusion. Since the extraction method in this study detects the Kuroshio axis based on the local maximum velocity along a subsidiary line, it is hard to reveal the KBC, i.e., the secondary Kuroshio surface axis. However, the ADT, SLA, and SST images (Figures 13-15) are able to show them obviously. Further work will be directed toward a more detailed analysis of the Kuroshio surface currents from more remote sensing data and the three-dimensional structure of the Kuroshio from numerical modelling results.

Impact of Mesoscale Eddies on the Kuroshio Surface Axis Variability
As mentioned in Section 1, the westward propagating open-ocean eddies can affect the variability of the Kuroshio surface axis. Figure 16 shows the monthly mean geostrophic velocity anomaly and the relative vorticity, ζ = ∂v ∂x − ∂u ∂y , in January and July, 2009 as an example. The reason why the monthly mean results in 2009 is selected is that the impacts of eddies are relatively obvious. January and July correspond to boreal winter and summer, respectively. The first-order principle of potential vorticity (PV) conservation, which is expressed as PV = f + ζ H = constant, is considered, where f is the Coriolis parameter, H is the water depth. In general, because of the difference of vertical stratification, the PV in the upper layer is larger on the shelf than in the open ocean and decreases gradually offshore, the strong PV gradient is generated across the slope and often blocks most of the Kuroshio water intrusion [20,52]. Over the small geographic region being considered, f changes little. Under the impact of anticyclonic eddies, the monthly mean ζ is negative during January 2009 in the region where the interaction between the eddies and Kuroshio often occurs (black square in Figure 16a), it indicates that the net relative vorticity decreases, then the water depth H must decrease in order to conserve the PV. Therefore, the water moves onshore because the water depth on the ECS shelf is shallower and the Kuroshio surface axis presents BMU M pattern (green line in Figure 16a). In contrast, in July 2009, the monthly mean ζ is positive in the same region because of the cyclonic eddies, then the water moves offshore and the Kuroshio surface axis presents BMU S pattern (Figure 16b). The impacts of the anticyclonic and cyclonic eddies are opposite. region where the interaction between the eddies and Kuroshio often occurs (black square in Figure  16a), it indicates that the net relative vorticity decreases, then the water depth H must decrease in order to conserve the PV. Therefore, the water moves onshore because the water depth on the ECS shelf is shallower and the Kuroshio surface axis presents BMUM pattern (green line in Figure 16a). In contrast, in July 2009, the monthly mean ζ is positive in the same region because of the cyclonic eddies, then the water moves offshore and the Kuroshio surface axis presents BMUS pattern (Figure 16b). The impacts of the anticyclonic and cyclonic eddies are opposite. The correlation coefficient between the monthly mean ζ and the IOFs of BMUS is 0.82 during the study period which is from 1993 to 2016, implying that the mean relative vorticity increases in summer months, and the corresponding IOF of BMUS is higher, i.e., the straight-path pattern more likely occurs when the relative vorticity increases especially in summer. In contrast, from Figure 17b one can see that the meandering-path pattern of the Kuroshio axis more likely occurs when the relative vorticity decreases. In conclusion, for the multi-year monthly mean scale, the anticyclonic eddies can carry negative PV flux and promote the Kuroshio axis to present the meandering-path pattern, while the cyclonic eddies can carry positive PV flux and promote the Kuroshio axis to present the straight-path pattern. The results are similar to the analysis of Figures 16 and the dynamical interpretations in [18][19][20].   Figure 17 indicate that the multi-year monthly mean ζ in the region shown as the black square in Figure 16 is positively correlated with the monthly IOFs of BMU S and negatively correlated with BMU M . The correlation coefficient between the monthly mean ζ and the IOFs of BMU S is 0.82 during the study period which is from 1993 to 2016, implying that the mean relative vorticity increases in summer months, and the corresponding IOF of BMU S is higher, i.e., the straight-path pattern more likely occurs when the relative vorticity increases especially in summer. In contrast, from Figure 17b one can see that the meandering-path pattern of the Kuroshio axis more likely occurs when the relative vorticity decreases. In conclusion, for the multi-year monthly mean scale, the anticyclonic eddies can carry negative PV flux and promote the Kuroshio axis to present the meandering-path pattern, while the cyclonic eddies can carry positive PV flux and promote the Kuroshio axis to present the straight-path pattern. The results are similar to the analysis of Figure 16 and the dynamical interpretations in [18-20]. region where the interaction between the eddies and Kuroshio often occurs (black square in Figure  16a), it indicates that the net relative vorticity decreases, then the water depth H must decrease in order to conserve the PV. Therefore, the water moves onshore because the water depth on the ECS shelf is shallower and the Kuroshio surface axis presents BMUM pattern (green line in Figure 16a). In contrast, in July 2009, the monthly mean ζ is positive in the same region because of the cyclonic eddies, then the water moves offshore and the Kuroshio surface axis presents BMUS pattern (Figure 16b). The impacts of the anticyclonic and cyclonic eddies are opposite. The correlation coefficient between the monthly mean ζ and the IOFs of BMUS is 0.82 during the study period which is from 1993 to 2016, implying that the mean relative vorticity increases in summer months, and the corresponding IOF of BMUS is higher, i.e., the straight-path pattern more likely occurs when the relative vorticity increases especially in summer. In contrast, from Figure 17b one can see that the meandering-path pattern of the Kuroshio axis more likely occurs when the relative vorticity decreases. In conclusion, for the multi-year monthly mean scale, the anticyclonic eddies can carry negative PV flux and promote the Kuroshio axis to present the meandering-path pattern, while the cyclonic eddies can carry positive PV flux and promote the Kuroshio axis to present the straight-path pattern. The results are similar to the analysis of Figures 16 and the dynamical interpretations in [18][19][20].

Conclusions
This study examines the variability of the Kuroshio surface axis northeast of Taiwan Island using satellite altimeter data. The 24-year time series of the Kuroshio axis are derived from the AGV data from 1993 to 2016 using the detection method established by [23] with particularly selected parameters, including the length of the subsidiary line, the intervals between two adjacent points, and the distance that the subsidiary line moves downstream once. The results indicate that the Kuroshio mean axis generally moves between 1.0 m and 1.2 m contours of the mean ADT, and large fluctuations on the region northeast of Taiwan Island appear because of the variation of the Kuroshio intrusion into the ECS shelf.
Based on the day-to-day Kuroshio axis detection, the daily distances away from the location of the mean axis along the Transect A are obtained. The EMD method is used to analyze the daily distance series and quantify the variability of the Kuroshio axis. The mean period and the energy of each of the IMFs are calculated using the discrete Hilbert transform. The results indicate that the intra-seasonal and inter-annual variability of the Kuroshio surface axis with a mean period of about 3.2 months and 1.3 year are the dominant components. To examine the intra-seasonal variation of the surface Kuroshio axis, 24-year mean daily series of IMF4-IMF9 is calculated. The fluctuation of IMF4 is the strongest in winter, followed by spring and autumn. IMF5 and IMF6 correspond to the monthly variations, the variations of the Kuroshio axis between shoreward and seaward are distinct, but the regularity remains non-obvious. IMF7-9 correspond to the intra-seasonal variations, the intra-seasonal variability of the Kuroshio axis occurs quasi-periodically, and is consistent with the variations of the Kuroshio intrusion. Transects B and C are used to analyze the relationship between the relatively upstream Kuroshio axis east of Taiwan Island and the Kuroshio axis northeast of Taiwan Island. The comparison results indicate that the KC flows closely to Taiwan Island when the Kuroshio intrusion into the ECS shelf occurs.
The SOM analysis identifies four patterns of the Kuroshio axis variation. BMU M and BMU S correspond to meandering-path and straight-path states, while BMU T1 and BMU T2 transition stages. The occurrence percentage of BMU M is only about 1/6 of that of BMU S , implying that the strong deflection of the Kuroshio axis occurs not frequently. The monthly IOFs of BMU S show strong seasonality: higher value in summer, and lower value in spring and autumn. The correlation analysis with the upstream surface current transport along Transect B indicates that the straight-path pattern more likely occurs when the upstream surface transport increases.
To examine the transient evolution of the four Kuroshio axis BMU patterns, the ADT and geostrophic currents, SLA and geostrophic velocity anomaly, and AVHRR SST during the typical period from November 26, 2012 to January 27, 2013 are used. The results indicate that during this period the patterns of the Kuroshio axis change successively in order of BMU T1 →BMU M →BMU T2 →BMU S , implying that the Kuroshio axis migrates from the meandering-path to the straight-path pattern. Meanwhile, the warm water intrusion and a mesoscale eddy occur at the second stage and migrate northwestward gradually at the last two stages. It should be noted that the evolution order appears only for the typical period and is not common during the whole study period. More detailed analysis to find out generalized characteristics of the transient evolution need further works in future.
The monthly mean relative vorticity is calculated and analyzed to evaluate the impact of the eddies on the Kuroshio surface axis variability, the results show that the anticyclonic eddies can make the Kuroshio water move onto the ECS shelf and promote the Kuroshio surface axis to present the meandering-path pattern because of the PV conservation, while the cyclonic eddies can make the Kuroshio water move offshore and promote the axis to present the straight-path pattern. It is worth noting that the PV conservation equation is the first-order approximation to the PV in the real and stratified ocean, where the baroclinic effect also contributes to the variation of vorticity [33,53]. The mechanistic explanation is simplified in this first-order form. Therefore, more quantitative vorticity analysis and more detailed study for longer period need to be performed in future.