Decadal Variability in Salinity of the Indian Ocean Subtropical Underwater During the Argo Period

Salinity of the Indian Ocean Subtropical Underwater (IOSTUW) showed robust decadal variability during the Argo period, with a freshening trend between 2005 and 2013 and a salinification trend between 2014 and 2019. The salinity variability of the IOSTUW originated from both the sea surface salinity maximum region in the subtropical South Indian Ocean and the southeastern tropical Indian Ocean. The mixed layer salinity variabilities in these two regions are dominated by horizontal advection, which is mainly influenced by anomalous local surface currents and salinity anomalies from the western tropical Pacific and Maritime Continent. Both of these two factors have a close relationship with the El Niño‐Southern Oscillation (ENSO) signal. As a result, the mixed layer salinity in those two regions changed dramatically during the strong and long‐lasting ENSO events and thus led to the decadal variability in salinity of the IOSTUW.


Introduction
Sea surface salinity (SSS) is high in the subtropical South Indian Ocean due to an excess of evaporation over precipitation. The location of the SSS maximum, with a highest value exceeding 35.9 practical salinity unit (psu), is further poleward and eastward in comparison with those in other subtropical gyres (e.g., Nie et al., 2016;O'Connor et al., 2005;Qu et al., 2019). Some of the high-salinity surface water enters the upper thermocline during early austral spring when the mixed layer shoals rapidly and forms the Indian Ocean Subtropical Underwater (IOSTUW) which is usually identified as a subsurface salinity maximum (Grundlingh et al., 1991;O'Connor et al., 2005;Wijffels et al., 2002). Spreading in the South Indian Ocean subtropical gyre, the IOSTUW contributes directly to the meridional overturning circulation and interocean water exchange between the Indian Ocean and the South Atlantic (e.g., Miyama et al., 2003;Schott et al., 2009). Variation of the IOSTUW is closely related to the global hydrological cycle through changes in surface freshwater fluxes (e.g., Du et al., 2015;Hu et al., 2019). It also influences regional sea level, especially in the southeast Indian Ocean, where halosteric sea-level change is dominant (Llovel & Lee, 2015). Thus, analyses about the variability of the IOSTUW and underlying dynamics are necessities for a better understanding of the climate system on both regional and global scales.
The salinity variability of the IOSTUW has a close relationship with anomalies formed in the upper layer of the southeast Indian Ocean, which is dominated by interannual to decadal signals (Hu et al., 2019;Zhang et al., 2018). Recent studies indicated that both local precipitation and advection anomalies are responsible for the interannual variability (Hu et al., 2019;Zhang et al., 2016). While the decadal signal originates from the Maritime Continent (MC) and is conveyed into the southeast Indian Ocean via advection processes and eddy propagations (Du et al., 2015;Llovel & Lee, 2015). Changes in transport of ocean currents such as the Indonesian Throughflow (ITF) also contribute to the decadal variability (Du et al., 2015;Hu & Sprintall, 2017;Llovel & Lee, 2015).
By far, the variability of the IOSUTW has not been well understood due to a lack of subsurface salinity observations before the Argo era. Since the early 2000s, Argo profiles have been sharply accumulated (Roemmich et al., 2004) and allow for detailed investigations of interannual to decadal variabilities in subsurface salinity and underlying dynamics. Using a gridded Argo data set for the period between 2005 and 2019, this study provides the first observation-based estimate of the decadal variability of the IOSTUW. Relevant physical processes that contribute to this variability were also discussed, which might be valuable for future studies that focus on regional climate changes in the south Indian Ocean.

