Spatiotemporal Variability of Mesoscale Eddies in the Indonesian Seas

: Mesoscale eddies are ubiquitous in the world ocean and well researched both globally and regionally, while their properties and distributions across the whole Indonesian Seas are not yet fully understood. This study investigates for the first time the spatiotemporal variations and generation mechanisms of mesoscale eddies across the whole Indonesian Seas. Eddies are detected from altimetry sea level anomalies by an automatic identification algorithm. The Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea are the main eddy generation regions. More than 80% of eddies are short-lived with a lifetime below 30 days. The properties of eddies exhibit high spatial inhomogeneity, with the typical amplitudes and radiuses of 2–6 cm and 50–160 km, respectively. The most energetic eddies are observed in the Sulawesi Sea and Seram Sea. Eddies feature different seasonal cycles between anticyclonic and cyclonic eddies in each basin, especially given that the average latitude of the eddy centroid has inverse seasonal variations. About 48% of eddies in the Sulawesi Sea are highly nonlinear, which is the case for less than 30% in the Sulu Sea and Banda Sea. Instability analysis is performed using high-resolution model outputs from Bluelink Reanalysis to assess mechanisms of eddy generation. Barotropic instability of the mean flow dominates eddy generation in the Sulu Sea and Sulawesi Sea, while baroclinic instability is slightly more in the Maluku Sea and Banda Sea.


Introduction
The progress of numerical models and remote sensing techniques, especially satellite altimeter, in the past few decades has greatly advanced our understanding of mesoscale eddies in the ocean, which contain most of the oceanic kinematic energy [1]. Mesoscale eddies have characteristic spatial and temporal scales ranging from tens to hundreds of kilometers and from several days to years, respectively [2]. Both cyclonic (CE) and anticyclonic (AE) eddies are ubiquitous in the ocean [1], serving as a key bridge of energy cascade between large-scale and sub-mesoscale oceanic dynamics [3][4][5]. Because of their nonlinearity, mesoscale eddies play a vital role in the transport of momentum, mass, heat and biogeochemical properties and further impact tracer budgets and primary production [2,6]. Mesoscale eddies in the Indonesian Seas exhibit multiscale temporal variations associated with the ITF. For the intraseasonal scale, a 50-day oscillation of horizontal velocity was observed from moorings at the entrance of Sulawesi Sea [30]. Based on a 1.5-layer reduced-gravity model, Qiu et al. [28] pointed out that this intense 50-day oscillation in the Sulawesi basin is a result of baroclinic Rossby wave resonance. The 50-day oscillation signals also exhibit intense interannual variability modulated by active eddy shedding, and enhanced 50-day oscillation can freshen the upper-ocean water mass in the Sulawesi Sea and Makassar Strait [31]. Through an ocean general circulation model (OGCM), three eddies in the Flores Sea were simulated in austral summer when ITF transport is low, and these eddies vary synchronously at an interannual scale, thus named as "Lombok Eddy Train" [32]. In terms of the seasonal cycle, the variation of eddy kinematic energy (EKE) in the Sulawesi Sea has different periods at different depths, and the annual and semiannual peaks are in the 0-100 m layer and the 100-300 m layer, respectively [33]. In the Sulu Sea, AEs and CEs display an opposite seasonal variability, with more and larger AEs in boreal winter and CEs in boreal summer [34]. However, a comprehensive description of properties, seasonal variability and generation mechanism for eddies in the Indonesian Seas is still lacking.
This study provides a comprehensive statistical description of the spatiotemporal variability of eddies in the Indonesian Seas using long term altimeter data. The mechanisms of eddy generation are investigated through instability analysis based on high-resolution model outputs. We believe that an overall presentation of mesoscale eddies in the Indonesian Seas will not only improve local ocean forecasts but also facilitate a better understanding of their roles in climate and ecosystems. The remainder of this paper is organized as follows: Section 2 describes the details of data and methods used in this study. Section 3 provides a statistical description of eddy properties. The mechanisms of eddy generation are discussed in Section 4. Section 5 provides a summary of the results.

Altimeter Data
The delayed-time altimeter data used in this study is a Ssalto/Ducas gridded product of SSH, sea level anomaly (SLA) and geostrophic current from January 1992 to December 2019 provided by Archiving, Validation, and Interpolation of Satellite Oceanographic (AVISO) data and distributed by the Copernicus Marine Environment Monitoring Service (CMEMS). This product has been widely used to detect mesoscale eddies globally [17,35] and regionally [12,19,34,36]. The SLA product was constructed by merging multi-mission satellites since 1992 optimally interpolated onto 0.25° × 0.25° grid with a daily resolution. The geostrophic velocity was computed by the Lagerloef methodology [37] introducing the β-plane approximation in the equatorial band (5°S~5°N), and by the 9-point stencil width ("stencil width" means the number of grid points utilized to estimate the finite difference approximation to the derivative on a grid) methodology outside the equatorial band [38]. Following Chelton et al. [2], the SLA fields are therefore spatially high-pass filtered with half-power filter cutoffs of 20° longitude × 10° latitude to remove steric heating and cooling effects, as well as other large-scale variability of SLA. The readers should be reminded that altimeter data become less accurate over the shelf area due to some aliases from the tides, coastal wave signals and sea-land transition [39]. Thus, the altimeter data in areas shallower than 200 m are masked in this study.

