Impact of Indian Ocean surface temperature gradient reversals on the Indian Summer Monsoon

a Department of Earth Science, University of California, Santa Barbara, CA 93106-9630, USA b Bundesanstalt für Geowissenschaften und Rohstoffe, Stilleweg 2, 30655 Hannover, Germany c Department of Geography, University of California, Santa Barbara, CA 93106-4060, USA d School of Energy, Geoscience, Infrastructure and Society, Lyell Centre, Heriot-Watt University, Edinburgh, UK e A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, 119017, Moscow, Russia f Institut für Geowissenschaften, Universität Kiel, D-24118 Kiel, Germany g Laboratoire des Sciences du Climat et de l’Environnement (LSCE/IPSL), Université Paris-Saclay, Gif-sur-Yvette, France

Indian Summer Monsoon (ISM) precipitation is the main determinant of livelihood in a densely populated world region. The interannual variability of the ISM is influenced by several modes of climate variability, including anomalous seasonal sea surface temperature (SST) gradient reversals between the eastern, western, and northeastern Indian Ocean. With global warming, the frequency of zonal and meridional Indian Ocean's SST gradient changes is projected to increase but its impact on the ISM is debated. Here we present a 25,000-year proxy record of SST and inferred Ganges-Brahmaputra-Meghna (GBM) River runoff that provides a spatially integrated measure of ISM precipitation changes. This record indicates a monotonic deglacial strengthening of the ISM system when the SST gradient between the Bay of Bengal surface water and the eastern equatorial Indian Ocean was reversed. We posit that the reversal in the meridional SST gradient reduced the impact of Heinrich Event 1 and Younger Dryas on the low elevation part of the ISM domain. Furthermore, the proxy record shows that the strongest Holocene ISM strengthening occurred between 7900±470 and 5700±360 years before present, coinciding with and causally linked to the reversal of the Indian Ocean zonal SST gradient and ensuing changes in the wind fields, a sequence of events that is inferred from and supported by the results of our climate simulation. Our study demonstrates that changes in the Indian Ocean's zonal and meridional thermal gradient strongly shaped the timing of Holocene monsoon strengthening and the response of ISM to the last deglacial freshwater forcing.

Indian Summer Monsoon (ISM) precipitation
The onset and northward progression of the Indian Summer Monsoon (ISM) precipitation begins, on average, in middle May/early June and is triggered by reversal in the thermal gradient in response to a warmer Indian subcontinent relative to the adjacent ocean surface water (Noska and Misra, 2016). In late September/early October warmer ocean surface water than the Indian subcontinent initiates the southward retreat and major decrease in ISM precipitation (Noska and Misra, 2016). Interannual * Corresponding author.
E-mail address: sweldeab@ucsb.edu (S. Weldeab). variability of ISM precipitation is strongly influenced by the El Nino/Southern Oscillation (ENSO) phenomenon, anomalies in the eastern tropical Indian Ocean's north-south SST gradient, and the development of a positive Indian Ocean Dipole (pIOD) Noska and Misra, 2016;Saji et al., 1999;Weller and Cai, 2014;Yadav et al., 2018). The pIOD is a coupled ocean-atmosphere mode of variability that is marked by lower sea surface temperature (SST) than normal in the eastern equatorial Indian Ocean (EEIO) and warmer than normal in the western tropical Indian Ocean (WEIO) coupled with a reversal of zonal wind direction (Saji et al., 1999). The SST gradient between the eastern equatorial Indian Ocean and the Bay of Bengal is not only important for the seasonal climatic transition, but also for changes on year-to-year variability . Anomalously heavy precipitation in the southern parts of the ISM region coupled with severe https://doi.org/10.1016/j.epsl.2021.117327 0012-821X/© 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). droughts and wildfires in western Indonesia are linked to the development of anomalous zonal and meridional Indian Ocean SST gradient changes (Abram et al., 2003;Weller and Cai, 2014). In a warming climate, the frequency of the Indian Ocean's zonal and meridional SST gradient changes is projected to increase DiNezio et al., 2020). Assessing the influence of the current and future thermal gradient changes on ISM precipitation is, however, complicated due to concurrent aerosol increases and differences in the rate with which the Indian subcontinent and the Indian Ocean are warming (Roxy et al., 2015). Paleoclimate records can inform about the long-term relationships between changes in the Indian Ocean's zonal and meridional SST structure and changes in the ISM precipitation.