Data and Methods
A gridded Argo data set referred to as "MOAA GPV" (Grid Point Value of the Monthly Objective Analysis using the Argo data) is used in the present study. This data set is produced by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) and it shows monthly global temperature and salinity from January 2001 in quasi-real time (Hosoda et al., 2008). The Argo float, TRITON mooring, and available CTD data, in addition to a 2-D optimal interpolation method are applied in the creation of this gridded data set. This data set has a horizontal resolution of 1 × 1°with 25 standard pressure levels from 10 to 2,000 dbar. The gridded temperature and salinity fields for the period 2005-2019 are employed in this study. We abandoned the data between 2001 and 2004 for the rarity of the Argo floats during that period, which may lead to unrealistic results. Geostrophic velocity is derived from the gridded temperature and salinity fields by assuming a level of no-motion at 2,000 dbar.
Monthly salinity, potential temperature, and current velocity model outputs from the ECMWF Ocean Reanalysis System 5 (ORAS5) during January 2005 to December 2018 are used to perform the mixed layer salinity (MLS) budget analyses. Produced by the OCEAN5 system in its Behind-Real-Time (BRT) stream, ORAS5 is a global ocean ensemble reanalysis and covers the period from 1979 onwards, with a backward extension until 1958. The model is forced by surface fluxes based on the ERA-40 reanalysis until December 1978, the ERA-Interim reanalysis from January 1979 to December 2014, and the ECMWF operational NWP from January 2015 onward (Zuo et al., 2019). ORAS5 assimilates the in-site salinity and temperature profiles from the quality-controlled data set EN4 (Good et al., 2013), and sea-level anomaly observations produced by Archiving Validation and Interpretation of Satellite Oceanographic data (AVISO) Data Unification and Altimeter Combination System (DUACS; version DT2014; Pujol et al., 2016).
The time evolution of MLS can be written as follows: where S, E, and P represent the MLS, evaporation, and precipitation, respectively. h m is the MLD, defined as the depth where the density is 0.125 kg m −3 heavier than that on the surface. u and w e represent horizontal and entrainment velocity, respectively. ΔS is the salinity jump crossing the bottom of the mixed layer. The terms from left to right are named as the MLS tendency, surface freshwater forcing (SFF), horizontal advection, and vertical entrainment. For more details, please refer to Ren and Riser (2009) and Ren et al. (2011).
The SFF is derived from estimates of evaporation and precipitation from the ERA-Interim atmospheric reanalysis, produced by the European Centre for Medium-Range Weather Forecasts (Dee et al., 2011). These model outputs are produced by a forecast model, based on temperature and humidity information derived from the assimilated observations. For comparison, we also calculated the SFF using evaporation from the Objectively Analyzed Air-sea Fluxes (OAFlux; Yu, 2007;Yu et al., 2008) and precipitation from the Global Precipitation Climatology Project (GPCP; Adler et al., 2003). The GPCP product is based on satellite data and surface rain gauge observations, while the OAFlux product is constructed from an optimal blending of satellite observations and model reanalyses.

Decadal Variability in Salinity of the IOSTUW and Its Origins
Former studies have shown that the IOSTUW originates in the subtropical South Indian Ocean northeast of 35°S and 58°E, and resides between 24.7 and 26.3 σ θ , with a mean density of 25.5 σ θ and a salinity higher than 35.41 psu (Grundlingh et al., 1991;O'Connor et al., 2005). This study focuses on the characteristics of the IOSTUW between 25 and 26 σ θ , where the high-salinity core exists. As shown in Figure 1a, the IOSTUW is mainly confined in the region south of 15°S, which has also been noticed by Wijffels et al. (2002). After subduction, the IOSTUW first flows anticyclonic with the subtropical gyre and then joins the westward South Equatorial Current (SEC). Mean salinity of the IOSTUW (averaged within the region enclosed by solid black lines in Figure 1a) shows robust decadal variability, with a freshening trend during 2005-2013 and a salinification trend during 2014-2019 ( Figure 1b). The first empirical orthogonal function (EOF) modes of salinity over the IOSTUW region (Figures 1c-1d), which explains 63% of the total variation, show that there are two centers for this decadal variability. One resides in the winter outcropping region of the IOSTUW, while the other locates at around 75°E and 17°S. It seems like that the salinity signal originated from the outcropping region gradually weakened along its northwestward pathway in the upstream region (enclosed by green lines in Figure 1c) and get strengthened again after joining the SEC in the downstream region (enclosed by blue lines in Figure 1c). To investigate the relationship between the salinity signals in the upstream and downstream regions, regional mean salinity anomalies were calculated in these two regions, respectively ( Figure 1b). Results show that the salinity variations in these two regions are almost identical with each other; however, the signal in the downstream region leads the signal in the upstream region by 10 months. This suggests that the origin of the salinity variability in the downstream region should be different from that in the upstream region.
As shown in Figure 1c, the salinity signal in the upstream region of the IOSTUW originates from winter outcropping region in the SSS maximum region (black boxes in Figure 2a). The correlation coefficient between the regional mean MLS in the SSS maximum region and the regional mean salinity in the upstream region of the IOSTUW is 0.81. The salinity variability in the mixed layer enters the subsurface through the subduction processes. Following Qiu and Huang (1995), the long-term mean annual subduction rate and its two main components, vertical pumping at the bottom of the mixed layer and lateral induction due to the uneven MLD, were calculated (Figures 2d-2f ). Results show that both vertical pumping and lateral induction  Figure 1a), the upstream region (enclosed by green lines in Figure 1c) and the downstream region (enclosed by blue lines in Figure 1c). (c) The first empirical orthogonal function (EOF) mode for salinity over the IOSTUW region and (d) its time series. Black contours in Figures 1a and 1c indicate mean winter outcropping lines of 25 and 26 σ θ . A 13-month running mean filter has been adopted to remove the seasonal cycle before the analyses.

