Zonal Indian Ocean Variability Drives Millennial‐Scale Precipitation Changes in Northern Madagascar

The low latitude Indian Ocean is warming faster than other tropical basins, and its interannual climate variability is projected to become more extreme under future emissions scenarios with substantial impacts on developing Indian Ocean rim countries. Therefore, it has become increasingly important to understand the drivers of regional precipitation in a changing climate. Here we present a new speleothem record from Anjohibe, a cave in northwest (NW) Madagascar well situated to record past changes in the Intertropical Convergence Zone (ITCZ). U‐Th ages date speleothem growth from 27 to 14 ka. δ18O, δ13C, and trace metal proxies reconstruct drier conditions during Heinrich Stadials 1 and 2, and wetter conditions during the Last Glacial Maximum and Bølling–Allerød. This is surprising considering hypotheses arguing for southward (northward) ITCZ shifts during North Atlantic cooling (warming) events, which would be expected to result in wetter (drier) conditions at Anjohibe in the Southern Hemisphere tropics. The reconstructed Indian Ocean zonal (west‐east) sea surface temperature (SST) gradient is in close agreement with hydroclimate proxies in NW Madagascar, with periods of increased precipitation correlating with relatively warmer conditions in the western Indian Ocean and cooler conditions in the eastern Indian Ocean. Such gradients could drive long‐term shifts in the strength of the Walker circulation with widespread effects on hydroclimate across East Africa. These results suggest that during abrupt millennial‐scale climate changes, it is not meridional ITCZ shifts, but the tropical Indian Ocean SST gradient and Walker circulation driving East African hydroclimate variability.

carbon stable isotope ratios as well as trace metal (i.e., Mg, Sr, Ba, U) concentrations (Fairchild & Treble, 2009;Lachniet, 2009;McDermott, 2004).One of the key benefits of speleothems is that they can be precisely and accurately dated using uranium-thorium (U-Th) geochronology, an absolute dating method that takes advantage of U-series disequilibrium in carbonate precipitates.
Speleothem records in opposite hemispheres tend to show an anti-phasing in precipitation intensity over millennial timescales that is strongly correlated with interhemispheric temperature gradients (McGee et al., 2014;Mohtadi et al., 2016;T. Schneider et al., 2014).Particularly notable are the many tropical records showing pronounced changes through abrupt North Atlantic warming (e.g., Dansgaard-Oeschger) and cooling (e.g., Heinrich) events during the last glacial period and the most recent deglaciation (Cheng et al., 2012;Kanner et al., 2012).For instance, through the Dansgaard-Oeschger (D-O) events of MIS 3, speleothem δ 18 O from Hulu Cave in China and Paraíso Cave in Brazil are almost perfectly in anti-phase (Wang et al., 2017).This anti-phase behavior is also remarkably pronounced between sites in tropical Brazil and East Asia during Heinrich stadials (Wang et al., 2004).Such findings have led to the prevalence of the concept of the Global Paleomonsoon (Cheng et al., 2012) in which meridional (i.e., north-south) shifts in the ITCZ are the cause for this precipitation anti-phasing.
In contrast, there is mounting evidence for a zonal (i.e., east-west) mechanism operating in the Indian Ocean over nearly all timescales beyond the annual seasonal cycle (see Section 2.2).The Indian Ocean Dipole (IOD) is known to dominate interannual hydroclimate variability across the basin (Saji et al., 1999).Over multidecadal timescales, there is observational, proxy, and model evidence for low frequency zonal variability in the tropical Indian Ocean (Tierney et al., 2013;Ummenhofer et al., 2017).While meridional ITCZ shifts would result in anti-phased behavior between hemispheres with more rainfall in the warmer hemisphere, a zonal mechanism instead predicts hemispherically in-phase hydroclimate variability linked to zonal tropical sea surface temperature (SST) gradients.Speleothem reconstructions of the Common Era from Oman and Madagascar (in opposite hemispheres of the western Indian Ocean) show coherent centennial-scale hydroclimate variability which suggests the possibility of a zonal mechanism over these timescales (Scroxton et al., 2017).On orbital timescales, there is evidence Northern Hemisphere summer insolation drives changes in the reconstructed W-E tropical Indian Ocean gradient with notable effects on East African rainfall (Johnson et al., 2016).
Although a zonal mechanism is prominent in driving tropical Indian Ocean variability over interannual, multidecadal, and perhaps even centennial and orbital timescales, it remains unclear whether meridional ITCZ shifts or zonal mechanisms dominate over millennial timescales.Previous studies have disagreed on the drivers of millennial-scale hydroclimate variability in East Africa, with some seeing regional changes as reflecting global meridional ITCZ shifts (Johnson et al., 2002;Mohtadi et al., 2014), while others have emphasized more local control by Indian Ocean SSTs (Niedermeyer et al., 2014;Stager et al., 2011;Tierney et al., 2008;Weldeab et al., 2014) or a combination of both (Scroxton et al., 2019).Others posit that the ITCZ expands and contracts rather than shifts latitudinally (Collins et al., 2011).Additionally, the location of the Congo Air Boundary can affect the amount of moisture transport from the Indian Ocean to East Africa's monsoon regions (Tierney et al., 2011).
Northern Madagascar sits under the climatological ITCZ position in December, January, and February (DJF) and is affected by interannual Indian Ocean zonal variability, making it an ideal study site for exploring the nature of millennial-scale hydroclimate variability in East Africa.Like many areas of the Global South, Madagascar is relatively under-sampled and its climate history largely unresolved through the last glacial period and deglaciation.Exceptions include a speleothem record from Tsimanampesotse in Madagascar's arid southwest (Scroxton et al., 2019), diatom records from Lake Tritrivakely (Gasse & Van Campo, 1998) and Lake Maudit (Teixeira et al., 2021), and a runoff record northwest of the island (Ma et al., 2021).There are also radiocarbon-dated pre-Holocene subfossils, some of which have been analyzed for δ 13 C and δ 15 N to reconstruct diet and habitat moisture (Crowley, 2010;Crowley et al., 2021;Godfrey et al., 2016;Hixon et al., 2021;Samonds et al., 2019), though they are also spatially limited.
Here, we reconstruct glacial and deglacial tropical Indian Ocean hydroclimate using AB12, a new speleothem record from Anjohibe ("big cave" in Malagasy) in northwest Madagascar, with the goal of understanding millennial-scale drivers of tropical Indian Ocean hydroclimate variability.We begin with an overview of Madagascar's modern hydroclimate (Section 2) and a description of speleothem proxy systems (Section 3), followed by an overview of the methods used (Section 4), age model and proxy results (Section 5), and a discussion comparing the new record with others from the region as well as further discussion regarding the underlying dynamics driving millennial-scale hydroclimate variability (Section 6).

Seasonal Variability
Anjohibe is located in northwest Madagascar (15.54°S, 46.89°E, 131 masl), a region featuring a tropical savanna climate (Aw Köppen classification) with highly seasonal precipitation.From observations in Mahajanga, a city 73 km southwest of Anjohibe, the wet season peaks during austral summer and makes up 85% of the mean annual total (Scroxton et al., 2017;Voarintsoa et al., 2021).This region is dominated by Madagascar's northwestern monsoon which is largely driven by northwesterly winds advecting moisture to the island from the tropical western Indian Ocean due to the sensible heating of Madagascar's land surface (Jury, 2016).The ITCZ is located over northern Madagascar during austral summer (Jury et al., 1995) as indicated by outgoing longwave radiation (OLR) minima (Burns et al., 2022), making Anjohibe an ideal study site for reconstructing past changes in the local monsoon and ITCZ over the southwest Indian Ocean (Figure 1).In austral winter, dry conditions prevail at Anjohibe due to the rain shadow cast by Madagascar's central mountains blocking prevailing easterlies and causing orographic rainfall along the island's eastern coast (Scroxton et al., 2017).

