Tropical vegetation productivity and atmospheric methane over the last 40,000 years from model simulations and stalagmites in Sulawesi, Indonesia

Abstract Recent research has shown the potential of speleothem δ13C to record a range of environmental processes. Here, we report on 230Th-dated stalagmite δ13C records for southwest Sulawesi, Indonesia, over the last 40,000 yr to investigate the relationship between tropical vegetation productivity and atmospheric methane concentrations. We demonstrate that the Sulawesi stalagmite δ13C record is driven by changes in vegetation productivity and soil respiration and explore the link between soil respiration and tropical methane emissions using HadCM3 and the Sheffield Dynamic Global Vegetation Model. The model indicates that changes in soil respiration are primarily driven by changes in temperature and CO2, in line with our interpretation of stalagmite δ13C. In turn, modelled methane emissions are driven by soil respiration, providing a mechanism that links methane to stalagmite δ13C. This relationship is particularly strong during the last glaciation, indicating a key role for the tropics in controlling atmospheric methane when emissions from high-latitude boreal wetlands were suppressed. With further investigation, the link between δ13C in stalagmites and tropical methane could provide a low-latitude proxy complementary to polar ice core records to improve our understanding of the glacial–interglacial methane budget.


INTRODUCTION
Direct measurements of atmospheric methane concentrations in air trapped within layers of ice have provided a high-quality record of methane variability over the last 800,000 yr (e.g., Brook et al., 1996Brook et al., , 2000;;Schaefer et al., 2006;Grachev et al., 2007;Fischer et al., 2008;Loulergue et al., 2008;Baumgartner et al., 2012;Möller et al., 2013;Rhodes et al., 2015).The ice core records show that atmospheric methane is sensitive to climate variations, with glacial-interglacial amplitudes of around 300-400 parts per billion by volume (ppbv), which broadly change alongside temperature.Over the last 40,000 yr, methane concentrations reached a minimum of ∼370 ppbv during the last glacial maximum (LGM; ∼21 ka, where ka is 1000 yr before 1950 CE), before increasing to ∼675 ppbv at ∼10 ka.While the concentration of atmospheric methane is well constrained, quantifying the changing sources and sinks of methane through time is hindered by the limited availability of information about methane sources (e.g., Brook et al., 2000;Fischer et al., 2008;Levine et al., 2011;Hopcroft et al., 2017).
Methane is relatively well mixed within the atmosphere; however, the modern-day dominance of methane sources in the Northern Hemisphere creates a gradient in methane between the Northern and Southern Hemispheres (Chappellaz et al., 1997;Brook et al., 2000;Dällenbach et al., 2000;Baumgartner et al., 2012).The difference in methane concentrations between Greenland and Antarctica (referred to as the methane "gradient") has therefore been used to infer the hemispheric contribution to methane sources through time.The methane gradient has been relatively stable over the last 32 ka, suggesting that northern methane sources were not completely shut off during the LGM, when large areas of the high latitudes were frozen (Baumgartner et al., 2012).Quantifying changes in the methane gradient, however, is less useful for attributing sources to the tropics, which contribute methane to both the Northern and Southern Hemisphere budgets throughout the year with the seasonal migration of the Intertropical Convergence Zone.
Model simulations of methane emissions since the LGM similarly suggest that wetlands were the predominant source of atmospheric methane.Kleinen et al. (2020) found that wetland emissions make up 93-96% of the net methane flux during the LGM.General circulation models coupled with vegetation models also tend to suggest a more dominant role for the tropics during the LGM, when large Northern Hemisphere ice sheets and cooler climates reduced boreal wetland areas (Valdes et al., 2005;Kaplan et al., 2006), and changes since the LGM are largely thought to be driven by wetlands (Kleinen et al., 2023).The LGM wetland reduction is supported by isotopic analyses of atmospheric methane in ice cores (Fischer et al., 2008;Bock et al., 2017;Hopcroft et al., 2018); however, attribution of methane isotopes to a specific source is not definitive (Möller et al., 2013).Large ice sheets and the associated lowering of sea level also influence methane emissions with the exposure and enlargement of low-lying tropical wetland areas, such as shallow maritime continents (Ridgwell et al., 2012;Kleinen et al., 2020Kleinen et al., , 2023)).
At present, there is no proxy for past tropical wetland methane emissions.It has been suggested, however, that carbon isotope ratios (δ 13 C) in speleothems may provide insights into tropical methane emissions through time by recording changes in vegetation productivity, which is closely related to methane production in wetlands (Burns, 2011;Griffiths et al., 2013;Fohlmeister et al., 2020).In wet tropical settings, speleothem δ 13 C primarily reflects C 3 vegetation productivity, with most of the carbon in infiltrating waters originating from CO 2 in the soil atmosphere, produced by vegetation root respiration and microbial activity (e.g., Genty et al., 2001;Wong and Breecker, 2015;Fohlmeister et al., 2020).Light carbon from soil CO 2 and heavier carbon derived from bedrock combine to influence the δ 13 C of dissolved inorganic carbon in cave seepage waters (Vogel and Kronfeld, 1997;Genty et al., 2001;Hou et al., 2003;Griffiths et al., 2012;Wong and Breecker 2015;Burns et al., 2016).Cave seepage waters then precipitate as speleothems, preserving a δ 13 C signature driven primarily by changes in vegetation productivity (Dorale and Liu, 2009;Breecker et al., 2012;Hartmann et al., 2013;Burns et al., 2016;Fohlmeister et al., 2020).
In this study, we use stalagmite δ 13 C from southwest Sulawesi, Indonesia (Fig. 1) as a record of vegetation productivity to explore the contribution of Indonesia and the broader tropics to the atmospheric methane budget over the last 40 ka.The robust age control of the Sulawesi stalagmite δ 13 C record, afforded by precise 230 Th dating, enables us to examine the link between tropical vegetation productivity and atmospheric methane concentrations recorded in ice cores over this period.Hypotheses drawn from the stalagmite record are then tested against model output from the Hadley Centre Coupled Model version 3 (HadCM3) and the Sheffield Dynamic Global Vegetation Model (SDGVM).The model results are used to explore possible relationships between Indonesian vegetation productivity and the global methane budget of the last 40 ka.