Bluelink ReANalysis
To study the mechanisms of eddy generation, we use the latest model outputs from Bluelink ReANalysis (BRAN, version of 3p5) to perform an instability analysis. The BRAN model is a multi-year integration of Ocean Forecasting Australian Model (OFAM) assimilating observations of SLA and SST from satellite, sea level data from tide gauges and in situ temperature and salinity profiles by using Bluelink Ocean Data Assimilation System (BODAS) [40]. This model domain covers the Indonesian Seas and surrounding oceans with a horizontal resolution of 0.1° and a vertical resolution of 10 m in the upper 300 m [40]. Compared with existing observations (e.g., moorings, drifters and fields survey) which are not assimilated into the model, BRAN performs well in the Indonesian Seas and surrounding ocean and can especially capture the realistic details of seasonal circulation and its interactions with various topographic features [41][42][43]. The entire dataset covers the period from January 1994 to August 2016. The daily outputs of sea level, zonal velocity, meridional velocity, potential temperature and salinity in the upper 300 m from January 1994 to December 2015 are used for the instability analysis.
To compare the time-varying signals of SSH between BRAN and AVISO, an empirical orthogonal function (EOF) analysis is performed to both daily BRAN and AVISO SSH data with a daily resolution from 1994 to 2015. The first two modes of BRAN and AVISO reveal good resemblance in both spatial patterns and principal components (PC), and their cumulative explained variances are 64.42% and 63.89%, respectively, which can describe the main characteristics of SSH variations in the Indonesian Seas and surrounding oceans. The spatial patterns and PCs of EOF mode 1 from BRAN and AVISO are displayed in Figure 2. Mode 1, respectively, explains 45.55% and 42.89% of the total variance in the BRAN and AVISO. BRAN and AVISO share a similar spatial pattern for EOF mode 1 (Figure 2a,b). The time series of two PCs also match favorably with a high correlation coefficient of 0.98. The EOF mode 1 physically characterizes the large-scale variations of SSH in response to wind forcing dominantly modulated by the El Niño-Southern Oscillation (ENSO) [31]. The spatial patterns and PCs for EOF mode 2 of SSH from BRAN and AVISO are shown in Figure 3. BRAN captures the EOF mode 2 of observed SSH field faithfully, not only in spatial patterns (Figure 3a,b) but also in PCs (Figure 3c). The variances explained by the EOF mode 2 are 18.87% and 21% for BRAN and AVISO, respectively, and the correlation coefficient between two PCs is 0.98. From the PCs (Figure 3c), it is clear that the EOF mode 2 represents the seasonal SSH variations in the Indonesian Seas and surrounding oceans. For example, Figure 3a,b indicate the surface circulations in the Sulu Sea and southern SCS exhibit a clear seasonal cycle, with an anticyclonic circulation during summer and a cyclonic circulation during winter, which has been proved based on both numerical models and observations [44,45]. The fact that BRAN simulates successfully the first two EOF modes of observed SSH signals is important since the large-scale circulation patterns tend to modulate mesoscale eddy activities.