Interannual Variability
Although meridional ITCZ shifts dominate the seasonal cycle of precipitation at Anjohibe, there are other modes of variability that can affect local hydroclimate at interannual timescales.The Indian Ocean Dipole (IOD) is a mode of coupled atmosphere-ocean variability corresponding to changes in the zonal tropical Indian Ocean SST gradient and surface winds with notable effects on tropical hydroclimate across the basin (Saji et al., 1999).The IOD peaks in austral spring during September, October, and November (SON) when southeasterly winds in the eastern Indian Ocean drive upwelling along the Javan and Sumatran coasts (Abram, Hargreaves, et al., 2020;Abram, Wright, et al., 2020).Much like the El Niño Southern Oscillation (ENSO), eastern upwelling is capable of spurring ocean-atmosphere coupling via the Bjerknes feedback which modulates the strength of the basinwide Walker circulation (X.Wang & Wang, 2014).This also causes deviations from the climatological SST gradient, which tends to be warmer in the east within the Indo-Pacific Warm Pool and cooler in the west near the Seychelles-Chagos Thermocline Ridge.A positive IOD (pIOD) event is associated with a weakening of the climatological zonal SST gradient due to enhanced eastern upwelling which results in warmer SSTs, enhanced moisture convergence, and increased precipitation across the western Indian Ocean and East Africa (Ummenhofer et al., 2009).The opposite is true for negative IOD (nIOD) events, which are characterized by a strengthening of  (Kummerow et al., 1998).Yellow contour represents the 4 mm/day isohyet which has been used in previous work to designate the boundary of the Indian Ocean tropical rain belt (Li et al., 2020).The red circle indicates the location of Anjohibe (this study), and green circles indicate the location of hydroclimate records discussed in this paper.(1): Moomi Cave (Shakun et al., 2007).( 2): Gulf of Aden (Tierney & DeMenocal, 2013).( 3): Lake Rutundu (Garelick et al., 2021).( 4): Lake Tanganyika (Tierney et al., 2008).( 5): Lake Maudit (Teixeira et al., 2021).( 6): Anjohibe (this study).( 7): Mitoho Cave (Scroxton et al., 2019).
the climatological zonal SST gradient and Indian Ocean Walker circulation.While previous work has explored the influence of ENSO on Madagascar's arid southwest (Randriamahefasoa & Reason, 2017), the effects of ENSO and IOD variability on precipitation over northern Madagascar have not yet been explored in much detail.
To explore the relationship between precipitation at Anjohibe and modes of tropical interannual climate variability, gridded precipitation data from the Global Precipitation Climatology Centre (GPCC, 0.25°) were correlated with the Optimum Interpolation Sea Surface Temperature (OISST, 0.25°) data set.All data were pre-processed by removing seasonality and long-term trends.The GPCC cell closest to Anjohibe as well as the eight adjacent cells were averaged using a queens contiguity distance-weighted method (see Diem et al., 2019).This method results in data primarily reflecting rainfall in the grid cell containing the Mahajanga rain gauge station, the city closest to Anjohibe, with some influence from other rain gauge stations in the region.The resulting Anjohibe precipitation time series was then correlated with Indian Ocean SST seasonal means.For determining the significance of these correlations as well as the rest shown in this paper, a nonparametric isospectral method implemented in Pyleoclim was used to account for reduced degrees of freedom in autocorrelated data (Ebisuzaki, 1997;Khider et al., 2022).
Correlations between Indian Ocean SSTs in SON and Anjohibe precipitation in SON reveal a statistically significant relationship between the zonal SST variability characteristic of Indian Ocean Dipole events and precipitation amount at Anjohibe at the onset of the austral summer rainy season (Figure 2a).Increased (decreased) precipitation at Anjohibe is highly correlated with warmer (cooler) SSTs in the western Indian Ocean and cooler (warmer) SSTs in the eastern Indian Ocean, a pattern reminiscent of pIOD (nIOD) events.Remarkably, Anjohibe precipitation shows a greater correlation with the eastern region (ρ = −0.431,p = 0.011) and DMI (ρ = 0.381, p = 0.001) than with the more proximal western region (ρ = 0.264, p = 0.055).Altogether, these findings indicate that the Indian Ocean zonal SST gradient and the eastern upwelling zone are more important controls on interannual hydroclimate variability at our study site rather than local SSTs alone.
While there is a relationship between site-specific precipitation and the W-E SST gradient in SON, it is important to note that shoulder season precipitation is a small percentage of the annual total.However, the impact of the seasonally locked IOD may extend into the subsequent DJF rainy season at Anjohibe, more likely by sea surface temperature changes rather than atmospheric circulation changes.A seasonal composite analysis shows that years with positive IOD events tend to have greater SON and DJF precipitation than years with negative IOD events (Figure S4 in Supporting Information S1).It is also worth noting that only some of the SST grid cells of highest correlation in the western Indian Ocean are located within the classic western IOD region, with some of the highest correlated cells further south.This suggests more proximal sea surface temperatures play a role in precipitation variability at Anjohibe in addition to basin-wide dynamics (Figure 2a).
There is strong evidence for coupled IOD and ENSO variability throughout much of the last millennium (Abram, Wright, et al., 2020).ENSO affects the tropics and subtropics differently in the western Indian Ocean region, with tropical regions receiving more rain and subtropical regions receiving less rain during El Niño years (and the opposite for La Niña years) (Zhang et al., 2015).Southwest Madagascan hydroclimate is impacted by ENSO and mirrors the subtropical south African response where El Niño years are drier and La Niña years are wetter (Cook, 2000;Randriamahefasoa & Reason, 2017).However, there is no significant correlation between Anjohibe precipitation and the Niño 3.4 index in SON (ρ = 0.002, p = 0.607) or DJF (ρ = 0.019, p = 0.923), suggesting that internal Indian Ocean variability related to the IOD is more important for Anjohibe precipitation than the pantropical inter-basin interactions characteristic of ENSO variability.This is in agreement with model evidence for multidecadal internal Indian Ocean variability controlling East African rainfall (Tierney et al., 2013).
Although the underlying dynamics require further investigation, the IOD plays a role in driving NW Madagascar's interannual hydroclimate variability.With the southern migration of the ITCZ bringing the annual rainy season, and zonal Indian Ocean SSTs gradients affecting interannual precipitation variability, it remains unclear whether meridional or zonal processes dominate Anjohibe's hydroclimate variability over longer timescales.