Geophysical Research Letters
contribute to the total subduction rate in the SSS maximum region. Closely related to the downward Ekman pumping, vertical pumping contributes 66% of the total subduction rate while the contribution from lateral induction is 34%.
The salinity signal in the downstream region of the IOSTUW is confined within the latitudes between 23°S and 15°S (Figure 1c), where is dominated by the strong westward SEC originated from the southeastern tropical Indian Ocean (SETIO) (red box in Figure 2a). The first EOF mode of salinity along the pathway of the SEC in Figure 2b, which explain 61% of the total variance, shows that the salinity variability formed in the upper layer of the SETIO is transported westward along the SEC pathway and gradually subducts into the depth range of the IOSTUW between 25 and 26 σ θ in west of 90°E. The subduction rate in the region between 15°S and 23°S is dominated by vertical pumping, while the contribution from the lateral induction is negligible (Figures 2d-2f ). The correlation coefficient between the regional mean MLS in the SETIO and the regional mean salinity in the downstream region of the IOSTUW is 0.67. Further analyses show that the variations of MLS in the SETIO and the SSS maximum region are closely related (Figures 3a and 3e). The lead-lag correlation coefficient between them is 0.91, and the signal in the SETIO leads the signal in the SSS maximum region by 12 months. Associating with the former result that the salinity variability in the downstream region of the IOSTUW leads the variability in the upstream region by 10 months, we can speculate that the SETIO should be the major origin of salinity signal in the downstream region, instead of the SSS maximum region.

MLS Budget Analyses
To better understand the underlying dynamics that contribute to the decadal variability in salinity of the IOSTUW, MLS budget analyses were performed in the SETIO and the SSS maximum region. The  (Kalney, 1996) and winter mixed layer depth (m), respectively.
variations of the MLS in these two regions are well simulated by the ORAS5 (Figures 3a and 3e); therefore, the products of this model are used in the MLS budge analyses.

Salinity Budget in the SETIO
In the SETIO, the sums of monthly integrated SFF, horizontal advection and vertical entrainment, both show a good consistency with the MLS tendency (∂S/∂t) (Figure 3b). The correlation coefficient is 0.61 between ∂S/∂t and the sum that contains the SFF derived from the ERA-Interim data and 0.80 between ∂S/∂t and the sum that contains the SFF based on OAFlux and GPCP products. However, the sum that contains the SFF derived from the ERA-Interim data shows unrealistic positive values before 2009, which is inconsistent with the freshening trend in the MLS during that period. The sum that contains the SFF based on OAFlux and GPCP products shows unrealistic negative values after 2017 which is inconsistent with the salinification trend of the MLS. One possible explanation for these biases could be that the degrees of correction for these products depend on the coverages of the observations on which they were based, which change in time. Nevertheless, ocean dynamic terms may also contribute to these discrepancies, which is not certain for lack of comparison.
For the long-term mean, the positive evaporation minus precipitation (E-P) tends to increase the MLS in the SETIO (Figure 3c). The vertical entrainment forcing is weak but also contributes to increase the MLS, due to the higher salinity below the mixed layer. The horizontal freshwater advection tends to lower the salinity for its close relationship with the ITF, which transports large amount of freshwater from the western tropical Pacific into the southeast Indian Ocean. The variation of the ∂S/∂t in the SETIO is dominated by the horizontal advection, for the correlation coefficient is 0.76 between ∂S/∂t and the horizontal advection, 0.33 between ∂S/∂t and the SFF based on the OAFlux and GPCP products, and 0.15 between ∂S/∂t and the SFF derived from the ERA-Interim data. The variability of the horizontal advection is mainly modulated by its zonal component, although the meridional component has a large contribution to the mean value (Figure 3d).