Eddy Detection and Tracking Algorithm
Numerous automatic eddy detection algorithms have been developed based on the physical or geometric criteria, and they can be divided into three categories: (1) the physical parameter method, such as the Okubo-Weiss parameter method [46]; (2) the flow geometry method, including the winding-angle method [18,47] and the vector geometry method [48]; (3) the SSH-based method [2,17,35]. However, not all algorithms are suitable for identifying mesoscale eddies in the Indonesian Seas according to three reasons. Firstly, the SSH-based method performs better than the Okubo-Weiss parameter method because of its ability to avoid noise and excess eddy detections [2]. Secondly, the flow geometry algorithms require higher resolution data to get an accurate flow field to identify eddies, and the existing observational data in the Indonesian Seas cannot satisfy this demand [35]. Therefore, we adopted the SSH-based method developed by Faghmous et al. [17] (hereafter JHF15) which has been applied in the Kuroshio Extension Region [10], the Bay of Bengal [36] and the Southeastern Indian Ocean [12]. JHF15 identifies an eddy as a closed SLA contour with a single extreme. It is considered as a parameter-free method in which no empirical parameters are applied and the identified eddies' edges depend only on the single extreme approximation. In this study, both altimeter data and model outputs are used for eddy detection. For each dataset, all identified eddies larger than 9 corresponding grid cells are kept, to avoid some spurious features, and more details of JHF15 refer to Faghmous et al. [17]. The reader should note that the results in Section 3 are only from the altimeter data, while eddy detections from the model outputs are used to assess the capability of BRAN in terms of capturing mesoscale eddies.
After eddies are detected in each daily SLA snapshot, they are tracked by an algorithm developed by Penven et al. [49]. The nondimensional distance of an eddy pair with the same polarity from two consecutive maps is defined as: where ∆ , ∆ and ∆ are, respectively, the differences in eddy centroid location, surface area and amplitude between and in two consecutive maps; and the characteristic length scale , the characteristic surface area , and the characteristic amplitude are 25 km, 60 km 2 and 2 cm, respectively. The smaller , , the higher similarity of the eddy pair. Due to the sampling errors and measurement noise of satellite altimeter, some eddies may disappear at some time and reappear several time steps later. If an eddy moves into the gap among satellites' orbits, it will also vanish in the SLA map. To solve this problem, we repeat this tracking algorithm for a longer time step, from 2 to 7 days, for those unpaired eddies. We also set the maximum search distance ∆ to 150 km per week to avoid eddies jumping to another track. Following He et al. [34] and Zhang et al. [12], the number of eddy tracks and the number of eddies are different. The former is counted when an eddy track is once identified, while the latter is the total number of eddies along the eddy tracks.

Definition of Eddy Properties
The amplitude of eddy a is defined as: where 〈 〉 is the mean value of SLA at the eddy edge.
The surface area of eddy A is the area delimited by the outermost contour of SLA, the eddy edge, whereas its apparent radius R corresponds to the radius of an equivalent circular eddy with the same area. Thus, R is calculated as: The mean EKE of eddy is calculated as: where and are the zonal and meridional components of geostrophic velocity anomaly, respectively, calculated from SLA.
The mean vorticity of eddy is calculated as:

Instability Analysis
The EKE is generally converted from mean-flow kinematic energy (MKE) through barotropic instability and Kelvin-Helmholtz instability or from eddy potential energy (EPE) via baroclinic instability [3][4][5]50]. Therefore, the energy conversion rates for the above three instabilities are calculated as follows [3,12,50]: Barotropic conversion rate (BTR) from MKE to EKE via barotropic instability, depending on horizontal shears of mean flow, is defined as: Kelvin-Helmholtz conversion rate (KHR) from MKE to EKE via Kelvin-Helmholtz instability depending on vertical shears of mean flow and Reynolds stresses, is defined as: Baroclinic conversion rate (BCR) from EPE to EKE via baroclinic instability is defined as: where is the background density set to 1030 kg/m 3 . The overbars and primes represent time mean (1 month) and anomalies from time mean, respectively.

Results
In the Indonesian Seas, a total of 46,676 AEs ( Figure 4a) and 47,004 CEs ( Figure 4b) are identified from January 1993 to December 2018. Most identified eddies are concentrated on four larger and deeper basins (i.e., the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea). Due to the shallower depth (<200 m) and smaller basin scale, fewer eddies are detected in the rest of the Indonesian Seas. To obtain robust statistical description for eddy properties, eddies with their amplitudes larger than 2 cm are analyzed in our study, with the consideration that the accuracy of SLA gridded data from AVISO is about 2 cm [15,51].

Eddy Geneis and Decay
The spatial distributions of the numbers of eddy genesis and decay events for AE and CE are shown in Figure 5. The locations of eddy for generation and decay are the first and the last point in each eddy track, respectively. Because the minimum resolvable temporal scale of gridded SLA product from AVISO is about 10 days [52], we excluded those eddies with a lifespan shorter than 10 days. Consequently, a total of 469 AE tracks and 500 CE tracks are identified in the Indonesian Seas. AE tracks in the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea, respectively, account for 15%, 53%, 7% and 25% of the total AE tracks, while the proportion of CE tracks are 15%, 52%, 7% and 27%, respectively. Hence, no apparent regional preference for AE and CE is found in the Indonesian Seas as is also observed in the Tropical Atlantic Ocean [19]. Due to the removal of SLA data in the area shallower than 200 m, the eddy generation and decay are less frequent near the coast. It is interesting that both AE and CE preferentially generate in the areas north of Ombai Strait ( Figure 5a,b), which may result from the interactions of intense mean flow with topographic features [53,54]. However, there are some differences of the geographical pattern for the formation and decay between AE and CE. We find that the generation of AE is concentrated on the southeastern portion of the Sulu Sea (Figure 5a) while CE tend to form in northeastern portion of the Sulu Sea (Figure 5b). Additionally, the larger value of CE formation southwest of the Mindanao Island may partially result from the interaction between eastward currents and the coastline [53,54]. In the Sulawesi Sea, AE usually decays in the central and western part (Figure 5c) while the death of CE concentrates on the western and northeastern part (Figure 5d).