Setting: imprint of runoff on northern Bay of Bengal surface water
The Ganges-Brahmaputra-Meghna (GBM) river basin covers a substantial part (1.7 million km 2 ) of the Indian Summer Monsoon (ISM) domain, resulting in an annual GBM runoff of 940±127 km 3 /yr ( Figs. 1 and 2). Precipitation over the northern Bay of Bengal is estimated at 2200 mm/yr whose freshening effect is reduced by an evaporation rate of 880 mm/yr (Middlehun et al., 2013;Sengupta and Sarkar, 2006;Singh et al., 2014). While the δ 18 O of ISM precipitation over low-to-middle elevations (< 2,200 m) and the Bay of Bengal varies between −7.4 and −5 , the δ 18 O values of snow and glacier meltwater varies between −21 and −13 (Bowen and Revenaugh, 2003;Breitenbach et al., 2010;Sengupta and Sarkar, 2006;Bowen and Revenaugh, 2003;Hren et al., 2009) (Fig. 2B). The GBM runoff and precipitation over the northern Bay of Bengal significantly lowers the surface salinity and δ 18 O values of the northern Bay of Bengal and leads to a stratification of the upper 11-40 m of the water column (Figs. 1 and 2) (Achyuthan et al., 2013;Durand et al., 2011;Li et al., 2017;Middlehun et al., 2013;Schmidt et al., 1999).
The GBM runoff is enriched in dissolved barium (Ba 2+ ) (Carroll et al., 1993;Moore, 1997;Singh et al., 2013). As a result, and qualitatively similar to several estuaries and bays ( Figure S1), dissolved Ba 2+ in northern Bay of Bengal surface waters is highly enriched and varies between 750 nmol/kg at a salinity of 10 psu (practical salinity) and 200 nmol/kg at a salinity of 30 psu compared to 42.4 nmol/kg at salinity of 35.4 psu in the open tropical Indian Ocean (Carroll et al., 1993;Moore, 1997;Singh et al., 2013). Although incursion of saline, Ba-rich groundwater release during a period of low river discharge is suggested (Moore, 1997), a compilation of global data clearly shows that sea surface water enrichment in dissolved Ba 2+ primarily results from runoff enriched in dissolved Ba 2+ and desorption of Ba 2+ from riverine suspended sediments at very low salinity ( Figure S1). Ba 2+ desorption linearly decreases as salinity increases ( Figure S1). The load of suspended particulate matter in the GBM water is highly correlated (r=0.99) with the flux of fresh water (Samanta and Dalai, 2016). Thus, dissolved Ba 2+ and runoff-induced mixed layer salinity are tightly coupled and directly linked to changes in the GBM runoff and, hence, to changes in ISM precipitation. Taken together, GBM runoff and by extension ISM precipitation, leaves a strong isotopic and trace element imprint on the northern Bay of Bengal surface water. The manifestation of this runoff imprint in the calcite tests of mixed layer-dwelling planktonic foraminifers makes the northern Bay of Bengal an excellent site for a spatially integrated hydroclimate reconstruction.

