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 and 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 relating to the salinity variability. It is found that the subsurface salinity anomaly of the BoB positively correlates to the IOD, with a lag of 3 months on average, while the subsurface temperature anomaly correlates negatively. The model results suggest the remote forcing from the equatorial 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 diffusion dominates the salinity changes near the surface, while 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.

1 Introduction

The Bay of Bengal (BoB) is a monsoon-controlled tropical ocean located in the northeast of the Indian Ocean. The robust monsoon significantly influences the ocean circulation, vertical water exchange, and water characteristics in the BoB (Shetye et al., 1991, 1996; Vecchi and Harrison, 2002; Li et al., 2017). During the summer monsoon, the Southwest Monsoon Current brings saltier Arabian Sea water into the BoB, whereas the Northeast Monsoon Current brings fresher water from the BoB to the Arabian Sea during the winter monsoon (Vinayachandran et al., 1999; Jensen, 2001; Sanchez‐Franks et al., 2019). Water from the Arabian Sea also enters the BoB as a subsurface flow during the northeast monsoon, which is proven by observations and model works (Wijesekera et al., 2015; Gordon et al., 2016). The salinity exchanges between the BoB and the equatorial Indian Ocean also show a seasonality associated with the monsoon (Jensen et al., 2016; Trott et al., 2019). In addition to the local monsoon, remote forcing from the Equator also affects the ocean circulation and thermocline in the BoB, by which equatorial signals pass through the Andaman Sea (Potemra et al., 1991; Yu et al., 1991; McCreary et al., 1993, 1996; Girishkumar et al., 2013). The Andaman and Nicobar Islands, as well as the eastern border of the Andaman Sea, significantly alter the circulation in the BoB (Chatterjee et al., 2017). A recent numerical study suggests that the equatorial forcing plays a dominant role in interannual variations of sea surface height and thermocline in the BoB during Indian Ocean Dipole (IOD) and El Niño–Southern Oscillation (ENSO) events, especially relating to their spatial pattern (Pramanik et al., 2019).

The IOD is an east–west dipole mode that dominates the interannual sea surface temperature (SST) variability in the tropical Indian Ocean (Saji et al., 1999; Webster et al., 1999; Schott et al., 2009; Deser et al., 2010), and it is a physical entity independent of the ENSO (Ashok et al., 2003; Fischer et al., 2005). The spatiotemporal coupling among ocean dynamics, SST, winds, and rainfall revealed by the IOD have inspired many studies regarding the relationship and processes between the IOD and variations of surface and subsurface temperature and salinity in the tropical Indian Ocean (Rao et al., 2002; Shinoda et al., 2004; Thompson et al., 2006; Grunseich et al., 2011; Du et al., 2012; Zhang et al., 2013; Sayantani and Gnanaseelan, 2015; Kido and Tozuka, 2017; Kido et al., 2019a). The research on both sea level and annual mean subsurface temperature anomalies revealed a seesaw in the thermocline that related to the IOD (Saji et al., 1999). During the positive IOD (pIOD) phase, westerly winds weaken, allowing cold water to rise in the eastern equatorial Indian Ocean and warm water to move toward to the west, which therefore lifts the equatorial thermocline in the east; during the negative IOD (nIOD) phase, this process is reversed and thus lifts the equatorial thermocline in the west.

