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 δ¹³C to record a range of environmental processes. Here, we report on ²³⁰Th-dated stalagmite δ¹³C 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 δ¹³C 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 CO₂, in line with our interpretation of stalagmite δ¹³C. In turn, modelled methane emissions are driven by soil respiration, providing a mechanism that links methane to stalagmite δ¹³C. 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 δ¹³C 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., 1996, 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).

It is generally agreed that wetland methane emissions were the most important contributor to the atmospheric methane budget in the past (e.g., Brook et al., 2000; Valdes et al., 2005; Kaplan et al., 2006; Fischer et al., 2008; Korhola et al., 2010; Weber et al., 2010; Baumgartner et al., 2012; Guo et al., 2012; Ringeval et al., 2013; Rhodes et al., 2015, 2017; Hopcroft et al., 2017, 2018, 2020; Kleinen et al., 2020, 2023; Dyonisius et al., 2020). The modern global methane cycle is dominated by methane from natural wetlands, accounting for ~60–80% of natural methane emissions (Kirschke et al., 2013; Rosentreter et al., 2021). About 60% of modern wetland methane emissions are from tropical sources, with ~40% from boreal sources (Aselmann and Crutzen, 1989; Cao et al., 1996; Guo et al., 2012).

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., 2020, 2023).

At present, there is no proxy for past tropical wetland methane emissions. It has been suggested, however, that carbon isotope ratios (δ¹³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 δ¹³C primarily reflects C₃ vegetation productivity, with most of the carbon in infiltrating waters originating from CO₂ 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₂ and heavier carbon derived from bedrock combine to influence the δ¹³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 δ¹³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 δ¹³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 δ¹³C record, afforded by precise ²³⁰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.

Figure 1. Map of the study region. Star indicates location of Gempa Bumi Cave, Sulawesi (5°S, 120°E, ~140 m above sea level). Locations of other paleoclimate reconstructions referenced in this study include: marine sediment cores (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., 2015, 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).

MATERIALS AND METHODS

Stalagmites

We present new δ¹³C, ²³⁴U/²³⁸U, and Mg/Ca records for two ²³⁰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 δ¹⁸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 δ¹⁸O and δ¹³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).

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 ²³⁰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 ²³⁰Th dates. All ages are in stratigraphic sequence, within error. Details of the ²³⁰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.

