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., Reference Brook, Sowers and Orchardo1996, Reference Brook, Harder, Severinghaus, Steig and Sucher2000; Schaefer et al., Reference Schaefer, Whiticar, Brook, Petrenko, Ferretti and Severinghaus2006; Grachev et al., Reference Grachev, Brook and Severinghaus2007; Fischer et al., Reference Fischer, Behrens, Bock, Richter, Schmitt, Loulergue and Chappellaz2008; Loulergue et al., Reference Loulergue, Schilt, Spahni, Masson-Delmotte, Blunier, Lemieux, Barnola, Raynaud, Stocker and Chappellaz2008; Baumgartner et al., Reference Baumgartner, Schilt, Eicher, Schmitt, Schwander, Spahni, Fischer and Stocker2012; Möller et al., Reference Möller, Sowers, Bock, Spahni, Behrens, Schmitt, Miller and Fischer2013; Rhodes et al., Reference Rhodes, Brook, Chiang, Blunier, Maselli, McConnell, Romanini and Severinghaus2015). 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., Reference Brook, Harder, Severinghaus, Steig and Sucher2000; Fischer et al., Reference Fischer, Behrens, Bock, Richter, Schmitt, Loulergue and Chappellaz2008; Levine et al., Reference Levine, Wolff, Jones, Hutterli, Wild, Carver and Pyle2011; Hopcroft et al., Reference Hopcroft, Valdes, O'Connor, Kaplan and Beerling2017).
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., Reference Brook, Harder, Severinghaus, Steig and Sucher2000; Valdes et al., Reference Valdes, Beerling and Johnson2005; Kaplan et al., Reference Kaplan, Folberth and Hauglustaine2006; Fischer et al., Reference Fischer, Behrens, Bock, Richter, Schmitt, Loulergue and Chappellaz2008; Korhola et al., Reference Korhola, Ruppel, Seppä, Väliranta, Virtanen and Weckström2010; Weber et al., Reference Weber, Drury, Toonen and van Weele2010; Baumgartner et al., Reference Baumgartner, Schilt, Eicher, Schmitt, Schwander, Spahni, Fischer and Stocker2012; Guo et al., Reference Guo, Zhou and Wu2012; Ringeval et al., Reference Ringeval, Hopcroft, Valdes, Ciais, Ramstein, Dolman and Kageyama2013; Rhodes et al., Reference Rhodes, Brook, Chiang, Blunier, Maselli, McConnell, Romanini and Severinghaus2015, Reference Rhodes, Brook, McConnell, Blunier, Sime, Xavier and Mulvaney2017; Hopcroft et al., Reference Hopcroft, Valdes, O'Connor, Kaplan and Beerling2017, Reference Hopcroft, Valdes and Kaplan2018, Reference Hopcroft, Ramstein, Pugh, Hunter, Murguia-Flores, Quiquet, Sun, Tan and Valdes2020; Kleinen et al., Reference Kleinen, Mikolajewicz and Brovkin2020, Reference Kleinen, Gromov, Steil and Brovkin2023; Dyonisius et al., Reference Dyonisius, Petrenko, Smith, Hua, Yang, Schmitt and Beck2020). The modern global methane cycle is dominated by methane from natural wetlands, accounting for ~60–80% of natural methane emissions (Kirschke et al., Reference Kirschke, Bousquet, Ciais, Saunois, Canadell, Dlugokencky and Bergamaschi2013; Rosentreter et al., Reference Rosentreter, Borges, Deemer, Holgerson, Liu, Song and Melack2021). About 60% of modern wetland methane emissions are from tropical sources, with ~40% from boreal sources (Aselmann and Crutzen, Reference Aselmann and Crutzen1989; Cao et al., Reference Cao, Marshall and Gregson1996; Guo et al., Reference Guo, Zhou and Wu2012).
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., Reference Chappellaz, Blunier, Kints, Dällenbach, Barnola, Schwander, Raynaud and Stauffer1997; Brook et al., Reference Brook, Harder, Severinghaus, Steig and Sucher2000; Dällenbach et al., Reference Dällenbach, Blunier, Flückiger, Stauffer, Chappellaz and Raynaud2000; Baumgartner et al., Reference Baumgartner, Schilt, Eicher, Schmitt, Schwander, Spahni, Fischer and Stocker2012). 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., Reference Baumgartner, Schilt, Eicher, Schmitt, Schwander, Spahni, Fischer and Stocker2012). 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. (Reference Kleinen, Mikolajewicz and Brovkin2020) 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., Reference Valdes, Beerling and Johnson2005; Kaplan et al., Reference Kaplan, Folberth and Hauglustaine2006), and changes since the LGM are largely thought to be driven by wetlands (Kleinen et al., Reference Kleinen, Gromov, Steil and Brovkin2023). The LGM wetland reduction is supported by isotopic analyses of atmospheric methane in ice cores (Fischer et al., Reference Fischer, Behrens, Bock, Richter, Schmitt, Loulergue and Chappellaz2008; Bock et al., Reference Bock, Schmitt, Beck, Seth, Chappellaz and Fischer2017; Hopcroft et al., Reference Hopcroft, Valdes and Kaplan2018); however, attribution of methane isotopes to a specific source is not definitive (Möller et al., Reference Möller, Sowers, Bock, Spahni, Behrens, Schmitt, Miller and Fischer2013). 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., Reference Ridgwell, Maslin and Kaplan2012; Kleinen et al., Reference Kleinen, Mikolajewicz and Brovkin2020, Reference Kleinen, Gromov, Steil and Brovkin2023).
At present, there is no proxy for past tropical wetland methane emissions. It has been suggested, however, that carbon isotope ratios (δ13C) 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, Reference Burns2011; Griffiths et al., Reference Griffiths, Drysdale, Gagan, Hellstrom, Couchoud, Ayliffe, Vonhof and Hantoro2013; Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020). In wet tropical settings, speleothem δ13C primarily reflects C3 vegetation productivity, with most of the carbon in infiltrating waters originating from CO2 in the soil atmosphere, produced by vegetation root respiration and microbial activity (e.g., Genty et al., Reference Genty, Baker, Massault, Proctor, Gilmour, Pons-Branchu and Hamelin2001; Wong and Breecker, Reference Wong and Breecker2015; Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020). Light carbon from soil CO2 and heavier carbon derived from bedrock combine to influence the δ13C of dissolved inorganic carbon in cave seepage waters (Vogel and Kronfeld, Reference Vogel and Kronfeld1997; Genty et al., Reference Genty, Baker, Massault, Proctor, Gilmour, Pons-Branchu and Hamelin2001; Hou et al., Reference Hou, Tan, Cheng and Liu2003; Griffiths et al., Reference Griffiths, Fohlmeister, Drysdale, Hua, Johnson, Hellstrom, Gagan and Zhao2012; Wong and Breecker Reference Wong and Breecker2015; Burns et al., Reference Burns, Godfrey, Faina, McGee, Hardt, Ranivoharimanana and Randrianasy2016). Cave seepage waters then precipitate as speleothems, preserving a δ13C signature driven primarily by changes in vegetation productivity (Dorale and Liu, Reference Dorale and Liu2009; Breecker et al., Reference Breecker, Payne, Quade, Banner, Ball, Meyer and Cowan2012; Hartmann et al., Reference Hartmann, Eiche, Neumann, Fohlmeister, Schröder-Ritzrau, Mangini and Haryono2013; Burns et al., Reference Burns, Godfrey, Faina, McGee, Hardt, Ranivoharimanana and Randrianasy2016; Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020).
In this study, we use stalagmite δ13C 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 δ13C record, afforded by precise 230Th 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., Reference Stott, Poulsen, Lund and Thunell2002; Linsley et al., Reference Linsley, Rosenthal and Oppo2010 and references therein), cave temperature record for Gunung Mulu National Park, northern Borneo (Løland et al., Reference Løland, Krüger, Fernandez, Buckingham, Carolin, Sodemann, Adkins, Cobb and Meckler2022), and leaf wax records for Sulawesi (Russell et al., Reference Russell, Vogel, Konecky, Bijaksana, Huang, Melles, Wattrus, Costa and King2014; Wicaksono et al., Reference Wicaksono, Russell and Bijaksana2015, Reference Wicaksono, Russell, Holbourn and Kuhnt2017). 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 δ13C, 234U/238U, and Mg/Ca records for two 230Th-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 δ18O by Krause et al. (Reference Krause, Gagan, Dunbar, Hantoro, Hellstrom, Cheng, Edwards, Suwargadi, Abram and Rifai2019). 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 δ18O and δ13C 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., Reference Wang, Cheng, Edwards, An, Wu, Shen and Dorale2001; Dorale and Liu, Reference Dorale and Liu2009).
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 230Th 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 230Th dates. All ages are in stratigraphic sequence, within error. Details of the 230Th age data are given in supplementary table 1 in Krause et al. (Reference Krause, Gagan, Dunbar, Hantoro, Hellstrom, Cheng, Edwards, Suwargadi, Abram and Rifai2019). 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.
230Th/234U chronology
The chronologies for stalagmites GB09-3 and GB11-9 are based on 44 U-Th dates described in Krause et al. (Reference Krause, Gagan, Dunbar, Hantoro, Hellstrom, Cheng, Edwards, Suwargadi, Abram and Rifai2019) 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, Reference Hellstrom2003) and the University of Minnesota (Cheng et al., Reference Cheng, Edwards, Shen, Polyak, Asmerom, Woodhead and Hellstrom2013). All samples were corrected for small amounts of detrital thorium using an initial [230Th/232Th] ratio of 3.0 ± 0.75 determined by stratigraphic constraint analysis of the measured U-Th dates (Hellstrom, Reference Hellstrom2006). 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 234U/238U values ([234U/238U]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 δ18O and δ13C 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. CO2 was liberated from the carbonate by reaction with 105% H3PO4 under vacuum at 90°C. Measurements of δ18O and δ13C were corrected using the NBS19 (δ18O = −2.20‰, δ13C = 1.95‰) and NBS18 (δ18O = −23.0‰, δ13C = 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 δ18O and ±0.02‰ for δ13C (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 δ18O) were reanalysed to minimise errors related to the mass spectrometry. Additionally, samples giving δ13C 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 δ13C 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 δ13C. 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, Reference Fairchild and Treble2009). During drier conditions, PCP increases 13CO2, Mg2+, and Sr2+ (relative to Ca2+) in drip waters and raises stalagmite δ13C, Mg/Ca and Sr/Ca simultaneously (Baker et al., Reference Baker, Ito, Smart and McEwan1997).
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 (Reference Schrag1999). 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 HNO3. 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. (Reference de Villiers, Greaves and Elderfield2002). 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., Reference Gordon, Cooper, Senior, Banks, Gregory, Johns, Mitchell and Wood2000). The resolution of the atmospheric model is 2.5° latitude by 3.25° longitude with 19 unequally spaced vertical levels (Gordon et al., Reference Gordon, Cooper, Senior, Banks, Gregory, Johns, Mitchell and Wood2000). 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., Reference Cattle, Crossley and Drewry1995). Coupling between the model components occurs daily (Gordon et al., Reference Gordon, Cooper, Senior, Banks, Gregory, Johns, Mitchell and Wood2000). The HadCM3 simulations include dynamic vegetation simulated with the TRIFFID vegetation scheme (Cox, Reference Cox2001), 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. (Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix and Davies-Barnard2017).
Climate simulations with HadCM3 have been evaluated against observations (Gordon et al., Reference Gordon, Cooper, Senior, Banks, Gregory, Johns, Mitchell and Wood2000; Collins et al., Reference Collins, Tett and Cooper2001), proxy records (Singarayer and Valdes, Reference Singarayer and Valdes2010; DiNezio and Tierney, Reference DiNezio and Tierney2013), and other general circulation models (Braconnot et al., Reference Braconnot, Otto-Bliesner, Harrison, Joussaume, Peterchmitt, Abe-Ouchi and Crucifix2007a, Reference Braconnot, Otto-Bliesner, Harrison, Joussaume, Peterchmitt, Abe-Ouchi and Crucifix2007b; Flato et al., Reference Flato, Marotzke, Abiodun, Braconnot, Chou, Collins, Cox, Stocker, Qin, Plattner, Tignor, Allen, Boschung and Nauels2013). HadCM3 represents LGM climate conditions relatively well when compared with reconstructions and other Paleoclimate Modelling Intercomparison Project models (Braconnot et al., Reference Braconnot, Otto-Bliesner, Harrison, Joussaume, Peterchmitt, Abe-Ouchi and Crucifix2007a, Reference Braconnot, Otto-Bliesner, Harrison, Joussaume, Peterchmitt, Abe-Ouchi and Crucifix2007b). 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 Reference DiNezio and Tierney2013; DiNezio et al., Reference DiNezio, Timmermann, Tierney, Jin, Otto-Bliesner, Rosenbloom, Mapes, Neale, Ivanovic and Montenegro2016). 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., Reference Hopcroft, Valdes and Beerling2011, Reference Hopcroft, Valdes, Wania and Beerling2014; Singarayer et al., Reference Singarayer, Valdes, Friedlingstein, Nelson and Beerling2011).
SDGVM vegetation and wetland model
The SDGVM is a global primary productivity and phytogeography model (Woodward et al., Reference Woodward, Smith and Emanuel1995; Beerling and Woodward, Reference Beerling and Woodward2001). 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 CO2 concentrations. SDGVM accounts for the main factors driving vegetation productivity, including climate (surface temperature, precipitation, relative humidity), atmospheric CO2 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), CO2, 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 CH4 Inter-comparison of Models Project (WETCHIMP) project in 2013, which aimed to compare and validate available methane models (Melton et al., Reference Melton, Wania, Hodson, Poulter, Ringeval, Spahni and Bohn2013; Wania et al., Reference Wania, Melton, Hodson, Poulter, Ringeval, Spahni and Bohn2013) and the Global Carbon Project methane budget (Saunois et al., Reference Saunois, Bousquet, Poulter, Peregon, Ciais, Canadell and Dlugokencky2016). 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 CO2 concentrations, air temperature, and precipitation (Melton et al., Reference Melton, Wania, Hodson, Poulter, Ringeval, Spahni and Bohn2013). SDGVM showed a 40% increase in global methane emissions due to a 557 ppm step increase in CO2 (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., Reference Melton, Wania, Hodson, Poulter, Ringeval, Spahni and Bohn2013). The tropics proved most sensitive to an increase in CO2 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 km2 (50% more expansive) compared with the present (Sathiamurthy and Voris, Reference Sathiamurthy and Voris2006). 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., Reference Bird, Taylor and Hunt2005; Wurster et al., Reference Wurster, Bird, Bull, Creed, Bryant, Dungait and Paz2010, Reference Wurster, Rifai, Zhou, Haig and Bird2019; Nguyen et al., Reference Nguyen, Moss, Wasson, Stewart and Ziegler2022; Cheng et al., Reference Cheng, Wu, Luo, Liu, Huang, Zhao, Dai and Weng2023); however, the spatial extent of savanna versus forest is debated (e.g., Bird et al., Reference Bird, Taylor and Hunt2005; Wurster et al., Reference Wurster, Bird, Bull, Creed, Bryant, Dungait and Paz2010). 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, Reference Singarayer and Valdes2010; Singarayer et al., Reference Singarayer, Valdes, Friedlingstein, Nelson and Beerling2011), 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 (Reference Singarayer and Valdes2010) and Singarayer et al. (Reference Singarayer, Valdes, Friedlingstein, Nelson and Beerling2011). 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., Reference Singarayer, Valdes, Friedlingstein, Nelson and Beerling2011), 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 δ13C records are generally in good agreement across their interval of overlap (40–26 ka). The millennial-scale δ13C 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 δ13C 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., Reference Partin, Cobb, Adkins, Tuen and Clark2013; Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020).
Figure 3. Stalagmite δ13C, δ18O, initial 234U/238U, and Mg/Ca records for Sulawesi over the last 40 ka. (A) δ13C for GB09-3 and GB11-9 corrected for the effect of atmospheric CO2 on carbon isotope fractionation in C3 plants (Breecker, Reference Breecker2017). Uncorrected δ13C is shown in grey. The large deglacial δ13C 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) δ18O for GB09-3 and GB11-9 corrected for the effect of ice volume (Krause et al., Reference Krause, Gagan, Dunbar, Hantoro, Hellstrom, Cheng, Edwards, Suwargadi, Abram and Rifai2019). Uncorrected δ18O shown in grey. (C) Initial 234U/238U 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 234U/238U 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.
It is important to note that the δ13C records have been corrected for the effect of atmospheric pCO2 on the δ13C of C3 plants (Schubert and Jahren, Reference Schubert and Jahren2012). The transfer of this effect on carbon isotope fractionation in C3 plants above a cave to stalagmites growing within the cave was identified by Breecker (Reference Breecker2017) in a study assessing globally averaged speleothem δ13C records over the past 90 ka. They found that, after accounting for other processes, the effect of atmospheric CO2 is best explained by a C3 plant δ13C sensitivity of −1.6‰ for every 100 parts per million by volume (ppmv) increase in pCO2 from the LGM to the Holocene. Therefore, it is important to correct for the change in stalagmite δ13C that occurs as a result of glacial–interglacial atmospheric pCO2 before investigating δ13C as a recorder of glacial–interglacial vegetation productivity. The Sulawesi stalagmite δ13C values were adjusted by −1.6‰/100 ppmv (Breecker, Reference Breecker2017) relative to modern atmospheric pCO2 (190 ppmv) using the Antarctic ice core composite pCO2 record (Bazin et al., Reference Bazin, Landais, Lemieux-Dudon, Toyé Mahamadou Kele, Veres, Parrenin and Martinerie2013; Bereiter et al., Reference Bereiter, Eggleston, Schmitt, Nehrbass-Ahles, Stocker, Fischer, Kipfstuhl and Chappellaz2015). The corrected records are shown in Figure 3 and the Supplementary Material and are used throughout the analysis.
The δ13C 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 δ13C values. The deglacial interval contains a prominent ~4.2‰ shift in corrected δ13C from the maximum δ13C value of −4.8‰ at 17.7 ka, near the onset of deglaciation (Pedro et al., Reference Pedro, van Ommen, Rasmussen, Morgan, Chappellaz, Moy, Masson-Delmotte and Delmotte2011), to −9‰ at 11.3 ka. This transition from highest to lowest δ13C 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 δ13C.
The Holocene section of the record shows a surprisingly high degree of δ13C variability, most notably a prominent V-like pattern in the Early to Middle Holocene. During this time, the δ13C 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: δ18O, initial uranium isotope activity [234U/238U]i), and Mg/Ca. The Sulawesi stalagmite δ18O data are explored in detail in Krause et al. (Reference Krause, Gagan, Dunbar, Hantoro, Hellstrom, Cheng, Edwards, Suwargadi, Abram and Rifai2019) and interpreted to reflect changes in rainfall and deep atmospheric convection over the IPWP. The deglacial transition towards wetter conditions, signified by lower δ18O values, occurs at ~11.5 ka.
[234U/238U]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., Reference Fairchild, Borsato, Tooth, Frisia, Hawkesworth, Huang, McDermott and Spiro2000, Reference Fairchild, Smith, Baker, Fuller, Spötl, Mattey and McDermott2006; Hellstrom and McCulloch, Reference Hellstrom and McCulloch2000; Fairchild and Treble, Reference Fairchild and Treble2009). Because uranium isotopes are not thought to be fractionated by natural processes such as calcite precipitation, [234U/238U]i is expected to reflect the activity ratio in seepage waters forming speleothems. The [234U/238U]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, [234U/238U]i can increase as a result of preferential leaching of 234U from alpha recoil–weakened crystal lattice sites in limestone bedrock (Hellstrom and McCulloch, Reference Hellstrom and McCulloch2000). 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, Reference Hellstrom and McCulloch2000). During wetter periods, when water is moving quickly through the epikarst and bedrock dissolution is more uniform, [234U/238U]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, Reference Fairchild and Treble2009); 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., Reference Fairchild, Borsato, Tooth, Frisia, Hawkesworth, Huang, McDermott and Spiro2000; Fairchild and Treble, Reference Fairchild and Treble2009). PCP tends to occur when infiltration rates are low, drip intervals are long, and seepage waters encounter an air-filled space with a pCO2 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, Reference Fairchild and Treble2009).
Sulawesi stalagmite [234U/238U]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 δ18O, 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 [234U/238U]i, is interpreted to reflect wetter conditions.
Previous studies have shown that similar decreases in stalagmite δ18O during the mid-late stages of the last deglaciation are related to increased rainfall in the Indonesian region (Partin et al., Reference Partin, Cobb, Adkins, Clark and Fernandez2007; Griffiths et al., Reference Griffiths, Drysdale, Gagan, Zhao, Ayliffe, Hellstrom and Hantoro2009; Ayliffe et al., Reference Ayliffe, Gagan, Zhao, Drysdale, Hellstrom, Hantoro and Griffiths2013; Carolin et al., Reference Carolin, Cobb, Adkins, Clark, Conroy, Lejau, Malang and Tuen2013). The multiproxy agreement between the Sulawesi δ18O, [234U/238U]i, and Mg/Ca records, alongside other regional rainfall δ18O 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 ([234U/238U]i, Mg/Ca) stabilize at interglacial values of ~10 ka, whereas δ18O 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 δ18O) over the IPWP continued to evolve.
Comparison of the Sulawesi stalagmite δ13C and the three stalagmite rainfall proxies ([234U/238U]i, Mg/Ca, δ18O) provides information about the potential relationship between rainfall and stalagmite δ13C via the influence of rainfall on vegetation (Fig. 3). The Sulawesi δ13C record shows little similarity in large-scale trends with the Sulawesi stalagmite rainfall proxies. The initiation of the deglacial transition in δ13C (~17.5 ka) leads the shift in rainfall (~11.5 ka), on average, by ~6 ka. The early deglacial decrease in stalagmite δ13C, therefore, cannot be driven by the deglacial increase in rainfall.
Stalagmite δ13C as a proxy for vegetation productivity
Previous studies have shown that large shifts in stalagmite δ13C, as observed in the Sulawesi record, can be produced by natural and anthropogenic changes in tropical vegetation cover (e.g., Cruz et al., Reference Cruz, Burns, Karmann, Sharp, Vuille and Ferrari2006; Griffiths et al., Reference Griffiths, Drysdale, Gagan, Hellstrom, Couchoud, Ayliffe, Vonhof and Hantoro2013; Hartmann et al., Reference Hartmann, Eiche, Neumann, Fohlmeister, Schröder-Ritzrau, Mangini and Haryono2013; Burns et al., Reference Burns, Godfrey, Faina, McGee, Hardt, Ranivoharimanana and Randrianasy2016; Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020). In tropical landscapes dominated by C3 plants, the δ13C of dissolved inorganic carbon (DIC) in cave drip waters is primarily set by the mass balance of carbon derived from plant root–respired CO2 and soil microbial activity (~80–90% of the carbon) and carbonate from limestone dissolution (Vogel and Kronfeld, Reference Vogel and Kronfeld1997; Genty et al., Reference Genty, Baker, Massault, Proctor, Gilmour, Pons-Branchu and Hamelin2001; Hou et al., Reference Hou, Tan, Cheng and Liu2003; Griffiths et al., Reference Griffiths, Fohlmeister, Drysdale, Hua, Johnson, Hellstrom, Gagan and Zhao2012; Meyer et al., Reference Meyer, Feng, Breecker, Banner and Guilfoyle2014; Wong and Breecker Reference Wong and Breecker2015; Burns et al., Reference Burns, Godfrey, Faina, McGee, Hardt, Ranivoharimanana and Randrianasy2016). The δ13C value of DIC (and speleothems) in such settings is generally around −8 to −12‰ (McDermott, Reference McDermott2004; Fairchild et al., Reference Fairchild, Smith, Baker, Fuller, Spötl, Mattey and McDermott2006). By contrast, in the absence of vegetation, the δ13C of drip water would reflect the mixing of carbon from atmospheric CO2 (e.g., around −6 to −7‰ at the LGM) and local bedrock (~ +1‰), with stalagmite δ13C 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 δ13C.
It is likely, therefore, that most of the Sulawesi δ13C signal reflects changes in temperature and atmospheric CO2 concentrations, through their combined influence on vegetation type, plant root respiration, and soil microbial activity over Gempa Bumi Cave (e.g., Cosford et al., Reference Cosford, Qing, Mattey, Eglington and Zhang2009; Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020; Lechleitner et al., Reference Lechleitner, Day, Kost, Wilhelm, Haghipour, Henderson and Stoll2021). Temperature and CO2 covary on glacial−interglacial timescales (e.g., Petit et al., Reference Petit, Jouzel, Raynaud, Barkov, Barnola, Basile and Bender1999; NGRIP, 2004; EPICA, 2006), and their individual effects on vegetation productivity (and stalagmite δ13C) are not easily separated. However, model studies designed to look at the relative influence of temperature and CO2 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, Reference Harrison and Prentice2003). Other studies lend support to CO2 as the dominant determinant of vegetation productivity in the tropics (Bennett and Willis, Reference Bennett and Willis2000; Bragg et al., Reference Bragg, Prentice, Harrison, Eglinton, Foster, Rommerskirchen and Rullkötter2013; Claussen et al., Reference Claussen, Selent, Brovkin, Raddatz and Gayler2013; Zhu et al., Reference Zhu, Piao, Myneni, Huang, Zeng, Canadell and Ciais2016; Chen et al., Reference Chen, Zhu, Ciais, Huang, Viovy and Kageyama2019), particularly during the LGM, when atmospheric CO2 is relatively low (Cowling and Field, Reference Cowling and Field2003).
Comparison of Sulawesi stalagmite δ13C with leaf wax δ13C records from Lake Towuti (Russell et al., Reference Russell, Vogel, Konecky, Bijaksana, Huang, Melles, Wattrus, Costa and King2014), Lake Matano (Wicaksono et al., Reference Wicaksono, Russell and Bijaksana2015), and Mandar Bay (Wicaksono et al., Reference Wicaksono, Russell, Holbourn and Kuhnt2017) 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 δ13C corresponds with the relative abundance of C3:C4 plants and/or changes in water and carbon use efficiency by C3 plants, often related to factors such as soil moisture, precipitation, temperature, and humidity (Diefendorf et al., Reference Diefendorf, Mueller, Wing, Koch and Freeman2010). The similarity of the Gempa Bumi stalagmite δ13C 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 δ13C, indicating that vegetation changes above the cave site are less sensitive to rainfall. On the other hand, the Sulawesi stalagmite δ13C and Borneo cave temperature record (Løland et al., Reference Løland, Krüger, Fernandez, Buckingham, Carolin, Sodemann, Adkins, Cobb and Meckler2022) 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 δ13C from Sulawesi, wherein the authors attribute changes in local rainfall as the main driver influencing vegetation type (Russell et al., Reference Russell, Vogel, Konecky, Bijaksana, Huang, Melles, Wattrus, Costa and King2014; Wicaksono et al., Reference Wicaksono, Russell and Bijaksana2015). 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 δ13C of Sulawesi leaf wax. (A) δ13C for stalagmite GB09-3, reflecting changes in vegetation productivity above Gempa Bumi Cave. (B) 230Th-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., Reference Løland, Krüger, Fernandez, Buckingham, Carolin, Sodemann, Adkins, Cobb and Meckler2022). (C) Leaf wax δ13C records for Lake Matano (Wicaksono et al., Reference Wicaksono, Russell and Bijaksana2015), Lake Towuti (Russell et al., Reference Russell, Vogel, Konecky, Bijaksana, Huang, Melles, Wattrus, Costa and King2014), and Mandar Bay (Wicaksono et al., Reference Wicaksono, Russell, Holbourn and Kuhnt2017). The figure is adapted from Wicaksono et al. (Reference Wicaksono, Russell, Holbourn and Kuhnt2017). Leaf wax δ13C corresponds with the relative abundance of C3:C4 plants and/or changes in water and carbon use efficiency by C3 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 δ13C, regional sea-surface temperatures (SSTs), global temperature, and atmospheric CO2 over the last 40 ka supports our interpretation that δ13C is recording changes in vegetation productivity, driven primarily by temperature and CO2 (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., Reference Linsley, Rosenthal and Oppo2010). SSTs then rise concurrently with atmospheric CO2 during the last deglaciation, starting at ~18.5–17.5 ka (Lea et al., Reference Lea, Pak and Spero2000; Stott et al., Reference Stott, Poulsen, Lund and Thunell2002; Visser et al., Reference Visser, Thunell and Stott2003; Linsley et al., Reference Linsley, Rosenthal and Oppo2010), completing the transition by ~11.5 ka. The timing of the late-glacial and deglacial trends in SST and atmospheric CO2 is mirrored in the Gempa Bumi Cave stalagmite δ13C record (Fig. 5). The three records diverge during the Holocene, suggesting that neither temperature nor CO2 is the dominant driver of Sulawesi δ13C at this time.
Figure 5. Relationship between Sulawesi stalagmite δ13C, temperature, atmospheric CO2, and CH4 over the last 40 ka. (A) δ13C 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., Reference Stott, Poulsen, Lund and Thunell2002) and composite SST anomalies for the western IPWP (Linsley et al., Reference Linsley, Rosenthal and Oppo2010 and references therein). (C) Antarctic temperature inferred from ice core δD (Jouzel et al., Reference Jouzel, Masson-Delmotte, Cattani, Dreyfus, Falourd, Hoffmann and Minster2007). (D) Composite Antarctic ice core CO2 concentrations (Bereiter et al., Reference Bereiter, Eggleston, Schmitt, Nehrbass-Ahles, Stocker, Fischer, Kipfstuhl and Chappellaz2015 and references therein). (E) Antarctic ice core CH4 concentrations (Loulergue et al., Reference Loulergue, Schilt, Spahni, Masson-Delmotte, Blunier, Lemieux, Barnola, Raynaud, Stocker and Chappellaz2008). Ice core records are plotted on the AICC2012 chronology (Bazin et al., Reference Bazin, Landais, Lemieux-Dudon, Toyé Mahamadou Kele, Veres, Parrenin and Martinerie2013). Heinrich stadial 1 (HS1), Bølling-Allerød (B-A), and Younger Dryas (YD) are indicated by shaded bars. The close association between Sulawesi δ13C, regional SSTs and air temperature, and atmospheric CO2, particularly during abrupt deglacial climate events, supports the interpretation that Sulawesi δ13C is recording changes in vegetation and soil productivity, driven by changes in temperature and CO2.
Pollen records from marine sediment cores around Sulawesi provide a basis for evaluating the potential influence of shifts in C3:C4 vegetation cover on the Gempa Bumi Cave δ13C record. The pollen assemblages in some sediment cores throughout the IPWP region suggest that C4 grasslands became more common at the LGM (Hope, Reference Hope2001; Bird et al., Reference Bird, Taylor and Hunt2005; Russell et al., Reference Russell, Vogel, Konecky, Bijaksana, Huang, Melles, Wattrus, Costa and King2014; Wicaksono et al., Reference Wicaksono, Russell, Holbourn and Kuhnt2017). 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., Reference Visser, Thunell and Goñi2004). Thus, we cannot rule out the possibility that the balance between C3:C4 vegetation types varied substantially throughout the IPWP region over the last glacial cycle.
In summary, the agreement between Sulawesi δ13C, regional SSTs, temperature, and atmospheric CO2 supports our conclusion that δ13C is recording changes in vegetation productivity. Through this mechanism, we explore the use of the Sulawesi stalagmite δ13C 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 δ13C and climate conditions driving vegetation productivity above the cave system would also affect regional terrestrial sources of atmospheric methane. For example, as rising CO2 and temperature increase vegetation productivity above the cave, warmer conditions in the tropics may also enhance biochemical processes in wetlands (Salimi et al., Reference Salimi, Almuktar and Scholz2021), prompting an increase in methane emissions (Cao et al., Reference Cao, Marshall and Gregson1996; Kleinen et al., Reference Kleinen, Mikolajewicz and Brovkin2020). Thus, while stalagmite δ13C 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 δ13C (vegetation productivity) shows good correspondence with EPICA ice core methane (Loulergue et al., Reference Loulergue, Schilt, Spahni, Masson-Delmotte, Blunier, Lemieux, Barnola, Raynaud, Stocker and Chappellaz2008) 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., Reference Kaplan, Folberth and Hauglustaine2006), 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 CO2 and temperature rise at 17.5 ka, marking deglacial onset in Sulawesi vegetation. The highest δ13C 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 CO2, 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., Reference Kienast, Hanebuth, Pelejero and Steinke2003; Weaver et al., Reference Weaver, Saenko, Clark and Mitrovica2003; Rosen et al., Reference Rosen, Brook, Severinghaus, Blunier, Mitchell, Lee, Edwards and Gkinis2014) and at ~11.5 ka during the Younger Dryas (YD) (Fairbanks, Reference Fairbanks1989; McManus et al., Reference McManus, Francois, Gherardi, Keigwin and Brown-Leger2004; Cheng et al., Reference Cheng, Zhang, Spötl, Baker, Sinha, Li and Bartolomé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 δ13C) occur at the end of HS1 and the YD and correspond with times of abrupt increases in atmospheric CO2 and methane. Previous studies have suggested that the rapid shifts in global methane are driven by tropical wetlands (Schaefer et al., Reference Schaefer, Whiticar, Brook, Petrenko, Ferretti and Severinghaus2006; Rosen et al., Reference Rosen, Brook, Severinghaus, Blunier, Mitchell, Lee, Edwards and Gkinis2014). Thus, the tropics may be a key contributor to the global methane budget during times of increasing CO2 and/or large-scale heat exchange across hemispheres.
The agreement between the Sulawesi stalagmite δ13C and ice core methane becomes decoupled after 10 ka, when stalagmite δ13C 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 δ13C. 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, CO2, and moisture are readily available (Cowling and Field, Reference Cowling and Field2003). Strengthening of the summer monsoon and strong seasonality could also influence productivity patterns (Vargas-Terminel et al., Reference Vargas-Terminel, Flores-Rentería, Sánchez-Mejía, Rojas-Robles, Sandoval-Aguilar, Chávez-Vergara, Robles-Morua, Garatuza-Payan and Yépez2022).
Comparison with global vegetation model simulations
To put these findings into broader context, we investigate the relationship between atmospheric methane concentrations and Sulawesi stalagmite δ13C 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., Reference Melton, Wania, Hodson, Poulter, Ringeval, Spahni and Bohn2013).
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, Reference Kaplan2002; Kleinen et al., Reference Kleinen, Mikolajewicz and Brovkin2020). 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., Reference Wurster, Rifai, Zhou, Haig and Bird2019; Nguyen et al., Reference Nguyen, Moss, Wasson, Stewart and Ziegler2022; Cheng et al., Reference Cheng, Wu, Luo, Liu, Huang, Zhao, Dai and Weng2023). This study, however, investigates the Sulawesi stalagmite δ13C 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., Reference Lea, Pak and Spero2000; Linsley et al., Reference Linsley, Rosenthal and Oppo2010; Gagan et al., Reference Gagan, Hendy, Haberle and Hantoro2004; Løland et al., Reference Løland, Krüger, Fernandez, Buckingham, Carolin, Sodemann, Adkins, Cobb and Meckler2022; 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, Reference Aselmann and Crutzen1989; Cao et al., Reference Cao, Marshall and Gregson1996; Guo et al., Reference Guo, Zhou and Wu2012). The relative source changes simulated by the SDGVM agree well with previous studies (Chappellaz et al., Reference Chappellaz, Blunier, Kints, Dällenbach, Barnola, Schwander, Raynaud and Stauffer1997; Dällenbach et al., Reference Dällenbach, Blunier, Flückiger, Stauffer, Chappellaz and Raynaud2000; Valdes et al., Reference Valdes, Beerling and Johnson2005; Kaplan et al., Reference Kaplan, Folberth and Hauglustaine2006; Fischer et al., Reference Fischer, Behrens, Bock, Richter, Schmitt, Loulergue and Chappellaz2008; Hopcroft et al., Reference Hopcroft, Valdes, O'Connor, Kaplan and Beerling2017; Kleinen et al., Reference Kleinen, Mikolajewicz and Brovkin2020).
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 δ13C record with the SDGVM model output, we identify soil respiration as the model parameter closest to stalagmite δ13C and use this parameter as a proxy for our record within the model. Soil respiration is the emission of CO2 from the soil surface (Schlesinger and Andrews, Reference Schlesinger and Andrews2000) that is produced within the soil profile by roots and soil organisms (Raich and Schlesinger, Reference Raich and Schlesinger1992). The predominant climatic driver of soil respiration rates is debated, but it is generally agreed that temperature, CO2, and soil moisture all play important roles in driving soil respiration rates (Raich and Schlesinger, Reference Raich and Schlesinger1992; Bragg et al., Reference Bragg, Prentice, Harrison, Eglinton, Foster, Rommerskirchen and Rullkötter2013; Hursh et al., Reference Hursh, Ballantyne, Cooper, Maneta, Kimball and Watts2017), with seasonality and forest structure also exerting control (Vargas-Terminel et al., Reference Vargas-Terminel, Flores-Rentería, Sánchez-Mejía, Rojas-Robles, Sandoval-Aguilar, Chávez-Vergara, Robles-Morua, Garatuza-Payan and Yépez2022). It also has been found that wetland drying significantly increases the temperature sensitivity of soil respiration rates (Chen et al., Reference Chen, Zou, Cui, Nie and Fang2018).
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, Reference Schlesinger and Andrews2000), and within the SDGVM, it correlates strongly with net primary productivity (r = 0.98). The rate of soil respiration sets the concentration of CO2 within the soil profile (Raich and Schlesinger, Reference Raich and Schlesinger1992), which is the most likely primary source for carbon in the Sulawesi stalagmites. Therefore, we use soil respiration as a qualitative proxy for stalagmite δ13C 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 CO2 concentration during and since the glacial period. Increasing atmospheric CO2 (and its fertilising influence on vegetation) accounts for half of the total LGM to preindustrial amplitude increase in soil respiration. Thus, atmospheric CO2 is a primary driver of vegetation productivity for modelled soil respiration rates throughout the LGM. This is consistent with the underlying hypothesis for atmospheric CO2 and temperature as external factors driving Sulawesi stalagmite δ13C.
The relationship between Sulawesi stalagmite δ13C 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 δ13C record (Figure 8). Stalagmite δ13C 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 δ13C 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 δ13C 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 δ13C. (A–C) Time series of modelled mean soil respiration for the grid points corresponding to Sulawesi, Indonesia, and tropics (±30°). Sulawesi stalagmite δ13C 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 δ13C, 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 δ13C as a reliable indicator of local-to-regional methane emissions, we examine the correlation between Sulawesi stalagmite δ13C 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 δ13C time series. The timing of deglacial increases in methane emissions across all three regions coincides with Sulawesi stalagmite δ13C (Fig. 9). Each time series is correlated with stalagmite δ13C, 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 δ13C. (A–C) Time series of modelled methane emissions totals for Sulawesi, Indonesia, and tropics (±30°). Sulawesi stalagmite δ13C 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 δ13C, 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 δ13C 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. (Reference Singarayer, Valdes, Friedlingstein, Nelson and Beerling2011) and Burns (Reference Burns2011) 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 δ13C record appears to support this. For example, increased convective rainfall in the Holocene is supported by Sulawesi stalagmite δ18O (Krause et al., Reference Krause, Gagan, Dunbar, Hantoro, Hellstrom, Cheng, Edwards, Suwargadi, Abram and Rifai2019), [234U/238U]i, and Mg/Ca, and by other IPWP records (e.g., Partin et al., Reference Partin, Cobb, Adkins, Clark and Fernandez2007; Griffiths et al., Reference Griffiths, Drysdale, Gagan, Zhao, Ayliffe, Hellstrom and Hantoro2009; Ayliffe et al., Reference Ayliffe, Gagan, Zhao, Drysdale, Hellstrom, Hantoro and Griffiths2013; Scroxton et al., Reference Scroxton, Gagan, Ayliffe, Hantoro, Hellstrom, Cheng, Edwards, Zhao, Suwargadi and Rifai2022). The disconnect between stalagmite δ13C V-like pattern and Holocene temperature–atmospheric CO2 (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 δ13C 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 δ13C and simulated tropical methane correspond well over the glacial period (Fig. 10). Departures of ice core methane from simulated tropical methane and the Sulawesi δ13C 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 δ13C (Fig. 10). The plateau in stalagmite δ13C 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., Reference Singarayer, Valdes, Friedlingstein, Nelson and Beerling2011); thus, the step change in the deglacial pattern in the model is likely occurring due to step changes in the corresponding atmospheric CO2 supplied to the model (Singarayer and Valdes, Reference Singarayer and Valdes2010).
Figure 10. Sulawesi δ13C as a potential indicator of the contribution of tropical methane to global atmospheric methane. Comparison of Sulawesi stalagmite δ13C, ice core methane concentrations (plotted on the AICC2012 chronology; Bazin et al., Reference Bazin, Landais, Lemieux-Dudon, Toyé Mahamadou Kele, Veres, Parrenin and Martinerie2013) and modelled total methane emissions for the tropics. The Sulawesi δ13C 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 δ13C 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 CO2 and temperature. This soil respiration signature was, in turn, recorded by stalagmite δ13C via seepage waters, which retain the carbon isotope signature of the plant matter and soil CO2 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., Reference Brook, Harder, Severinghaus, Steig and Sucher2000; Fischer et al., Reference Fischer, Behrens, Bock, Richter, Schmitt, Loulergue and Chappellaz2008; Weber et al., Reference Weber, Drury, Toonen and van Weele2010; Baumgartner et al., Reference Baumgartner, Schilt, Eicher, Schmitt, Schwander, Spahni, Fischer and Stocker2012; Guo et al., Reference Guo, Zhou and Wu2012; Rhodes et al., Reference Rhodes, Brook, Chiang, Blunier, Maselli, McConnell, Romanini and Severinghaus2015, Reference Rhodes, Brook, McConnell, Blunier, Sime, Xavier and Mulvaney2017; Kleinen et al., Reference Kleinen, Mikolajewicz and Brovkin2020). In the SDGVM model simulations, tropical wetland methane emissions are largely controlled by changing soil respiration rates, raising the possibility that the Sulawesi stalagmite δ13C 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 δ13C 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 δ13C 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 δ13C show some evidence for tropical methane sources contributing to Late Holocene methane variability. A disconnect between stalagmite δ13C, temperature, global atmospheric CO2, 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 δ13C 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 δ13C 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.
Acknowledgments
The fieldwork in Indonesia was carried out under Kementerian Negara Riset dan Teknologi (RISTEK) research permit numbers 04/TKPIPA/FRP/SM/IV/2009 and 1b/TKPIPA/FRP/SM/I/ 2011 with the support of the Research Center for Geotechnology, Indonesian Institute of Sciences (LIPI). We are grateful for the invaluable field assistance provided by Neil Anderson, Dan Zwartz, Garry Smith, Linda Ayliffe, Nick Scroxton, Engkos Kosasih, Djupriono, and the staff of Bantimurung-Bulusaraung National Park (with special thanks to Syaiful Fajrin). We also thank Heather Scott-Gagan, Joan Cowley, Joe Cali, Linda McMorrow, Chris Vardanega, and Daniel Becker for laboratory assistance, and Joy Singarayer and David Beerling for providing HadCM3 and SDGVM simulations for analysis. The work was funded by an Australian Postgraduate Award to CEK; Australian Research Council (ARC) Discovery grants DP0663274, DP1095673, DP110101161, and DP180103762 to MKG, WSH, JCH, RLE, and HC; ARC Future Fellowship FT130100801 to JCH; NERC UK projects NE/I010912/1 and NE/P002536/1 to POH; U.S. National Science Foundation grant 2202913 to RLE and HC; and National Natural Science Foundation of China grants NSFC 41731174 and 41888101 to HC. The authors declare no competing interests.
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.
Open access