Materials and methods
Sediment core SO188-KL342 (19 • 58.41'N/90 • 02.03'E, seafloor depth: 1258 m) was recovered from the northern Bay of Bengal where the surface and mixed layer water is strongly influenced by the GBM river runoff (Figs. 1 and 2). Results of a parasound survey over the SO188-KL342 site indicate a gently inclining slope that consists of a relatively homogeneous sediment package (Fig. 3). The sediment sequence of SO188-KL342 consists of fine-grained and relatively dark gray sediment with a lithic fraction varying between 49.5% and 48% ( Figure S2). The absence of turbidites in the investigated section of SO188-KL342 is indicated by the lack of sudden or large changes in the color and lithic fraction, as well as the absence of sudden increase of elements, such as Ti and Fe, in the bulk sediment ( Figure S2). The age model of SO188-KL342 is based on 10 radiocarbon-dates of G. Globigerinoides ruber and Globigerinoides tests (Table S1). Radiocarbon data were analyzed at the Center for Accelerator Mass Spectrometer at Lawrence Livermore National Laboratory and the National Ocean Science Accelerator Mass Spectrometer, Woods Hole Oceanographic Institution. Radiocarbon dates were converted to calendar ages using calibration software CALIB (version 7.10), Marine13 data set, and a constant reservoir age correction of R= 76±120 years for the Bay of Bengal (Stuiver et al., 2019). The final age model is based on linear interpolation between age control points.
For oxygen isotope and trace element analyses, 30-to-35 individual tests of G. ruber sensu stricto (white variety) were picked from the 250-300 μm size fraction. Shells were gently crushed and homogenized. A split of the homogenized samples was cleaned using a foraminifera cleaning procedure that includes oxidative and reductive steps, as described in Martin and Lea (Martin and Lea, 2002). Dissolved samples were analyzed by the isotope dilution/internal standard method described in Martin and Lea (2002) using a Thermo Finnigan Element2 sector field ICP-MS. Analytical reproducibility of Mg/Ca and Ba/Ca measurements, assessed by analyzing consistency standards matched in concentration to dissolved foraminifera solutions and analyzed over the course of entire study, is estimated at ±0.04 mmol/mol and ±0.017 μmol/mol, respectively. Four samples were rejected due anomalously high Ba/Ca (> 3 μmol/mol). Al, Fe, and Mn were also analyzed and used to check cleaning quality of the samples. The range of Al/Ca, Fe/Ca, and Mn/Ca is 3-355, 3-200, 51-210 μmol/mol, respectively. None of these element/Ca ratios are significantly correlated with Mg/Ca and Ba/Ca.
For the down core estimates of runoff-induced mixed layer salinity changes, we analyzed Ba/Ca in G. ruber sensu stricto (G. ruber ss) white variety (test size: 250-300 μm) from northern Bay of Bengal core top samples that were collected from an area with an annually averaged SSS gradient of 3.5 (30.3-33.8 psu). The core top Ba/Ca and mixed layer salinity (S ml ) ( Fig. 3) are linearly correlated with S ml =[(−2.81±0.60)*Ba/Ca G. ruber +36.59±0.41], r 2 = 0.95 and 2σ = ±1.1 salinity unit. Mg/Ca was converted to SST estimate using the following equation (Gray and Evans, 2019;Gray et al., 2018): Mg/Ca (mmol/mol)= exp[0.036*(S-35)+0.064*T-0.87*(pH-8)-0.03] with a residual standard error of 0.65 • C (1σ ). Salinity and pH correction for the SO188-KL342 record is carried out using global salinity and global pH (calculated from pCO2) changes following the protocol provided by Gray and Evans (2019) and Gray et al. (2018). We opt to correct the Mg/Ca using the global change in salinity only, rather than using the reconstructed local salinity estimates, as the reconstructed salinities decrease below the currently calibrated salinity range for Mg/Ca in G. ruber during certain intervals. However, applying the local salinity correction has a negligible impact on the reconstructed temperatures and δ 18 O of the mixed layer and increases the reduction in the N-S temperature gradient observed at ∼7 ka ( Figure S4). The use of the global salinity correction is thus a more conservative approach.
For the reconstruction of the meridional and zonal SST gradient in the tropical Indian Ocean, we focused on and compiled only Mg/Ca data analyzed in G. ruber ss from sediment cores GeoB12605-3 (5.5775 S/39.1136E) (Kuhnert et al., 2014),  (Huffman, 2017). Blue and gray lines indicate tributaries of GBM rivers and selected topographic elevation, respectively. Color shading and contours in the Bay of Bengal shows sea surface salinity (SSS) in September (4-year average) (Fore et al., 2016). Inset (left) depicts monthly SSS changes (Fore et al., 2016) (blue envelope: 2 sigma uncertainty) above SO188-342KL site which is indicated by a red star in the main figure. The black stars in the GBM Basin indicate cave locations of stalagmites discussed below. Lower panel: pIOD conditions in 1997. SST anomaly of November 1997 relative to a 16-year average SST of November (Behringer and Xue, 2004)). Sites of cores whose data are used to reconstruct SST gradients are indicated by purple stars. Rectangles show areas whose SON SST is used to define the pIOD (Saji et al., 1999) and to create the modeled pIOD occurrence. Within the Bay of Bengal, the large and small symbols indicate the location of SO188-342KL and core top samples used for the analysis of Ba/Ca and its relation to mixed layer salinity. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) GeoB12615-4 (7.1416S/39.8458E) (Romahn et al., 2014), and SO189-39KL (0.78S/99.9E) . The cleaning procedure of the samples from SO188-342KL (northern Bay of Bengal, this study) and these of GeoB12605 (Kuhnert et al., 2014) includes a reductive step and therefore the Mg/Ca data were corrected for an estimated 10% Mg loss due to dissolution during the reductive cleaning step (Rosenthal et al., 2004). All Mg/Ca data were converted to SST estimates using the equation and protocol of Gray and Evans (2019), which account for the whole ocean change in salinity and pH. The conversion algorithm requires an input of modern values of pCO 2 (the CO 2 partial pressure difference between sea water and air), salinity, and total alkalinity (A T ) of the water in which G. ruber calcifies (mixed layer: 0-30 m). For SO188-342KL (this study) we used pCO 2 = −10 μatm, A T = 2252 μmol/kg, mixed layer T=28 • C and mixed layer S=31.7. For SO189-39KL  we used pCO 2 = 0.1 μatm, A T = 2245 μmol/kg, mixed layer T=28.96 • C and mixed layer S=33.96. For GeoB12605 (Kuhnert et al., 2014) and GeoB12615-4 (Romahn et al., 2014) we used pCO 2 = 10 μatm, A T = 2298 μmol/kg, mixed layer T=26.5 • C and mixed layer S=35.1. Total alkalinity values were calculated using the salinity and temperature values from the mixed layers of each site and the equation of Lee et al. (2006). The SST gradient between northern Bay of Bengal and eastern equatorial Indian Ocean was established by resampling the highly resolved SST record from eastern equatorial Indian Ocean  using the resolution of SO118-342KL (this study) and subtracting the SST record of Bay of Bengal (this study) from the SST record of eastern equatorial Indian Ocean . Similarly, the SST gradient between the eastern equatorial Indian Ocean and western equatorial Indian Ocean was obtained by resampling the eastern equatorial Indian Ocean SST record  using the resolution of GeoB12605 (Kuhnert et al., 2014) (1.3-9.7 kyr BP) and GeoB12615-4 (9.8-25.8 kyr BP) and subtracting the SST record of western equatorial Indian Ocean from the SST record of eastern equatorial Indian Ocean .
We analyzed the oxygen isotope composition in tests of G.  Precipitation data set shown here is a merged satellite-gauge precipitation estimate with a monthly and 0.1 • x.0.1 • resolution (Huffman, 2017). C) δ 18 O of precipitation (June-to-September) over the GBM Basin (Bowen and Revenaugh, 2003) and monthly averaged δ 18 O of Bay of Bengal surface water (0-1 m) in September (Schmidt et al., 1999). The red star and dots indicate the location of SO188-342KL and core top samples used for the analysis of Ba/Ca and its relation to mixed layer salinity. D) Monthly record of GBM runoff (Durand et al., 2011), sea surface salinity (SSS), and mixed layer salinity (Zweng et al., 2018) over the core site. E) Monthly salinity and temperature changes within the upper 100 m of the water column over the core site. The depth of the mixed layer is estimated at 30 m (0-30 m) (gray areas) and is inferred from the temperature profile.
ice volume corrected δ 18 O of mixed layer seawater (δ 18 O sw ). We use the Mg/Ca-based calcification temperature estimate, the tem- (Bemis et al., 1998), and an estimate of δ 18 O seawater changes related to secular changes in the ice volume (Lambeck et al., 2014). The error estimate of the temperatureand ice volume-corrected δ 18 O sw values is ±0.21 (1σ ). The propagation of uncertainty error is estimated using the follow- 14 , σ I=±0.07 , σ CW=±0.147 , and σ IV=±0.05 is the uncertainty related to the uncertainty of temperature reconstruction (±0.65 • C, 1σ ), uncertainty of analysis of δ 18 O calcite , error associated with the calculation of δ 18 O seawater using calcification temperature and δ 18 O calcite , and reconstruction of δ 18 O seawater changes related to ice volume changes, respectively.
The Kiel Climate Model (KCM) (Park et al., 2009) was employed to perform a transient climate simulation over the last 9500 years. The KCM consists of the atmospheric general circulation model ECHAM5 (Roeckner et al., 2003) with a horizontal resolution of T31 (3.75 • × 3.75 • ) and 19 vertical levels coupled to the ocean-sea ice model NEMO (Madec and Team, 2008) with a horizontal resolution of 2 • × 2 • and enhanced meridional resolution of 0.5 • close to the Equator. The model has been shown to realistically simulate the present-day climate (Park et al., 2009). The KCM has previously been used to analyze transient Holocene simulations with accelerated orbital forcing (Khon et al., 2018). Here, we use nonaccelerated forcing so that the Holocene is represented by 9,500 model years starting from 9.5 kyr BP. Greenhouse gas concentrations are held constant at pre-industrial levels. Changes in the ice-sheet configuration and high latitude meltwater discharge are not varied. The modeled SST record used to calculate the zonal gradients was obtained by averaging the monthly SST of the regions defined to characterize the occurrence of pIOD as shown in Fig. 1.
We employ the MPI (ECHAM5) AGCM (Roeckner et al., 2003) with a module for water isotopes (Werner et al., 2011). The model is run at T106 (about 100 km) resolution and is coupled to a 50 m slab ocean. A cyclostationary climatological heat flux is added to the ocean temperature tendency equation to maintain a seasonal cycle of ocean temperature and sea ice condition that is close to that observed in the present day. Two key experiments are performed to represent the model's response to the middle and early Holocene insolation during 7 kyr BP and 10 kyr BP, respectively, with other boundary conditions derived from modern conditions: 360 ppmv CO 2 and continental geometry, orography, and ice sheets. All of the experiments are run for at least 40 years, with output from the last 30 years used to construct climatologies and climatological differences.