Oxygen Isotopes
Measuring and modeling isotopes in modern precipitation is important for understanding the drivers of δ 18 O variability in proxy reconstructions.A recent study using the Experimental Climate Prediction Center's Isotope-incorporated Global Spectral Model (IsoGSM) found a strong inverse correlation between DJF δ 18 O and precipitation amount over northwest Madagascar and the surrounding region (Duan et al., 2021).This phenomenon is known as the amount effect, which is the preferential depletion of heavier precipitation isotopes with increasing rainout.The amount effect in interannual DJF δ 18 O is also observed in Antananarivo as well as Rodrigues and Mauritius based on the Global Network of Isotopes in Precipitation (Li et al., 2020), suggesting that this process is dominant across the tropical southwest Indian Ocean.While source changes can influence δ 18 O in other tropical regions, such as the Maritime Continent (Atwood et al., 2021), NW Madagascar's summer rainfall is sourced almost exclusively from the equatorial western Indian Ocean (Scroxton et al., 2017).For these reasons, we interpret Anjohibe δ 18 O records to reflect precipitation amount and the strength of the regional summer monsoon.Recent cave monitoring work further supports the viability of using speleothem δ 18 O at Anjohibe for monsoon reconstructions (Voarintsoa et al., 2021), though additional measurements of precipitation isotopes through multiple annual cycles would help support this interpretation.Anjohibe correlated with SON-mean Indian Ocean SST anomalies from 1982 to 2020.Anjohibe's location is shown as a green dot.Red indicates a positive correlation, and blue a negative correlation.The Spearman correlation statistic (ρ) was used to account for precipitation's non-normal distribution.Areas within the black contours are significant at the 95% confidence level based on a nonparametric isospectral method accounting for autocorrelation (Ebisuzaki, 1997;Khider et al., 2022).The black dashed boxes represent the regions used to calculate the Dipole Mode Index (DMI) (western: 50-70°E, 10°N-10°S; eastern: 90-110°E, 0-10°S) (Saji et al., 1999).(b)-(d): Scatterplots of Anjohibe precipitation in SON versus western Indian Ocean SSTs in SON (b), eastern Indian Ocean SSTs in SON (c), the Dipole Mode Index (DMI), the difference between the western and eastern Indian Ocean regions in SON (d).The solid black lines represent the lines of best fit, and gray shading represents the 95% confidence intervals based on 1,000 bootstrapped resamples.Note: the IOD is known to peak in SON.SST data are from OISST, and site-specific precipitation data are from GPCC.

Carbon Isotopes
The carbon isotopic composition of speleothem calcite depends on the initial δ 13 C of multiple carbon sources (McDermott, 2004) as well as the degree of CO 2 degassing in the epikarst (Johnson et al., 2006).Carbon in cave waters reflects contributions from atmospheric CO 2 , respired CO 2 in soils, and CO 2 from bedrock dissolution (McDermott, 2004).As detailed in a previous paper on a speleothem reconstruction from Anjohibe, plant root respired CO 2 is likely the greatest contributor to δ 13 C variability at this site (Burns et al., 2016).Depending on the photosynthetic pathway used, different plant types will generate distinct δ 13 C signals recorded in speleothems.C 4 plant δ 13 C values are about 8-10 permil more enriched compared to C 3 plant δ 13 C (McDermott, 2004).Speleothem δ 13 C can also be influenced by the strong Rayleigh fractionation effects of CO 2 degassing from solution before reaching a drip site (Fairchild et al., 2006).The lighter 12 C preferentially degasses, leading to enrichment of 13 C and an increase in δ 13 C of the host solution.Arid conditions reduce the amount of water flowing through the vadose zone which aerates the aquifer and allows greater degassing, altogether increasing δ 13 C (Fairchild et al., 2006;Fairchild & Treble, 2009).
Though plant coverage as well as degassing can drive changes in δ 13 C their combined effect is often confounded by the fact that they tend to drive δ 13 C in the same direction under similar hydroclimatic conditions.For instance, drier conditions might result in reduced microbial and root respiration in the soil and enhanced degassing, both of which increase δ 13 C.For this reason, a multiproxy approach using speleothem trace elements alongside δ 13 C is often beneficial in confirming paleoclimatic and paleoenvironmental proxy interpretations.

Trace Elements
Trace elements in speleothems such as Mg, Sr, and U can act as proxies for local water balance because they are preferentially excluded from the calcite lattice during precipitation.As a result, increased prior calcite precipitation (PCP) results in higher X/Ca in speleothem source waters, where "X" here refers to Mg, Sr, or other elements with distribution coefficients relative to calcium that are less than 1 (Day & Henderson, 2013).With PCP increasing under more arid conditions, speleothem X/Ca ratios can act as a proxy for local paleo-aridity.While Mg/Ca and Sr/Ca tend to be in good agreement due to the effects of PCP, U/Ca ratios often show trends opposite those seen in Mg/Ca and Sr/Ca (Asmerom et al., 2020).This is attributed to processes besides PCP, such as the formation of mobile uranyl-phosphate complexes (Treble et al., 2003), U scavenging onto mineral surfaces during PCP (Johnson et al., 2006), and prior aragonite precipitation (PAP) (Jamieson et al., 2016).It is worth noting that calcite and aragonite have distinct trace element distribution coefficients due to differences in cation site sizes (Finch et al., 2001).
Recent cave monitoring work at Anjohibe spanning almost 2 years reveals that with decreasing drip water Ca, Mg/Ca and Sr/Ca increases, a trend indicative of PCP (Voarintsoa et al., 2021).This study also found that drip water Mg/Ca and Sr/Ca decrease after the peak of the DJF rainy season with up to a one season delay due to the "epikarst storage effect" (Voarintsoa et al., 2021).Additional monitoring work analyzing drip water U/Ca could shed light on the drivers of this less understood proxy system.These preliminary findings demonstrate the viability of using trace elements proxies at Anjohibe for climate reconstruction.With other non-climatic processes also having an effect on trace element to calcium ratios, it is important to check whether these trace element proxies agree with one another as well as with stable isotope proxies to verify whether they are recording climatic and environmental signals.

Sample Collection
AB12 is a 1.7m-long calcite stalagmite collected from the main chamber of Anjohibe in 2019.It was broken into four pieces to assist in transport from the study site.These pieces were each cut lengthwise and polished prior to sub-sampling for U-Th dating, stable isotope analysis, and trace element analysis.

U-Th Geochronology
Thirty samples weighing ∼200 mg were drilled from AB12 with a vertical mill.Once the powders were dissolved and spiked with a 229 Th-233 U-236 U tracer, U and Th were isolated by Fe-coprecipitation and purified with ion-exchange chromatography using columns containing AG1-X8 resin following (Edwards et al., 1987).The subsequent fractions were analyzed alongside a total procedural blank using a Nu Plasma II-ES MC-ICP-MS at MIT with the methods detailed in a prior study on a speleothem from Anjohibe (Burns et al., 2016).An initial 230 Th/ 232 Th ratio of 17 ± 8.5 × 10 −6 was used based on a prior isochron analysis (Burns et al., 2016).The 230 Th and 234 U half-lives from Cheng et al. (2013) and the 238 U half-life from Jaffey et al. (1971) were used in age calculations.COPRA (Breitenbach et al., 2012) was used to build an age model and develop proxy-age relationships by using a Monte Carlo approach.