Recently it has been reported that the Wyrtki Jets (Wyrtki, 1973) affect the salinity balance of the BoB by means of their northward bifurcation (Wang, 2017). Furthermore, it could be shown that on an interannual timescale the Wyrtki Jets are highly associated with the IOD. This holds particularly true for the fall jet, which develops in the peak phase of the IOD (Nyadjro and McPhaden, 2014; McPhaden et al., 2015). Previous studies have discussed the impact of IOD on the subsurface dynamics in the equatorial Indian Ocean. However, mechanisms and quantitative understandings of the impact of IOD on subsurface dynamics in the BoB have not been established yet, especially for the impact on subsurface salinity. Subsurface salinity is of great importance for determining the ocean barrier layer and mixed-layer depth (Lukas and Lindstrom, 1991; Montégut et al., 2007; Li et al., 2018; Kido et al., 2019b). Understanding the variations and dynamics of subsurface salinity is helpful for understanding the evolution of stratification and upper-ocean properties, in order to further understand the response of the ocean to the atmosphere and its role on climate. Nevertheless, it is not clear how the subsurface salinity in the BoB varies and whether it is affected by the IOD. The questions addressed here are as follows: first, is there an identifiable correlation between the subsurface salinity in the BoB and the surface temperature in the tropical Indian Ocean on the interannual scale? Second, how are these two variabilities related? Third, how does the subsurface salinity in the BoB respond to the IOD?

To answer the above questions, the subsurface salinity variability in the BoB and its relation with the IOD and the corresponding mechanisms are investigated in this paper. Unless otherwise specified, anomalies used in this paper are residuals subtracting monthly climatology from monthly data. The rest of this paper is organized as follows. In Sect. 2, we introduce four data sets and a regional ocean model used in this study, and the model validation is also presented in this section. In Sect. 3, we examine the correlation between the subsurface temperature and salinity anomaly of the BoB and the IOD by analyzing the four independent data sets and the model results. Connecting mechanisms and contributions of advection and diffusion are discussed in Sect. 4. Section 5 gives the summary and discussion.

2 Data and model

2.1 Data sets

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 an interannual scale, four independent data sets (Table 1) are used in this study. 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), hereafter 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 the mixed resolution of MPI-ESM (Jungclaus et al., 2013) under historical conditions, hereafter MPI-ESM-MR. The fourth is the Roemmich–Gilson Argo Climatology (Roemmich and Gilson, 2009), hereafter 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 using its monthly climatology from 2004 to 2016. Monthly anomalies of the other three data sets are defined using their monthly climatology from 1971 to 2000.

Table 1 Summary of data sets.

2.2 Model setting

For the purpose of discussing the relevant processes and mechanisms, a regional ocean model is performed. The Hamburg 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, 1996, 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 using z coordinates and an Arakawa C grid under hydrostatic and Boussinesq assumptions. For temperature and salinity, a 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, 2006).

In principle, we perform a dynamic downscaling simulation on the model domain using HAMSOM, with the external forcing derived from the MPI-ESM-MR historical scenario. The model domain (Fig. 1b) covers the Bay of Bengal and the Andaman Sea, ranging meridionally from 0 to 22.83^∘ N and zonally from 77.4 to 103.5^∘ E, with bathymetry derived from SRTM30_PLUS (Becker et al., 2009). The horizontal model resolution is set to ${\mathrm{5}}^{\prime }\times {\mathrm{5}}^{\prime }$. A total of 58 model layers are specified in vertical, of which there are 26 layers over upper 200 m and 33 layers over the upper 400 m. To stabilize the inner model domain, a sponge layer is implemented along the lateral open boundaries in order to damp disturbances arising from inconsistencies within the prescribed boundary conditions extracted from the MPI-ESM-MR. Hence, equatorial processes that are not directly resolved in our model domain are still able to enter the inner domain through the prescription of open boundary conditions. A correlation analysis (not shown) could demonstrate that the sponge layer does not block the propagation of low-frequency signals (seasonal scale and below). Therefore, only the BoB region (marked in Fig. 1b) is analyzed in our HAMSOM simulation.

Figure 1a offers an overview of our research domain and the tropical Indian Ocean and also shows a composite of sea surface temperature anomalies (SSTa) during August–October (ASO, for ease of presentation, the months in the following are simplified to initials) of pIOD years from MPI-ESM-MR. The pIOD years (1974, 1978, 1993, 1997, 2000) are identified by the normalized dipole mode index (DMI) calculated from the data set itself in Fig. 7d. This distribution of composited SSTa (Fig. 1a) shows a significant dipole mode in the tropical Indian Ocean, which is in good agreement with previous studies (Webster et al., 1999; Deser et al., 2010). The corresponding DMI time series (Fig. 7d) exhibits reasonable interannual variation characteristics. These indicate that the global model used in this downscaling study can realistically reproduce IOD events.