Eddy Propagation
The propagation velocities are calculated by a forward difference scheme for eddy centroid displacement every day. The velocities of each eddy track were further filtered using a central moving average of 7 days to reduce the random noise and then averaged in each 0.25° × 0.25° grid, following He et al. [34]. Figure 6 displays the mean propagation velocity fields and the standard deviation (STD) of azimuth of propagation direction relatively to the west for AE and CE. In the Sulu Sea, the AE ( Figure 6a) and CE (Figure 6b) translation speed are approximately 2.3 ± 1.0 cm/s and 0.9 ± 0.5 cm/s, respectively. AEs move southeastward at the northeastern half basin because of the advection of mean flow from the Mindoro Strait during boreal spring and summer when AEs are more frequent (Figure 6a). In the Sulawesi Sea, the AE translation speed is approximately 2.5 ± 2.3 cm/s, with a westward propagation in the central basin (Figure 6a). The CE migrations across the entire basin are significantly affected by the intense intrusion of MC (Figure 6b), with a speed of 6.3 ± 2.7 cm/s. In the Maluku Sea, all eddies propagate southwestward, which is parallel to the long axis of basin, at a velocity of 1.9 ± 0.7 cm/s for AE and 1.7 ± 0.6 cm/s for CE (Figure 6a,b). In the Banda Sea, the CEs are advected by the eastward mean flow during austral summer with more CEs (Figure 6b). More frequent AEs move southeastward with mean flow in the austral spring and autumn (Figure 6a). The corresponding propagation speeds of AE and CE in the Banda Sea are 5.7 ± 2.5 cm/s and 3.9 ± 2.5 cm/s (More details refer to Section 3.5 and Appendix A).

Eddy Lifespan
In this study, the lifespan of an eddy is defined as the number of days between the genesis and decay of the eddy. The upper-tail cumulative histograms of the eddy lifespan for the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea are presented in Figure 7. The short-lived (10~30 days, referenced in [15]) AEs (CEs) account for 73% (84%), 85% (86%), 100% (94%) and 86% (80%) in the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea, respectively. This suggests that most mesoscale eddies in the Indonesian Seas are shortlived, as also observed in the tropical oceans including the equatorial current systems and tropical oceanic warm pools [15]. The lifespans of AE and CE have a generally similar distribution in all these seas with a rapid decrease within 40 days and a slow decrease beyond 40 days except for the Maluku Sea where the eddy lifespan is mostly shorter than 30 days. However, there are some differences between AE and CE in each sea. In the Sulu Sea (Figure 7a), the portion of AE with lifespan larger than 30 days, 27%, is higher than CE with a percentage of 16%. The distributions of lifespan for eddies are more skewed to high values for AE than CE in the Sulawesi Sea (Figure 7b). For the relatively longevous eddies, CEs are more abundant than AEs in the Banda Sea, with a percentage of 20% (Figure 7d).

Distribution of Eddy Properties
The spatial distributions of amplitude and radius for AEs and CEs averaged in 1° × 1° grid are shown in Figure 8. Eddies have typical amplitudes of 2-6 cm and radiuses of 50-160 km in the Indonesian Seas. The patterns of amplitude (Figure 8a,b) are generally similar to those of radius (Figure 8c,d), indicating that the amplitude is positively correlated with the radius. Both amplitude (Figure 8a,b) and radius (Figure 8c,d) are larger in the interior basin and smaller near the coast for the Sulu Sea, Sulawesi Sea and Banda Sea. This may be due to the strong energy dissipation near the coast [34], and eddy growth in the central basin resulted from the eddy interactions [55]. Meanwhile the amplitude and radius for AE and CE in the Maluku Sea and Seram Sea are smaller due to small basin scale, as also observed in the Red Sea by Zhan et al. [47]. The large values of eddy radius are observed in the Banda Sea, where the maximum value of 157 km is observed in the central basin and the majority is more than 100 km. Additionally, the amplitude and radius of CE (Figure 8b,d) are stronger than those of AE (Figure 8a,c) in the Sulawesi Sea, where MC intrudes into Sulawesi Sea as a cyclonic loop structure [22].  Figure 9 displays the spatial distributions of mean EKE and mean relative vorticity for AE and CE averaged in 1° × 1° grid. There is no significant difference between AE and CE in the spatial pattern of mean EKE in the Indonesian Seas. The spatial pattern of mean EKE is similar to that of amplitude, which indicates that the tangential speed of eddy is largely determined by amplitude [12,34]. Both the magnitude and variation of mean EKE in the Sulu Sea agree well with the estimate of He et al. [34]. The most energetic eddies are observed in the Sulawesi Sea and Seram Sea with high MKE [42,56]. MKE is converted to EKE through the barotropic pathway [3][4][5], which will be discussed in Section 4. Although high eddy amplitude is also observed in the central Banda Sea (Figure 8a,b), the value of mean EKE (Figure 9a,b) is smaller than those in the Sulawesi Sea due to larger eddy area (Figure 8c,d). The spatial distribution of mean relative vorticity for AE and CE are generally similar (Figure 9c,d). Eddies with high mean relative vorticity were found in the Sulawesi Sea and Seram Sea, which is similar to the distribution of mean EKE. However, compared with other eddy properties, high value of mean vorticity is located near the coast instead of in the interior basin. On the one hand, the sampling number of eddies is very small in the coastal area. On the other hand, strong lateral friction leads to strong horizontal shear of velocities in the coastal region [55].