Stable Isotope Analysis
The stable oxygen and carbon isotope ratio measurements were performed at the University of Massachusetts.Subsamples were microdrilled from a cut and polished slab at a 2 mm resolution along the central growth axis of the stalagmite (n = 866).The samples were analyzed in an on-line carbonate preparation system, GasBench II, linked to a Thermo Scientific Delta V ratio mass spectrometer.Results are reported as the permil (‰) difference between sample and the Vienna Pee Dee Belemnite (VPDB) standard in delta notation where δ 18 O = (R sample / R standard − 1) × 1,000, and R is the ratio of the minor to the major isotope.Reproducibility of standard materials is 0.08‰ for δ 18 O and 0.05‰ for δ 13 C.
For the comparison between AB12 δ 18 O and the δ 18 O from Anjohibe Common Era records discussed in Section 6.1, calcite δ 18 O was converted to precipitation δ 18 O based on the known temperature dependence of oxygen isotope fractionation during calcite formation (Tremaine et al., 2011).Six records from Anjohibe spanning parts of the Common Era were used: AB2, AB3, MA2, MA3, ANJB-2, and ANJ94-5 (Burns et al., 2016;Scroxton et al., 2017;Voarintsoa, Railsback, et al., 2017;Voarintsoa, Wang, et al., 2017;Wang et al., 2019).Monitoring work at Anjohibe found a modern mean cave temperature of 25.0°C (Voarintsoa et al., 2021).Using the Last Glacial Maximum Reanalysis (LGMR), the LGM was 3.7°C cooler than the Common Era which resulted in a mean local temperature of 21.3°C (Osman et al., 2021).After using these temperatures to calculate precipitation δ 18 O and converting to VSMOW, the glacial record was corrected by 1.1 permil to account for the effect of ice volume on seawater δ 18 O (Schrag et al., 2002).The same methods to correct speleothem calcite δ 18 O for ice volume and temperature effects and convert to precipitation δ 18 O have been used in previous work (Comas-Bru et al., 2019;Osman et al., 2021;Tierney et al., 2020).

Trace Element Analysis
The powders remaining after stable isotope analysis were down-sampled to an 8 mm resolution (i.e., one out of every four of the stable isotope samples) and dissolved in 2% nitric acid.These samples were analyzed on an Agilent 7900 ICP-MS at MIT for Mg, Sr, U, and Ca and calibrated using a multi-element standard with similar elemental ratios at variable concentrations.Samples with excess powder were run in replicate to ensure measurement reproducibility.Standards were run every 10 samples to check signal stability through the course of each run.Samples were identified as outliers using Tukey's outlier criterion (i.e., falling above the third quartile 1.5X the interquartile range or below the first quartile −1.5X the interquartile range).Outliers might be attributed to a lack of instrumental signal and/or sampling of small zones of aragonite within the sample (Figure S1 in Supporting Information S1).All replicates including outliers (n = 72) show a 7.76%, 19.75%, and 6.89% reproducibility for Mg/Ca, Sr/Ca, and U/Ca.Non-outlier replicates (n = 54) better reflect the data used later on in this paper and show an improved reproducibility of 1.67%, 1.55%, and 3.53% for Mg/Ca, Sr/ Ca, and U/Ca.

Event Timing
Millennial-scale event timing was based on the INTIMATE event stratigraphy for Heinrich Stadial 1 (HS-1, 17.5-14.7 ka), the Bølling Oscillation (14.7-14.1 ka), the Allerød Oscillation (14.1-12.9ka) (together, the B-A), and the Younger Dryas (YD, 12.9-11.7 ka) (Rasmussen et al., 2014).Heinrich Stadial 2 (HS-2) timing is not as well constrained as HS-1, though there is evidence from Greenland dust records it occurred in two distinct pulses 26-25.1 and 24.1-23.4ka (Rasmussen et al., 2014).This is further supported by a compilation of North Atlantic SST records (Heath et al., 2018).HS-2 was also divided into two events in a recent paper on sedimentary runoff records from the Mozambique Channel northwest of Madagascar (Ma et al., 2021).

Age Model
AB12 appears to have grown continually and with minimal changes in growth rate over 12,000 years of the last glacial period and deglaciation, including HS-2, the LGM, HS-1, and the onset of the B-A (Figure 3).The oldest and youngest U-Th age determinations for sample AB12 are 26,486 ± 133 years before present (yr BP) and 14,260 ± 75 yr BP, respectively.With 2σ errors in U-Th ages averaging 125 years through the record, AB12 provides tight age constraints on millennial-scale climate variability during the last glacial period and deglaciation.There are no obvious hiatuses across growth axis changes or other parts of the sample.While most of AB12's crystal fabric is continuous and dense, some parts were rather porous.High porosity can indicate dissolution and recrystallization which can result in spurious age determinations (Bajo et al., 2016).For this reason, these areas were avoided in the drilling process where possible.Three of the 30 dates were not included because they were too old for their stratigraphic position, and one was excluded because it was too young (Figure 3).The old outliers were drilled from areas with high porosity, which suggests that they could have been affected by diagenetic U-loss.U-loss due to post-depositional dissolution and open-system behavior can result in anomalously old U-Th determinations (Bajo et al., 2016).The anomalously old outliers were resampled from less-porous layers and reanalyzed which resulted in better stratigraphic agreement (Figure 3).

Stable Isotope Proxies
Based on the COPRA-derived age model, the temporal resolution of the δ 18 O and δ 13 C records averages 14 years with a range of 5-30 years.Speleothem δ 18 O varies between −6.77 and −2.95‰ and averages −4.67‰, with the highest sustained values between 15 and 18 ka (Figure 4).Speleothem δ 13 C varies between −10.65 and −5.13‰ and averages −9.07‰.δ 18 O and δ 13 C are strongly correlated (r = 0.55, p < 0.001).Though the δ 18 O and δ 13 C records have many features in common, δ 13 C appears to have more pronounced variability than δ 18 O over centennial and millennial time scales.(Breitenbach et al., 2012).A photograph of AB12 is shown to the right of the plot.

Trace Metal Proxies
Trace metals were sampled at lower resolution than stable isotope measurements, resulting in a ∼60 years average resolution.Some of the extreme Sr/Ca and U/Ca values were excluded as outliers based on thresholds determined by Tukey's criterion (see Materials and Methods) (Figure S1 in Supporting Information S1).The Sr/Ca and U/ Ca data in the Figure S2 in Supporting Information S1 scatterplots and Figure 5 time series have these outliers removed.Most of these outliers in Sr/Ca and U/Ca are located near the top of AB12.Though the cause of these heightened values remains unknown, it is possible that mm-scale shifts to aragonite increased Sr/Ca and U/Ca because of the high compatibility of Sr and U in aragonite compared to calcite.
Mg/Ca versus Sr/Ca have a positive correlation, while U/Ca versus Mg/Ca and U/Ca versus Sr/Ca both have negative correlations.Mg/Ca and Sr/Ca show a statistically significant correlation at the p = 0.05 level after outlier removal (Figure S2 in Supporting Information S1).Although the U/Ca versus Mg/Ca relationship is significant, the U/Ca versus Sr/Ca relationship is not.Correlations with oxygen and carbon stable isotope proxies are also all significant beyond the p = 0.05 level.δ 13 C and δ 18 O are both positively correlated with Mg/Ca and Sr/Ca, and negatively correlated with U/Ca.There is a particularly strong relationship between δ 13 C and Mg/Ca.The δ 13 C-trace metal correlations are overall stronger than those with δ 18 O.
The slope of the line of best fit between ln(Sr/Ca) and ln(Mg/Ca) is diagnostic of PCP activity based on a model of PCP's influence on drip water and speleothem trace metal concentrations (Wassenburg et al., 2020).According to this model, ln(Sr/Ca) versus ln(Mg/Ca) slopes between 0.709 and 1.45 are indicative of PCP.AB12 ln(Sr/ Ca) versus ln(Mg/Ca) has a slope of 0.575, which is slightly shallower than what is predicted for PCP.This Figure 5. AB12 stable isotope and trace metal proxy time series comparison.From top to bottom: δ 18 O, δ 13 C, Mg/Ca, Sr/Ca, U/Ca.Trace metals are presented with outliers removed based on Tukey's method.Trace metal data are reported in molar ratio units.Note that the y-axes for all variables except U/Ca are flipped so that inferred dry conditions are oriented down.Black lines show the 12-point and 3-point running means for stable isotope and trace metal records, respectively.Different running means are applied to these records because they have different temporal resolution.Gray lines show unaveraged data.
suggests the influence of calcite recrystallization in combination with PCP in driving changes to Mg/Ca and Sr/ Ca (Wassenburg et al., 2020).
Strongly positive correlations between Mg/Ca and Sr/Ca with δ 13 C would suggest the influence of PCP-related degassing on δ 13 C, as detailed in previous work (Johnson et al., 2006).While it is clear that δ 13 C and trace metal ratios are in agreement with respect to the large changes reconstructed, approximately half of the variance in δ 13 C remains unexplained by PCP-sensitive trace element ratios.It is probable that δ 13 C is also driven by changes in soil respiration rates.
The U/Ca record is most likely driven by the effects of U scavenging onto mineral surfaces during PCP rather than PAP (Johnson et al., 2006).If PAP were occurring at Anjohibe, we would expect to find less of a relationship between Sr/Ca and U/Ca because Sr has a distribution coefficient of ∼1 in aragonite and would therefore be relatively unaffected by increasing aragonite precipitation (Fairchild & Treble, 2009).
The AB12 proxy time series show a general agreement in millennial-scale variability as well as some agreement in centennial-scale variability (Figure 5).One notable feature is the coherent, punctuated shift to higher values 16-17 ka in δ 13 C and trace metal ratios which coincides with a maximum in δ 18 O values.The δ 18 O record overall has more muted variability in comparison to δ 13 C and the trace metal records.It is worth noting that the abundance of trace metal outliers around 15 ka results in a lack of data there which makes this part of the record difficult to interpret from trace metals alone.