Figure 1 Composite of sea surface temperature anomalies during ASO of pIOD years from MPI-ESM-MR (a). The contour intervals are 0.2 ^∘C. Anomalies significant at the 95 % confidence level using a two-tailed Welch's t test are hatched with grey dots. The western tropical Indian Ocean (WTIO) and southeastern tropical Indian Ocean (SETIO), the two areas related to the dipole mode index (DMI), are marked with black boxes. The bathymetry used in the downscaling simulation is shown in (b). Our research domain, the Bay of Bengal (BoB), is marked with black borders in (a) and (b).

Sea level height, temperature, and salinity at lateral boundaries are prescribed monthly and derived from the oceanic part MPI-OM of MPI-ESM-MR. Atmospheric forcing, such as air temperature, cloud cover, precipitation, specific humidity, air pressure, wind stress, and wind speed at the open boundary, are prescribed every 6 h and derived from the atmospheric part of MPI-ESM-MR (ECHAM6). Under the consideration of a large amount of freshwater input through river discharge to our research domain, the river discharge is also prescribed every 6 h and 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 reference data via this bias-correction procedure. The reference data for atmospheric forcing except air pressure are extracted from ERA5 (Hersbach et al., 2020). The air pressure is 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. Although we applied the bias correction, the interannual signals from MPI-ESM-MR have not been changed, and thus the IOD signal input to our regional model is consistent with MPI-ESM-MR. Nevertheless, it is noteworthy that HAMSOM is a regional uncoupled ocean model, which means that 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 min time step and daily average output. In addition to typical output like temperature, salinity, and velocity, six terms concerning salinity change rate related to the contribution of advection and diffusion in the 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 residuals when directly using monthly data and allows us to precisely close the salinity budget and track relevant exchange processes of salinity.

2.3 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 its 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 and upwelling or downwelling systems may be involved. Both at 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.

Figure 2 Spatial pattern of climatological surface (a, b, c, d) and subsurface (100 m; e, f, g, h) salinity from HAMSOM, MPI-ESM-MR, GECCO2, and EN4, respectively. The contour intervals are 0.3 psu for the surface but 0.1 psu for the subsurface.