Results
The sediment sequence of SO188-342KL continuously covers the last 25 thousand years (kyr) (Figs. 3A-D). The sediment accumulation rate is variable with the lowest value during early deglacial and glacial (4-5 cm/kyr), highest values (10-20 cm/kyr) between the early-middle Holocene and deglacial, and a declining trend toward the late Holocene (Fig. 3D). The Mg/Ca values vary between 3.6 (glacial) and 4.5 mmol/Ca (middle Holocene). Similarly, Ba/Ca record exhibits lowest values (1 μmol/mol) during the glacial and a highest value of 1.9 μmol/mol during the middle Holocene (Fig. 4A). The range of δ 18 O G. ruber is 3.15 with the highest (−0.61 ) and lowest (−3.76 ) value centered at the early deglacial and middle Holocene, respectively. The Ba/Ca record and Ba/Ca-based salinity estimates indicate that the late deglacial (14.5-11.6 kyr BP) is marked by a monotonic runoff increase. In contrast, the δ 18 O ml reconstruction shows, on average,

Glacial and deglacial hydroclimate
In runoff-influenced marine settings, both planktonic foraminiferal Ba/Ca and δ 18 O are used as proxies for spatially integrated precipitation changes (Gebregiorgis et al., 2016;Saraswat et al., 2013;Weldeab et al., 2007;Weldeab, 2012;Weldeab et al., 2014a,b). In this study, our hydroclimate interpretation is based on the Ba/Ca record for the intervals that exhibit substantial differences in the trend and amplitude between the Ba/Ca and δ 18 O time-series. Based on a number of analyses discussed farther below and above, we conclude that the Ba/Ca record is not biased by changes in the chemical composition of terrigenous sediments, diagenetic imprints, and productivity changes in the Bay of Bengal, respectively ( Figure S2). δ 18 O analysis of rainfall samples collected over the core region of ISM shows that the δ 18 O harbors a strong imprint that is unrelated to the rainfall amount effect (Breitenbach et al., 2010) (Weldeab et al., 2007). The latter is interpreted to contain a significant δ 18 O signature unrelated to the runoff amount (Leduc et al. 2013). The deglacial Ba/Ca-based estimate of runoff-induced S ml changes reveals two outstanding features centered at the middle Holocene and deglacial segments of the record (Fig. 4A). With the exception of a centennial-scale increase in runoff that has been attributed to an increase in Himalayan glacier meltwater input (Weldeab et al., 2019), the Ba/Ca-based salinity reconstruction shows no substantial difference in runoff between the early deglacial (14.6-18 kyr BP) and the Last Glacial Maximum (LGM: 23-19 kyr BP). Furthermore, the Ba/Ca record indicates a monotonic ISM strengthening during the late deglacial (14.5-11.6 kyr BP). Change in the chemical composition of riverine sediment load is not indicated by the results of major element analysis in bulk sediment ( Figure S2). Similarly, an increase of ocean productivity during Heinrich Event 1 and Younger Dryas is not supported by the concentration of chlorins ( Figure S2), degradation products of chlorophyll, that is used as proxy for paleoproductivity (Harris et al., 1996;Schubert et al., 2005). The lowest Ba/Ca value within the early deglacial is 0.9 μmol/mol and is well above the lowest Ba/Ca value (0.68 μmol/mol) recorded during the early deglacial and LGM in the eastern Arabian Sea (Saraswat et al., 2013). This indicates that Ba/Ca values and trend are not affect by the lower limit of Ba incorporation into the test of G. ruber white. Taken together, the impact of atmospheric circulation related to the North Atlantic surface water cooling during Heinrich Event 1 and the Younger Dryas on the basin-wide integrated precipitation imprint is not discernable in the Ba/Ca-based runoff reconstruction, a finding that is consistent with results of a climate simulation (Pausata et al., 2011). Also consistent with our Ba/Ca record, foraminiferal Ba/Ca records from the eastern Arabian Sea and Andaman Sea indicate no distinctive imprint of Heinrich Event 1 ( Figure S6) (Gebregiorgis et al., 2016;Saraswat et al., 2013).
In contrast to the Ba/Ca, the δ 18 O ml record shows an increase δ 18 O ml value between 17.1 and 14.1 kyr BP relative to the LGM values and a pause at 12.8-11.3 kyr BP, interrupting a steady decrease of δ 18 O ml during the late deglacial (Fig. 4B). The timing and duration (including transitions) of these two episodes is synchronous, within age model uncertainties, with the timing of Heinrich Event 1 (17.1-14.6 kyr BP) and the Younger Dryas (12.9.-11.6 kyr BP), prominent intervals of northern high latitude melt-water induced weakening of the Atlantic meridional overturning circulation ( As demonstrated in δ 18 O values of drip water from a cave within the GBM Basin (Breitenbach et al., 2010), we hypothesize that the spatially integrated imprint of δ 18 O precip and consequently the δ 18 O runoff contained a strong component unrelated to changes in the amount of precipitation during Heinrich Event 1 and to lesser degree during the Younger Dryas. As detailed below, we hypothesize that an interplay between the thermal mean state of the eastern tropical Indian Ocean and Northern Hemisphere climate changes resulted in a spatially heterogenous ISM response during the Younger Dryas and Heinrich Event 1 events.

North-South SST gradient reversal in the Eastern tropical Indian Ocean
We create time-series of meridional and zonal SST gradients using Mg/Ca-SST records from the eastern equatorial Indian Ocean , western equatorial Indian Ocean (Kuhnert et al., 2014;Romahn et al., 2014), and northern Bay of Bengal (this study) ( Figures, 1 and 4D-F) (see also Method section and Figures S3). During the deglacial and LGM, the meridional SST gradient is anomalous in exhibiting warmer northern Bay of Bengal surface waters than the eastern equatorial Indian Ocean (Fig. 4E). Our North-South SST gradient change is not accompanied by a sustained reversal of the low resolved, zonal SST gradient (Fig. 4F), indicating the importance of the northern Bay of Bengal's thermal state for the establishment of the north-south SST gradient changes. Similar to previous findings based on two independent SST proxies (Kudrass et al., 2001;Marzin et al., 2013) ( Figure S4), our northern Bay of Bengal SST record exhibits a weak Holocene-LGM temperature contrast. The consistency of three independently reconstructed SST records ( Figure S4) indicates that the weak Holocene-LGM temperature contrast is a robust feature of the northern Bay of Bengal mixed layer. The magnitude and trend of our northern Bay of Bengal SST record is, within method uncertainties, indistinguishable from the planktonic foraminiferal based summer SST reconstructions carried out at an adjacent site (Marzin et al., 2013) (Figure S4). Under modern conditions with surface salinity comparable to our deglacial-glacial salinity reconstruction, 65% of the annual flux of G. ruber in the Bay of Bengal (at 13.14 • N/84.41 • E) occurs during the ISM season (Guptha et al., 1997). Acknowledging that past oceanographic changes can alter seasonal G. ruber flux not resolved in our record, the agreement between two independent foraminiferal SST reconstructions indicates likely a dominated summer and early autumn thermal imprint. The exposure of the Sunda and Sahul shelves during the LGM has been invoked as a crucial player for SST and hydroclimate changes across the tropical Indo-Pacific realm (DiNezio et al., 2018). However, our record indicates that major changes in the eastern Indian Ocean meridional SST gradient occurred during the LGM and the entire deglacial period (Fig. 4E). The latter is marked by rapid sea level rise (Lambeck et al., 2014), suggesting that shelf exposure as a single driver for our observation is unlikely. To summarize, our SST gradient reconstruction reveals a warmer northern Bay of Bengal surface water compared to the Eastern Equatorial Indian Ocean. As detailed below, we propose that the SST gradient reversal is crucial toward understanding the climate dynamics responsible for the deglacial hydroclimate of the ISM region.
The SST gradient change between the northern Bay of Bengal and the Eastern Equatorial Indian Ocean likely played a role in co-shaping the ISM precipitation amount and its δ 18 O signature observed during the Younger Dryas and Heinrich Event 1. We posit that an absence of tropical Indian Ocean SST gradient changes specific to the Younger Dryas and Heinrich Event 1 episodes requires an interplay between the climate dynamics of the northern high latitudes and the thermal mean state of the eastern tropical Indian Ocean. In response to large-scale atmospheric changes associated with the Younger Dryas and Heinrich Event 1events (Pausata et al., 2011), we suggest that the seasonal northernmost extension of the intertropical convergence zone (ITCZ) was limited over the middle-to-low elevation range of the ISM region. Using published stalagmite data (Dutt et al., 2015;Kathayat et al., 2016;Sinha et al., 2005), we calculated the elevation (latitudinal) gradient of δ 18 O during the Younger Dryas and Heinrich Event 1 relative to the Bølling-Allerød (BA) and LGM, respectively (Fig. 5) (see Figure S5

65). E) Mg/Ca-based SST gradient reconstruction between Eastern Equatorial Indian Ocean (EEIO) and northern
Bay of Bengal using Mg/Ca data from the EEIO  and Bay of Bengal as shown in Fig. 3E). F) Mg/Ca-based SST gradient between EEIO and Western Equatorial Indian Ocean (WEIO) using published Mg/Ca data (Kuhnert et al., 2014;Mohtadi et al., 2014;Romahn et al., 2014). G) δD alk-ic : ice volume-corrected hydrogen isotoped analyzed in plant waxes in SO188-KL342 sediments (Contreras-Rosales et al., 2014). I) δ 18 O time-series analyzed in stalagmites from Mawmluh Cave (Dutt et al., 2015) (for location see Fig. 1). Triangles along the x-axis indicate position of age model control points and their 2σ uncertainties. The age range of YD and HE1 is indicated. Gray and yellow shaded areas indicate episodes of meridional and zonal reversal of SST gradient, respectively. for sampling approach). The results reveal a strong δ 18 O increase in the high elevation, the northern part of the ISM domain, during Younger Dryas and Heinrich Event 1 compared to the low elevation southern ISM region (Fig. 5). While we don't rule out the possibility of changes in winter precipitation in the high elevation part of the GBM Basin, we argue that δ 18 O changes in the stalagmites from the south-facing Himalayan foreland (Fig. 2) contain substantial imprints of changes in the amount of summer precipitation.
Hence, the δ 18 O gradient supports (Fig. 5) our hypothesis of latitudinal rainfall zone contraction. In other words, during the Younger Dryas and Heinrich Event 1 monsoonal precipitation in the high elevation ISM region was strongly reduced compared to the low elevation ISM domain.
Noting that the ITCZ follows regions of elevated surface temperatures (Noska and Misra, 2016), we posit that a warmer northern Bay of Bengal surface water than Eastern Equatorial Indian Ocean during the deglacial (Fig. 3F) leads to a delay in the southward retreat of ITCZ toward the Eastern Equatorial Indian Ocean during late boreal summer and early boreal autumn and, consequently, results in seasonally prolonged ISM precipitation over the low elevation part of the ISM region. Both the weakening of precipitation over the high elevation (Fig. 5) and the effect of north-south SST  (Dutt et al., 2015;Kathayat et al., 2016;Sinha et al., 2005) (see Fig. 1 for site locations and Figures S5 for sampling details). Hydroclimate changes within a cave during the Younger Dryas (YD) and Heinrich Event 1 (HE1) are assessed relative to the conditions during the Bølling-Allerød (BA) and Last Glacial Maximum (LGM).
Dots and vertical bars indicate difference of the mean δ 18 O stalagmites values and the variance of the δ 18 O stalagmites data within a given time window. gradient reversal on precipitation over the low elevation ISM region during the Younger Dryas and Heinrich Event 1 provide a scenario that explains two outstanding observations. Firstly, there is no discernable decline in the basin-wide integrated runoff record, as inferred from the Ba/Ca record. Secondly, there is an increase of the δ 18 O values of the spatially integrated runoff due to a reduced contribution from high elevation rainfall that is characterized be low δ 18 O values (Figs. 2B and 2C).