Seasonal Variability
The mean seasonal cycle of eddy properties in the Sulu Sea is shown in Figure 10. For AE, the mean number (Figure 10a), amplitude ( Figure 10b) and radius (Figure 10c) display similar patterns, which are large in boreal summer and small in boreal winter, with peak values of 0.4, 3.7 cm and 123.5 km in August, respectively. The mean number, amplitude and radius of CE display similar patterns but in an opposite phase with AE, with large values in boreal winter and small values in boreal summer (Figure 10a-c). The latitude of eddy centroid presents opposite seasonal variability between AE and CE (Figure 10d). This indicates more AEs than CEs in the northern Sulu Sea from November to February, while more CEs populate in the southern Sulu Sea from August to October. The mean EKE of CE have a seasonal cycle with a bimodal structure (Figure 10e), which reaches the annual maximum of 352.5 cm 2 /s 2 in June, and then drops rapidly to the annual minimum of 239.6 cm 2 /s 2 in October. From then it increases dramatically to the local maximum of 317.4 cm 2 /s 2 in December and decreases to the local minimum of 248.6 cm 2 /s 2 in January. While the mean EKE of AE shows a single peak in May with a value of 334.5 cm 2 /s 2 , and a single trough with a value of 204.2 cm 2 /s 2 in March. The mean relative vorticity of CE (Figure 10f) also has a bimodal structure, in which two peaks exist in September and December, with corresponding values of 5.1 × 10 −6 s −1 and 4.4 × 10 −6 s −1 , and two troughs appear in February and October with corresponding values of 3.5 × 10 −6 s −1 and 3.8 × 10 −6 s −1 . No significant seasonal variation is observed for the mean relative vorticity of AE ( Figure  10f). Although the magnitudes of monthly average eddy properties are different from He et al. [34] due to different limitations of amplitude in eddy identification algorithms, the seasonal cycle of eddy characteristics displays almost identical tendency not only for AE but also for CE.  corresponding values of 0.7, 3.5 cm and 99.9 km, and reach local maximum in July, with values of 0.8, 3.6 cm and 105.2 km, respectively. While the mean number, amplitude and radius of CE vary in an opposite phase compared with those of AE (Figure 11a-c). It is noted that CEs have larger amplitude than AEs through the whole year in the Sulawesi Sea. This further confirms the fact found above from the spatial perspective, the CE amplitude is larger than AE amplitude. The latitude of eddy centroid between AE and CE presents the opposite tendency in the Sulawesi Sea, with more AEs in the southern Sulawesi Sea from March to September while more CEs in the northern Sulawesi Sea from October to February (Figure 11d). Both the seasonal cycle of mean EKE and mean relative vorticity have bimodal structure for CE (Figure 11e,f). Two peaks appear in January with values of 1816 cm 2 /s 2 and 12.9 × 10 −6 s −1 and May with values of 1688 cm 2 /s 2 and 12.3 × 10 −6 s −1 , respectively, while two troughs exist in April with values of 1480 cm 2 /s 2 and 11.3 × 10 −6 s −1 and September with values of 1493 cm 2 /s 2 and 11.2 × 10 −6 s −1 , respectively. The mean EKE of AE have a similar pattern to that of CE, with two peaks in December and May and two troughs in April and September (Figure 11e). No apparent seasonal variations are observed for the mean relative vorticity of AE in the Sulawesi Sea, with a maximum value of 13.2 × 10 −6 s −1 and a minimum value of 11.5 × 10 −6 s −1 (Figure 11f). The mean seasonal cycle of eddy properties in the Maluku Sea is shown in Figure 12. The monthly mean number of AE and CE are less than 0.2 in the Maluku Sea (Figure 12a).
The amplitude of CE has a bimodal curve with annual maximum of 3.1 cm in February, annual minimum of 2.6 cm in July, local maximum of 2.8 cm in October and local minimum of 2.7 cm in December (Figure 12b). While there is no clear seasonal cycle of AE amplitude with the maximum of 3 cm in November and the minimum of 2.7 cm in July (Figure 12b). The seasonal variation of CE radius is similar to that of CE amplitude. Two peaks appear in March and October, with corresponding values of 98.7 km and 103.7 km, and two troughs are in February and May, with corresponding values of 88.51 km and 91.4 km (Figure 12c). Like AE amplitude, the AE radius has no apparent seasonal cycle, with the maximum of 99.8 km in October and the minimum of 87.9 km in January ( Figure  12c). AEs in the Maluku Sea prefer to locate south of the equator from July to November and appear north of equator from December to June (Figure 12d). On the contrary, CEs in the Maluku Sea tend to populate south of the equator in the first half of the year and north of the equator in the second half of the year (Figure 12d).   Figure 13a,b, the mean number and amplitude of CE hit peak in austral summer with corresponding values of 0.7 and 4.8cm. From then, they decrease gradually to the trough in austral winter with values of 0.3 and 2.9 cm, respectively. Meanwhile, the CE radius reaches the maximum of 143.1 km in January and the minimum of 124.0 km in March (Figure 13c). The latitude of AE and CE shows opposite patterns in the Banda Sea ( Figure  13d). AEs tend to seat in the southern Banda Sea during austral winter and in the northern Banda Sea during austral summer. Meanwhile, CEs prefer to populate in the southern Banda Sea for austral summer and in the northern Banda Sea for austral winter. More energetic AEs exist during austral summer with a peak value of 736 cm 2 /s 2 , and the mean EKE of AE reaches the minimum of 335.2 cm 2 /s 2 in austral winter (Figure 13e). The seasonal cycle of mean EKE for CE is a bimodal curve with two peaks in February and June with corresponding values of 672.5 cm 2 /s 2 and 585.6 cm 2 /s 2 , and two troughs in April and October with the values of 509.2 cm 2 /s 2 and 488.5 cm 2 /s 2 , respectively (Figure 13e). No significant seasonal variability of mean relative vorticity for CE is observed, with the maximum of 9.5 × 10 −6 s −1 in March and the minimum of 8.0 × 10 −6 s −1 in October (Figure 13f). While the seasonal cycle of mean relative vorticity for AE displays a bimodal structure in the Banda Sea: two peaks exist in March and July with corresponding values of 9.4 × 10 −6 s −1 and 7.9 × 10 −6 s −1 and two troughs appear in June and September with values of 7.1 × 10 −6 s −1 and 6.9 × 10 −6 s −1 , respectively (Figure 13f).