Salinity Budget in the SSS Maximum Region
In the SSS maximum region, there is a larger discrepancy between the ∂S/∂t and the sum that contains the SFF derived from the ERA-Interim data (Figure 3f). This may due to the difference between the data coverages that used to produce these reanalysis products. Besides, the neglection of contribution from eddy-induced meridional salinity flux could also contributes to this large discrepancy, which has been found to play an essential role in the MLS budget in the subtropical South Indian (Qu et al., 2019).
The salt gain from the positive E-P in this region is mostly taken away by horizontal advection, while vertical entrainment plays a minor role (Figure 3g). The contribution from vertical entrainment increased when the mixed layer was greatly deepened during 2005-2008 and 2016-2017 (supporting information Figure S1), and larger amount of fresher subsurface water merged into the mixed layer. The horizontal advection in this region is dominated by its meridional component (Figure 3h), due to the large meridional MLS gradient. The correlation coefficient is 0.19 between ∂S/∂t and the SFF based on the ERA-Interim data, 0.29 between ∂S/∂t and the SFF derived from the OAFlux and GPCP products, and 0.45 between ∂S/∂t and the horizontal advection. Therefore, the variability of the ∂S/∂t in the SSS maximum region is also dominated by horizontal advection.

Link to the ENSO Events
Then how did the horizontal advection generate the decadal variability in the MLS in the SSS maximum region and the SETIO during the Argo period? Previous studies suggested that the upper layer salinity variability in the southeast Indian Ocean is strongly affected by signals from the western tropical Pacific and the Maritime Continent (WTP-MC) (Du et al., 2015;Hu & Sprintall, 2017;Llovel & Lee, 2015;Zhang et al., 2018). The E-P variability in the WTP-MC is closely related to the ENSO signal ( Figure 4) and thus induces positive (negative) salinity anomalies during El Niño (La Niña) events. These salinity anomalies are conveyed into the southeast Indian Ocean via the ITF and then spread out through the westward SEC and the southward East Gyral Current (EGC) and the Leeuwin Current (LC) (Domingues et al., 2007;Gordon et al., 1997;Hu et al., 2019;Meyers et al., 1995;Wijffels et al., 2002), as illustrated in Figure 2a.
The variations of surface currents in the southeast Indian Ocean also show an ENSO-related signal. The mixed layer of the SETIO is dominated by a westward transport by the SEC and a southward Ekman transport generated by the background southeast trade winds. While in eastern part of the SSS maximum region (between 100°and 115°E) where the most robust decadal variability is found, the mixed layer is dominated by a southward transport associated with the EGC and the LC. These surface currents show similar variabilities that get strengthened (weakened) during La Niña (El Niño) ( Figure 4) and contribute to lower (increase) the MLS in corresponding regions. The variation of the LC was found to be connected with the anomalous precipitation in the MC and Indonesian-Australian Basin .
Under the influences of salinity anomalies from the WTP-MC and the anomalous local surface currents, the MLS in the SSS maximum region and the SETIO changed dramatically during the strong and long-lasting ENSO events between 2005 and 2019 (Figures 3a and 3e), including two La Niña events (2007-2009 and 2010-2013) and one El Niño event (2014)(2015)(2016) (Figure 4). During those two La Niña events, intensified precipitation in the WTP-MC and strengthened local surface currents transports contributed to negative anomalies in the horizontal advection (Figures 3c and 3g). This in turn led to negative anomalies in the MLS tendency (Figures 3b and 3f ), corresponding to declines in the MLS (Figures 3a and 3e). As a result, the MLS in these two regions showed a significant freshening trend during 2005-2013, although the MLS increased sharply during the relatively short El Niño in 2009. Notably, the decreased local E-P in the

Conclusions
Based on a gridded Argo data set named "MOAA GPV," we identified the robust decadal variability in salinity of the IOSTUW, with a freshening trend during 2005-2013 and a salinification trend during 2014-2019. The salinity variability in the upstream region of the IOSTUW originated from the SSS maximum region in the subtropical south Indian Ocean, while that in the downstream region originated from the SETIO. MLS budget analyses show that the variations of MLS in the SSS maximum region and the SETIO are both dominated by horizontal advection, while the SFF plays a secondary role. This result is consistent with former studies (Hu et al., 2019;Zhang et al., 2016Zhang et al., , 2018. The relevant physical processes contributing to variability in the horizontal advection include the advection of salinity anomalies from the WTP-MC region and anomalous local surface currents. Both of these two contributors have a strong ENSO-related signal and thus induce dramatic changes in the MLS in the SSS maximum region and the SETIO during strong and long-lasting ENSO events. During the Argo period, the MLS in those two regions decreased dramatically