The Last Glacial Period Compared to the Common Era
After correcting for temperature and ice volume effects and converting to precipitation δ 18 O, δ 18 O from the Last Glacial Period (27-19 ka) is more depleted than δ 18 O from Common Era (2 ka-present) records at Anjohibe (Figure S5 in Supporting Information S1).The difference in the medians of these populations is 1.43‰, and they are significantly distinct (p < 0.001 for both a t-test and Kolmogorov-Smirnov test).The depleted δ 18 O during the Last Glacial Period might signal wetter conditions and a stronger monsoon at Anjohibe compared to present-day conditions.Simulations under LGM boundary conditions support this claim, showing greater western Indian Ocean precipitation and reduced eastern Indian Ocean precipitation compared to preindustrial runs (DiNezio et al., 2018;Thirumalai et al., 2019).This pattern is likely due to sea level effects associated with Sunda and Sahul Shelf exposure in the Maritime Continent region (DiNezio et al., 2018).Future work examining influences on precipitation isotopes under LGM boundary conditions will be required to better understand the relationship between depleted precipitation δ 18 O values and NW Madagascar precipitation during the Last Glacial Period.

Interpreting Millennial-Scale Variability
AB12 δ 18 O shows relatively unchanging conditions from the start of the record to ∼18 ka, suggesting a stable monsoon regime with minor fluctuations during much of the period 27-18 ka (Figure 5).A slight δ 18 O increase ∼25-26 ka suggests moderately dry conditions through the first HS-2 pulse, though these conditions don't exceed the range of glacial hydroclimate variability.AB12 δ 18 O decreases at the peak of the LGM ∼21 ka, signaling moderately wet conditions.At the start of HS-1 ∼18 ka, δ 18 O values begin increasing until they reach a maximum ∼16 ka, indicating a shift to drier conditions.After this interval of maximum δ 18 O and aridity during the stadial, δ 18 O begins to decrease slowly until the start of the B-A ∼14.7 ka when δ 18 O sharply decreases, marking a punctuated shift to wetter conditions consistent with the onset of wetter conditions throughout eastern and northern tropical Africa (Otto-Bliesner et al., 2014).Finally, δ 18 O sharply increases and growth stops at the end of the Bølling Oscillation.Overall, the δ 18 O hydroclimate record shows stable conditions during the latter part of the last glacial period with the possible exception of a drier HS-2, drier conditions over HS-1, and an abrupt shift to wetter conditions at the onset of the B-A.δ 13 C and trace elements similarly feature a shift to dry conditions 25.5ka, coinciding with the first pulse of HS-2 (Figure 5).A wet LGM at ∼21ka and a dry HS-1 starting ∼18 ka are also supported by δ 13 C and trace elements, with the wettest parts of these records often occurring during the LGM and the driest parts of these records occurring during HS-1 (Figure 5).While δ 18 O appears to respond to the broader Greenland stadial, δ 13 C and trace metals indicate a shorter-lived dry period at 16.0-16.5ka.This shift to dry conditions is in good agreement with Heinrich Event 1.1 occurring 15.5-17.1 ka spurred by massive discharge of the Laurentide Ice Sheet as evidenced by a substantial increase in North Atlantic ice-rafted debris (IRD) (Hodell et al., 2017).Heinrich Event 1.1 was found to correspond to drier conditions recorded in speleothems from southwest Madagascar and the Southeast Asian Monsoon domain, and wetter conditions in the South American Monsoon domain (Scroxton et al., 2019;Stríkis et al., 2015;Wu et al., 2009).This timing is also consistent with the timing of Heinrich-related δ 18 O changes in the Hulu record (Wang et al., 2001).Though both Sr/Ca and U/Ca show an initial shift to drier conditions coeval with the onset of this event, the Sr/Ca and U/Ca records are lacking in non-outlier data points that might signal a recovery from Heinrich Event 1.1.
The δ 18 O, δ 13 C, and U/Ca records show an abrupt shift to wetter conditions at the onset of the B-A.However, Mg/ Ca does not appear to support this interpretation.Sr/Ca might also signal disagreement with a wetter B-A, but this is difficult to interpret given the record's lack of data points from the abundance of outliers through this interval.In summary, most of the records presented here agree on a wetter B-A, though some of the trace element data through this interval are an unresolved exception to the pattern suggested by δ 18 O, δ 13 C, and U/Ca.