It is noteworthy that seasonality is one of the most crucial characteristics in the research region. For the overall monthly climatology of salinity, HAMSOM results also show a reliable seasonal variability (Fig. 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 averaged Pearson correlation between each line shown in Fig. 3b is 0.63, 0.42, 0.60, 0.73, and 0.68 for EN4, GECCO2, RG_Clim, MPI-ESM-MR, and HAMSOM, respectively. The lack of subsurface observations and more complex subsurface thermodynamics and hydrodynamics can all be the reason for this difference.

Figure 3 Normalized monthly climatology of surface (a) and subsurface (b) salinity of the BoB from EN4, GECCO2, RG_Clim, MPI-ESM-MR, and HAMSOM, respectively. The standard deviation σ corresponding to each data set is labeled.

The upper-ocean circulation is also validated in two sections (Figs. 4 and 5). In general, circulations from HAMSOM are in good 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 (Fig. 4), which should be attributed to the higher resolution.

Figure 4 Depth–longitude section of climatological V velocity (averaged over 10 to 12^∘ N) during MJJ (a, b, c) and NDJ (d, e, f) from HAMSOM, MPI-ESM-MR, and GECCO2, respectively. The contour intervals are 0.03 m s⁻¹.

Figure 5 Depth–latitude section of climatological U velocity (averaged over 88 to 90^∘ E) during MJJ (a, b, c) and NDJ (d, e, f) from HAMSOM, MPI-ESM-MR, and GECCO2, respectively. The contour intervals are 0.03 m s⁻¹.

Figure 6 shows the Taylor diagram (Taylor, 2001) of the surface and subsurface salinity from HAMSOM and other data sets. In this Taylor diagram, the standard deviation reflects both the temporal variability and the spatial variability. The surface salinity standard deviation of HAMSOM is consistent with the observation-based EN4, indicating the realistic extent of HAMSOM in realistically simulating the amplitude of variations. Even though HAMSOM has a finer grid than EN4, the sea surface feature simulated by HAMSOM is largely determined by the coarser atmospheric forcing due to our simulation strategy, so a good agreement of surface salinity variability between HAMSOM and the reference data set is expected. For the subsurface salinity, the standard deviation of HAMSOM is larger than EN4. Considering that the high horizontal resolution of HAMSOM data allows more mesoscale features, while the low resolution of other data sets does not, the relatively large standard deviation of HAMSOM subsurface salinity is acceptable since it shows more spatial variabilities. The root-mean-square difference (RMSD) is often used to quantify differences in two fields. Compared to GECCO2 and MPI-ESM-MR, HAMSOM shows a relatively large difference to the reference data set EN4, and the resolution difference could be the reason. Overall, HAMSOM-simulated salinity variabilities are in a reasonable range.

Figure 6 The Taylor diagram of (a) the surface salinity and (b) the subsurface salinity from 1971 to 2000 for different data sets. The observation-based EN4 (pentagon) is chosen as the reference data set. Grey lines indicate the centered RMSD from the reference data set.

The above validation indicates that HAMSOM model can reproduce reasonable climatological fields and that it is reliable for use as a numerical approach to study the physical processes and their specific contributions to 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 physical processes simulated by HAMSOM in the following sections.

3 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 and salinity anomalies of the BoB and the SSTa of the tropical Indian Ocean (specifically the IOD) on an interannual scale. The DMI describes the difference in SSTa between the western tropical Indian Ocean (WTIO) and the southeastern tropical Indian Ocean (SETIO; see Fig. 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 et al., 1999). The time series of DMI can indicate different phases of the IOD, and thus 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 the domain-averaged subsurface temperature and salinity anomalies of the BoB. Normalized time series from HAMSOM, MPI-ESM-MR, GECCO2, and EN4 are shown in the first column of Fig. 7. Subsurface is defined at a depth of 100 m, where wind-induced mixing is negligible, while upwelling and downwelling still plays a role. The DMI time series of HAMSOM is extracted from MPI-ESM-MR. DMI extracted from different data sets have similar interannual variation characteristics, but they do not precisely match. EN4 and GECCO2 both capture some typical pIOD events well, for example in the years 1994 and 1997. The free run of MPI-ESM-MR shows a reasonable amplitude and interannual variations and does not exactly repeat the positive event in 1994, which is to be expected, since in this historical run only the statistical features have to be consistent. Lead–lag running Pearson correlation coefficients with a window of 30 years between the DMI and the domain-averaged subsurface temperature and salinity anomaly of the BoB (as indicated in Fig. 1b) are calculated. The value of the Pearson correlation determines the extent of linearity between two variables. All of these four data sets show that the subsurface temperature anomaly of the BoB negatively correlates to the DMI with a notable lag of about 3 months on average. Results from HAMSOM, MPI-EMS-MR, and GECCO2 also show a similar but positive correlation between the subsurface salinity anomaly of the BoB and the DMI, but results from EN4 do not show a correlation.

Figure 7 Normalized 3-month running mean of DMI, temperature anomaly, and salinity anomaly at subsurface (100 m) of the BoB from HAMSOM (a), MPI-ESM-MR (d), GECCO2 (g), and EN4 (j), respectively, are shown in the first column. The standard deviation σ is labeled with the corresponding color. Lead–lag running Pearson correlation coefficient with a window of 30 years between the DMI and the subsurface temperature (salinity) anomaly from each data set is shown in the second (third) column, respectively. The contour intervals are 0.1. Only significant correlation coefficients with p value <0.05 are shaded.

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. The term “degree of freedom” is used here to describe the inherent complexity of a data set, and this complexity is mainly determined by the number of processes involved in the data set itself. For example, HAMSOM has a lower degree of freedom than MPI-EMS-MR because it is a regional ocean model that does not include the ocean–atmosphere feedback processes. GECCO2 has a higher degree of freedom because assimilation processes are included. EN4 is supposed to have the highest degree of freedom of these four data sets because it is based on observations. The difference in correlation magnitudes between them can be explained by the difference in their own degrees of freedom. Therefore, it can be reasonably inferred that this correlation does exist in the real ocean, but at the same time there are some processes in reality that obscure the correlation. In this sense, our HAMSOM simulation is more suitable to investigate the physical processes behind the presented correlation.

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 that states the total weighting given to the observation increments when forming this analyses, and the mean weight of subsurface temperature and salinity for the BoB is 0.57 and 0.22, respectively, which points out the lack of salinity observations in the BoB. For example, there are almost no observations from 1951 to 1956 for subsurface salinity of the BoB, and thus the subsurface salinity anomaly shows artificial oscillations during this period (Fig. 7j).

The situation of the lack of subsurface salinity observations in the BoB has improved since 2000, especially with the development of Argo. We present related time series and their lead–lag relation calculated from RG_Clim in Fig. 8. These time series support the correlations described by the other four data sets, and the subsurface salinity anomaly positively correlates with the DMI. The DMI leads the subsurface temperature anomaly, which can be seen for 100 and 150 m depth. Such a leading relationship corresponding to the subsurface salinity anomaly suggested by model-related results does not show for 100 m, but a broad, positive correlation with a peak value close to 0.4 is clear. Moreover, at 150 m depth, the DMI leads the salinity anomaly for 4 months with a peak correlation of over 0.4. Although the time length of RG_Clim is not as long as the other data sets, this Argo-based data clearly shows that the subsurface salinity anomaly of the BoB is correlated to the zonal gradient of SSTa in the tropical Indian Ocean.

Lead–lag Pearson correlation coefficients between the DMI and the salinity anomalies of the BoB at different depths are shown in Fig. 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 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 3-month delay. The local intense wind-induced mixing and other surface factors that are not closely related to the IOD are the reasons that the upper 50 m of the BoB does not reflect this correlation.

Figure 8 Normalized 3-month running mean of DMI, temperature anomaly, and salinity anomaly at the subsurface of the BoB from Argo-based RG_Clim are shown in (a). The standard deviation σ is labeled with the corresponding color. Their respective lead–lag relations described by the Pearson correlation coefficient between the DMI and subsurface anomalies are shown in (b) for 100 m and (c) of 150 m. Only significant correlation coefficients with p value <0.05 are shaded.

Figure 9 Lead–lag Pearson correlation coefficient between the DMI and salinity anomalies of the BoB at different depths from HAMSOM (a), MPI-ESM-MR (b), GECCO2 (c), and RG_Clim (d), respectively. The contour intervals are 0.1. The analysis period is from 1960 to 2005 for (a), from 1951 to 2005 for (b) and (c), and from 2004 to 2018 for (d), respectively. Only significant correlation coefficients with p value <0.05 are shaded.

By analyzing time series and their Pearson correlation coefficients, a time delay of about 3 months and positive correlation 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 these two variabilities are 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 factors may result in changes in salinity anomalies of the BoB: for example, 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.

4.1 Connecting mechanism

To figure out the general feature of the response of the subsurface salinity in the BoB to the IOD, we construct composites of subsurface salinity anomalies during ASO, NDJ, FMA, and MJJ of pIOD and nIOD years (Fig. 10). The nIOD years (1979, 1988, 1992, 1998, 2004) are defined similarly to pIOD years but with valleys below −2 (see Fig. 7d). The pattern of subsurface salinity anomalies is the opposite during pIOD and nIOD events in ASO and in NDJ when IOD events end with the DMI returning to 0. This opposite feature becomes 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 first show large and statistically significant anomalies (Fig. 10a, e). Next, these large and statistically significant anomalies show in most areas of the eastern basin but are limited to the western boundary areas (Fig. 10b, f). This developing process is consistent with the characteristics of the coastal Kelvin wave and westward Rossby waves. The Welch's t test is designed for working with small samples. Of course more IOD events would strengthen the power of the test. However, from our results, the five IOD events are sufficient to show the significant differences between IOD years and climatological conditions and the propagation of these IOD-related waves in a statistical meaning.

We speculate the propagation process is as follows. First, the subsurface disturbance signals in the eastern equatorial Indian Ocean related to the IOD propagate counterclockwise along the BoB coasts in the form of coastal Kelvin waves. The estimated wave phase speed is about 2.65 m s⁻¹, and thus it takes approximately 2 weeks to propagate from the Equator to the northern BoB (Moore and McCreary, 1990; Cheng et al., 2013). These coastal Kelvin waves travel quickly, explaining why the related significant anomalies first show up near the coasts. Subsequently, these signals are reflected at the eastern boundary and propagate westward into the interior of the basin, since the phase speed of Rossby waves is predominantly moving westward. This also explains why the signal near the western boundary seems to be trapped there. From the significant area of influence in NDJ, these Rossby waves travel slower, which accounts for the domain-averaged subsurface salinity anomaly lags in the DMI.

Figure 10 Composite of subsurface (100 m) salinity anomalies during ASO (a, e), NDJ (b, f), FMA (c, g), and MJJ (d, h) of pIOD and nIOD years, respectively, from HAMSOM. The contour intervals are 0.1 psu. Anomalies significant at the 95 % confidence level using a two-tailed Welch's t test are hatched with grey dots. Selected subareas are marked with magenta boxes in (a) and (e).

Five subareas are selected in order to further investigate the response of different areas to the IOD (Fig. 10a, e). Figure 11 shows the typical salinity profiles and the equivalent vertical displacement, yielding a 0.6 psu salinity anomaly in the BoB and five subareas. These salinity profiles of subareas indicate the distribution of salinity stratification and their seasonal changes in our research domain. This distribution is consistent with our model validation shown in Fig. 2. Due to the influence of river discharge and distribution, as well as the saline water from the open lateral boundary, the salinity stratification shows a weakening gradient from northeast to southwest. As shown in Fig. 10, averaged subsurface salinity anomaly of extreme IOD events can reach to about 0.6 psu in areas close to the coast. As indicated in Fig. 11, a 0.6 psu salinity anomaly is equivalent to about 20–50 m vertical displacement. Apparently, these displacements are already significant for vertical water motions and would be larger for some extreme IOD events since the results shown here represent a mean state.

Figure 11 Domain-averaged climatological vertical salinity profiles of the BoB (a) and subareas (b, c, d, e, f) during ASO (red thick line) and NDJ (blue thick line). Solid thin lines with corresponding colors indicate the salinity of 100 m. Dashed thin lines with corresponding colors indicate the equivalent vertical displacement, yielding a 0.6 psu salinity anomaly.

A similar plot to Fig. 9 but for subareas is presented in Fig. 12. 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 western boundary (Fig. 12d). Results from the subarea closest to the Equator also show weaker correlations (Fig. 12e), while results from a subarea that is far away from the Equator but closer to the eastern boundary show stronger correlations (Fig. 12c), suggesting that the signal propagates along the boundary rather than directly going north. The lags for these subareas also indicate that the subarea SAS is affected first, followed by NBB and EBB (subareas are define in Fig. 10). 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.

Figure 12 Lead–lag Pearson correlation coefficient between the DMI and the salinity anomaly of subareas SAS (a), EBB (b), NBB (c), WBB (d), and SBB (e) at different depths from HAMSOM. The contour intervals are 0.1. The analysis period is from 1960 to 2005. Only significant correlation coefficients with p value <0.05 are shaded. Solid purple lines indicate the depth-averaged correlation over 50–150 m. Dashed purple lines indicate the highest correlation and the corresponding lag.

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 are used for tracking Rossby waves. As the Hovmöller diagram of daily climatological subsurface salinity averaged over 10 to 12^∘ N shows in Fig. 13a, a westward Rossby wave signal can be seen by the westward low-salinity water. This signal takes approximately 4 months to cross the basin zonally, and from this it can be estimated that its propagation speed is about 0.16 m s⁻¹. Low-salinity water already appears at the western boundary before the westward Rossby wave has arrived there. 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 that coastal Kelvin waves travel faster and dominate the coastal zone in the BoB.

Figure 13 Hovmöller diagram of daily climatological subsurface (100 m) salinity (a; averaged over 10 to 12^∘ N) from HAMSOM. The intervals are 0.1 psu. Panels (b) and (c) are the same as (a) but for the composite of subsurface salinity anomalies for pIOD (b) and nIOD (c), respectively. The intervals are also 0.1 psu. Anomalies significant at the 95 % confidence level by a two-tailed Welch's t test are hatched with grey dots.

In pIOD and nIOD years, the propagation characteristics of coastal Kelvin waves and westward Rossby waves are essentially the same as the climatology but carry positive and negative anomalies, respectively (Fig. 13b, c). These statistically significant anomalies first appear at the eastern boundary, then at the western boundary, and finally 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 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 speeds.

We also calculated the correlation between the local wind and salinity in the BoB and all other 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 interannual scale. This indicates that interannual signals are not primarily induced by coastal Kelvin waves forced by the local wind along the east coast of the Andaman Sea but instead by far-field signals originating from the Equator. 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 Sumatran–Javan coasts (Du et al., 2012; Chen et al., 2015). It has been shown that enhanced upwelling occurs in the eastern Indian Ocean during pIOD years (Chen et al., 2016). 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.

4.2 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:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial S}{\partial t}}& =-u{\displaystyle \frac{\partial S}{\partial x}}-v{\displaystyle \frac{\partial S}{\partial y}}-w{\displaystyle \frac{\partial S}{\partial z}}+{\displaystyle \frac{\partial }{\partial x}}\left({\mathit{\kappa }}_{\mathrm{H}}{\displaystyle \frac{\partial S}{\partial x}}\right)\\ & +{\displaystyle \frac{\partial }{\partial y}}\left({\mathit{\kappa }}_{\mathrm{H}}{\displaystyle \frac{\partial S}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{V}}{\displaystyle \frac{\partial S}{\partial z}}\right),\end{array}\end{array}$$