Eddy Nonlinearity
The advective nonlinearity parameter is calculated to assess the degree of eddy nonlinearity in the Indonesian Seas. The advective nonlinearity parameter is a nondimensional ratio of U/c, where U is the maximum rotation speed and c is the propagation speed of each eddy [2]. The definition of c is shown in Section 3.2. Figure 14 displays the uppertail histogram of U/c for the Sulu Sea, Sulawesi Sea and Banda Sea. Most eddies in these three seas are nonlinear by this measure, among which U/c exceeding 1 account for 89% (88%), 100% (100%) and 94% (96%) for AE (CE) in the Sulu Sea, Sulawesi Sea and Banda Sea, respectively. Most highly nonlinear eddies are observed in the Sulawesi Sea, with 48% of the U/c values exceeding 5 and 17% exceeding 10, as observed in major unstable and meandering currents of global ocean [2]. Fewer eddies with high nonlinearity are detected in the Sulu Sea and Banda Sea (Figure 14a,c), where the U/c values above 10 is less than 10%. The distributions of U/c for eddies are more skewed to higher values for AE than for CE in the Sulu Sea, but the opposite is found in the Banda Sea.

Discussion
To investigate the generation mechanism of eddies in the Indonesian Seas, we calculated the BTR, KHR and BCR, respectively, using BRAN outputs from 1994 to 2015. Due to the complex current patterns and the different spatial distribution of eddy genesis in each basin, we select four subregions in the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea (Figure 1, dashed boxes). These four subregions are defined in the Appendix A. We also find that BRAN can reasonably reproduce mesoscale eddies in the Indonesian Seas, and more details of evaluation refer to Appendix B. The seasonal cycles of volume integral BTR, KHR and BCR in each box are displayed in Figure 15. In these four regions, the volumetric integral BCR is positive throughout the year, which means EPE converts to EKE through baroclinic instability. Meanwhile, the volume integration of BTR for each box is positive in most seasons, except for some seasons in the Sulu Box and Maluku Box (Figure  15a,c). The weak negative BTR may result from inverse energy cascade from EKE to MKE [4,5], and the detailed discussions for this beyond the scope of present paper. However, the BTR almost dominates the energy conversion in the Sulawesi box throughout the year (Figure 15b), as in the Gulf Stream [4] and Kuroshio [5]. This result is similar to the conclusions of Yang et al. [33] that the EKE is governed by barotropic instability of mean flow. Compared with BTR and BCR, KHR are much smaller and make little net contribution to EKE in all four boxes. Therefore, barotropic and baroclinic instability dominate eddy generation in the Indonesian Seas. There are some seasonal variations of BTR and BCR, likely associated with the seasonal variations if the mean flow. To further determine which instability dominates for individual eddy generation over the 1994-2015 period, the monthly BTR and BCR integrated vertically in the upper 300 m are firstly calculated. Secondly, the corresponding values of BTR and BCR for each eddy genesis event were extracted according to its generation position and month. Thirdly, a comparison between BTR and BCR is performed to assess an eddy genesis event dominated by either barotropic instability (BTR > BCR > 0) or baroclinic instability (BCR > BTR > 0). Figure 16 presents the spatial distribution of the number of eddy genesis events dominated by barotropic and baroclinic instability. In the Sulu Sea, eddy genesis events dominated by barotropic instability account for 48% while those dominated by baroclinic instability are 35% of total eddy generation. This is similar to He et al. [34], who concluded that the instability of mean flow is one of the major eddy generation mechanisms. Barotropic instability dominates eddy generation in the Sulawesi Sea with a percentage of 55%, as in Yang et al [33], while the corresponding value of baroclinic dominance is 26% in the Sulawesi Sea. It is clear that barotropic instability is more frequent in the central Sulawesi Sea where the most intrusion of MC is located (Figure 16a), and eddy generation is dominated by baroclinic instability in the northern Sulawesi Sea (Figure 16b). In the Maluku Sea, 47% of eddy genesis events are dominated by baroclinic instability, while 43% are dominated by barotropic instability. Baroclinic instability in the Banda Sea dominates by 49% of eddy generation (Figure 16b). Additionally, the percentage of eddy generation dominated by barotropic instability is 40% in the Banda Sea (Figure 16a). For the Indonesian Seas as a whole, barotropic instability dominates the eddy generation with a percentage of 49%, which is higher than that dominated by baroclinic instability with a percentage of 35%. This is similar to the previous conclusions that barotropic instability is the dominant mechanisms in EKE generation along intense boundary currents, such as the Gulf Stream and Kuroshio, while baroclinic instability dominates in the broader open ocean [3][4][5]. A total of 16% of eddies in the Indonesian Seas are not dominated by either barotropic or baroclinic instability, so that costal Kelvin waves and complex topography may also affect eddy generation [57]. In our analysis, eddies with amplitude less than 2 cm are excluded according to the fact that the accuracy of gridded SLA product provided by AVISO is about 2 cm [15,51]. Generally, the variability of eddy amplitude during its lifetime undergoes growing stage, mature stage and decaying stage [15,19]. Inevitably, the eddy amplitudes less than 2 cm exist in the beginning and end of eddy lifespan, which are below the resolving ability of altimeter data [2,51]. There are two possible effects: one is underestimate of eddy lifespan; the other is that locations of eddy generation and decay tracked by eddy tracking algorithm may have some offset. However, our results in the Sulu Sea are similar to that of He et al. [34], especially with same spatial and temporal variability for eddy properties. The eddy amplitudes dramatically increase and decrease in growing stage and decaying stage, respectively, which account for less than 20% of the eddy lifespan [15,19]. Besides, the propagation speed of eddies in the Indonesian Seas is really small, with a mean value below 6 cm/s and identified eddies whose lifespan are less 30 days account for more than 80%. Thus, there is likely no significant difference between the positions detected by eddy tracking algorithm and the actual positions for eddy generation and decay in the region.
Additionally, the reader should be reminded that a free running model cannot reproduce real particular eddy events, while data assimilation inevitably interferes the dynamical or energy balance [40]. The new version of BRAN have used globally balanced forcing and low update cycle of assimilation (4 days) which improved the dynamic imbalance in the previous version [40]. Zhang et al. [12] have proved that BRAN performs well for the instability analysis and the robustness of their conclusions is further proved by a free running model, the eddy-resolving OGCM for the Earth Simulator (OFES). Besides, the magnitude of BTR averaged vertically in the Sulawesi Sea using BRAN outputs in this study is equal to the magnitude, 10 −3 cm 2 /s 3 , estimated based on OFES outputs [33]. Although the detailed assessment of the difference of EKE budget analysis between BRAN and a free running model beyond the scope of this paper, future studies about this assessment are on the way. This paper reports the statistical characteristics and spatiotemporal variations for the eddies across the whole Indonesian Seas for the first time. We believe that our results can be used to assess the performance of numerical models in terms of reproducing mesoscale eddies and to improve local ocean forecasts. Moreover, the significant meridional migrations and the high nonlinearity of eddies revealed in this study indicate considerable eddy transports of heat and salt [29,31], which will be the focus of our ongoing research. In addition, because eddies play a key role in the circulation of semi-closed and enclosed basins [34,47], our conclusions can also be used for the further investigations of eddy-ITF interaction.