Comparison With Other Sites From the Region
When comparing our record with others, we use δ 18 O because it is a tracer of regional monsoon strength.δ 13 C and trace metals record local water availability which is correlated with regional monsoon strength, but with greater variability than δ 18 O.At present there are not enough regional records with the centennial-scale resolution necessary to compare to this variability, future work developing new records and monitoring drip water chemistry will be required.
The deglacial millennial-scale hydroclimatic and environmental changes reconstructed in other records from Madagascar and the western Indian Ocean region are in good agreement with those reconstructed in AB12 (Figure 6).A pollen record from Lake Maudit in northernmost Madagascar, a diatom assemblage from Lake Tritrivakely in central Madagascar, and a speleothem record from Tsimanampesotse in southwest Madagascar all feature drier conditions during Heinrich 1 and wetter conditions during the B-A (Gasse & Van Campo, 1998;Scroxton et al., 2019;Teixeira et al., 2021).This suggests that millennial-scale hydroclimate variability was coherent across much of Madagascar during the last deglaciation.
Additionally, there are multiple hydroclimate records from the tropical East African mainland that also show close agreement with those from Madagascar during the deglaciation.Leaf wax δD reconstructions from Lake Malawi, Lake Tanganyika, and Lake Rutundu in the African Great Lakes region feature stability during the LGM, shifts to drier conditions during HS-1, and a recovery during the B-A (Garelick et al., 2021;Konecky et al., 2011;Tierney et al., 2008).
Paleoclimate modeling studies focusing on the tropical Indian Ocean and East Africa during the LGM and deglaciation show agreement with proxy reconstructions.A previous study using the Community Earth System Model 1.2 found an increase in SON and DJF precipitation over Madagascar and the western Indian Ocean during the LGM compared to preindustrial simulations (DiNezio et al., 2018;Thirumalai et al., 2019).Additionally, a transient model of the deglaciation and Holocene (TRaCE) supports drier conditions during HS-1 and wetter conditions during the B-A in southeast Africa (Otto-Bliesner et al., 2014).This broad agreement across proxy records and models indicates that AB12 δ 18 O is representative of the tropical western Indian Ocean and East Africa, which is consistent with results from an isotope-enabled model (Duan et al., 2021).
Surprisingly, this coherence is not restricted to the Southern Hemisphere, with a leaf wax δD reconstruction from the Gulf of Aden and a speleothem record from Socotra Island also showing a drier HS-1 and wetter B-A (Figure 6) (Shakun et al., 2007;Tierney & DeMenocal, 2013).This coherence goes beyond the African continent, with same-sign millennial-scale variability observed in deglacial Indian Summer Monsoon and East Asian Monsoon reconstructions (Tierney et al., 2016;Wang et al., 2001).This feature of widespread drought in the Afro-Asian monsoon domain during HS-1 has been noted in previous work (Stager et al., 2011).
Though AB12 shows good agreement with many records from the Afro-Asian monsoon domain, there are inconsistencies with others from the region.Sedimentary runoff records from the northeastern Mozambique Channel have been interpreted to reflect southern African monsoon strengthening during HS-1 (Ma et al., 2021).While these records show marked agreement through some intervals, such as the Younger Dryas, only one out of three  (Shakun et al., 2007), Gulf of Aden core P178-15P leaf wax δD (Tierney & DeMenocal, 2013), Lake Rutundu leaf-wax derived precipitation δD (Garelick et al., 2021), Lake Tanganyika cores NP04-KH04-3A-1K and NP04-KH04-4A-1K leaf wax δD (Tierney et al., 2008), Lake Maudit pollen record, north Madagascar (Teixeira et al., 2021), Anjohibe speleothem AB12 δ 18 O, northwest Madagascar (this study), Tsimanampesotse speleothem MT1 (dark red) and AD4 (light red) δ 18 O, southwest Madagascar (Scroxton et al., 2019).Gray and red shading indicate cold and warm NH intervals, respectively.cores reconstruct a notably wetter HS-1.Additional work improving consensus among sedimentary archives from this under-sampled region will be essential prior to cross-validation with other records.
Paleoclimate data suggests conditions in subtropical southern Africa were opposite those recorded in tropical East Africa.A precipitation stack derived from South African pollen records reconstructs wetter conditions during HS-1, which stands in contrast to AB12 and other tropical African records (Chevalier & Chase, 2015).Similarly, a runoff record from the Zambezi basin suggests wetter conditions during HS-1 and drier conditions during the B-A (Schefuß et al., 2011).It is known that rainfall in tropical East Africa and subtropical southern Africa are anti-correlated over interannual and decadal timescales due to ENSO and the IOD (Zhang et al., 2015).Earlier work noted this phenomenon as well, finding negative correlations in OLR, SSTs, and rainfall between northern Madagascar and southern Africa (Jury & Pathack, 1991).This is likely because southern African sites are influenced by austral westerlies and southern air masses more than sites closer to the tropics (Weldeab et al., 2014).The southern African rainfall dipole might explain some of the reconstructed differences between AB12 in northwest Madagascar, which takes the tropical East African signal, and other African sites further south.

The Role of Zonal Indian Ocean Variability
The interhemispheric temperature gradient drives changes in the global zonal-mean ITCZ position, which in turn affects cross-equatorial atmospheric heat transport (Donohoe et al., 2013;McGee et al., 2014;T. Schneider et al., 2014).During NH stadials, the NH is cooler relative to the SH, driving the ITCZ south and increasing northward cross-equatorial atmospheric heat transport to compensate for hemispheric energy imbalance likely spurred by a weakening Atlantic Meridional Overturning Circulation (AMOC) (McGee et al., 2014;McManus et al., 2004;Otto-Bliesner et al., 2014;T. Schneider et al., 2014).The opposite is true for instances when the NH is warm relative to the SH, such as the B-A, in which the ITCZ shifts north.It is surprising that the AB12 multiproxy record reconstructs drier conditions at Anjohibe during Heinrich Stadials 1 and 2 (i.e., a warmer SH relative to NH) and wetter conditions during the Bølling-Allerød (i.e., a cooler SH relative to NH).This is the opposite response to that expected from a meridional ITCZ shift into the warmer hemisphere.Furthermore, agreement between many tropical East African hydroclimate records north and south of the equator through abrupt millennial-scale changes during the deglaciation cannot be explained by interhemispheric ITCZ shifts alone, which would result in an anti-phasing on either side of the equator rather than the coherence observed.
Here, we argue zonal SST gradients in the tropical Indian Ocean drive millennial-scale hydroclimate variability in tropical East Africa.To calculate the tropical Indian Ocean zonal SST gradient, we use the Last Glacial Maximum Reanalysis (LGMR), a product which integrates global SST records into the water isotope-enabled Community Earth System Model versions 1.2 and 1.3 (Osman et al., 2021).See Figure 1 in Osman et al. (2021) for details regarding the spatial and temporal distribution of SST records used as well as the proxies used in each record.Prior to calculating the gradient, a 10ka highpass filter was applied to global SSTs to isolate centennial-to-millennial scale variability.SSTs over the western and eastern IOD regions were then averaged and subtracted to generate a gradient (region outlines shown in Figure 2).Uncertainties associated with LGMR SSTs were added in quadrature to propagate errors through the calculation.Although a gradient calculated from proxies alone shows the same timing and direction of change as a gradient calculated from the assimilated product, their magnitudes are different (not shown), indicating a possible model bias which constitutes an important direction for future research.It is also worth noting that the unfiltered LGMR-derived gradient shows a strong trend from 20 to 10 ka which closely tracks deglacial warming.It is unclear why AB12 δ 18 O does not show this deglacial trend, but it is possible that AB12's growth termination before the end of the deglaciation precludes analysis of these longer-term signals.
Variability in the tropical Indian Ocean zonal SST gradient coincides with instances of abrupt millennial-scale variability during the last glacial period and deglaciation (Figure 7).HS-1 and the YD are associated with negative values (i.e., relatively warmer eastern SSTs), and the B-A is associated with positive values (i.e., relatively warmer western SSTs).Remarkably, AB12 δ 18 O closely tracks the tropical Indian Ocean zonal SST gradient, particularly through HS-1 and the B-A, with higher (lower) δ 18 O associated with lower (higher) values of the SST gradient (Figures 7a and 7b).In other words, a stronger (weaker) southeast African monsoon is associated with warmer (cooler) SSTs in the western Indian Ocean relative to the eastern Indian Ocean.
A correlation analysis similar to that of Figure 2 was performed using LGMR δ 18 O of precipitation in NW Madagascar and Indian Ocean SSTs (Figure S6 in Supporting Information S1).This correlation demonstrates that local δ 18 O values are significantly correlated with the tropical Indian Ocean zonal SST gradient further highlighting these zonal relationships at the millennial scale.Additional support for a zonal mechanism during HS-1 comes  (Osman et al., 2021).The dark red line is the LGMR ensemble mean, and the darker and lighter shading represents 1σ and 2σ uncertainty, respectively.A 10 ka highpass filter was applied to LGMR SSTs to isolate millennial-scale variability and remove long-term trends.See Figure S3 in Supporting Information S1 for the calculation of the gradient.from recent work using iTRACE, a transient, isotope-enabled model.During HS-1, iTRACE shows an east-west dipole in SST anomalies, precipitation, and δ 18 O across the tropical Indian Ocean, all consistent with the findings presented here (Du et al., 2023).The calculation of an Indian Ocean zonal gradient was similarly performed in a previous study which hypothesized the importance of these gradients in influencing precipitation changes over the Zambezi basin (Weldeab et al., 2014).
The connection between tropical SST gradients and East African hydroclimate is forged by the Indian Ocean Walker circulation which is sensitive to zonal tropical SST gradients and thermocline depth via the Bjerknes feedback (X.Wang & Wang, 2014).A synthesis of modeled and proxy thermocline depth reconstructions during HS-1 show a shallower thermocline in the western Indian Ocean and a deeper thermocline in the eastern Indian Ocean, a configuration favorable for the strengthening of the Indian Ocean Walker circulation (Du et al., 2023).This would weaken vertical ascent and convection over East Africa, thereby reducing the strength of regional monsoons.This ocean-atmosphere connection is also supported by evidence from the GFDL CM2.1 simulation in which tropical Indian Ocean SSTs drive changes in the strength of the Walker circulation and East African hydroclimate over multidecadal timescales independent of tropical Pacific SSTs (Tierney et al., 2013).Here we hypothesize that this mechanism also extends to millennial timescales and offers an explanation for the coherent interhemispheric hydroclimate variability reconstructed across tropical East Africa during the deglaciation.
The effects of variable Walker circulation strength on δ 18 O in the modern have been previously explored, lending further support to the zonal hypothesis.An analysis of precipitation isotopes from a network of observations revealed that global precipitation δ 18 O is sensitive to changes in the strength of the Indo-Pacific Walker circulation (Falster et al., 2021).More locally, the modern IOD was found to drive changes in precipitation δ 18 O over East Africa in observations and an isotope-enabled model (ECHAM-4 AGCM) (Vuille et al., 2005).It is plausible Indian Ocean Walker circulation strength influences precipitation amount locally, which would explain the relationship between large-scale atmospheric circulation regimes and the δ 18 O of regional precipitation over Anjohibe and tropical East Africa.This is supported by a correlation between the δ 18 O of precipitation in NW Madagascar and Indian Ocean SSTs in the LGMR (Figure S6 in Supporting Information S1).