where S is salinity; u, v, and w are zonal, meridional, and vertical velocity, respectively; and κ_H and κ_V are horizontal and vertical diffusion coefficients, respectively. The left side represents the salinity tendency (ST), while the right side, from left to right, represents the salinity change rate of zonal (UADV), meridional (VADV), and vertical (WADV) advection and of zonal (UDIF), meridional (VDIF), vertical (WDIF) diffusion.

A salinity budget of the BoB at 100 m during ASO is presented in Fig. 14. For terms on the right side of the Eq. (1) at this depth, the advection term is much larger than the diffusion term, and the vertical diffusion term is larger than the 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 (Fig. 14b). All advection terms show significant differences in pIOD and nIOD events compared to the climatology. The salinity changes caused by advection at each direction increase, especially during nIOD, which is believed to be the result of increased zonal subsurface salinity gradients associated with nIOD. Meanwhile, it can be seen that at this depth, on 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 Domain-averaged subsurface (100 m) salinity tendency and related salinity change rate terms of the BoB during ASO of pIOD years, nIOD years, and the climatological period, respectively, from HAMSOM (a). The values of ST and WDIF in different cases are labeled with the corresponding color. The sum of all advection terms and the sum of all diffusion terms are shown in (b). Dots with corresponding color indicate that they are significantly different at the 95 % confidence level using a two-tailed Welch's t test compared to the climatology.