Timing of Holocene ISM strengthening and changes in the zonal SST gradient
Both the Ba/Ca and the δ 18 O time-series indicate that the ISM peak strengthening was centered between at 10 and 5.7 kyr BP and was interrupted by an ISM weakening at 8.4-7.9 kyr BP. Contrasting the Ba/Ca and the δ 18 O records, the latter suggests a comparable magnitude of ISM strengthening between the early Holocene (10-8.5 kyr BP) and middle Holocene (8.4-5.7 kyr BP).
Furthermore, the Ba/Ca and δD alkane records (Figs. 4A, G) indicate a significantly weakened ISM during the late Holocene (5-0 kyr BP) relative to the middle-to-early Holocene. In contrast, the δ 18 O time-series suggest a weak contrast between the late and middleto-early Holocene. The divergence of the late-middle Holocene ISM trend suggested by the δ 18 O record from the Ba/Ca and δD alkane records is consistent with observations in the West African and East Asian monsoon records (Leduc et al., 2013;Wang et al., 2010;Weldeab et al., 2007). These observations, compiled by Leduc et al. (2013), demonstrate that while istopologue-, Ba/Ca-and pollenbased proxies indicate a strong monsoon weakening trend, the δ 18 O records suggest a moderate late monsoon weakening relative to the middle Holocene (Leduc et al., 2013;Wang et al., 2010;Weldeab et al., 2007). As discussed in the first paragraph of the Discussion section, we argue the Holocene δ 18 O record contains a substantial component that is unrelated to the rainfall amount effect. Hence, our interpretation of ISM changes during the Holocene is based on the Ba/Ca record.
The most outstanding feature is that the Ba/Ca record indicates that the Holocene ISM occurred between 7.9±0.47 and 5.7±0.36 (2σ ) kyr BP, significantly later than the timing of Holocene precipitation maximum over the West African, East African, and East Asian monsoon domains (Cheng et al., 2012;Weldeab et al., 2007;Weldeab et al., 2014a). Based on the results of chlorins and major element analyses ( Figure S2), we suggest that the mid Holocene Ba/Ca record is not biased by changes in the chemical composition of the riverine terrigenous sediment input and primary productivity changes in northern Bay of Bengal. Hence, the middle Ba/Ca record most likely indicates increase in runoff and, by extension, precipitation over the GBM Basin. With 31.5±0.2 (n=22), the average mixed layer salinity (S ml ) of the middle Holocene is significantly lower than the early Holocene (32.1±0.1, n=24) and the late Holocene (32. 7±0.1, n=16). The magnitude of this change is substantial given that under modern conditions a 0.5 mixed layer salinity (S ml ) change in the northern Bay of Bengal is associated with a 1-to-2 unit change in sea surface salinity and up to a 55-62 km 3 /month GBM river discharge difference during the monsoon season (Fig. 2D). Though low resolved, a foraminiferal Ba/Ca-based runoff reconstruction from Andaman Sea suggests a Holocene ISM precipitation maximum at 5.7 kyr BP (Gebregiorgis et al., 2016) ( Figure S6), a timing that overlaps with the timing of the peak GBM river runoff and corroborates our finding.
The ISM is sensitive to changes in northern high latitude climate, El Niño activity, and pIOD (Saji et al., 1999). Results of climate simulations also suggest that in a warming climate the ISM will be and was sensitive during the LGM to a pIOD-like mode of interannual variability of the Indian Ocean that is assumed to be triggered by an atmospheric circulation anomaly over the western Indian Ocean and Arabian Sea (DiNezio et al., 2020;Thirumalai et al., 2019). The available paleoclimate records from northern high latitudes do not exhibit any distinct feature coinciding with the timing of peak ISM strengthening (Moossen et al., 2015), though the reconstruction of global mean temperature trend indicates a warming centered at ∼6.5 kyr (Kaufman et al., 2020). Unraveling changes in a regional climate mode variability requires, however, a focus on highly resolved climate and welldated archives and identical methods/proxies of SST reconstruction. Proxy time-series and results of climate model simulations indicate suppressed El Niño activity during the middle Holocene (Pausata et al., 2017). Results of several climate simulations consistently suggest an increase of pIOD occurrence during the middle Holocene (Abram et al., 2007;Iwakiri and Watanabe, 2019;Zhao et al., 2005). Our proxy-based temperature gradient reconstruction reveals a middle-Holocene episode of a warmer western equatorial Indian Ocean than eastern equatorial Indian Ocean, exhibiting pIOD-like conditions (Fig. 4F). Indicating that it is a robust feature, the warming of western equatorial Indian Ocean is paralleled by air warming of the adjacent land mass and surface water cooling of eastern equatorial Indian Ocean which is accompanied by cooling of the westernmost part of the equatorial Pacific Thompson et al., 2002) (see Figure S3 for compilation). Emphasizing that pIOD is a mode of interannual climate variability not resolved by our zonal SST gradient reconstruction, we consider that the SST gradient reversal evident in our reconstruction reflects a mixed signature of strong and frequent pIOD-like occurrences, alternating with conditions of warmer eastern equatorial Indian Ocean surface water compared with the western equatorial Indian Ocean. Thus, our SST gradient reconstruction most likely underestimates the strength of pIOD-like events because of its low resolution relative to the climate anomaly and because of sediment mixing. Nonetheless, the zonal SST gradient reconstruction captures a major change in the middle Holocene thermal structure of the equatorial Indian Ocean. Our interpretation of middle Holocene increase of pIOD occurrence is independently supported by a monthly-resolved coral-based SST reconstruction that exhibits increased middle Holocene occurrences of strong (4.3±0.6 • C below the average) and prolonged (approximately 5 months) cooling of eastern equatorial Indian Ocean surface water compared with the late Holocene (Abram et al., 2007). In summary, the consistency of climate simulations and proxy results strongly indicates that a strong and frequent pIOD occurrence is a robust feature of middle Holocene climate. The timing of the Holocene maximum ISM strengthening (centered between 7.9 and 5.7 kyr BP) is synchronous, within age model uncertainties, with the zonal SST gradient reversal across the Equatorial Indian Ocean. Below we explore causal link between the middle Holocene thermal and hydroclimate changes.
To examine the underlying climate dynamics behind the hydrological and thermal changes observed in the Holocene proxy records (Figs. 4A-F), we analyze the results of the δ 18 O-enabled high resolution ECHAM5 atmospheric general circulation model coupled with a global 50 m slab ocean (see Method). In addition, we use results from a transient simulation of the Kiel Climate Model (KCM) (see Method). The results of the KCM transient simulation indicate that during boreal autumn (Septemberto-November: SON) the middle Holocene was marked by frequent and strong pIOD occurrences (Figs. 6A-C). With 2.4 pIODs per decade, the middle Holocene experienced 33% more pIOD occurrences than the late Holocene (1.6 pIODs per decade). Although the number of the simulated pIOD occurrences per decade was similar for the middle and early Holocene (Figs. 6B-C), there is a difference in that the east-west SST gradient reversal during the early Holocene is weak ( Figure S7). The KCM simulation qualitatively reproduces the increased occurrence of zonal SST gradient reversal indicated in the proxy record. Hence, it provides insights into the seasonality, frequency and strength of the increased middle Holocene pIOD-like occurrences that stand out relative to those of the late and early Holocene.
The high resolution ECHAM5 simulation results suggest that relative to the early (10 kyr BP) and late (0 kyr BP) Holocene the middle Holocene (7 kyr BP) autumn (SON) precipitation over the ISM region was substantially higher and had relatively positive δ 18 O values (Figs. 6D-E). Consistent with the temperature results of the KCM simulation results (Figs. 6A-C), the results of the ECHAM5 simulation indicate a middle Holocene (7 kyr BP) reversal of the equatorial Indian Ocean wind system during SON (Figs. 6D-E). At the same time, strong westerly and northwesterly winds developed over the Arabian Sea and Bay of Bengal (Figs. 6D-E), carrying moisture into the ISM domain. The increase of middle Holocene SON precipitation is thus directly linked to the strong and frequent pIOD-like development and associated atmospheric circulation changes. Relative to boreal summer precipitation (JJA) of the early Holocene, summer precipitation during the middle Holocene was reduced and its δ 18 O signature was significantly higher (Figs. 6F-G). The sum of simulated middle Holocene boreal summer and autumn precipitation over the GBM Basin was slightly higher than that of the early Holocene. The model results, however, suggest a reduction of precipitation from March to May (MAM) during the middle Holocene and, as a result, no significant difference in the annual precipitation between 7 kyr and 10 kyr BP ( Figure S7). Notably, the simulated δ 18 O of middle Holocene precipitation is significantly heavier than the early Holocene (Figs. 6E-G), in line with the reconstructed δ 18 O ml record (Fig. 4B). Overall, though it does not reproduce all proxy-based observations, the result of the climate simulations stands out in highlighting a strengthening of pIOD-like conditions and enhanced SON precipitation during the middle Holocene as the most distinctive changes compared with the early and late Holocene, consistent with the proxy records.
Modern pIOD occurrences are sometimes coincident with the development of El Niño conditions (Saji et al., 1999). However, available middle Holocene proxy data and climate simulation results indicate reduced El Niño activity at this time (Pausata et al., 2017). Hence, the middle Holocene pIOD-like thermal mode could have either evolved independent of El Niño activity or it reflects changes in the threshold sensitivity of the Indian Ocean's response to changes in the thermal state of the tropical Pacific. Some studies have suggested that the increased occurrence of middle Holocene pIODs was caused by the strengthening of the Asian Monsoon (Abram et al., 2007). However, Asian Monsoon strengthening occurred in the early Holocene (8-10 kyr BP) (Cheng et al., 2012) and there is no evidence for a strong and sustained development of pIOD-like conditions during the early Holocence. We point to an increased latitudinal insolation gradient during late summer of the middle Holocene relative to the early Holocene (Fig. 7) and emphasize a positive feedback between insolation gradient, SST, and wind speeds in the Indian Ocean as the primary mechanism for the amplification of middle Holocene pIOD, as also invoked by other simulation studies (Iwakiri and Watanabe, 2019;Zhao et al., 2005). Compared to the early Holocene, it stands out that the insolation gradient between 5 • S-0 • and 20 • N-15 • N during the middle Holocene is positive and stronger starting around mid August (Day 225 in Fig. 7) with the northern Indian Ocean receiving more insolation (Fig. 7). A difference of about two months between the substantial increase of insolation gradients in August and the development of frequent pIOD in SON, as suggest by our simulation results (Figs. 5A-D), is consistent with a delayed response of the mixed layer (Timmermann et al., 2014). Late summer and early autumn warming in the northern Indian Ocean, for instance in the Arabian Sea ( Figure S7), would cause a delay of the ITCZ southward retreat. Consequently, the alongshore winds off Somalia become weaker, resulting in reduced upwelling off Somalia and hence stronger warming of the western Indian Ocean. At the same time, the delay of the ITCZ retreat also allows southeasterly winds to develop in the southeastern equatorial Indian Ocean. The alongshore component of the surface winds off Sumatra likely caused coastal upwelling and surface cooling (Fig. 6). In summary, although we cannot rule out remote oceanic-atmospheric influences, we suggest that the occurrence of enhanced middle Holocene pIOD-like conditions and the associated ISM intensification were driven by changes within the tropical Indian Ocean and triggered by increased middle Holocene insolation gradient (Fig. 7).