High Latitude Teleconnections
Deglacial variability in the zonal SST gradient also closely tracks North Atlantic sedimentary 231 Pa/ 230 Th, a proxy for AMOC strength (Figures 7b and 7c) (Böhm et al., 2015;McManus et al., 2004).Reduced AMOC strength during HS-1 and the YD was synchronous with a cooler western Indian Ocean and warmer eastern Indian Ocean, and strong Atlantic overturning during the B-A resulted in the opposite pattern (Figures 7b and 7c).Although the LGMR does not cover both HS-2 pulses, the 231 Pa/ 230 Th records do, and they also might show reduced AMOC strength synchronous with reduced monsoon intensity in NW Madagascar within marine sedimentary archive age errors (Figures 7a and 7c).This suggests the strength of the tropical Indian Ocean SST gradient and southeast African monsoon were tightly linked with North Atlantic conditions, particularly during the deglaciation.
While Walker Cell variability and ocean-atmosphere coupling offer a compelling explanation for the zonal structure of millennial-scale phenomena in the Indian Ocean, the dynamics of North Atlantic-Indian Ocean teleconnection over these timescales remains largely unknown.Why would abrupt millennial-scale events altering interhemispheric temperature gradients result in east-west SST gradient changes in the Indian Ocean?Recent work hypothesizes that a weakening of the SH Hadley circulation in response to warming of the SH during HS-1 weakens southeasterly trades in the eastern Indian Ocean, reducing upwelling and driving warming (Du et al., 2023).A mechanism involving atmospheric teleconnections between the North Atlantic and Arabian Sea driving western Indian Ocean cooling during HS-1 has also been proposed in previous work (Tierney et al., 2016).A different possibility is that Africa cools more than the Indian Ocean which intensifies westerlies across the basin and strengthens the zonal SST gradient.Over interannual timescales, Plinean eruptions are capable of altering the Indian Ocean Walker circulation and zonal SST gradient by changing land-sea temperature gradients between Africa and the Indian Ocean, resulting in an nIOD-like state in simulations (Izumo et al., 2018).This leads to reduced rainfall over tropical East Africa similar to what is reconstructed through HS-1.Though volcanism is an imperfect analog to Heinrich stadials, it is plausible they might share a similar mechanism involving land-sea gradients between Africa and the Indian Ocean.Another mechanism is AMOC's influence on Indo-Pacific climate variability via tropical South Atlantic warming, with AMOC shutdown resulting in a La Niña mean state (Orihuela-Pinto et al., 2022) which often leads to drier conditions in East Africa (Ummenhofer et al., 2018).This is consistent with previous work which finds that the mean state of Pacific multidecadal variability modulates Indian Ocean Dipole event phasing (Ummenhofer et al., 2017).

Caveats to the Zonal Hypothesis
If millennial-scale changes in zonal SST gradients do in fact drive changes to Indian Ocean Walker circulation strength, then it follows that opposite changes in precipitation and δ 18 O would be recorded through these intervals at sites in the tropical eastern Indian Ocean, such as those within the Maritime Continent.However, deglacial hydroclimate reconstructions from the tropical eastern Indian Ocean instead show anti-phased anomalies between hemispheres and lack the interhemispheric agreement observed in western Indian Ocean records (cf.Ayliffe et al., 2013;Denniston et al., 2013;Krause et al., 2019;Wurtzel et al., 2018).Records from China, the Philippines, Borneo, and Sumatra show an opposite response to sites further south in Flores and Australia, and areas in-between, such as Sulawesi, show little to no change during the deglaciation (see Figure 4 in Krause et al., 2019;Figure 8 in Wurtzel et al., 2018).A possible explanation for this discrepancy is that during HS-1, ITCZ shifts overprinted the effects of zonal Indian Ocean dynamics in the eastern Indian Ocean (Du et al., 2023).Increased aridity over the Maritime Continent during glacial periods due to Sunda and Sahul shelf exposure could also further complicate these interpretations (DiNezio et al., 2018;Windler et al., 2019).
Finally, it is important to mention these findings only provide evidence for a zonal mechanism during HS-1 the B-A, and possibly HS-2.There is some disagreement between the reconstructed SST gradient and the AB12 speleothem record during the LGM ∼22-20ka; additional work will be required to determine whether this disagreement results from age model uncertainties in the SST records assimilated in the LGMR or other factors.It also remains unknown whether the zonal SST gradient drove hydroclimate changes during other stadials and interstadials or variability over even longer timescales.For instance, there is evidence a different mechanism can operate under interglacial boundary conditions.Several speleothem records from northwest Madagascar show wetter conditions during the 8.2 ka event, possibly due to a southward ITCZ shift (Duan et al., 2021;Voarintsoa, Railsback, et al., 2017;L. Wang et al., 2019).It is also unclear whether it is the IOD operating over millennial timescales driving the reconstructed changes or a lower frequency zonal mode.Recent evidence suggests that Pacific decadal variability can be attributed to the residual of long-term ENSO mean state changes rather than a distinct multidecadal mode (Power et al., 2021).Additional work developing interannually resolved proxy records and paleoclimate simulation ensembles could shed light on this issue.

