Atmospheric greenhouse gas concentrations are thought to have synchronized global temperatures during Pleistocene glacial–interglacial cycles, yet their impact relative to changes in high-latitude insolation and ice-sheet extent remains poorly constrained. Here, we use tropical glacial fluctuations to assess the timing of low-latitude temperature changes relative to global climate forcings. We report 10Be ages of moraines in tropical East Africa and South America and show that glaciers reached their maxima at ~29 to 20 ka, during the global Last Glacial Maximum. Tropical glacial recession was underway by 20 ka, before the rapid CO2 rise at ~18.2 ka. This “early” tropical warming was influenced by rising high-latitude insolation and coincident ice-sheet recession in both polar regions, which lowered the meridional thermal gradient and reduced tropical heat export to the high latitudes.
The interhemispheric synchrony of the ice age cycles is one of the greatest questions in paleoclimate research. During much of the Pleistocene, Northern Hemisphere (NH) high-latitude summer insolation appears to have paced global glaciation, yet the mechanisms by which this climate response was propagated to the Southern Hemisphere (SH) are uncertain (1, 2). The rapid terminations of ice ages complicate the problem, as these occurred during periods of both high- and low-amplitude NH high-latitude summer insolation changes (3).
Atmospheric greenhouse gases (GHGs), particularly CO2, are commonly invoked to explain the global synchrony of glacial-interglacial cycles (4). This hypothesis is supported by near-concurrent changes in atmospheric CO2 levels and Antarctic temperatures over the last ~800,000 years (5), and numerous studies point to atmospheric GHGs as a key factor in warming Earth during the last deglaciation (6–8). However, other climate processes may have also played a role in unifying global temperatures. For example, low NH high-latitude summer insolation and coincident high ice-sheet albedo altered the meridional thermal gradient, which likely affected the strength and position of wind belts (9). Such changes in atmospheric circulation may have altered the strength of atmospheric and oceanic heat transport during glacial periods and induced cooling separate from that caused by GHG changes (9). The interconnected nature of high-latitude and GHG forcings complicates efforts to assess the sensitivity of Earth’s climate system to changes in these boundary conditions.
The tropics are an ideal region in which to investigate the impact of global forcings on glacial-interglacial temperature change as they are far from the direct forcing from high-latitude insolation and large ice sheets (10, 11). Tropical glaciers, in particular, provide a valuable record of past change in the low latitudes as these glaciers are highly sensitive to changes in temperature (see the Supplementary Materials) (12), and their past fluctuations reflect changes in mid-tropospheric temperature (13). Prior work using cosmogenic beryllium-10 (10Be) surface exposure dating and analyses of glacially influenced lake sediments indicates that some tropical glaciers achieved their maximum extents either before or early during the Last Glacial Maximum (LGM; ~26.5 to 19.0 ka) (14–17) and that deglaciation from their LGM maxima was underway by ~20 ka before the rapid CO2 rise at ~18.2 ka (18, 19). Critically, these existing data are limited to tropical South America. It is therefore unknown whether these records reflect temperature across the wider tropics or more regional conditions. In addition, the recent determination of low-latitude, high-elevation cosmogenic nuclide production rates (20, 21) requires recalculation of 10Be ages from tropical South America to allow for accurate comparison with global records.
To assess the timing of glacial fluctuations across the tropics during the LGM, we determined a chronology of past glacial extents using 10Be dating of moraines in tropical East Africa, far from the South American tropics. Our 10Be chronology includes 17 new and 8 previously published ages (22) that constrain the timing of glacial fluctuations in the equatorial Rwenzori Mountains (~0.3°N, 30.0°E), located on the border between Uganda and the Democratic Republic of the Congo, during the LGM and the onset of deglaciation (Fig. 1). We also recalculated 177 10Be ages of 48 LGM moraines in tropical South America from 10 prior studies (nine sites) (Fig. 2 and the Supplementary Materials). All ages are calculated using a low-latitude, high-altitude 10Be production rate (21) and time-independent (“St”) scaling (see Materials and Methods) (23–25). There are notable uncertainties in tropical cosmogenic nuclide production rate scaling during the LGM (26). Although these uncertainties affect comparisons of tropical 10Be glacial chronologies with other paleoclimate records, they do not affect the larger conclusions of this work. All 10Be ages reported here are of boulders on the crests of moraines, which we infer as representing the final sedimentation on a moraine. Therefore, we interpret a 10Be age as the time of glacial recession from a given moraine and, thus, the onset of warming.
During the last ice age, glaciers from the central peaks in the Rwenzori Mountains flowed down the Bujuku and Mubuku valleys and merged to form a single glacier in the lower Mubuku valley that terminated at ~2000-m elevation (Fig. 1). At its maximum extent, the glacier overtopped the right-lateral (south) valley wall, flowed ~500 m southward, and formed a series of moraines that enclose Lake Mahoma at ~3000 meters above sea level (m.a.s.l.) (Fig. 1 and tables S1 to S3). We term these moraines, from the outermost to innermost, the Mahoma 4, 3, 2, and 1 moraines. 10Be ages of the two outermost moraines yield mean ages of 28,960 ± 1020 years ago (Mahoma 4; n = 6) and 28,740 ± 620 years ago (Mahoma 3; n = 3; one outlier omitted). Inboard of the Mahoma 3 moraine, the Mahoma 2 moraine has a mean age of 24,940 ± 500 years ago (n = 4) (22). The Mahoma 2 moraine is cross cut by the largest of these moraines, Mahoma 1, which has a mean age of 21,500 ± 650 years ago (n = 4) (22) and extends down the Mubuku valley to ~2300 m.a.s.l. The Mahoma 1 moraine has the greatest relief of any moraine within the series (~150 m relief above the Mubuku river valley floor) and marks the farthest down-valley extent of ice in the Mubuku valley during the last ice age. Approximately 5 km up the Mubuku valley from the Mahoma 1 terminus, the Mahoma 0 moraine has a mean age of 17,860 ± 830 years ago (n = 3).
We also dated a single large (30 to 50 m relief above the valley floor), sharp-crested moraine in the Moulyambouli valley, ~10 km south of the Mubuku valley (Fig. 1). While there are smaller, partially preserved moraine segments on the steep bedrock valley wall opposite the large Moulyambouli moraine, there is no series of moraines comparable to that in the lower Mubuku valley. 10Be ages of the Moulyambouli moraine yield a mean age of 18,750 ± 340 years ago (n = 4). Because there is only a single moraine dated in the Moulyambouli valley, we rely more heavily on the moraine chronology from the lower Mubuku valley where we can track glacial fluctuations throughout the LGM. Together, the Mahoma and Moulyambouli moraine chronologies indicate that Rwenzori glaciers reached their maximum extents by ~28 ka. The glaciers fluctuated near their maxima until ~21.5 ka in the Mubuku valley and ~18.8 ka in the Moulyambouli valley, at which point deglaciation was underway in both catchments.
A comparison of the Rwenzori moraine chronology with 48 recalculated moraine ages from nine sites in tropical South America shows a broad similarity in the timing and structure of tropical glacial fluctuations during the LGM (Fig. 2 and the Supplementary Materials). South American glaciers fluctuated throughout the LGM, with some glaciers achieving their maximum extents by ~28 to 29 ka. Similar to the Rwenzori glaciers, the tropical South American glaciers retreated from or from near their LGM maxima by ~20 to 19 ka. Together, the Rwenzori and South American moraine chronologies show a coherent signal of cool tropical temperatures during the LGM until the onset of glacial recession, indicating warming at ~20 to 19 ka.
The coincidence of tropical glacial maxima with a minimum in NH high-latitude summer insolation, high global ice-sheet volume, and low atmospheric GHG concentrations shows that tropical temperatures were sensitive to these global climate forcings during the LGM (Fig. 3). Milankovitch theory holds that low NH high-latitude summer insolation enabled the growth of NH ice sheets, which, in turn, cooled the planet via increased albedo (8, 27). Because late Pleistocene sea-level changes largely reflect changes in NH ice sheets (28), we use sea level as a proxy for NH ice-sheet volume. By this metric, the period of most extensive NH ice sheets (~29 to 19 ka), and likely the period of greatest ice-sheet albedo forcing, is coincident with the time of expanded tropical glaciers. GHG radiative forcing was also low from ~34 to 18 ka (Fig. 3) (29–31), which presumably reinforced tropical cooling during the LGM (32), particularly at high altitudes (33).
Although atmospheric GHG concentrations influenced tropical cooling during the LGM, the glacial chronologies presented here indicate that tropical warming began “early” (i.e., at ~20 to 19 ka) before the rapid rise of CO2 at ~18.2 ka (34) that is hypothesized to have influenced global deglaciation [e.g., (3)]. The early warming registered by tropical glaciers is supported by temperature reconstructions from East African lake sediments that document the onset of deglacial warming at ~20 ka (Fig. 4) (35, 36). Moreover, synthesized tropical sea surface temperature (SST) records (7) indicate that warming was underway across the tropical ocean (30°N-30°S) by at least ~19 ka (Fig. 4). This early onset of tropical warming does not align with changes in atmospheric GHG levels or local insolation. GHG radiative forcing rose by ~0.5 W/m2 between ~25 and 20 ka but stagnated or fell between ~20 and 18 ka, a time when tropical glaciers receded and temperatures increased (Figs. 3 and 4). In addition, mean annual equatorial insolation decreased by ~2.0 W/m2 between ~30 and 16 ka. Although the rate and magnitude of change were low, decreasing insolation would have encouraged glacial advance throughout this period, including at ~20 to 19 ka. An alternative mechanism is required to explain the timing of tropical warming at the end of the LGM.
Coincident with, or possibly before, the tropical warming at ~20 to 19 ka, NH and SH high-latitude regions warmed. The Laurentide Ice Sheet began to retreat from or from near its maximum extents before ~20 ka (37, 38), driven by rising NH high-latitude summer insolation after ~24 ka (38). Deglaciation on the Antarctic Peninsula (39) and the onset of warming in West Antarctica (40) are dated to ~20 ka and are attributed to increasing SH summer insolation duration (41). The rise in global sea level at ~21 to 19 ka (28) likely reflects the recession of the large polar ice sheets. Compilations of other NH high- and mid-latitude terrestrial and oceanic records likewise suggest that warming was underway by ~20 ka; SH mid-latitude oceanic records (7) show the onset of warming by ~21 to 20 ka.
We suggest that the deglaciation and warming in NH and SH high-latitude regions by at least ~20 ka led to a reduction of Earth’s meridional thermal gradient (40) and that this reduced gradient influenced warming in the tropics (Fig. 4). The global atmosphere can be modeled as a “heat engine,” with net flow from warm “source” regions at low latitudes to cool heat “sinks” at high latitudes (42). The efficiency of heat flow is greatest when the temperature difference between source and sink regions is large and decreases when the temperature gradient between these regions is reduced. The export of heat from the tropics is dominated by the Hadley circulation, the strength of which is positively related to the meridional thermal gradient (43). This relationship is apparent in simulations of anthropogenic global warming that show decreased Hadley circulation strength as the thermal gradient between the tropics and the poles is reduced (44–46). By analogy, we infer that the insolation-induced warming of both polar regions marked by the recession of NH and SH ice sheets at ~21 to 20 ka influenced tropical warming through reducing the Hadley circulation strength and, thus, net tropical heat export. Slackened Hadley circulation also would have likely decreased poleward oceanic heat export by tropical surface currents (47, 48), further reinforcing tropical warming. Tropical SSTs likewise began to rise at or just before ~19 ka (Fig. 4) (7).
One hypothesis arising from this mechanism is that tropical glacial recession should have occurred nearly coincident with, or immediately after, increases in high-latitude temperatures. Tropical glacial recession was underway at ~20 to 19 ka. In contrast, compilations of higher-latitude temperature records suggest high-latitude warming initiated at ~21 to 20 ka (Fig. 4). Future work testing this hypothesis should focus on further high-resolution reconstructions of tropical temperatures, including dating of tropical glacial extents in Africa and elsewhere in the global tropics, which will reduce uncertainties in the precise phasing of high- and low-latitude glacial fluctuations during and subsequent to the LGM. Nevertheless, our results show that tropical temperatures respond to remote, high-latitude climate forcing and that climate models of the last deglaciation should account for tropical warming before ~18 ka absent the influence of GHG forcing.
MATERIALS AND METHODS
We conducted three field seasons in the Mubuku, Bujuku, and Moulyambouli valleys of the Rwenzori Mountains (~0.3°N, 30°E) between 2012 and 2016. We sampled large (2 to 4 m in diameter) boulders on the crests of moraines that showed little or no sign of postdepositional movement or modification. We obtained ~0.5- to 1-kg pieces of the upper 3 to 5 cm of boulders using a hammer and chisel or a battery-powered hammer drill. We measured the surface dip of each sample using a handheld Suunto compass and any potential shielding by topography using a Suunto clinometer. Where possible, we sampled flat rock surfaces with no dip to minimize potential uncertainties in the shielding correction. We recorded the location and elevation of each sample at the time of collection using a handheld Garmin Global Positioning System (GPS). We averaged GPS measurements over the course of 2 to 5 min, sufficient to record ≥100 measurements. We then recorded the average location and elevation values for each sample (vertical error, ±3 m; horizontal error, ±1 m).
We measured the thickness of each sample to the nearest millimeter (±1 mm) using calipers before crushing and milling the sample to a grain size of 210 to 750 μm. We then put the 210- to750-μm grain-size fraction through a series of o-phosphoric acid, sodium hydroxide, and dilute hydrofluoric-nitric acid leaches to etch quartz and dissolve other minerals. Once we obtained pure quartz, we spiked the samples and a process blank with a known mass of 9Be carrier and digested them in hydrofluoric acid. We used a 9Be carrier developed from a deeply buried beryl crystal. For carrier concentrations, see table S1. We then used a modified version of the procedures described in (49) to isolate beryllium from each sample and the process blank. 10Be/9Be ratios of samples were measured at the Center for Accelerator Mass Spectrometry at Lawrence Livermore National Laboratory (LLNL) and normalized to the 07KNST3110 standard (table S1) (50). Once obtained, we subtracted the 10Be/9Be process blank ratio from the 10Be/9Be ratios of measured samples and used the blank-corrected sample ratios to determine the concentration of 10Be in each sample per gram of quartz. For process blank ratios, see table S1.
10Be age calculation
We calculated the Rwenzori 10Be ages using version 3 of the online exposure age calculator described by Balco et al. [(51) and subsequently updated] using a low-latitude, high-altitude 10Be production rate (21) and the time-independent St scaling scheme (tables S1 and S2 and fig. S1) (23–25). Two 10Be production rate calibrations from tropical high-altitude sites (both in South America) (20, 21) yield the most robust calibration results when paired with St scaling, which suggests that this method may be most appropriate for sites such as the Rwenzori. We note that these 10Be production rate calibrations use glacial deposits that date to ~13 to 11 ka BP, younger than the LGM.
Assuming that changes in the geomagnetic field influenced the cosmic-ray flux and, thus, the 10Be production rate during and since the LGM, a time-dependent scaling scheme (e.g., “Lm” or “LSDn”) (52) may be more appropriate for calculating the Rwenzori 10Be ages. Time-dependent scaling schemes such as Lm and LSDn (51) yield less accurate calibrations of the Kelly et al. (21) dataset than St, with Lm scaling producing the greatest scatter in the resulting calibration data. Moreover, Lm scaling is based on paleomagnetic models (51); because the theorized impacts of changes in the magnetic field on nuclide production are greatest in the low-latitude regions (53), we suggest that Lm scaling is less appropriate than other scaling schemes for sites in the tropics. LSDn scaling yields more accurate calibration results than the Lm scaling and is based on the cosmic-ray flux rather than magnetic field models (52, 53). We show the Rwenzori 10Be ages calculated using Lm and LSDn scaling in table S2. However, because of uncertainties in the 10Be production rate and scaling schemes in the low latitudes (26), we choose to base our interpretations on 10Be ages calculated using St scaling. Our interpretations regarding the timing of tropical glaciation (during the global LGM) and onset of deglaciation (at ~20 to 19 ka) do not change with LSDn-based calculations. Our interpretations would change using Lm scaling, but because of the relative inaccuracy of Lm when used to calibrate tropical exposure age data, we consider Lm less appropriate to use in this study.
We did not correct the Rwenzori 10Be ages for the potential influence of boulder surface erosion. Although many samples showed evidence of surface erosion (e.g., exfoliation, raised quartz veins), 10Be ages obtained from raised quartz veins are indistinguishable from those obtained from boulder surfaces. For example, on the Mahoma 4 moraine, samples RZ-16-49 (28,650 ± 570 years) and RZ-16-50 (30,100 ± 790 years) are from raised quartz veins (1 to 3 cm above the boulder surface) on gneissic boulders (table S2). The other four samples from the Mahoma 4 moraine are from gneissic boulder surfaces and yielded ages between ~29.9 and 27.5 ka. The mean age of the Mahoma 4 moraine including these quartz vein samples (28,960 ± 1020 years) is statistically indistinguishable from the mean age of the moraine with the quartz vein samples excluded (28,750 ± 1100 years).
In addition, we did not correct the Rwenzori 10Be ages for the influence of cover by snow and/or ice, soil, or vegetation. We never observed snow or ice cover at the sample sites and suggest that snow and ice do not persist on the ground for any length of time at the sample elevations (~2635 to 2990 m.a.s.l.) owing to the relatively warm temperatures and intense equatorial solar radiation. All samples were covered by soil, moss, and various amounts of other vegetation. However, on the basis of the consistency of 10Be ages on individual moraines, we believe that shielding by soil and vegetation was small or negligible.
The χ2 value of each moraine, except for the Mahoma 3 moraine, is lower than the expected value (table S3). This indicates that any age scatter in the 10Be ages on individual moraines can be explained by analytical uncertainty alone rather than by postdepositional processes or modification. We identified outliers in the dataset using Chauvenet’s criterion. After removing a single outlier (RZ-13-63) from the Mahoma 3 moraine ages, the χ2 value of the Mahoma 3 moraine is lower than the expected χ2 value. Sample RZ-13-63 from the Mahoma 3 moraine is the only sample identified as an outlier and is not included in our interpretations or discussion. All “mean ages” described in the text are arithmetic means of the sample ages discussed.
Certain samples were measured twice (i.e., two separate quartz aliquots from the same rock sample were processed and measured) to check for internal sample age consistency. Secondary aliquot measurements are marked as aliquot “a” in tables S1 and S2, and all 10Be ages are shown in tables S1 and S2. We used the secondary aliquot (aliquot “a”) measurements for our interpretations and discussion of 10Be ages from moraine Mahoma 0, because the primary sample measurements returned less scatter between ages. We do not base our larger interpretations on these samples, and we note that the aliquot “a” measurements for all other samples are within error of primary measurements.
Recalculation of preexisting tropical 10Be datasets
We recalculated 10Be ages from prior studies on tropical glacier fluctuations using the same methodology described above for the Rwenzori 10Be ages. We provide a brief review of the sites and samples from prior studies within the Supplementary Materials of this manuscript. We used the St scaling scheme (23–25) for calculating the 10Be ages from prior work, but note that the differences in 10Be ages calculated using the Lm (51) and LSDn (52) scaling schemes are similar to those for the Rwenzori (see tables S1 and S4).
For the purpose of making the most direct comparisons possible in the discussion of the paper, we chose to include only prior work that used 10Be dating to develop glacial chronologies. Therefore, we did not include prior studies that used radiocarbon dating of organic material associated with moraines or surface exposure dating of moraines using other cosmogenic nuclides such as 36Cl and 3He. However, we note that these additional data neither alter nor counter our broader interpretations.
With the recalculated 10Be ages from prior studies in the tropics, we assigned moraine ages at sites using the following methodology. First, we excluded 10Be ages identified as outliers by the original authors. We also excluded 10Be ages that were not from moraines (i.e., 10Be ages of bedrock surfaces or boulders on bedrock surfaces). We then plotted the locations of 10Be ages in Google Earth and evaluated their glacial geologic contexts. We checked these locations using the original published maps of 10Be ages and glacial geomorphology to ensure that samples plotted correctly as reported in the prior studies. We assigned each sample to a moraine based on the original authors’ interpretations. However, where authors grouped 10Be ages of multiple landforms into a single mean age or moraine group, we attempted to assign 10Be ages to individual moraines using Google Earth. Once each 10Be age was assigned to a moraine, we calculated the mean age of the moraine. Descriptions of the recalculated 10Be datasets and assigned moraine ages are presented in the Supplementary Materials (see the “Tropical 10Be Site Descriptions” section). Because our study focuses on the LGM and the onset of the last deglaciation, we excluded 10Be ages >40 ka and mean moraine ages >30 ka, as well as 10Be ages <17 ka from prior studies.
Where possible for prior studies, we also assessed whether individual moraines mark the onset of recession from or from near the LGM maximum position of a given glacier. This was not always possible to determine because, at some sites, only lateral, composite moraines are dated. We detail these classifications in the Supplementary Materials. All recalculated 10Be ages and associated data are shown in table S4.
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/12/eaaw2610/DC1
Fig. S1. Camel plots showing probability-distribution curves for individual moraine ages with sample age statistics.
Table S1. Sample information for Moulyambouli and Mahoma moraines.
Table S2. Calculated 10Be ages from the Moulyambouli and Mahoma moraines.
Table S3. Distribution of 10Be ages from the Moulyambouli and Mahoma moraines as presented in fig. S1.
Table S4. Recalculated 10Be ages of tropical South American moraines.
Table S5. Calibration dataset for Kelly et al. (21) as provided for use with version 3 of the online exposure age calculator described by Balco et al. (51) and subsequently updated.
This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Acknowledgments: We thank the Ugandan Wildlife Authority and the Uganda National Council for Science and Technology for their support of this project and permission to conduct research and collect samples in the Rwenzori Mountains. We are particularly grateful for the essential logistical and safety support from Rwenzori Mountaineering Services and Rwenzori Trekking Services and for the concerted efforts of their guides and porters in sample collection. M. B. Baber and S. E. Loomis helped with field logistics and sample collection. This is LLNL-JRNL-759653. Funding: This research was supported by grants from the NSF (EAR-1222705, EAR-1702293, and GSS-1558358), the National Geographic Society, and the Comer Family Foundation. Additional student support was provided by Sigma Xi and Dartmouth College. Satellite imagery was granted by the Digital Globe Foundation. Author contributions: M.S.J., M.A.K., J.M.R., and A.M.D. collected samples, analyzed results, and wrote and edited the paper and figures. M.S.J., J.A.H., and S.R.H.Z. processed samples for analysis. D.C. collected and processed samples for analysis. S.R.H.Z. analyzed samples. J.W.C. and S.R.H.Z. edited the paper and figures. M.A.K., J.M.R., and B.N. conceptualized the project. B.N. coordinated the project within Uganda. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).
Credit: Source link