Stalagmites
We present new δ 13 C, 234 U/ 238 U, and Mg/Ca records for two 230 Th-dated stalagmite specimens collected from Gempa Bumi Cave in tower karst overlain by tropical rainforest in the Maros limestone district of southwest Sulawesi, Indonesia (Fig. 1), which were previously analysed for δ 18 O by Krause et al. (2019).Stalagmite GB09-3 (Fig. 2) was collected ∼300 m from the entrance of the cave (5°S, 120°E, ∼140 m above sea level).Isotopic equilibrium deposition of GB09-3 calcite was tested using a second stalagmite (GB11-9; Fig. 2) collected ∼10 m away from GB09-3.GB11-9 overlaps GB09-3 for the 40-26 ka period, providing a reasonable length of time over which to compare the isotopic behaviour of the two stalagmites.Good replication of δ 18 O and δ 13 C within this common growth interval provides a valuable test for isotopic equilibrium, because it is unlikely that different stalagmites could produce similar disequilibrium fractionated signals (e.g., Wang et al., 2001;Dorale and Liu, 2009). 230Th/ 234 U chronology The chronologies for stalagmites GB09-3 and GB11-9 are based on 44 U-Th dates described in Krause et al. (2019) and presented in supplementary table 1 therein.Briefly, 36 samples were collected for uranium and thorium isotope analysis to develop the chronology for GB09-3, and eight samples were analysed for GB11-9 (Fig. 2).The dating samples, with an average weight of 65 mg, were cut adjacent to the stable isotope sampling tracks on stalagmite slabs cut parallel to central growth axes.Samples were analysed using multi-collector inductively coupled plasma mass spectrometry at the University of Melbourne (Hellstrom, 2003) and the University of Minnesota (Cheng et al., 2013).All samples were corrected for small amounts of detrital thorium using an initial [ 230 Th/ 232 Th] ratio of 3.0 ± 0.75 determined by stratigraphic constraint analysis of the measured U-Th dates (Hellstrom, 2006).Two age outliers were not included in the final age model for GB09-3.All corrected dates are in stratigraphic order, within error, and have a median 2-sigma age uncertainty of ±1.6%.In the present study, the initial 234 U/ 238 U values ([ 234 U/ 238 U] i ) for GB09-3 and GB11-9 are used as an indicator of seepage water infiltration rates.