Summary and conclusion
The results of this study are distinct from previous hydroclimate reconstructions of the ISM in a number of aspects. The site of the climate archive, the application of a multiproxy approach, and reconstruction of SST gradient changes across the tropical Indian Ocean allowed us to unravel climate changes and their potential association with changes in the thermal states of the Indian Ocean. Our finding has significant implications for the interpretation of oxygen isotope records within the Indian Summer Monsoon region. Furthermore, the results of this study highlight the influence of the Indian Ocean's thermal mean state on the ISM response in conjunction with the deglacial meltwater-induced climate perturbation of the northern high latitude and changes in the middle Holocene insolation forcing. Our record suggests that the Younger Dryas and Heinrich Event1 impact was spatially heterogeneous with a relatively strong ISM precipitation weakening over the high elevation region, while the low elevation ISM region seems to be less affected. We hypothesize that the north-south reversal of SST gradient within the eastern tropical Indian Ocean during the deglacial played a role in co-shaping the hydroclimate gradient within the ISM domain.
Changes in the insolation-driven mean climate state during the middle Holocene produced conditions that led to a reversal of zonal SST gradient across the equatorial Indian Ocean. We hypothesize that the zonal SST gradient change was critical for the timing and duration of Holocene peak ISM intensification by strengthening boreal autumn precipitation, as suggested by the simulation results. Assessing the response of the ISM to an increase of pIOD events over the last 50-60 years is complicated by increasing aerosols and greenhouse gases and their effects on ISM precipitation (Roxy et al., 2015). Our study indicates that changes in the zonal and meridional thermal structure of the tropical Indian Ocean play a crucial role in shaping the hydroclimatic conditions of the GBM Basin.

CRediT authorship contribution statement
Weldeab conceived the research idea, picked the foraminifera, conducted foraminifera cleaning for trace element analysis, performed the stable isotope analysis, wrote the draft of the manuscript. Rühlemann provided access to the SO188-KL342 samples.
Gray assisted with data analysis. Ding performed the δ 18 O-enabled ECHAM simulation. Khon and Schneider performed the KCM simulation. All authors discussed data interpretation and contributed to the final text.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
the Collaborative Research Centre Climate-Biogeochemistry Interactions in the Tropical Ocean (SFB754) and the Excellence Cluster Future Ocean (grant FO EXC 80/2). We thank Dr. Luisa Palamenghi for sharing the parasound image with us. We thank Georges Paradis and Tom Guilderson for the operation of the ICP-MS and radiocarbon datings, respectively.