Conclusions
Ice volume and temperature-corrected δ 18 O from the Last Glacial Period is more depleted than δ 18 O from the Common Era at Anjohibe, likely indicating wetter conditions and a stronger monsoon in NW Madagascar during the glacial compared to the present.Over millennial timescales, the AB12 multiproxy hydroclimate record suggests moderately wetter conditions in NW Madagascar during the LGM, drier conditions during HS-1 and possibly during HS-2, and wetter conditions during the B-A.As the opposite would be expected from ITCZ shifts during Heinrich stadials and the B-A, this record stands in opposition to hypotheses espousing meridional ITCZ shifts driving millennial-scale variability in the tropical Indian Ocean during the last glacial period and deglaciation.Interhemispheric agreement between AB12 and other paleorecords from the western Indian Ocean through the millennial-scale events of the last deglaciation further undermines the meridional ITCZ shift hypothesis.We instead propose that tropical Indian Ocean SST gradients drive millennial-scale hydroclimate variability in NW Madagascar and the rest of East Africa by affecting the strength of the regional Walker circulation.This is supported by the close agreement between AB12 δ 18 O and a LGMR-derived tropical Indian Ocean zonal SST gradient.
With many abrupt deglacial climate changes hypothesized to be ultimately driven by North Atlantic variability and AMOC strength, it is prudent to understand the mechanisms linking this region to the tropical Indian Ocean.This is especially important considering model projections of a slowing AMOC under high emissions scenarios (Liu et al., 2020) and the strong influence of rapidly warming tropical Indian Ocean SSTs on AMOC strength and pantropical hydroclimate variability (Hu & Fedorov, 2019).It is therefore critical to continue exploring the dynamics of these teleconnections and their influence on past and future hydroclimate for vulnerable Indian Ocean rim communities.data from https://www.dwd.de/EN/ourservices/gpcc/gpcc.html(U.Schneider et al., 2014), and TRMM data from https://gpm.nasa.gov/data/directory(Kummerow et al., 1998)

Figure 2 .
Figure2.Correlations between Indian Ocean SSTs and site-specific precipitation at Anjohibe, northwest Madagascar.(a) SON-mean precipitation anomalies over Anjohibe correlated with SON-mean Indian Ocean SST anomalies from 1982 to 2020.Anjohibe's location is shown as a green dot.Red indicates a positive correlation, and blue a negative correlation.The Spearman correlation statistic (ρ) was used to account for precipitation's non-normal distribution.Areas within the black contours are significant at the 95% confidence level based on a nonparametric isospectral method accounting for autocorrelation(Ebisuzaki, 1997;Khider et al., 2022).The black dashed boxes represent the regions used to calculate the Dipole Mode Index (DMI) (western: 50-70°E, 10°N-10°S; eastern: 90-110°E, 0-10°S)(Saji et al., 1999).(b)-(d): Scatterplots of Anjohibe precipitation in SON versus western Indian Ocean SSTs in SON (b), eastern Indian Ocean SSTs in SON (c), the Dipole Mode Index (DMI), the difference between the western and eastern Indian Ocean regions inSON (d).The solid black lines represent the lines of best fit, and gray shading represents the 95% confidence intervals based on 1,000 bootstrapped resamples.Note: the IOD is known to peak in SON.SST data are from OISST, and site-specific precipitation data are from GPCC.

Figure 3 .
Figure 3. Age-depth relationship for speleothem AB12 based on Monte-Carlo modeling of U-Th ages.The median age is shown as a dashed black line, and the 95% confidence interval is shown in light purple.Ages used for the model are shown in black, outliers are shown in red.Error bars on these points represent the 2σ uncertainties for each age determination.One thousand Monte-Carlo simulations were run in COPRA to generate median ages and confidence intervals(Breitenbach et al., 2012).A photograph of AB12 is shown to the right of the plot.

Figure 4 .
Figure 4. AB12 oxygen (top) and carbon (bottom) stable isotope proxy time series.VPDB = the Vienna Pee Dee Belemnite standard.Black dots above the δ 18 O curve represent U-Th ages and their associated 2σ errors.Gray and red shading indicates cold and warm NH intervals, respectively.

Figure 7 .
Figure 7. Drivers of abrupt deglacial hydroclimate variability in northwest Madagascar.(a): AB12 δ 18 O record (this study).The black line shows the 12-point running mean, and the gray line the raw data.(b): Zonal Indian Ocean SST gradient calculated from the Last Glacial Maximum Reanalysis (LGMR)(Osman et al., 2021).The dark red line is the LGMR ensemble mean, and the darker and lighter shading represents 1σ and 2σ uncertainty, respectively.A 10 ka highpass filter was applied to LGMR SSTs to isolate millennial-scale variability and remove long-term trends.See FigureS3in Supporting Information S1 for the calculation of the gradient.(c): North Atlantic 231 Pa/ 230 Th records, a proxy for AMOC strength.Dark and light blue curves indicate the McManus et al. (2004) and Bohm et al. (2015) 231 Pa/ 230 Th records, respectively, and shading represents 2σ uncertainty.Gray shading represents HS-1, HS-2, and the YD; red shading the B-A.
Figure 7. Drivers of abrupt deglacial hydroclimate variability in northwest Madagascar.(a): AB12 δ 18 O record (this study).The black line shows the 12-point running mean, and the gray line the raw data.(b): Zonal Indian Ocean SST gradient calculated from the Last Glacial Maximum Reanalysis (LGMR)(Osman et al., 2021).The dark red line is the LGMR ensemble mean, and the darker and lighter shading represents 1σ and 2σ uncertainty, respectively.A 10 ka highpass filter was applied to LGMR SSTs to isolate millennial-scale variability and remove long-term trends.See FigureS3in Supporting Information S1 for the calculation of the gradient.(c): North Atlantic 231 Pa/ 230 Th records, a proxy for AMOC strength.Dark and light blue curves indicate the McManus et al. (2004) and Bohm et al. (2015) 231 Pa/ 230 Th records, respectively, and shading represents 2σ uncertainty.Gray shading represents HS-1, HS-2, and the YD; red shading the B-A.
B.H.T., R.D., N. S., S.J.B., and D.M. acknowledge support from NSF awards AGS-1702891/1702691 and AGS-2102923/2102975, L.R.G. and S.J.B. from NSF award BCS-1750598, and D.M. from NSF award EAR-1439559.Support for this research was also provided by a core center Grant P30-ES002109 from the National Institute of Environmental Health Sciences, National Institutes of Health.This research was conducted under a collaborative accord between the University of Massachusetts Amherst (Department of Anthropology, Department of Geoscience) and Université d'Antananarivo, Madagascar (Bassins sédimentaires, Evolution, Conservation).We thank the Ministry of Higher Educational and Research, Ministry of Mines and Strategic Resources, and Ministry of Communication and Culture which provided permits for field research.This work was supported by NSF award #2102975.We thank Adam Jost for assistance with U-Th chemistry and analysis, Angela Douglass and Markey Freudenburg-Puricelli for assistance with trace metal sample preparation, Raspberry Williams for stable isotope sample preparation and analysis, and Caroline Ummenhofer for fruitful conversations on contemporary Indian Ocean climate dynamics.The authors would also like to thank Jessica Tierney and one anonymous reviewer for their helpful comments, which substantially improved the paper.