Stable isotope analysis
Stalagmites GB09-3 and GB11-9 were slabbed and micro-milled with a 1-mm-diameter drill bit along their central growth axes at intervals of 0.8-1.2mm and 0.7 mm, respectively, equating to an average sample resolution of ∼55 yr.Measurements of δ 18 O and δ 13 C were conducted on 755 samples for GB09-3 and 323 samples for GB11-9.Sample powders (∼200 μg) were analysed on a Finnigan MAT 251 mass spectrometer equipped with an automated Kiel carbonate reaction device.CO 2 was liberated from the carbonate by reaction with 105% H 3 PO 4 under vacuum at 90°C.Measurements of δ 18 O and δ 13 C were corrected using the NBS19 (δ 18 O = −2.20‰,δ 13 C = 1.95‰) and NBS18 (δ 18 O = −23.0‰,δ 13 C = 5.0‰) standards and are reported in delta notation relative to Vienna Pee Dee Belemnite (VPDB).The analytical precision of the measurements for aliquots of NBS19 run in parallel with the stalagmite samples was ± 0.05‰ for δ 18 O and ±0.02‰ for δ 13 C (n = 270, 1σ).
The reproducibility for replicate aliquots of the stalagmite samples was determined to check for sample homogeneity.In the first instance, samples with a mass spectrometer measurement cycle standard deviation greater than 0.05‰ (for δ 18 O) were reanalysed to minimise errors related to the mass spectrometry.Additionally, samples giving δ 13 C values that deviated significantly from adjacent values in the time series were reanalysed to ensure that abrupt variations in the data set were not analytical artifacts.The mean standard error for duplicate/triplicate analyses of δ 13 C was 0.02‰ for GB09-3 (n = 126), and 0.03‰ for GB11-9 (n = 46).
Mg/Ca and Sr/Ca analysis Analysis of Mg/Ca and Sr/Ca in stalagmite GB09-3 was conducted to check for the occurrence of prior calcite precipitation (PCP) and its potential effect on δ 13 C. PCP is driven by the process of seepage waters degassing along flow pathways, resulting in "upstream" precipitation of calcite before the seepage waters reach the stalagmite surface (Fairchild and Treble, 2009).During drier conditions, PCP increases 13 CO 2 , Mg 2+ , and Sr 2+ (relative to Ca 2+ ) in drip waters and raises stalagmite δ 13 C, Mg/Ca and Sr/Ca simultaneously (Baker et al., 1997).
Measurements of Mg/Ca and Sr/Ca were made on aliquots of the same samples analysed for stable isotopes in GB09-3.Every second sample of GB09-3 (n = 378) was analysed at the Australian National University, Research School of Earth Sciences (RSES) by inductively coupled plasma atomic emission spectroscopy using methods based on Schrag (1999).These samples were measured on 0.5 mg (n = 189) and 1.5-2 mg (n = 189) aliquots dissolved in 5 mL of 2% v/v HNO 3 .Analytical precision was determined by repeat analyses of an in-house laboratory  (Stott et al., 2002;Linsley et al., 2010 and references therein), cave temperature record for Gunung Mulu National Park, northern Borneo (Løland et al., 2022), and leaf wax records for Sulawesi (Russell et al., 2014;Wicaksono et al., 2015Wicaksono et al., , 2017)).Base maps were created in QGIS 3.20 (https://qgis.org/en/site) using Shuttle Radar Topography Mission 1 Arc-Second Global by NASA/NGS/USGS (2015-01-01 EPSG4326_31m).
(coral) standard.Standards bracketed each stalagmite sample to correct for any instrument drift occurring within the runs.The analytical precision (relative standard deviation [RSD]) for repeat measurements on the laboratory standard was 0.70% for Mg/Ca and 0.64% for Sr/Ca (n = 376).Approximately every fourth sample of  was analysed at the Australian Nuclear Science and Technology Organisation (ANSTO) using methods based on de Villiers et al. (2002).These measurements were made on 1 mg aliquots dissolved in 5 mL of 3% v/v HNO3.The analytical precision for repeat measurements on the laboratory standard was 0.98% for Mg/Ca and 0.94% for Sr/Ca (RSD, n = 11).There is no significant offset between RSES and ANSTO measurements for Mg/Ca; however, there is a relatively consistent Sr/Ca offset of ∼29% between the two facilities (0.0052 mmol/mol from the two-record average of 0.0178 mmol/ mol).This offset is likely due to the low stalagmite Sr concentrations in solution being near instrument detection limits.For these reasons, Sr/Ca is not included in the results.

HadCM3 general circulation model
The HadCM3 is a coupled ocean-atmosphere-sea ice general circulation model (Gordon et al., 2000).The resolution of the atmospheric model is 2.5°latitude by 3.25°longitude with 19 unequally spaced vertical levels (Gordon et al., 2000).The ocean model resolution is 1.25°latitude by 1.25°longitude with 20 unequally spaced layers extending to a depth of 5200 m.The sea-ice model uses a simple thermodynamic scheme (Cattle et al., 1995).Coupling between the model components occurs daily (Gordon et al., 2000).The HadCM3 simulations include dynamic vegetation simulated with the TRIFFID vegetation scheme (Cox, 2001), which allows feedbacks to the atmosphere from changes in the distribution and structure of vegetation over time.The precise configuration of the model is called HadCM3BM2.1Dand is fully described by Valdes et al. (2017).
Climate simulations with HadCM3 have been evaluated against observations (Gordon et al., 2000;Collins et al., 2001), proxy records (Singarayer and Valdes, 2010;DiNezio and Tierney, 2013), and other general circulation models (Braconnot et al., 2007a(Braconnot et al., , 2007b;;Flato et al., 2013).HadCM3 represents LGM climate conditions relatively well when compared with reconstructions and other Paleoclimate Modelling Intercomparison Project models (Braconnot et al., 2007a(Braconnot et al., , 2007b)).Moreover, when compared with LGM proxy data from the Indo-Pacific Warm Pool (IPWP), HadCM3 emerges as one of the few models to successfully capture the climate conditions recorded by both terrestrial and marine proxies in the IPWP region (DiNezio and Tierney 2013;DiNezio et al., 2016).For the purposes of this study, the HadCM3 climate model is used to drive the SDGVM vegetation model to investigate the drivers of regional methane emissions (Hopcroft et al., 2011(Hopcroft et al., , 2014;;Singarayer et al., 2011).

SDGVM vegetation and wetland model
The SDGVM is a global primary productivity and phytogeography model (Woodward et al., 1995;Beerling and Woodward, 2001).SDGVM is driven with inputs from HadCM3 to simulate dynamic changes in vegetation distribution and leaf area index and productivity in response to changing climate and atmospheric CO 2 concentrations.SDGVM accounts for the main factors driving vegetation productivity, including climate (surface temperature, precipitation, relative humidity), atmospheric CO 2 concentration, and soil characteristics.Plant species are broadly categorised into "plant functional types" (PFTs), allowing tractable calculations of global vegetation distribution and facilitating simulation of their dynamic response to other model variables.The response of PFTs is driven by sensitivities to temperature, net precipitation (precipitation minus evapotranspiration), CO 2 , and inter-PFT competition.Vegetation response to changing climate and environment is not instantaneous, but is dependent on the cycle of mortality and establishment of PFTs.
The SDGVM does not model methane emissions in its standard configuration, and an additional methane module is used  (2019).The average growth rates are 1.74 mm/100 yr for GB09-3 and 1.40 mm/ 100 yr for GB11-9 (for 40-26 ka), with no detectable hiatuses.Data from the bottom of GB11-9 are not included in this study.
to simulate wetland extent and methane emissions.The methane module uses topography, surface air temperature, soil moisture, soil type, and soil respiration outputs from HadCM3 and the SDGVM to calculate methane emissions (see Singarayer et al. [2011] and Wania et al. [2013] for a detailed discussion of the methods used to calculate methane emissions.) The SDGVM participated in the Wetland and Wetland CH 4 Inter-comparison of Models Project (WETCHIMP) project in 2013, which aimed to compare and validate available methane models (Melton et al., 2013;Wania et al., 2013) and the Global Carbon Project methane budget (Saunois et al., 2016).The spatial distribution and absolute amount of methane emissions from SDGVM compare well with those of other models, showing a maximum in emissions across the tropics, driven largely by emissions from the Amazon.Sensitivity experiments demonstrate that the SDGVM is responsive to changing CO 2 concentrations, air temperature, and precipitation (Melton et al., 2013).SDGVM showed a 40% increase in global methane emissions due to a 557 ppm step increase in CO 2 (from ∼300 to 857 ppm), compared with a multi-model mean of 73 ± 49%, and a 2.4% increase in methane due to a temperature increase of 3.4°C, versus a multimodel mean of −2.5 ± 21% (Melton et al., 2013).The tropics proved most sensitive to an increase in CO 2 concentrations via fertilisation of tropical vegetation, while the extratropics were most sensitive to an increase in temperature.
Exposed land area across Indonesia varied significantly over the last 40 ka, due to fluctuations in glacial-interglacial sea level, which dropped by up to ∼130 m relative to modern.For example, the LGM sea-level low stand resulted in Sunda Shelf exposure of ∼2.4 million km 2 (50% more expansive) compared with the present (Sathiamurthy and Voris, 2006).Terrestrial and marine paleoenvironmental studies show evidence for a substantial savannah corridor occupying the interior of the exposed Sunda Shelf during the LGM (Bird et al., 2005;Wurster et al., 2010Wurster et al., , 2019;;Nguyen et al., 2022;Cheng et al., 2023); however, the spatial extent of savanna versus forest is debated (e.g., Bird et al., 2005;Wurster et al., 2010).Modelled vegetation on exposed continental shelves during the LGM relies on the simulation of dynamic vegetation coverage within SDGVM.

Model outputs
HadCM3 and SDGVM were run at 1 ka resolution for the period 22-0 ka and 2 ka resolution for the period 40-22 ka (Singarayer and Valdes, 2010;Singarayer et al., 2011), resulting in a total of 32 time-slice simulations for the period 40 ka to present.HadCM3 was forced with orbital parameters, ice-sheet volume (and sea level), and greenhouse gases for each time slice, as described by Singarayer and Valdes (2010) and Singarayer et al. (2011).All time slices were run from an equilibrated preindustrial control run and forced with boundary conditions appropriate to the time slice being run.The model was then allowed to re-equilibrate under these new conditions for 500 model years.The results presented here represent the climatology of the last 30 yr of each model run.
Importantly, abrupt millennial-scale events were not simulated in this experiment.This is not a "transient" experiment; thus, a continually evolving climate was not simulated, but rather the time evolution of climate was simulated through the use of 1-2 ka snapshots.SDGVM was forced by the mean climatological outputs derived from each HadCM3 simulation to produce a dynamic vegetation response to the modelled climate time slices.An extended set of similar simulations back to 130 ka has been used to produce a simplified estimate of the changing contributions to atmospheric methane (Singarayer et al., 2011), and this shows good agreement with the observed changes in atmospheric methane over the last glacial-interglacial cycle.

RESULTS
The GB09-3 and GB11-9 stalagmite δ 13 C records are generally in good agreement across their interval of overlap (40-26 ka).The millennial-scale δ 13 C variability in GB09-3 is mostly reproduced in GB11-9, and the trends in the two records are similar (Fig. 3).There are periods when the records diverge (40-38 ka and around 33 ka), but fine-scale differences between records  (Breecker, 2017).Uncorrected δ 13 C is shown in grey.The large deglacial δ 13 C transition (green shading) encompasses two abrupt negative excursions at ∼14.7-14.5 ka and 11.7-11.6ka that mark the terminations of Heinrich stadial 1 (HS1) and the Younger Dryas (YD), respectively.The Bølling-Allerød (B-A) is also shown (yellow).(B) δ 18 O for GB09-3 and GB11-9 corrected for the effect of ice volume (Krause et al., 2019) with small ranges in δ 13 C are to be expected due to localised effects associated with the degree of water-gas exchange in the soil zone and different seepage water flow pathways (e.g., Partin et al., 2013;Fohlmeister et al., 2020).
It is important to note that the δ 13 C records have been corrected for the effect of atmospheric pCO 2 on the δ 13 C of C 3 plants (Schubert and Jahren, 2012).The transfer of this effect on carbon isotope fractionation in C 3 plants above a cave to stalagmites growing within the cave was identified by Breecker (2017) in a study assessing globally averaged speleothem δ 13 C records over the past 90 ka.They found that, after accounting for other processes, the effect of atmospheric CO 2 is best explained by a C 3 plant δ 13 C sensitivity of −1.6‰ for every 100 parts per million by volume (ppmv) increase in pCO 2 from the LGM to the Holocene.Therefore, it is important to correct for the change in stalagmite δ 13 C that occurs as a result of glacial-interglacial atmospheric pCO 2 before investigating δ 13 C as a recorder of glacial-interglacial vegetation productivity.The Sulawesi stalagmite δ 13 C values were adjusted by −1.6‰/100 ppmv (Breecker, 2017) relative to modern atmospheric pCO 2 (190 ppmv) using the Antarctic ice core composite pCO 2 record (Bazin et al., 2013;Bereiter et al., 2015).The corrected records are shown in Figure 3 and the Supplementary Material and are used throughout the analysis.
The δ 13 C time series for GB09-3 can be divided into three main sections: glacial (40-18 ka), deglacial (18-11 ka), and Holocene (11 ka to the present) (Fig. 3).The glacial state includes Marine Isotope Stage 3 and the LGM and is characterized by relatively high δ 13 C values.The deglacial interval contains a prominent ∼4.2‰ shift in corrected δ 13 C from the maximum δ 13 C value of −4.8‰ at 17.7 ka, near the onset of deglaciation (Pedro et al., 2011), to −9‰ at 11.3 ka.This transition from highest to lowest δ 13 C includes two abrupt negative excursions from 14.7 to 14.1 ka (1.3‰ decrease), and from 11.9 to 11.6 ka (1.4‰ decrease).Together, the magnitude of these two events is relatively large, accounting for about two-thirds of the total deglacial transition in δ 13 C.
The Holocene section of the record shows a surprisingly high degree of δ 13 C variability, most notably a prominent V-like pattern in the Early to Middle Holocene.During this time, the δ 13 C increases from about −8.5‰ at ∼11 ka to a brief maximum of −6.3‰ at 7.5 ka, before decreasing to around −8‰ in the Late Holocene.

Sulawesi stalagmite rainfall proxies
Three additional geochemical proxies are presented for stalagmites GB09-3 and GB11-9: δ 18 O, initial uranium isotope activity [ 234 U/ 238 U] i ), and Mg/Ca.The Sulawesi stalagmite δ 18 O data are explored in detail in Krause et al. (2019) and interpreted to reflect changes in rainfall and deep atmospheric convection over the IPWP.The deglacial transition towards wetter conditions, signified by lower δ 18 O values, occurs at ∼11.5 ka.
[ 234 U/ 238 U] i and Mg/Ca are sensitive to groundwater movement through the epikarst and along flow pathways leading to the stalagmites, thus serving as additional proxies for rainfall (e.g., Fairchild et al., 2000Fairchild et al., , 2006;;Hellstrom and McCulloch, 2000;Fairchild and Treble, 2009).Because uranium isotopes are not thought to be fractionated by natural processes such as calcite precipitation, [ 234 U/ 238 U] i is expected to reflect the activity ratio in seepage waters forming speleothems.The [ 234 U/ 238 U] i of seepage waters can be altered by groundwater residence time and water-rock interactions.During drier periods, when there is less water moving through the epikarst and longer residence times, [ 234 U/ 238 U] i can increase as a result of preferential leaching of 234 U from alpha recoil-weakened crystal lattice sites in limestone bedrock (Hellstrom and McCulloch, 2000).Because this effect is sensitive to the amount of surface area that seepage waters are exposed to, waters moving through capillaries and pore spaces may be more strongly influenced (Hellstrom and McCulloch, 2000).During wetter periods, when water is moving quickly through the epikarst and bedrock dissolution is more uniform, [ 234 U/ 238 U] i is expected to be relatively low.
Mg and Ca are sourced primarily from the bedrock during dissolution.The partition coefficient for Mg is less than one (Fairchild and Treble, 2009); thus, Ca is preferentially lost from solution during calcite precipitation.Therefore, Mg/Ca increases when precipitation occurs before seepage waters reach the surface of a stalagmite; this process is known as PCP (Fairchild et al., 2000;Fairchild and Treble, 2009).PCP tends to occur when infiltration rates are low, drip intervals are long, and seepage waters encounter an air-filled space with a pCO 2 lower than that in the seepage waters.Thus, PCP is more likely to be active during drier periods, resulting in higher Mg/Ca values.In contrast, during wetter periods, when the cave network is saturated and water moves continuously through the epikarst, PCP is reduced or absent, and Mg/Ca values are expected to be low (Fairchild and Treble, 2009).
Sulawesi stalagmite [ 234 U/ 238 U] i is relatively high throughout the glacial period before abruptly decreasing after ∼11.8 ka.This transition towards lower values coincides with the deglacial decrease in δ 18 O, and Mg/Ca also shifts to lower values (Fig. 3).The shift in Mg/Ca culminates with a marked stabilization of Mg/Ca variability from ∼10 ka to the present, with an average of 0.63 mmol/mol and variance of 0.01 mmol/mol (σ 2 ).Before the deglacial transition, from ∼40 to 11.5 ka, Mg/Ca swings between 0.78 and 2.68 mmol/mol, with an average of 1.5 mmol/mol and variance of 0.13 (σ 2 ).The stabilization of Mg/Ca at lower values following the deglacial transition, corroborated by [ 234 U/ 238 U] i , is interpreted to reflect wetter conditions.
Previous studies have shown that similar decreases in stalagmite δ 18 O during the mid-late stages of the last deglaciation are related to increased rainfall in the Indonesian region (Partin et al., 2007;Griffiths et al., 2009;Ayliffe et al., 2013;Carolin et al., 2013).The multiproxy agreement between the Sulawesi δ 18 O, [ 234 U/ 238 U] i , and Mg/Ca records, alongside other regional rainfall δ 18 O records, supports our interpretation of an increase in rainfall amount initiating at ∼11.5 ka.Although the onset of the increase in rainfall is shared by all three hydrological proxies, proxies for recharge in the Sulawesi cave system ([ 234 U/ 238 U] i , Mg/Ca) stabilize at interglacial values of ∼10 ka, whereas δ 18 O stabilizes about 2000 yr later (∼8 ka).Thus, it is likely that the increase in rainfall was sufficient to saturate the karst by ∼10 ka, but changes in deep atmospheric convection (rainfall δ 18 O) over the IPWP continued to evolve.
Comparison of the Sulawesi stalagmite δ 13 C and the three stalagmite rainfall proxies ([ 234 U/ 238 U] i , Mg/Ca, δ 18 O) provides information about the potential relationship between rainfall and stalagmite δ 13 C via the influence of rainfall on vegetation (Fig. 3).The Sulawesi δ 13 C record shows little similarity in largescale trends with the Sulawesi stalagmite rainfall proxies.The initiation of the deglacial transition in δ 13 C (∼17.5 ka) leads the shift in rainfall (∼11.5 ka), on average, by ∼6 ka.The early deglacial decrease in stalagmite δ 13 C, therefore, cannot be driven by the deglacial increase in rainfall.

Stalagmite δ 13 C as a proxy for vegetation productivity
Previous studies have shown that large shifts in stalagmite δ 13 C, as observed in the Sulawesi record, can be produced by natural and anthropogenic changes in tropical vegetation cover (e.g., Cruz et al., 2006;Griffiths et al., 2013;Hartmann et al., 2013;Burns et al., 2016;Fohlmeister et al., 2020).In tropical landscapes dominated by C 3 plants, the δ 13 C of dissolved inorganic carbon (DIC) in cave drip waters is primarily set by the mass balance of carbon derived from plant root-respired CO 2 and soil microbial activity (∼80-90% of the carbon) and carbonate from limestone dissolution (Vogel and Kronfeld, 1997;Genty et al., 2001;Hou et al., 2003;Griffiths et al., 2012;Meyer et al., 2014;Wong and Breecker 2015;Burns et al., 2016).The δ 13 C value of DIC (and speleothems) in such settings is generally around −8 to −12‰ (McDermott, 2004;Fairchild et al., 2006).By contrast, in the absence of vegetation, the δ 13 C of drip water would reflect the mixing of carbon from atmospheric CO 2 (e.g., around −6 to −7‰ at the LGM) and local bedrock (∼ +1‰), with stalagmite δ 13 C values approaching ∼0‰.The large isotopic contrast between the two endmember mixing scenarios provides considerable scope for changes in tropical vegetation productivity to alter stalagmite δ 13 C.
It is likely, therefore, that most of the Sulawesi δ 13 C signal reflects changes in temperature and atmospheric CO 2 concentrations, through their combined influence on vegetation type, plant root respiration, and soil microbial activity over Gempa Bumi Cave (e.g., Cosford et al., 2009;Fohlmeister et al., 2020;Lechleitner et al., 2021).Temperature and CO 2 covary on glacial −interglacial timescales (e.g., Petit et al., 1999;NGRIP, 2004;EPICA, 2006), and their individual effects on vegetation productivity (and stalagmite δ 13 C) are not easily separated.However, model studies designed to look at the relative influence of temperature and CO 2 show a 30% reduction in the net primary productivity of tropical forests at the LGM, compared with a 10% reduction when only temperature was changed (Harrison and Prentice, 2003).Other studies lend support to CO 2 as the dominant determinant of vegetation productivity in the tropics (Bennett and Willis, 2000;Bragg et al., 2013;Claussen et al., 2013;Zhu et al., 2016;Chen et al., 2019), particularly during the LGM, when atmospheric CO 2 is relatively low (Cowling and Field, 2003).
Comparison of Sulawesi stalagmite δ 13 C with leaf wax δ 13 C records from Lake Towuti (Russell et al., 2014), Lake Matano (Wicaksono et al., 2015), and Mandar Bay (Wicaksono et al., 2017) spanning the last glacial period reveals a similar deglacial transition towards lower values from ∼17 to 11.3 ka (Fig. 4).The proximity of these sites to the Gempa Bumi Cave stalagmite locality is shown in Figure 1.Leaf wax δ 13 C corresponds with the relative abundance of C 3 :C 4 plants and/or changes in water and carbon use efficiency by C 3 plants, often related to factors such as soil moisture, precipitation, temperature, and humidity (Diefendorf et al., 2010).The similarity of the Gempa Bumi stalagmite δ 13 C and leaf wax records from Sulawesi lakes supports a broad shift in vegetation productivity and/or type over the deglacial transition.However, the multiproxy record of glacialinterglacial rainfall at Gempa Bumi Cave does not correspond with Sulawesi stalagmite δ 13 C, indicating that vegetation changes above the cave site are less sensitive to rainfall.On the other hand, the Sulawesi stalagmite δ 13 C and Borneo cave temperature record (Løland et al., 2022) show similar timing and trends across the deglacial period, supporting a link between increased vegetation productivity and increasing temperature (Fig. 4).This link between vegetation and temperature is at odds with the interpretation of leaf wax δ 13 C from Sulawesi, wherein the authors attribute changes in local rainfall as the main driver influencing vegetation type (Russell et al., 2014;Wicaksono et al., 2015).Thus, it is possible that heterogeneity in Sulawesi hydroclimate is driving these differences, or a combination of factors, including temperature, are influencing vegetation type near the Sulawesi lake regions.
Agreement between Sulawesi δ 13 C, regional sea-surface temperatures (SSTs), global temperature, and atmospheric CO 2 over the last 40 ka supports our interpretation that δ 13 C is recording changes in vegetation productivity, driven primarily by temperature and CO 2 (Fig. 5).SSTs calculated from Globigerinoides ruber Mg/Ca ratios in a composite of cores from the western IPWP show a 3-4°C cooling during the LGM relative to the Holocene (Linsley et al., 2010).SSTs then rise concurrently with atmospheric CO 2 during the last deglaciation, starting at ∼18.5-17.5 ka (Lea et al., 2000;Stott et al., 2002;Visser et al., 2003;Linsley et al., 2010), completing the transition by ∼11.5 ka.The timing of the late-glacial and deglacial trends in SST and atmospheric  230 Th-dated temperature record (with 2 SEM) for Gunung Mulu Cave, northern Borneo, corrected for the effect of changing elevation due to rising sea level (Løland et al., 2022).(C) Leaf wax δ 13 C records for Lake Matano (Wicaksono et al., 2015), Lake Towuti (Russell et al., 2014), and Mandar Bay (Wicaksono et al., 2017).The figure is adapted from Wicaksono et al. (2017).Leaf wax δ 13 C corresponds with the relative abundance of CO 2 is mirrored in the Gempa Bumi Cave stalagmite δ 13 C record (Fig. 5).The three records diverge during the Holocene, suggesting that neither temperature nor CO 2 is the dominant driver of Sulawesi δ 13 C at this time.
Pollen records from marine sediment cores around Sulawesi provide a basis for evaluating the potential influence of shifts in C 3 :C 4 vegetation cover on the Gempa Bumi Cave δ 13 C record.The pollen assemblages in some sediment cores throughout the IPWP region suggest that C 4 grasslands became more common at the LGM (Hope, 2001;Bird et al., 2005;Russell et al., 2014;Wicaksono et al., 2017).However, analysis of lignin phenol ratios in a sediment core from the Makassar Strait (immediately to the west of Sulawesi) recorded no major vegetation change during the LGM (Visser et al., 2004).Thus, we cannot rule out the possibility that the balance between C 3 :C 4 vegetation types varied substantially throughout the IPWP region over the last glacial cycle.
In summary, the agreement between Sulawesi δ 13 C, regional SSTs, temperature, and atmospheric CO 2 supports our conclusion that δ 13 C is recording changes in vegetation productivity.Through this mechanism, we explore the use of the Sulawesi stalagmite δ 13 C record as a proxy for regional vegetation change and, in turn, methane emissions from tropical wetlands during cooler glacial and deglacial times.

An indicator of tropical sources of glacial atmospheric methane?
Our finding of a link between Sulawesi stalagmite δ 13 C and climate conditions driving vegetation productivity above the cave system would also affect regional terrestrial sources of atmospheric methane.For example, as rising CO 2 and temperature increase vegetation productivity above the cave, warmer conditions in the tropics may also enhance biochemical processes in wetlands (Salimi et al., 2021), prompting an increase in methane emissions (Cao et al., 1996;Kleinen et al., 2020).Thus, while stalagmite δ 13 C does not record a direct relationship with atmospheric methane concentrations, it can be seen as an indicator of when conditions in this tropical region are suitable for methane production.
Sulawesi δ 13 C (vegetation productivity) shows good correspondence with EPICA ice core methane (Loulergue et al., 2008) during the glacial period, particularly from 40 to 25 ka (Fig. 5).During glacial times, large areas of northern boreal wetlands were impacted by ice-sheet growth and permafrost, reducing their methane output (Kaplan et al., 2006), while tropical sources remained a dominant source.The transition to minimum productivity in the Sulawesi record initiates around 19 ka, and recovery begins alongside initial atmospheric CO 2 and temperature rise at 17.5 ka, marking deglacial onset in Sulawesi vegetation.The highest δ 13 C value in the Sulawesi stalagmite record (minimum vegetation productivity) occurs at 17.7 ka, just before the onset of Heinrich stadial 1 (HS1).Like atmospheric CO 2 , Sulawesi vegetation productivity continues to rise throughout the deglaciation, leveling out at ∼14.7 ka during the Bølling-Allerød (B-A) (Kienast et al., 2003;Weaver et al., 2003;Rosen et al., 2014) and at ∼11.5 ka during the Younger Dryas (YD) (Fairbanks, 1989;McManus et al., 2004;Cheng et al., 2020) before continuing its deglacial rise.The absence of a substantial decrease in vegetation during the YD is a marked difference between the Sulawesi record and global methane, suggesting that this cold event had little impact on Sulawesi.The largest increases in Sulawesi vegetation productivity (lower δ 13 C) occur at the end of HS1 and the YD and correspond with times of abrupt increases in atmospheric CO 2 and methane.Previous studies have suggested that the rapid shifts in global methane are driven by tropical wetlands (Schaefer et al., 2006;Rosen et al., 2014).Thus, the tropics may be a key contributor to the global methane budget during times of increasing CO 2 and/or large-scale heat exchange across hemispheres.
The agreement between the Sulawesi stalagmite δ 13 C and ice core methane becomes decoupled after 10 ka, when stalagmite δ 13 C increases in a V-like pattern.It is possible that changes in boreal methane emissions during the Early to Middle Holocene counteract tropical methane emission variability, resulting in a muted global methane signal that is decoupled from the Sulawesi stalagmite δ 13 C.The disconnect also corresponds with the reestablishment of the Indo-Australian summer monsoon and attainment of interglacial temperatures that could prompt a shift in Sulawesi vegetation sensitivity.For example, vegetation above the cave may become more nutrient limited when temperature, CO 2 , and moisture are readily available (Cowling and Field, 2003).Strengthening of the summer monsoon and  (Jouzel et al., 2007).(D) Composite Antarctic ice core CO 2 concentrations (Bereiter et al., 2015 and references therein).(E) Antarctic ice core CH 4 concentrations (Loulergue et al., 2008).Ice core records are plotted on the AICC2012 chronology (Bazin et al., 2013).Heinrich stadial 1 (HS1), Bølling-Allerød (B-A), and Younger Dryas (YD) are indicated by shaded bars.The close association between Sulawesi δ 13 C, regional SSTs and air temperature, and atmospheric CO 2 , particularly during abrupt deglacial climate events, supports the interpretation that Sulawesi δ 13 C is recording changes in vegetation and soil productivity, driven by changes in temperature and CO 2 .

Comparison with global vegetation model simulations
To put these findings into broader context, we investigate the relationship between atmospheric methane concentrations and Sulawesi stalagmite δ 13 C over the last 40 ka using model outputs from the SDGVM.Methane sources are divided into three categories: tropics (±30°), boreal (≥35°N), and other (≥30°S and 30-35°N).These definitions are consistent with the convention established by the WETCHIMP project (Melton et al., 2013).
Model results for methane emissions from the three source areas using a modern-day land mask verses a dynamic land-sea mask, which includes exposure of shallow continental shelves, are shown in Figure 6.For comparison, simulated global emissions are also shown.The simulated methane emissions that account for exposure of shallow continental shelves show almost no effect on Sulawesi emissions.However, shallow landmass exposure for all of the tropics results in a 16% increase in total tropical methane flux from 40 to 10 ka; most of this increase (12%) is from exposed landmasses in Indonesia.The modern landmass configuration and shelf exposure scenarios both show lower LGM methane emissions compared with preindustrial for Sulawesi, the tropics, and global regions.However, inclusion of the exposed shelves produces drastically different emissions for the whole of Indonesia, with emission levels equal to or higher than preindustrial throughout the glacial, in broad agreement with other studies using different models (e.g., Kaplan, 2002;Kleinen et al., 2020).This is largely due to the major increase in the maritime continent landmass, which, in the model, is ∼95% more expansive during the LGM compared with the modern.Additionally, the simulated vegetation type over the maritime continent landmass is dominated by evergreen broadleaf trees, which is likely an overestimate, given the marine and terrestrial proxy data (e.g., Wurster et al., 2019;Nguyen et al., 2022;Cheng et al., 2023).This study, however, investigates the Sulawesi stalagmite δ 13 C record as a possible indicator of local and regional methane emissions via the response of vegetation productivity to climate and environmental conditions.Therefore, because this work is not comparing Sulawesi vegetation to emissions resulting from exposed maritime continental shelves, we have elected to perform the following analyses using modern landmass configuration.
Simulated methane emissions from the tropics remain relatively high throughout the last 40 ka, with only a small reduction in total emissions, likely due in part to the relatively small 3-4°C cooling of the tropics during the LGM (Lea et al., 2000;Linsley et al., 2010;Gagan et al., 2004;Løland et al., 2022;Fig. 7).Methane emissions from boreal sources, however, decrease dramatically during the LGM because of much lower temperatures throughout most of the year.During the LGM (26 to 20 ka), the tropics account for ∼70% of total emissions, compared with ∼20% from boreal sources (Fig. 7).During the Holocene (10 to 0 ka), their relative contributions converge, with the tropics contributing on average ∼50% of total methane emissions, compared with ∼45% from boreal sources, in line with modern observations (Aselmann and Crutzen, 1989;Cao et al., 1996;Guo et al., 2012).The relative source changes simulated by the SDGVM agree well with previous studies (Chappellaz et al., 1997;Dällenbach et al., 2000;Valdes et al., 2005;Kaplan et al., 2006;Fischer et al., 2008;Hopcroft et al., 2017;Kleinen et al., 2020).
To compare the Sulawesi stalagmite δ 13 C record with the SDGVM model output, we identify soil respiration as the model parameter closest to stalagmite δ 13 C and use this parameter as a proxy for our record within the model.Soil respiration is the emission of CO 2 from the soil surface (Schlesinger and Andrews, 2000) that is produced within the soil profile by roots and soil organisms (Raich and Schlesinger, 1992).The predominant climatic driver of soil respiration rates is debated, but it is generally agreed that temperature, CO 2 , and soil moisture all play important roles in driving soil respiration rates (Raich and Schlesinger, 1992;Bragg et al., 2013;Hursh et al., 2017), with seasonality and forest structure also exerting control (Vargas-Terminel et al., 2022).It also has been found that wetland drying significantly increases the temperature sensitivity of soil respiration rates (Chen et al., 2018).
Soil respiration acts as an indicator of vegetation productivity, as increased vegetation growth leads to an increase in organic material available to decomposers (Schlesinger and Andrews, 2000), and within the SDGVM, it correlates strongly with net primary productivity (r = 0.98).The rate of soil respiration sets the concentration of CO 2 within the soil profile (Raich and Schlesinger, 1992), which is the most likely primary source for carbon in the Sulawesi stalagmites.Therefore, we use soil respiration as a qualitative proxy for stalagmite δ 13 C within the SDGVM, noting that further work is needed to identify the processes underlying this link, for example, isotope-enabled wetland modelling.In the model, soil respiration in Indonesia responds strongly to the changing atmospheric CO 2 concentration during and since the glacial period.Increasing atmospheric CO 2 (and its fertilising influence on vegetation) accounts for half of the total LGM to preindustrial amplitude increase in soil respiration.Thus, atmospheric CO 2 is a primary driver of vegetation productivity for modelled soil respiration rates throughout the LGM.This is consistent with the underlying hypothesis for atmospheric CO 2 and temperature as external factors driving Sulawesi stalagmite δ 13 C.
The relationship between Sulawesi stalagmite δ 13 C and modelled soil respiration for different regions (e.g., Sulawesi, Indonesia, tropics) is tested using correlation analysis on the time series of mean simulated soil respiration rates and the Sulawesi δ 13 C record (Figure 8).Stalagmite δ 13 C correlates strongly with soil respiration across all three regions (Sulawesi r = −0.87,Indonesia r = −0.88,tropics r = −0.88;P < 0.001 in all cases).When the Holocene (10-0 ka) is excluded, correlations for the glacial and deglacial periods rise (Sulawesi r = −0.94,A-C) Time series of modelled mean soil respiration for the grid points corresponding to Sulawesi, Indonesia, and tropics (±30°).Sulawesi stalagmite δ 13 C is plotted on each graph for reference (the bold green curve has been resampled to match the 1 ka model resolution).(D-F) Relationships between modelled mean soil respiration and stalagmite δ 13 C, with regression statistics.Results for the glacial and deglacial periods only (40-10 ka) are in red; those for the full record (40-0 ka) are in grey.
Indonesia r = −0.93,tropics r = −0.92;P < 0.001 in all cases).These correlations support the link between speleothem δ 13 C and the modelled changes in vegetation productivity and soil respiration across the last 40 ka.Additionally, the strong agreement in soil respiration trends across local, regional, and latitudinal scales suggests that vegetation across the tropics may have varied coherently over the last 40 ka.In sum, the close agreement between modelled soil respiration and stalagmite δ 13 C suggests it is possible that Sulawesi stalagmite carbon isotopes are being driven by changes in vegetation productivity above the cave.
To explore the potential of the Sulawesi stalagmite δ 13 C as a reliable indicator of local-to-regional methane emissions, we examine the correlation between Sulawesi stalagmite δ 13 C and the total modelled methane emissions for each of the three regions (Sulawesi, Indonesia, tropics; Fig. 9).To do this, modelled time series of total methane emissions for the three regions were regressed against the Sulawesi δ 13 C time series.The timing of deglacial increases in methane emissions across all three regions coincides with Sulawesi stalagmite δ 13 C (Fig. 9).Each time series is correlated with stalagmite δ 13 C, with total methane emissions from Sulawesi and the tropics showing the strongest correlations (r = −0.88 and r = −0.87,respectively; P < 0.001).When the Holocene is excluded, the correlation with the Sulawesi grid box increases to −0.93 (P < 0.001).
Interestingly, the V-like feature during the Mid-Holocene in the Sulawesi stalagmite δ 13 C time series is not evident in the simulated total methane emissions time series for Sulawesi or Indonesia (Fig. 9).The data-model mismatch indicates that the reduction in vegetation productivity in Sulawesi is due to factors not represented in the model.The more subtle V-like feature in the modelled methane emissions time series for the tropics as a whole appears to have been driven by changes in methane emissions beyond Indonesia.Singarayer et al. (2011) and Burns (2011) found that precession-induced modification of seasonal precipitation in the Late Holocene and associated increases in modelled methane emissions from the Southern Hemisphere tropics can explain much of the Late Holocene trend in tropical methane.The V-like pattern in the Sulawesi stalagmite δ 13 C record appears to support this.For example, increased convective rainfall in the Holocene is supported by Sulawesi stalagmite δ 18 O (Krause et al., 2019), [ 234 U/ 238 U] i , and Mg/Ca, and by other IPWP records (e.g., Partin et al., 2007;Griffiths et al., 2009;Ayliffe et al., 2013;Scroxton et al., 2022).The disconnect between stalagmite δ 13 C V-like pattern and Holocene temperature-atmospheric CO 2 (see Fig. 5) coincides with increased rainfall in Sulawesi after ∼10 ka.It is possible that vegetation productivity becomes more sensitive to seasonal rainfall and/or nutrient availability during this time.
Sulawesi stalagmite δ 13 C and simulated tropical methane emissions share a similar general trend over the last 40 ka.When compared with methane measured from the EPICA ice core, stalagmite δ 13 C and simulated tropical methane correspond well over the glacial period (Fig. 10).Departures of ice core methane from simulated tropical methane and the Sulawesi δ 13 C record likely reflect major changes in boreal methane sources at higher latitudes and/or changes in other regions of the tropics.The deglacial increases in atmospheric methane measured in the EPICA ice core (at the end of HS1 and YD) coincide with negative shifts in stalagmite δ 13 C (Fig. 10).The plateau in stalagmite δ 13 C at ∼14-12 ka, during the B-A, is mirrored in the model.Because the SDGVM is only forced by climate changes every 1 ka, it does not include millennial-scale variability (Singarayer et al., 2011); thus, the step change in the deglacial pattern in the model is likely occurring due to step changes in the corresponding atmospheric CO 2 supplied to the model (Singarayer and Valdes, 2010).(A-C) Time series of modelled methane emissions totals for Sulawesi, Indonesia, and tropics (±30°).Sulawesi stalagmite δ 13 C is plotted on each graph for reference (the bold green curve has been resampled to match the 1 ka model resolution).(D-F) Relationships between modelled total methane emissions and stalagmite δ 13 C, with regression statistics.Results for the glacial and deglacial periods only (40-10 ka) are in red; those for the full record (40-0 ka) are in grey.

CONCLUSIONS
The new stalagmite δ 13 C record from Sulawesi is interpreted as a record of changing soil respiration rates through the past 40,000 yr.We explore a link to the natural methane cycle using a series of global climate and biogeochemical model simulations.These simulations show that soil respiration in Indonesia was predominantly controlled by vegetation productivity, primarily through the influence of atmospheric CO 2 and temperature.This soil respiration signature was, in turn, recorded by stalagmite δ 13 C via seepage waters, which retain the carbon isotope signature of the plant matter and soil CO 2 above the cave.
Previous work has identified the tropics as a likely source of methane emissions during the last glacial period (e.g., Brook et al., 2000;Fischer et al., 2008;Weber et al., 2010;Baumgartner et al., 2012;Guo et al., 2012;Rhodes et al., 2015Rhodes et al., , 2017;;Kleinen et al., 2020).In the SDGVM model simulations, tropical wetland methane emissions are largely controlled by changing soil respiration rates, raising the possibility that the Sulawesi stalagmite δ 13 C record indirectly reflects methane emissions related to vegetation productivity.A similar pattern in modelled soil respiration rates emerges across the whole tropics, suggesting that inferences drawn from Sulawesi may be applicable across the broader tropics.However, this is contingent on the spatial expression of the glacial-interglacial climate transition in the climate model.The good agreement between the stalagmite δ 13 C record and SDGVM output indicates that tropical vegetation productivity, and hence organic matter decomposition and methanogenesis, were active during the glacial period despite moderate decreases in temperature and precipitation.Our findings support the predominance of tropical sources of methane emissions during the glacial period when boreal sources were mostly dormant.
The likely relationship between Sulawesi δ 13 C and ice core methane is masked during the Holocene, when boreal wetland methane emissions become more influential in the atmospheric methane budget.However, the model results and stalagmite δ 13 C show some evidence for tropical methane sources contributing to Late Holocene methane variability.A disconnect between stalagmite δ 13 C, temperature, global atmospheric CO 2 , and methane emissions coincides with increased rainfall in Sulawesi after ∼10 ka.It is possible that vegetation productivity becomes more sensitive to seasonal rainfall and/or nutrient availability during this time.
We have established Sulawesi stalagmite δ 13 C as a proxy for changes in vegetation productivity via soil respiration, which, in the model examined, is also strongly related to changes in tropical methane production.These changes in tropical methane production appear to have made a substantial contribution to the glacial atmospheric methane budget.Sulawesi stalagmite δ 13 C may therefore provide an indirect tropical proxy of glacial methane emissions, offering a unique non-polar constraint on the likely sources of past atmospheric methane.Data Availability Statement.The supplementary material for this article can be found on the NOAA Paleoclimate Data repository: https://www.ncei.noaa.gov/access/paleo-search/study/38940.

Figure 2 .
Figure 2. Stalagmites GB09-3 and GB11-9 with age-depth models.Photographs of (A) GB09-3 and (B) GB11-9 show sampling tracks used for stable isotope analysis.Coloured dots indicate the locations of 230 Th dates, expressed as ka (where present is defined as 1950 CE and errors are 2σ).Two dates shown in grey for GB09-3 were not used in the final age model.(C) Age-depth models for each stalagmite with 2σ age uncertainties on 230 Th dates.All ages are in stratigraphic sequence, within error.Details of the 230 Th age data are given in supplementary table 1 in Krause et al.(2019).The average growth rates are 1.74 mm/100 yr for GB09-3 and 1.40 mm/ 100 yr for GB11-9 (for 40-26 ka), with no detectable hiatuses.Data from the bottom of GB11-9 are not included in this study.

Figure 3 .
Figure 3. Stalagmite δ 13 C, δ 18 O, initial 234 U/ 238 U, and Mg/Ca records for Sulawesi over the last 40 ka.(A) δ 13 C for GB09-3 and GB11-9 corrected for the effect of atmospheric CO 2 on carbon isotope fractionation in C 3 plants(Breecker, 2017).Uncorrected δ 13 C is shown in grey.The large deglacial δ 13 C transition (green shading) encompasses two abrupt negative excursions at ∼14.7-14.5 ka and 11.7-11.6ka that mark the terminations of Heinrich stadial 1 (HS1) and the Younger Dryas (YD), respectively.The Bølling-Allerød (B-A) is also shown (yellow).(B) δ 18 O for GB09-3 and GB11-9 corrected for the effect of ice volume(Krause et al., 2019).Uncorrected δ 18 O shown in grey.(C) Initial 234 U/ 238 U records for GB09-3 and GB11-9.(D) Mg/Ca record for GB09-3.The late-deglacial transition to lower values in all three hydroclimate proxies is interpreted as an increase in rainfall amount and a strengthened Indo-Australian summer monsoon (IASM).Initial 234 U/ 238 U is influenced by drip-water flow pathways; thus, coeval stalagmites are unlikely to share the same values and are therefore plotted on separate scales. 230Th dates with 2σ errors are shown at the top of the figure.PCP, prior calcite precipitation.
Figure 3. Stalagmite δ 13 C, δ 18 O, initial 234 U/ 238 U, and Mg/Ca records for Sulawesi over the last 40 ka.(A) δ 13 C for GB09-3 and GB11-9 corrected for the effect of atmospheric CO 2 on carbon isotope fractionation in C 3 plants(Breecker, 2017).Uncorrected δ 13 C is shown in grey.The large deglacial δ 13 C transition (green shading) encompasses two abrupt negative excursions at ∼14.7-14.5 ka and 11.7-11.6ka that mark the terminations of Heinrich stadial 1 (HS1) and the Younger Dryas (YD), respectively.The Bølling-Allerød (B-A) is also shown (yellow).(B) δ 18 O for GB09-3 and GB11-9 corrected for the effect of ice volume(Krause et al., 2019).Uncorrected δ 18 O shown in grey.(C) Initial 234 U/ 238 U records for GB09-3 and GB11-9.(D) Mg/Ca record for GB09-3.The late-deglacial transition to lower values in all three hydroclimate proxies is interpreted as an increase in rainfall amount and a strengthened Indo-Australian summer monsoon (IASM).Initial 234 U/ 238 U is influenced by drip-water flow pathways; thus, coeval stalagmites are unlikely to share the same values and are therefore plotted on separate scales. 230Th dates with 2σ errors are shown at the top of the figure.PCP, prior calcite precipitation.

Figure 6 .
Figure 6.Influence of shallow landmass exposure on total methane emissions in the Sheffield Dynamic Global Vegetation Model (SDGVM).Modelled total methane emissions from present-day land areas (black) for (A) Sulawesi, (B) Indonesia, (C) Tropics (±30°), and (D) global.Red curves show increases in emissions due to exposure of new land at times of lowered sea levels.Although the amount of methane emitted increases with landmass exposure, the patterns of emissions during glacial times remain relatively constant.

Figure 7 .
Figure 7. Glacial-interglacial evolution of tropical and higher-latitude methane sources from the Sheffield Dynamic Global Vegetation Model (SDGVM).(A) Map showing the spatial distribution of regions used in this study: tropics (green), boreal (blue), and other (grey).Inset shows Indonesia (red box) and Sulawesi (pink grid cell) as represented for the present day in the SDGVM.(B) Total methane emissions by region.(C) Stacked regional emissions showing the relative contribution to the global total.(D) Regional emissions as a percentage of total emissions.

Figure 8 .
Figure 8.Comparison of modelled mean soil respiration and Sulawesi stalagmite δ 13 C. (A-C) Time series of modelled mean soil respiration for the grid points corresponding to Sulawesi, Indonesia, and tropics (±30°).Sulawesi stalagmite δ 13 C is plotted on each graph for reference (the bold green curve has been resampled to match the 1 ka model resolution).(D-F) Relationships between modelled mean soil respiration and stalagmite δ 13 C, with regression statistics.Results for the glacial and deglacial periods only (40-10 ka) are in red; those for the full record (40-0 ka) are in grey.

Figure 9 .
Figure 9.Comparison of modelled total methane emissions and Sulawesi stalagmite δ 13 C.(A-C) Time series of modelled methane emissions totals for Sulawesi, Indonesia, and tropics (±30°).Sulawesi stalagmite δ 13 C is plotted on each graph for reference (the bold green curve has been resampled to match the 1 ka model resolution).(D-F) Relationships between modelled total methane emissions and stalagmite δ 13 C, with regression statistics.Results for the glacial and deglacial periods only (40-10 ka) are in red; those for the full record (40-0 ka) are in grey.

Figure 10 .
Figure 10.Sulawesi δ 13 C as a potential indicator of the contribution of tropical methane to global atmospheric methane.Comparison of Sulawesi stalagmite δ 13 C, ice core methane concentrations (plotted on the AICC2012 chronology; Bazin et al., 2013) and modelled total methane emissions for the tropics.The Sulawesi δ 13 C values and modelled methane emissions are approximately scaled to the glacial section of the ice core methane record to reflect the tropical contribution to global methane.Deviations between these records likely reflect major changes in boreal methane sources at higher latitudes and/or variations in other parts of the tropics.Heinrich stadial 1 (HS1), Bølling-Allerød (B-A), and the Younger Dryas (YD) are indicated by shaded vertical bars.