Correlation between subsurface salinity anomalies in the Bay of Bengal and the Indian Ocean Dipole and governing mechanisms

Lead-lag correlations between the subsurface temperature/salinity anomalies in the Bay of Bengal (BoB) and the Indian Ocean Dipole (IOD) are revealed in model results, ocean synthesis, and observations. Mechanisms for such correlations are further investigated using the Hamburg Shelf Ocean Model (HAMSOM), mainly on the salinity variability. It is found that the subsurface salinity anomaly of the BoB positively correlates to the IOD with a lag of three months on average, while the subsurface temperature anomaly negatively correlates. The model results suggest the remote forcing from the equatorial 5 Indian Ocean dominates the interannual subsurface salinity variability in the BoB. The coastal Kelvin waves carry signals of positive (negative) salinity anomalies from the eastern equatorial Indian Ocean and propagate counterclockwise along the coasts of the BoB during positive (negative) IOD events. Subsequently westward Rossby waves propagate these signals to the basin at a relatively slow speed, which causes a considerable delay of the subsurface salinity anomalies in the correlation. By analyzing the salinity budget of the BoB, it is found that the diffusion dominates the salinity changes near the surface, while 10 the advection dominates the subsurface; the vertical advection of salinity contributes positively to this correlation, while the horizontal advection contributes negatively. These results suggest that the IOD plays a crucial role in the interannual subsurface salinity variability in the BoB.

In order to examine the potential correlation between the surface temperature pattern in the tropical Indian Ocean and the subsurface salinity variability in the BoB on the interannual scale, four independent data sets (Table 1) are used in this study.

60
The first is the global quality controlled monthly ocean temperature and salinity objective analyses of version 4.2.1 of the Met Office Hadley Centre 'EN' series (Good et al., 2013), named as EN4. The second is an ocean synthesis, which is the German contribution of the Estimating the Circulation and Climate of the Ocean project GECCO2 (Köhl, 2015). The third is the free run of mixed resolution of MPI-ESM (Jungclaus et al., 2013) under historical condition, named as MPI-ESM-MR. The fourth is the Roemmich-Gilson Argo Climatology (Roemmich and Gilson, 2009), named as RG_Clim, which offers a basic description of the modern upper ocean based entirely on Argo data. Due to the limitation of the data period, monthly anomalies of RG_Clim are defined on its monthly climatology from 2004 to 2016. Monthly anomalies of the other three data sets are defined on their monthly climatology from 1971 to 2000.