²³⁰Th/²³⁴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 [²³⁰Th/²³²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 ²³⁴U/²³⁸U values ([²³⁴U/²³⁸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.2 mm and 0.7 mm, respectively, equating to an average sample resolution of ~55 yr. Measurements of δ¹⁸O and δ¹³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₂ was liberated from the carbonate by reaction with 105% H₃PO₄ under vacuum at 90°C. Measurements of δ¹⁸O and δ¹³C were corrected using the NBS19 (δ¹⁸O = −2.20‰, δ¹³C = 1.95‰) and NBS18 (δ¹⁸O = −23.0‰, δ¹³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 δ¹⁸O and ±0.02‰ for δ¹³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 δ¹⁸O) were reanalysed to minimise errors related to the mass spectrometry. Additionally, samples giving δ¹³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 δ¹³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 δ¹³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 ¹³CO₂, Mg²⁺, and Sr²⁺ (relative to Ca²⁺) in drip waters and raises stalagmite δ¹³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₃. Analytical precision was determined by repeat analyses of an in-house laboratory (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 GB09-3 (n = 192) 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.

Model simulations

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.1D and 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, 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, 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, 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₂ concentrations. SDGVM accounts for the main factors driving vegetation productivity, including climate (surface temperature, precipitation, relative humidity), atmospheric CO₂ 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₂, 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 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₄ 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₂ 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₂ (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 multi-model mean of −2.5 ± 21% (Melton et al., 2013). The tropics proved most sensitive to an increase in CO₂ 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² (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., 2010, 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 δ¹³C records are generally in good agreement across their interval of overlap (40–26 ka). The millennial-scale δ¹³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 with small ranges in δ¹³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).

Figure 3. Stalagmite δ¹³C, δ¹⁸O, initial ²³⁴U/²³⁸U, and Mg/Ca records for Sulawesi over the last 40 ka. (A) δ¹³C for GB09-3 and GB11-9 corrected for the effect of atmospheric CO₂ on carbon isotope fractionation in C₃ plants (Breecker, 2017). Uncorrected δ¹³C is shown in grey. The large deglacial δ¹³C transition (green shading) encompasses two abrupt negative excursions at ~14.7–14.5 ka and 11.7–11.6 ka 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) δ¹⁸O for GB09-3 and GB11-9 corrected for the effect of ice volume (Krause et al., 2019). Uncorrected δ¹⁸O shown in grey. (C) Initial ²³⁴U/²³⁸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 ²³⁴U/²³⁸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. ²³⁰Th dates with 2σ errors are shown at the top of the figure. PCP, prior calcite precipitation.

It is important to note that the δ¹³C records have been corrected for the effect of atmospheric pCO₂ on the δ¹³C of C₃ plants (Schubert and Jahren, 2012). The transfer of this effect on carbon isotope fractionation in C₃ plants above a cave to stalagmites growing within the cave was identified by Breecker (2017) in a study assessing globally averaged speleothem δ¹³C records over the past 90 ka. They found that, after accounting for other processes, the effect of atmospheric CO₂ is best explained by a C₃ plant δ¹³C sensitivity of −1.6‰ for every 100 parts per million by volume (ppmv) increase in pCO₂ from the LGM to the Holocene. Therefore, it is important to correct for the change in stalagmite δ¹³C that occurs as a result of glacial–interglacial atmospheric pCO₂ before investigating δ¹³C as a recorder of glacial–interglacial vegetation productivity. The Sulawesi stalagmite δ¹³C values were adjusted by −1.6‰/100 ppmv (Breecker, 2017) relative to modern atmospheric pCO₂ (190 ppmv) using the Antarctic ice core composite pCO₂ 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 δ¹³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 δ¹³C values. The deglacial interval contains a prominent ~4.2‰ shift in corrected δ¹³C from the maximum δ¹³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 δ¹³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 δ¹³C.

The Holocene section of the record shows a surprisingly high degree of δ¹³C variability, most notably a prominent V-like pattern in the Early to Middle Holocene. During this time, the δ¹³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: δ¹⁸O, initial uranium isotope activity [²³⁴U/²³⁸U]_i), and Mg/Ca. The Sulawesi stalagmite δ¹⁸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 δ¹⁸O values, occurs at ~11.5 ka.

[²³⁴U/²³⁸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., 2000, 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, [²³⁴U/²³⁸U]_i is expected to reflect the activity ratio in seepage waters forming speleothems. The [²³⁴U/²³⁸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, [²³⁴U/²³⁸U]_i can increase as a result of preferential leaching of ²³⁴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, [²³⁴U/²³⁸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₂ 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 [²³⁴U/²³⁸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 δ¹⁸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 (σ²). 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 (σ²). The stabilization of Mg/Ca at lower values following the deglacial transition, corroborated by [²³⁴U/²³⁸U]_i, is interpreted to reflect wetter conditions.

Previous studies have shown that similar decreases in stalagmite δ¹⁸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 δ¹⁸O, [²³⁴U/²³⁸U]_i, and Mg/Ca records, alongside other regional rainfall δ¹⁸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 ([²³⁴U/²³⁸U]_i, Mg/Ca) stabilize at interglacial values of ~10 ka, whereas δ¹⁸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 δ¹⁸O) over the IPWP continued to evolve.

Comparison of the Sulawesi stalagmite δ¹³C and the three stalagmite rainfall proxies ([²³⁴U/²³⁸U]_i, Mg/Ca, δ¹⁸O) provides information about the potential relationship between rainfall and stalagmite δ¹³C via the influence of rainfall on vegetation (Fig. 3). The Sulawesi δ¹³C record shows little similarity in large-scale trends with the Sulawesi stalagmite rainfall proxies. The initiation of the deglacial transition in δ¹³C (~17.5 ka) leads the shift in rainfall (~11.5 ka), on average, by ~6 ka. The early deglacial decrease in stalagmite δ¹³C, therefore, cannot be driven by the deglacial increase in rainfall.

Stalagmite δ¹³C as a proxy for vegetation productivity

Previous studies have shown that large shifts in stalagmite δ¹³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₃ plants, the δ¹³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₂ 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 δ¹³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 δ¹³C of drip water would reflect the mixing of carbon from atmospheric CO₂ (e.g., around −6 to −7‰ at the LGM) and local bedrock (~ +1‰), with stalagmite δ¹³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 δ¹³C.

It is likely, therefore, that most of the Sulawesi δ¹³C signal reflects changes in temperature and atmospheric CO₂ 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₂ covary on glacial−interglacial timescales (e.g., Petit et al., 1999; NGRIP, 2004; EPICA, 2006), and their individual effects on vegetation productivity (and stalagmite δ¹³C) are not easily separated. However, model studies designed to look at the relative influence of temperature and CO₂ 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₂ 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₂ is relatively low (Cowling and Field, 2003).

Comparison of Sulawesi stalagmite δ¹³C with leaf wax δ¹³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 δ¹³C corresponds with the relative abundance of C₃:C₄ plants and/or changes in water and carbon use efficiency by C₃ plants, often related to factors such as soil moisture, precipitation, temperature, and humidity (Diefendorf et al., 2010). The similarity of the Gempa Bumi stalagmite δ¹³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 glacial–interglacial rainfall at Gempa Bumi Cave does not correspond with Sulawesi stalagmite δ¹³C, indicating that vegetation changes above the cave site are less sensitive to rainfall. On the other hand, the Sulawesi stalagmite δ¹³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 δ¹³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.

Figure 4. Sulawesi vegetation productivity compared with Borneo cave temperature and δ¹³C of Sulawesi leaf wax. (A) δ¹³C for stalagmite GB09-3, reflecting changes in vegetation productivity above Gempa Bumi Cave. (B) ²³⁰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 δ¹³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 δ¹³C corresponds with the relative abundance of C₃:C₄ plants and/or changes in water and carbon use efficiency by C₃ plants related to climate conditions. Heinrich stadial 1 (HS1), Bølling-Allerød (B-A), and Younger Dryas (YD) are indicated by shaded bars.

Agreement between Sulawesi δ¹³C, regional sea-surface temperatures (SSTs), global temperature, and atmospheric CO₂ over the last 40 ka supports our interpretation that δ¹³C is recording changes in vegetation productivity, driven primarily by temperature and CO₂ (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₂ 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 CO₂ is mirrored in the Gempa Bumi Cave stalagmite δ¹³C record (Fig. 5). The three records diverge during the Holocene, suggesting that neither temperature nor CO₂ is the dominant driver of Sulawesi δ¹³C at this time.

Figure 5. Relationship between Sulawesi stalagmite δ¹³C, temperature, atmospheric CO₂, and CH₄ over the last 40 ka. (A) δ¹³C for stalagmites GB09-3 and GB11-9. (B) Summer sea-surface temperature (SST) reconstruction from core MD98-2181in the northern Indo-Pacific Warm Pool (IPWP; Stott et al., 2002) and composite SST anomalies for the western IPWP (Linsley et al., 2010 and references therein). (C) Antarctic temperature inferred from ice core δD (Jouzel et al., 2007). (D) Composite Antarctic ice core CO₂ concentrations (Bereiter et al., 2015 and references therein). (E) Antarctic ice core CH₄ 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 δ¹³C, regional SSTs and air temperature, and atmospheric CO₂, particularly during abrupt deglacial climate events, supports the interpretation that Sulawesi δ¹³C is recording changes in vegetation and soil productivity, driven by changes in temperature and CO₂.

Pollen records from marine sediment cores around Sulawesi provide a basis for evaluating the potential influence of shifts in C₃:C₄ vegetation cover on the Gempa Bumi Cave δ¹³C record. The pollen assemblages in some sediment cores throughout the IPWP region suggest that C₄ 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₃:C₄ vegetation types varied substantially throughout the IPWP region over the last glacial cycle.

In summary, the agreement between Sulawesi δ¹³C, regional SSTs, temperature, and atmospheric CO₂ supports our conclusion that δ¹³C is recording changes in vegetation productivity. Through this mechanism, we explore the use of the Sulawesi stalagmite δ¹³C record as a proxy for regional vegetation change and, in turn, methane emissions from tropical wetlands during cooler glacial and deglacial times.

DISCUSSION

An indicator of tropical sources of glacial atmospheric methane?

Our finding of a link between Sulawesi stalagmite δ¹³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₂ 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 δ¹³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 δ¹³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₂ and temperature rise at 17.5 ka, marking deglacial onset in Sulawesi vegetation. The highest δ¹³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₂, 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 δ¹³C) occur at the end of HS1 and the YD and correspond with times of abrupt increases in atmospheric CO₂ 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₂ and/or large-scale heat exchange across hemispheres.

The agreement between the Sulawesi stalagmite δ¹³C and ice core methane becomes decoupled after 10 ka, when stalagmite δ¹³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 δ¹³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₂, and moisture are readily available (Cowling and Field, 2003). Strengthening of the summer monsoon and strong seasonality could also influence productivity patterns (Vargas-Terminel et al., 2022).

Comparison with global vegetation model simulations

To put these findings into broader context, we investigate the relationship between atmospheric methane concentrations and Sulawesi stalagmite δ¹³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 δ¹³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.

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.

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).

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.

To compare the Sulawesi stalagmite δ¹³C record with the SDGVM model output, we identify soil respiration as the model parameter closest to stalagmite δ¹³C and use this parameter as a proxy for our record within the model. Soil respiration is the emission of CO₂ 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₂, 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₂ 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 δ¹³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₂ concentration during and since the glacial period. Increasing atmospheric CO₂ (and its fertilising influence on vegetation) accounts for half of the total LGM to preindustrial amplitude increase in soil respiration. Thus, atmospheric CO₂ 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₂ and temperature as external factors driving Sulawesi stalagmite δ¹³C.

The relationship between Sulawesi stalagmite δ¹³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 δ¹³C record (Figure 8). Stalagmite δ¹³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, Indonesia r = −0.93, tropics r = −0.92; P < 0.001 in all cases). These correlations support the link between speleothem δ¹³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 δ¹³C suggests it is possible that Sulawesi stalagmite carbon isotopes are being driven by changes in vegetation productivity above the cave.

Figure 8. Comparison of modelled mean soil respiration and Sulawesi stalagmite δ¹³C. (A–C) Time series of modelled mean soil respiration for the grid points corresponding to Sulawesi, Indonesia, and tropics (±30°). Sulawesi stalagmite δ¹³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 δ¹³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.

To explore the potential of the Sulawesi stalagmite δ¹³C as a reliable indicator of local-to-regional methane emissions, we examine the correlation between Sulawesi stalagmite δ¹³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 δ¹³C time series. The timing of deglacial increases in methane emissions across all three regions coincides with Sulawesi stalagmite δ¹³C (Fig. 9). Each time series is correlated with stalagmite δ¹³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).

Figure 9. Comparison of modelled total methane emissions and Sulawesi stalagmite δ¹³C. (A–C) Time series of modelled methane emissions totals for Sulawesi, Indonesia, and tropics (±30°). Sulawesi stalagmite δ¹³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 δ¹³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.

Interestingly, the V-like feature during the Mid-Holocene in the Sulawesi stalagmite δ¹³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 δ¹³C record appears to support this. For example, increased convective rainfall in the Holocene is supported by Sulawesi stalagmite δ¹⁸O (Krause et al., 2019), [²³⁴U/²³⁸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 δ¹³C V-like pattern and Holocene temperature–atmospheric CO₂ (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 δ¹³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 δ¹³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 δ¹³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 δ¹³C (Fig. 10). The plateau in stalagmite δ¹³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₂ supplied to the model (Singarayer and Valdes, 2010).

Figure 10. Sulawesi δ¹³C as a potential indicator of the contribution of tropical methane to global atmospheric methane. Comparison of Sulawesi stalagmite δ¹³C, ice core methane concentrations (plotted on the AICC2012 chronology; Bazin et al., 2013) and modelled total methane emissions for the tropics. The Sulawesi δ¹³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.

CONCLUSIONS

The new stalagmite δ¹³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₂ and temperature. This soil respiration signature was, in turn, recorded by stalagmite δ¹³C via seepage waters, which retain the carbon isotope signature of the plant matter and soil CO₂ 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., 2015, 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 δ¹³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 δ¹³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 δ¹³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 δ¹³C show some evidence for tropical methane sources contributing to Late Holocene methane variability. A disconnect between stalagmite δ¹³C, temperature, global atmospheric CO₂, 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 δ¹³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 δ¹³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.