Conclusions
In this study, we provided a detailed description for both spatial and temporal variability of eddies in the Indonesian Seas using altimeter data and discussed their generation mechanisms using high-resolution model outputs. A total of 469 AE tracks and 500 CE tracks were identified from January 1993 to December 2018. Most of eddies in the Indonesian Seas are short-lived (below 30 days) with a percentage of 73% (84%), 85% (86%), 100% (94%) and 86% (80%) for AE (CE) in the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea, respectively. There is no geographical preference for AE and CE in above four basins. For the mechanisms of eddy generation, barotropic instability dominates over baroclinic instability in the Sulu Sea and Sulawesi Sea, while baroclinic instability is slightly more important in the Maluku Sea and Banda Sea. Most energetic eddies are observed in the Sulawesi Sea and Seram Sea. The spatial distribution of eddy properties shows high inhomogeneity with typical amplitudes and radiuses of 2-6 cm and 50-160 km, respectively. Because of the intrusion of MC in the Sulawesi Sea, CE has larger amplitude and radius than AE. Limited by basin scale, the eddies in the Maluku Sea show smaller amplitude and radius.
The climatological properties of eddies in the Indonesian Seas exhibit different seasonal variations between AE and CE not only in each basin but also among four major basins. In the Sulu Sea, AEs prevail with relatively large amplitude and radius in the boreal summer, while CEs display an opposite seasonal variability. The mean number, amplitude and radius of eddies in the Sulawesi Sea present bimodal structure in seasonal cycles, with opposite phases between AE and CE. In the Banda Sea, the seasonal cycle of AE properties, including mean number, amplitude and radius, shows a bimodal structure, while CEs are more abundant with larger amplitude during austral summer than the other three seasons. There is no consistent tendency in the seasonal cycle of the mean number, amplitude and radius for both AE and CE in the Maluku Sea. It is interesting that every basin in the Indonesian Seas displays a reversal meridional distribution of eddy polarity with season. Nevertheless, our results only focused on a comprehensive description of eddy properties at the surface due to the limitation of altimeter data, future studies should explore the three-dimensional structure of eddies in the Indonesian Seas and their contributions to the ITF dynamics and thermodynamics.
Supplementary Materials: The following are available online at www.mdpi.com/2072-4292/13/5/1017/s1, Figure S1: Spatial distribution of amplitude and radius of AE and CE over the 1994-2015 period from BRAN. Figure S2: Spatial distribution of mean EKE and mean relative vorticity of AE and CE over the 1994-2015 period from BRAN. Figure