Model setting
For the purpose of discussing the relevant processes and mechanisms, a regional ocean model is performed. The Hamburg 70 Shelf Ocean Model (HAMSOM) we applied in this study is a three-dimensional baroclinic primitive equation model based upon a semi-implicit numerical scheme (Backhaus, 1985;Pohlmann, 1996Pohlmann, , 2006. In contrast to explicit shelf sea models, the semi-implicit scheme proposed is faster and allows the simulation of the shelf and the deep ocean regions together without being limited by stability considerations for the free surface (Backhaus, 1985). The underlying primitive equations are defined in z-coordinates and Arakawa C-grid under the hydrostatic and Boussinesq assumption. For temperature and salinity, the 75 second order Lax-Wendroff scheme is applied for advection, the horizontal eddy viscosity is defined according to Smagorinsky diffusivity (Smagorinsky, 1963), and vertical viscosity is calculated using the Kochergin scheme (Pohlmann, 1996(Pohlmann, , 2006.   (1974,1978,1993,1997,2000) are identified by 85 the normalized dipole mode index (DMI) calculated from the data set itself, in the places with peaks above two and located around September (see Figure 7d). The distribution of composited SSTa presents a significant dipole mode in the tropical Indian Ocean, which indicates the MPI-ESM-MR can reproduce IOD events.
Sea level height, temperature, and salinity at lateral boundaries are monthly prescribed and derived from the oceanic part MPI-OM of MPI-ESM-MR. Atmospheric forcing, such as air temperature, cloud cover, precipitation, specific humidity, air 90 pressure, wind stress, and wind speed at the open boundary, are six-hourly prescribed and derived from the atmospheric part ECHAM6 of MPI-ESM-MR. Under the consideration of a large amount of freshwater input through river discharge to our research domain, the river discharge is also prescribed six-hourly derived from ECHAM6. We applied a bias correction for forcing parameters on the climatological scale in order to bring the climatology of our simulation closer to the reality, because they are extracted from a purely free run. Principally, the monthly climatology of the external forcing was corrected by 95 reference data by this bias correction procedure. The reference data for atmospheric forcing except air pressure are extracted from ERA5 (Hersbach et al., 2020). The air pressure kept unchanged since it only affects sea surface height due to the inverse barometer effect in HAMSOM, and its seasonal pattern matches well with the local monsoon system. The reference data for sea temperature and salinity are derived from the World Ocean Atlas 2018. The amplitude of river discharge was corrected by WaterGAP (Döll et al., 2003), and the location where the river discharge enters the ocean was also corrected. regional model is consistent with MPI-ESM-MR. Nevertheless, it is noteworthy that HAMSOM is a regional uncoupled ocean model, so no response of the atmosphere to the ocean is considered. Therefore, our model result must be treated as a pure response of the ocean to external signals rather than a two-way air-sea coupling simulation.
The simulation runs from 1951 to 2005 with 3 minutes time step and daily average output. In addition to typical output 105 like temperature, salinity, and velocity, six terms concerning salinity change rate related to the contribution of advection and diffusion at U, V, W directions are conducted separately. Terms estimated from monthly outputs may differ significantly from those directly outputted by an online calculation, especially when high-frequency changes occur (Hasson et al., 2013;Köhler et al., 2018). This so-called 'online analysis' avoids the problem of large residual when directly using monthly data and allows us to precisely close the salinity budget and track relevant exchange processes of salinity.

Model Validation
Before investigating the interannual variability and detailed mechanisms of subsurface salinity in the BoB, the HAMSOM result is validated on the climatological scale by comparing it with other data sets. Figure 2 shows their spatial pattern of climatological surface and subsurface salinity. River discharge and distribution affect this spatial pattern at the surface, as well as the saline water from the western boundary. Reasons for the subsurface salinity pattern are complicated, ocean circulation 115 and upwelling/downwelling systems may be involved. Both in the surface and in the subsurface, the climatological salinity from HAMSOM presents a gradient from southwest to northeast, which is more consistent with both individual data sets GECCO2 and EN4 than from MPI-ESM-MR. Hence it can be concluded that the bias correction we applied here improves our modeling by offering a more realistic climatological background. Noteworthy that the seasonality is one of the most crucial characteristics in the research region. For the overall monthly 120 climatology of salinity, HAMSOM results also show a reliable seasonal variability (Figure 3). At the surface, the significant seasonal salinity variability is supposed to be the consequence of freshwater flux variability caused by the monsoon. All five data sets show a consistent seasonality of the surface, which indicates the domination of monsoon in this region. The monthly climatological salinity of these data sets differs more for the subsurface than for the surface, the lack of subsurface observations and more complex subsurface thermodynamics and hydrodynamics can be the reasons.

125
The upper ocean circulation is also validated in two sections ( Figure 4 and Figure 5). In general, circulations from HAMSOM is in well agreement with those from MPI-ESM-MR and GECCO2. The direction of upper ocean currents is reversed in MJJ and NDJ, which indicates that the monsoon dominates the upper ocean flow field in the BoB. Given the higher model resolution and more accurate terrain, HAMSOM is expected to perform better in coastal areas. The western boundary current simulated by HAMSOM, also known as the East Indian Current, is stronger than that given by GECCO2 (Figure 4), which should be 130 attributed to the higher resolution.   The above validation indicates that HAMSOM model can reproduce reasonable climatological fields and is reliable to be used as a numerical approach to study the physical processes and their specific contributions for the BoB. The interannual variations simulated by HAMSOM are combinations of external signals from MPI-ESM-MR and internal variabilities produced by HAMSOM itself. Hence, it can be concluded that it is reasonable to discuss the interannual variability and corresponding 140 physical processes simulated by HAMSOM in the following sections.

Lead-lag correlation
Four individual data sets and the downscaling model results are used in this section to examine if there is a statistically significant relationship between the subsurface temperature/salinity anomalies of the BoB and the SSTa of the tropical Indian Ocean on the interannual scale, specifically, the IOD. The DMI describes the difference in SSTa between the western tropical 145 Indian Ocean (WTIO) and the southeastern tropical Indian Ocean (SETIO, see Figure 1). It has a strong correlation with the principal component of EOF2 in the tropical Indian Ocean and is considered to be a reliable representation of the IOD (Saji 8 https://doi.org/10.5194/os-2020-75 Preprint. Discussion started: 7 August 2020 c Author(s) 2020. CC BY 4.0 License. et al., 1999). The time series of DMI can indicate different phases of the IOD, so in this study, DMI also covers the meaning of IOD variability.
In order to focus on interannual variations, a 3-month running mean is applied on the monthly time series of DMI and At the same time, by comparing the correlation magnitudes obtained from different data sets, it is noticeable that the correlation is stronger when the data set shows a lower degree of freedom. For example, HAMSOM has less complexity than 165 MPI-EMS-MR because it is a regional ocean model not including the ocean-atmosphere feedback processes. GECCO2 has a higher degree of freedom because assimilation processes are included. EN4 shows the highest degree of freedom of these four data sets because it is based on observations. However, the lack of observations and the objective analysis method used in EN4 limits its capability to reproduce the interannual subsurface salinity variability in the BoB. There is a weight index from 0 to 1 in EN4 which states the total weighting given to the observation increments when forming this analyses, and the mean Lead-lag Pearson correlation coefficients between the DMI and the salinity anomalies of the BoB at different depths are shown in Figure 9. Three aspects shown by these data sets are noteworthy. First, the most significant positive correlation appears below 50 m, and it is possible to be as deep as 250 m. Second, the DMI is leading a few months. Third, no obvious 185 positive correlation is validated for the sea surface. These results suggest that the subsurface salinity anomalies of the BoB are indeed related to the IOD with a considerable delay. On average, their correlation reaches its maximum at a three-month delay.
The local intense wind-induced mixing and other surface factors that are not closely related to the IOD are the reasons for the upper 50 m of the BoB does not reflect this correlation.
relation between the subsurface salinity anomaly of the entire BoB and the IOD represented by the DMI is revealed by observations, ocean synthesis, and modeling. A similar but negative correlation is also revealed between the subsurface temperature anomaly of the BoB and the IOD. The correlation between them differs in different data sets and becomes smaller when the data set has a higher degree of freedom; this is because some other variations may obscure the correlation we are focusing on.
However, this correlation can still be detected in the observational data, and it is very significant in the model-related data.

4 Mechanisms
In the above analysis, we have determined and discussed the lead-lag correlation between the domain averaged subsurface salinity anomaly of the BoB and the IOD. In this section, we mainly study how are these two variabilities connected by analyzing the HAMSOM result. Besides the connecting mechanisms, related physical processes of BoB's responses to IOD events are also a subject of this section. Several reasons may result in changes in salinity anomalies of the BoB, for example, 200 the salinity redistribution within the BoB or the salinity exchange between the BoB and its surroundings. Whatever the reason is, it will eventually be reflected in the salinity advection and diffusion. In this manner, an online analysis of the salinity budget is used in this section.

Connecting mechanism
To figure out the general feature of the response of the subsurface salinity in the BoB to the IOD, we construct composites 205 of subsurface salinity anomalies during ASO, NDJ, FMA, and MJJ, of pIOD and nIOD years, respectively ( Figure 10). The nIOD years (1979,1988,1992,1998,2004) are defined similarly to pIOD years, but with valleys below minus two normalized anomaly(see Figure 7d). The pattern of subsurface salinity anomalies is opposite during pIOD and nIOD events in ASO, as well as in NDJ when IOD events end with the DMI is returning to 0. This opposite feature is getting weaker over time, which can be seen in FMA and MJJ. When a pIOD or nIOD event happened in the tropical Indian Ocean, areas near the BoB coasts  Five subareas are selected in order to further investigate the response of different areas to the IOD (Figure 10a, e). A similar plot as Figure 9 but for subareas is presented in Figure 11. The closer the subarea is to the eastern boundary, the stronger the Pearson correlation between the DMI and the local salinity anomaly. When the subarea is close to but not directly at the western boundary, the correlation is relatively weak, which indicates that the signal is trapped at the west boundary ( Figure   225 11d). Results from the subarea closest to the equator also show weaker correlations (Figure 11e), while results from the subarea that is far away from the equator but closer to the eastern boundary show stronger correlations (Figure 11c), suggesting that the signal propagates along the boundary rather than directly goes north. These features support our speculation that the interannual subsurface salinity variability in the BoB is connecting to the IOD through both coastal Kelvin waves and westward Rossby waves.

230
It is challenging to observe Rossby waves in the BoB if we only use monthly data because of the basin size. Therefore, daily data from HAMSOM is used for tracking Rossby waves. As the Hovmöller diagram of daily climatological subsurface salinity averaged over 10 • N to 12 • N shown in Figure 12a, a westward Rossby wave signal can be seen by the westward lowsalinity water. This signal takes approximately four months to cross the basin zonally, and from this, it can be estimated that its propagation speed is about 0.16 ms −1 . Low-salinity water already appears at the western boundary before the westward 235 Rossby wave had reached here. In May, even though the westward Rossby wave signal represented by the low-salinity water has reached the western boundary, the water at the western boundary is as salty as the water at the eastern boundary, demonstrating In pIOD and nIOD years, the propagation characteristics of coastal Kelvin waves and westward Rossby waves are essentially the same as the climatology, but carrying positive and negative anomalies, respectively (Figure 12b, c). These statistically 240 significant anomalies first appear at the eastern boundary, then at the western boundary, then in the basin interior, which indicates that the extreme IOD signal is propagated to the entire BoB by both coastal Kelvin waves and westward moving Rossby waves. Previous studies have demonstrated the dominant role of coastal Kelvin waves in sea level variability in the BoB, especially near the eastern and northern boundaries (Han and Webster, 2002;Cheng et al., 2013). Our analysis about subsurface salinity anomalies suggests that the coastal Kelvin waves also dominate the western boundary when extreme IOD 245 events occur. The positive anomalies associated with pIOD reduce the zonal gradients of subsurface salinity, while the negative anomalies associated with nIOD increase their zonal gradients and result in different baroclinic Rossby wave modes with different propagating speed.
We also calculated the correlation between the local wind curl and salinity in the BoB and all subareas (not shown). The results show that these two parameters are strongly correlated on the seasonal scale, but there is no significant correlation on the 250 interannual scale, which indicates that the interannual signal is not from the wind field. Therefore, the model result suggests that the propagation process through coastal Kelvin waves and westward Rossby waves is the primary connecting mechanism of the delayed positive correlation between the subsurface salinity anomaly of the BoB and the zonal SSTa gradient in the tropical Indian Ocean. The interannual variability of thermocline depth in the eastern Indian Ocean is dominated by equatorial Indian Ocean winds, which drive eastward moving equatorial Kelvin waves that are blocked at the Sumatra-Java coasts (Du et al.,255 2012; Chen et al., 2015). It has been shown that enhanced upwelling occurs in the eastern Indian Ocean during pIOD years . This enhanced upwelling signal is converted into coastal Kelvin waves, which propagate counterclockwise along the boundary of the BoB. Subsequently, this signal is reflected at the eastern boundary forming westward moving Rossby waves that keep propagating into the central basin. During nIOD events the related subsurface anomalies in the BoB are modulated in a similar way.

Contributions of advection and diffusion
By outputting terms concerning salinity change rate related to the contribution of advection and diffusion, we can precisely close the salinity budget and analyze changes of advection and diffusion in the BoB in different IOD phases. The salinity budget can be written as follows: horizontal diffusion term. The sum of these large advection terms becomes much smaller and about the same magnitude as the salinity tendency and the vertical diffusion ( Figure 13b). All advection terms show significant differences in pIOD and nIOD events comparing to the climatology. Especially during nIOD, the salinity changes caused by advection at each direction increase, which is believed to be the result of increased zonal subsurface salinity gradients associated with nIOD. Meanwhile,

275
it can be seen that at this depth, on the average for the entire BoB, the vertical advection contributes positively to the positive salinity tendency of pIOD and the negative salinity tendency of nIOD. In contrast, the summed horizontal advection contributes negatively. Figure 14 shows the sum of advection terms, the sum of diffusion terms, and the final salinity tendency at different depths of the BoB and other selected subareas during ASO. The salinity tendency shows a subsurface salinity increase (decrease), corresponding color indicate that they are significant different at the 95% confidence level by a two-tailed Welch's t-test comparing to the climatology. The second and third row as in the first row, but for sum of domain-averaged diffusion terms and salinity tendency, respectively.

Conclusions
In this study, we have investigated the subsurface salinity variability in the BoB on the interannual scale and its relation with the IOD through multiple data sets, and we have also investigated the corresponding mechanisms through a regional ocean model simulation. The regional downscaling model successfully reproduces the reasonable climatology of the salinity and flow field, 295 proving its capability for investigating the physical processes in the BoB. In order to further discuss advection and diffusion contributions to salinity, we have performed an online analysis of salinity budget. This approach can precisely close the salinity budget, and hence, reflects the response of salinity in the BoB to the IOD, also with respect to the driving mechanisms in a quantitative manner.
A delayed positive correlation between the subsurface salinity anomaly of the BoB and the IOD was revealed by analyzing 300 their Pearson correlation coefficient. This correlation is not only shown in the modeling data but also in the ocean synthesis and observations. On average, a lag of three months shows the strongest correlation. Meanwhile, the correlation is relatively weaker when the data set shows a higher degree of freedom, suggesting that some processes exist in reality, which are not well resolved by numerical simulations that may disturb the relation between the subsurface salinity variability of the BoB and the IOD on the interannual scale. From this perspective, the numerical simulation is a more suitable method for investigating the 305 physical processes behind this correlation.
The model results suggested that the interannual subsurface salinity variability in the BoB and the IOD variability in the tropical Indian Ocean are connected by both coastal Kelvin waves and westward moving Rossby waves. First, coastal Kelvin waves carry the disturbance signal in the eastern equatorial Indian Ocean that is related to the IOD, propagating counterclockwise along the BoB coasts. subsequently, the signal reflects at the eastern boundary and propagates westward to the basin interior 310 in the form of Rossby waves. The main reason that the domain averaged subsurface salinity anomaly lags the DMI several months is that the westward Rossby waves travel slowly. The analysis of the salinity budget revealed that the contribution of advection plays a dominant role in this correlation. Particularly, the vertical advection shows a positive contribution, while the horizontal advection shows a negative contribution.
For the eastern equatorial Indian Ocean, the weakening of Wyrtki jet and the strengthening of upwelling caused by the 315 easterly wind anomalies during pIOD result in freshening at the surface and saltening at the subsurface (Kido and Tozuka, 2017). Large-scale wind stress anomalies play the dominant role in the salinity anomalies of this area during IOD events through modulating salinity advection mainly (Kido et al., 2019a). For the BoB, the model results suggest that the remote forcing from the equatorial Indian Ocean converted into coastal Kelvin waves and westward moving Rossby waves is the principal mechanism, which is responsible for the interannual salinity variability in the subsurface. Correlation analysis shows 320 that the subsurface salinity anomaly positively correlates with the IOD, while the subsurface temperature anomaly negatively correlates, which implies that the IOD remotely modulates the vertical advection in the BoB subsurface. The salinity budget of HAMSOM results proves that the vertical advection positively contributes to the correlation between the subsurface salinity anomaly of the BoB and the IOD. The decomposition of advective anomalies (Zhang et al., 2013;Li et al., 2016;Kido and Tozuka, 2017) will be helpful to understand the specific contribution of each specific process. For instance, it will allow 325 separately diagnose the contribution of the anomalous vertical salinity gradient and the contribution of the anomalous vertical velocity. During the pIOD phase, intensified upwelling occurs in the eastern Indian Ocean (Nyadjro and McPhaden, 2014;Chen et al., 2016), inducing an uplift of cold, more saline water along the eastern BoB coasts. This anomaly in turn induces coastal Kelvin waves, which are reflected at topographic disturbances, inducing Rossby waves that move westward to the central basin. Through this chain of processes, the remote forcing from the equatorial Indian Ocean is able to dominate the 330 interannual subsurface temperature/salinity variability in the BoB. The BoB is known as a region with vigorous mesoscale eddy activity (Chen et al., 2012(Chen et al., , 2018. How these eddies affect the evolution of subsurface salinity anomalies requires future studies. Based on this discovered correlation and related mechanisms, one application is using the DMI to predict the subsurface ocean state in the BoB. The subsurface ocean state affects the barrier layer and mixed layer depth, as well as the near-surface 335 state and the air-sea energy transfer. However, how the subsurface parameters response to the IOD affects its local upper ocean still needs more studies. Previous studies have demonstrated the importance of forcing from the equator to the BoB, such as sea surface height, thermocline, and circulation structure (Girishkumar et al., 2013;Chatterjee et al., 2017;Pramanik et al., 2019), especially for the mechanisms forcing the East India Coastal Current (Yu et al., 1991;McCreary et al., 1996;Shankar et al., 1996). The response of the subsurface salinity field we discussed may also affect the local flow field and mesoscale eddies.

340
Coastal Kelvin waves and westward Rossby waves play a vital role in the process of receiving information from the equator in the BoB, and similarly, they also play their role on the interannual scale. Especially during the nIOD phase, the increase of the zonal subsurface salinity gradients makes it easier to excite high mode westward Rossby waves, suggesting the effect of IOD on the subsurface thermohaline circulation in the BoB.
In addition to these, in the future we will investigate how the relationship between the BoB subsurface parameters and the 345 tropical Indian Ocean surface parameters will be affected by the impact of climate change. The sea surface is more susceptible to global warming, and the BoB subsurface has a notable connection to the tropical Indian Ocean surface. The sea surface warming may also affect the subsurface and even deeper through dynamic mechanisms. Therefore, how the BoB subsurface responds to climate change is the next subject we are going to study.