Figure 15 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), indicating a positive (negative) anomaly for pIOD (nIOD) years, which in this season shows an obvious response to the IOD signal (Fig. 15m, n, o). The results in different regions show that the salinity tendency is dominated by diffusion near the surface, while it is dominated by advection for the subsurface. The dominant role of diffusion near the surface can be explained by the wind-induced mixing. The salinity change rate due to advection shows more obvious responses in the subsurface during both extreme IOD events, especially in the subareas SAS and EBB and the entire BoB, suggesting that the correlation we are discussing is mainly caused by advection processes.

Figure 15 Sum of domain-averaged advection terms (in psu per month) at different depths of the BoB (a) and its subareas (b, c, d, e, f) during ASO. The solid black line is for the climatology, and the dashed red and blue lines are for the composite of pIOD and nIOD years, respectively. Dots with corresponding color indicate that they are significantly different at the 95 % confidence level using a two-tailed Welch's t test compared to the climatology. The second and third row are the same as the first row but for sum of domain-averaged diffusion terms and salinity tendency, respectively.

As we analyzed through multiple data sets, there exists a delayed positive correlation between the subsurface salinity anomaly of the BoB and the zonal SSTa gradient in the tropical Indian Ocean represented by the DMI. Therefore, by analyzing the salinity budget of the BoB, the model results suggest that the contribution of advection plays a dominant role in this correlation. The vertical advection contributes positively, while the horizontal advection contributes negatively to the correlation stated above.