Appendix B
To evaluate whether BRAN reproduces realistic mesoscale eddies in the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea, we also detect eddies from the daily SLA snapshots of BRAN over the 1994-2015 period. The eddies with lifetime longer than 10 days and amplitude larger than 2 cm are considered for the analysis. The spatial distributions of eddy genesis and decay events from AVISO and BRAN over the 1994-2015 period are show in Figure A2. A total of 843 and 649 eddy tracks are detected from AVISO and BRAN over this period, respectively. The proportion of eddy tracks are 14% (22%), 53% (50%), 7% (1%) and 26% (27%) in the Sulu Sea, Sulawesi Sea, Maluku Sea and Banda Sea, respectively. The distributions of both eddy generation and eddy decay share a similar pattern between AVISO and BRAN in the Sulu Sea, Sulawesi Sea and Banda Sea. However, the numbers of eddy tracks in these four seas from BRAN are smaller than those from AVISO. This may result from the difference of data resolution, as is also found in the Subtropical Countercurrent regions from OFES outputs [59]. To further demonstrate the performance of BRAN in terms of mesoscale eddies, the comparison for the trajectories of eddies detected from AVISO and BRAN are performed. We find that BRAN can reproduce the realistic trajectories of mesoscale eddies in the regions we are interested, with the similar pathways and generation locations. Some examples in the Sulu Sea, Sulawesi Sea and Banda Sea are shown in Figure A3. Besides, BRAN can reasonably capture the spatial and temporal variations of eddy properties in the Sulu Sea, Sulawesi Sea and Banda Sea (see Supplementary Materials Figures S1-S5 for more details). Thus, we believe that BRAN is a credible reanalysis dataset to investigate the eddy generation mechanisms.