5 Conclusions

In this study, we have investigated the subsurface salinity variability in the BoB on an 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, 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 in a quantitative manner with respect to the driving mechanisms.

A delayed positive correlation between the subsurface salinity anomaly of the BoB and the IOD was revealed by analyzing 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 3 months shows the strongest correlation. Meanwhile, the correlation is relatively weak when the data set shows a higher degree of freedom, suggesting that some processes exist in reality that 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 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 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. 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 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 mainly through modulating salinity advection (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 responsible for the interannual salinity variability in the subsurface. Because of the unique topographic configuration, the BoB is more susceptible to equatorial signals than any other ocean region. Equatorial signals carried by equatorial Kelvin waves propagate eastward to the western coast of Sumatra, where coastal Kelvin waves are derived, and in turn influence the BoB (Cheng et al., 2013). Correlation analysis shows 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 the separate diagnosis of 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 interannual subsurface temperature and salinity variability in the BoB. The BoB is known as a region with vigorous mesoscale eddy activity (Chen et al., 2012, 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 state and the air–sea energy transfer. However, how the subsurface parameters response to the IOD affects its local upper ocean still needs more study. 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. Coastal Kelvin waves and westward Rossby waves play a vital role in the process of receiving information from the Equator in the BoB, and they also similarly 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. The biological processes are significantly affected by the salinity stratification and vertical mixing in the BoB (Prasanna Kumar et al., 2002). As our results show, the IOD significantly modulates the BoB subsurface salinity and further potentially affects the ocean barrier layer and mixed layer depth through the Kelvin and Rossby waves. Therefore, the correlation and the corresponding processes we discussed are expected to play an important role in the biology of the BoB. For example, the considerable IOD-related vertical displacement may transport nutrients across the halocline and then increase biological production, as the eddy pumping does (Prasanna Kumar et al., 2004). Furthermore, it is also expected that these waves affect the air–sea exchange processes in the BoB, which in turn influence the remote ocean feedback to the atmosphere.

In addition to these, in the future we will investigate how the relationship between the BoB subsurface parameters and the tropical Indian Ocean surface parameters will be affected by the impacts 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 areas through dynamic mechanisms. Therefore, how the BoB subsurface responds to climate change is the next subject we are going to study.