Research ArticleCLIMATOLOGY

Increasing Pleistocene permafrost persistence and carbon cycle conundrums inferred from Canadian speleothems

See allHide authors and affiliations

Science Advances  28 Apr 2021:
Vol. 7, no. 18, eabe5799
DOI: 10.1126/sciadv.abe5799


Permafrost carbon represents a potentially powerful amplifier of climate change, but little is known about permafrost sensitivity and associated carbon cycling during past warm intervals. We reconstruct permafrost history in western Canada during Pleistocene interglacials from 130 uranium-thorium ages on 72 speleothems, cave deposits that only accumulate with deep ground thaw. We infer that permafrost thaw extended to the high Arctic during one or more periods between ~1.5 million and 0.5 million years ago but has been limited to the sub-Arctic since 400,000 years ago. Our Canadian speleothem growth history closely parallels an analogous reconstruction from Siberia, suggesting that this shift toward more stable permafrost across the Pleistocene may have been Arctic-wide. In contrast, interglacial greenhouse gas concentrations were relatively stable throughout the Pleistocene, suggesting that either permafrost thaw did not trigger substantial carbon release to the atmosphere or it was offset by carbon uptake elsewhere on glacial-interglacial time scales.


Permafrost, which covers nearly a quarter of the Northern Hemisphere land surface, has warmed considerably over the past few decades and may be highly vulnerable to thaw with continued climate change (1). This thaw vulnerability has potential global implications because high-latitude permafrost is estimated to contain ~1500 gigatons of organic carbon (1), some of which may be released as CH4 and/or CO2 upon thawing. However, because the response of permafrost to climate change is mediated by a complex array of site-scale factors, models of 21st century near-surface permafrost degradation have high uncertainty (2). Furthermore, deep permafrost more than a few meters below the surface hosts a sizeable carbon pool (1), and even less is known about the response of this deep permafrost to climate change.

An alternative way to assess the vulnerability of permafrost to warming is to reconstruct its response during previous interglacial periods (35). The fate of permafrost and associated carbon stocks during marine isotope stage (MIS) 5e [the last interglacial; ~125 thousand years (ka) ago] and MIS 11 (~400 ka ago) is of particular interest because these were likely among the warmest interglacials of the past 1.5 million years (6, 7), with strong cryosphere response to warming expressed globally as substantial sea-level rise (8). Unexpectedly, Vaks et al. (9) recently found that ground thaw in northern Siberia was instead more common during interglacials before MIS 11 and restricted to southern Siberia since then based on episodes of cave deposition. They attributed this heightened vulnerability of permafrost before 400 ka ago to the absence of year-round Arctic sea ice (10, 11), which increased oceanic heat and moisture fluxes to surrounding landmasses, driving deeper ground thaw in summer and greater insulation under thicker snow packs in winter. However, these findings come from only a few parts of Siberia and seem to differ from patterns in sub-Arctic northwest Canada, where relict ice wedges at two sites indicate that permafrost has persisted for at least the past ~740 ka (3, 12).

We identify intervals of northern permafrost thaw through uranium-thorium (U-Th) dating of previously undated speleothems (130 ages on 72 speleothems) from three regions of Canada that today are in distinct permafrost zones: the southern Canadian Rockies (isolated zone; <10% permafrost coverage), Nahanni National Park in the Northwest Territories (discontinuous zone; 50 to 90% permafrost coverage), and northern Yukon (continuous zone; >90% permafrost coverage) (Fig. 1). Samples come from at least five caves in each region, representing varying topographic settings, cave depth and geometry, and glacial histories (Supplementary Text). The caves are on exposed slopes with minimal unconsolidated sediment cover and are thus likely among the first parts of the permafrost landscape to thaw with warming. The speleothems are predominantly small and commonly display visually sharp changes in coloration in the layered calcite that may indicate different growth intervals separated by hiatuses. Most samples are flowstones, with fewer stalactites and stalagmites.

Fig. 1 Current permafrost extent and sampled caves.

Permafrost zones in North America (purple shading), maximum Pleistocene glaciation extent (blue hatched line), and locations of caves and number of speleothems analyzed.

Although other factors can prevent speleothem growth locally despite thawed conditions (13, 14), permafrost is the first-order control on widespread presence or absence of speleothem growth in northern karst regions (5, 15). Calcite precipitation is suppressed in caves in permafrost regions because of lack of liquid water, reduced temperature, low groundwater CO2, and/or poor cave ventilation associated with ice plugs (fig. S1) (14). Thus, speleothems in contemporary permafrost regions are relicts of past thaw that enabled meteoric waters to seep into caves and deposit calcite, and have received increased attention as a paleo-permafrost proxy (5, 9, 1520). Speleothems can be precisely dated over the past 500 ka, spanning the last five glacial-interglacial cycles, with U-Th geochronology using inductively coupled plasma mass spectrometry (ICP-MS) (21).

Earlier U-Th dating (15, 16) of 64 speleothems from our caves using low-precision alpha spectrometry found growth associated with the Holocene and past three interglacials in the southern two areas, suggesting recurring interglacial permafrost thaw in the isolated and discontinuous permafrost zones (Fig. 2, C and E). In contrast, the continuous zone yielded only a few speleothems from recent interglacials, and these were collected from cave entrances subject to seasonal thaw, which thus cannot be attributed to deep ground thaw and are excluded in our discussion below. A large fraction of the speleothems from all these regions exceeded the ~350-ka limit of the dating technique at that time, which precluded assessment of MIS 11 permafrost history. In addition, a larger speleothem dataset is needed for a more robust characterization of the presence versus absence of growth in these caves during past interglacials. The 72 speleothems that we dated with ICP-MS were not analyzed in these previous studies (the speleothems from these earlier studies are unfortunately lost) and represent an essentially random sample, with speleothems collected where available throughout the caves. All speleothems are from at least several meters beneath the surface. We dated the outermost laminae from the speleothems to determine the youngest episode of growth and frequently obtained additional dates from stratigraphically older laminae to check for additional growth episodes on speleothems that yielded finite U-Th ages or to confirm the nonfinite ages of older speleothems beyond analytical limits of the method (figs. S2 to S4).

Fig. 2 Previously published and new Canadian speleothem ages.

(A and B) Benthic δ18O stack. (C and D) Histograms of Canadian speleothem growth based on U-Th dates, by MIS from isolated (light blue), discontinuous (blue), and continuous permafrost zones (dark blue). Results from previous studies using alpha spectrometry (15, 16) are shown in (C), and new results using ICP-MS are shown in (D). The vertical error bars reflect uncertainty only in the upper limit of the number of speleothems that grew during an interval because the age uncertainties of some samples overlapped multiple interglacials as shown in (E). Those samples are represented in the error bars for each interval they overlap. (E and F) Individual U-Th ages with 2σ uncertainties from previous studies (E) and this study (F). Nonfinite ages are shown by arrows.


Except for a few (n = 6) younger speleothems from the southern Canadian Rockies, all speleothems are ≥400 ka, and most (n = 59) exceed 500 ka, the limit of U-Th methods used in this study (Fig. 2, D and F). Earlier alpha-spectrometric dating obtained similar results to us in northern Yukon (15) but found a somewhat larger fraction of speleothems from recent interglacials than we did in Nahanni and the southern Canadian Rockies (Fig. 2) (16). This difference may be due to chance and the large uncertainties on many of the samples near the limit of alpha spectrometry at the time (Fig. 2). In any case, including these earlier results alongside ours does not affect our conclusions.

No ICP-MS–dated speleothem growth intervals from caves in the discontinuous (35 ages on 18 speleothems at six caves) or continuous (61 ages on 29 speleothems at seven caves) permafrost zones date to MIS 5e (Fig. 2, D and F), a period characterized by Arctic warming, sea ice retreat, and Greenland Ice Sheet mass loss (8). These results bolster the earlier alpha-spectrometric dating of 43 speleothems from these areas, which found no interior cave growth in the continuous permafrost zone and only a single sample in the discontinuous permafrost zone from MIS 5e (Fig. 2, C and E) (15, 16). The speleothem data are thus generally consistent with stratigraphic evidence for persistence of deeper permafrost through MIS 5e, despite widespread near-surface thaw, in the discontinuous zone of Alaska and Yukon (4, 22, 23).

The interglacial at MIS 11 stands out in our dataset, with speleothem deposition taking place in all three regions (Fig. 2, D and F, and table S1). Speleothem growth occurred in isolated permafrost of the southern Canadian Rockies, as it did during later interglacials (16). In the discontinuous zone, five speleothems from three caves have MIS 11 growth intervals, implying deep regional thaw of discontinuous permafrost at that time. One speleothem from the continuous permafrost zone dates to MIS 11 (419 ± 11 ka ago, error-weighted mean ± 2σ of four dates from the same growth layer), and notably, it was collected from the rear of the deepest cave in the area (fig. S5).

Speleothem 234U/238U ratios suggest that cave deposition was more common across our sites earlier in the Pleistocene than during the interglacials of the past 500 ka. The majority of speleothems from all three permafrost zones have 234U activities that are out of equilibrium (at 3σ) with the parent isotope 238U (fig. S6). Initial 234U/238U disequilibrium would generally decay away within 1500 ka, about six half-lives of 234U (t1/2 = 245 ka). Thus, the ubiquity of samples with residual disequilibrium represents substantial speleothem deposition during one or more intervals in the ~500- to 1500-ka window, in marked contrast to the paucity of MIS 11 and 5e speleothems from both the continuous and discontinuous permafrost zones (Fig. 3, A and B).

Fig. 3 Arctic speleothem growth and paleoclimate records.

(A and B) Histograms of Canadian (A) [this study; (15, 16)] and Siberian (B) (5, 9) speleothem growth based on U-Th and U-Pb dates and 234U/238U disequilibria by MIS and other age intervals (top) from isolated (light blue/purple), discontinuous (blue/purple), and continuous permafrost zones (dark blue/purple). Canadian histogram (A) shows alpha spectrometry results from previous studies [faint colors with no outlines (15, 16)] combined with our new ICP-MS results (solid colors with black outlines). 234U/238U uncertainties on nonfinite alpha spectrometry ages >350 ka from previous studies are generally too large to determine which interval these speleothems grew during (fig. S6), so they are represented by faint error bars in the MIS 11, 500 to 1500 ka, and >1500-ka bins to show the full range of possibilities. Note that some speleothems grew during multiple intervals. (C) Abundance of benthic ostracod Polycope in western Arctic, a genus signifying high local surface productivity in sea ice marginal areas and thus the presence of sea ice (10). (D and E) Atmospheric CO2 (D) and CH4 (E) concentrations from Antarctic ice (0- to 800-ka lines and open circles at 950 and 1500 ka) (31, 32, 46, 47) and the foraminifera δ11B proxy (circles and lines from 1050 to 1250 ka and 1375 to 1500 ka) (37, 38, 48). ppm, parts per million; ppb, parts per billion. (F) Global mean surface temperature (7).


Our Canadian speleothem growth history is similar to an analogous record from Siberia that is based on U-Th and U-Pb dating in a north-south transect of caves (5, 9). The northernmost caves in continuous permafrost from both areas are dominated by speleothem deposition in the ~500- to 1500-ka interval, a single speleothem from MIS 11, and none thereafter, while the southern caves in discontinuous and isolated permafrost contain speleothems from the most recent four interglacials (Fig. 3, A and B). It will be important to determine with U-Pb dating whether the northernmost speleothems in Canada were deposited throughout the 500- to 1500-ka interval, as was the case in northern Siberia, or whether they grew only during a few interglacials. In any event, these data are consistent with intervals of more widespread deep thaw across the high Arctic earlier in the Pleistocene, with subsequent long-term cooling and development of perennial Arctic sea ice stabilizing the continuous permafrost and restricting deep thaw to the sub-Arctic since MIS 11 (Fig. 3, C and F) (9). Other elements of the Arctic cryosphere, such as the Greenland Ice Sheet, also appear to have been more dynamic in the early Pleistocene before largely stabilizing in the past million years except during particularly warm interglacials (8, 24, 25). Notably, temperature (7) and sea ice (10, 11) conditions associated with these earlier intervals of widespread permafrost thaw may have been similar to those projected for later this century (26).

The amount of permafrost thaw represented by the single MIS 11 speleothems from northern Yukon and northern Siberia (5, 9) is open to interpretation. Speleothem deposition is, by nature, patchy within caves, even in tropical regions with contemporary speleothem formation (27). Thus, these two speleothems could imply deep regional-scale thaw in the continuous permafrost zone during MIS 11 (Fig. 2B). However, a more modest interpretation is that they reflect more limited, spatially restricted thaw, perhaps related to site-specific factors [e.g., fractures in carbonate bedrock and karstic dissolution features (28) and local taliks beneath depressions that concentrate surface water (5)] that could have enabled local cave seepage and speleothem formation, while deposition at other nearby caves, or even elsewhere within the same cave, was inhibited by permafrost. This view is supported by the lack of any other MIS 11 speleothem deposition in our seven sites from the continuous permafrost zone, representing 61 U-Th ages on 29 speleothems. Furthermore, a larger proportion of MIS 11 ages would be expected in northern Yukon if widespread speleothem growth last occurred then because of the preferential loss of older speleothems through overgrowth, frost shattering, clastic infilling, dissolution, collapse of cave roofs, and retreat of cave entrances (15, 16, 27). Speleothems from recent interglacials are more common in our southern two areas. Last, the limited deep thaw that we infer in the northern Yukon caves would imply that much of the continuous permafrost remained intact, given that the caves are on exposed hillslopes that are among the most sensitive parts of the landscape to warming.

Regardless of the degree of permafrost thaw in the modern continuous zone during MIS 11, the deep thaw of discontinuous permafrost at that time documented here, coupled with widespread shallow thaw thought to typify warm interglacials inferred by others (4, 23), would have presumably led to substantial release of greenhouse gases from respiration of thawed organic carbon. Models suggest that similar carbon release occurred with permafrost thaw during the last deglaciation (29) and is expected with projected warming over the next several centuries (1, 30). However, atmospheric CO2 and CH4 concentrations did not exceed Holocene levels during MIS 11 (31, 32) (Fig. 3, D and E), seemingly indicating that widespread ground thaw did not initiate a permafrost carbon feedback where increased atmospheric greenhouse gases and subsequent warming thaw permafrost that releases more greenhouse gases. Similarly, isotopic measurements on Antarctic ice core CH4 suggest minimal deglacial CH4 emissions from old carbon reservoirs, such as thawing permafrost (33).

This apparent conundrum could be explained in several ways. Perhaps there was a smaller stock of stored carbon in the permafrost at the start of MIS 11 relative to today. It is also possible that permafrost carbon release was slow enough to be buffered by oceanic uptake and, thus, not expressed in ice core records of atmospheric composition, a scenario supported by evidence from marine sediments for substantial carbonate dissolution during MIS 11 (34), or permafrost carbon release could have been balanced by lower emissions from other sources and/or increased uptake by the terrestrial biosphere, peatlands, and soils in newly deglaciated regions (35, 36). Alternatively, the thaw of deep permafrost inferred here from our exposed upland caves may not have been associated with substantial carbon release. Carbon-rich valley-bottom permafrost insulated under thick organic deposits may have instead remained perennially frozen, leading, for example, to persistence of fine-grained ice-rich syngenetic permafrost (i.e., yedoma) in central Yukon through MIS 11 and subsequent interglacials (3, 12).

More broadly, the general similarity of atmospheric CO2 and CH4 concentrations during warm interglacials of, respectively, the past ~1500 ka (37, 38) and ~800 ka (31, 32), despite differences in the spatial extent of deep permafrost thaw that we infer for the last 1500 ka, suggests that greenhouse gas concentrations during these interglacials were relatively insensitive to permafrost thaw (Fig. 3, D and E). Interglacial greenhouse gas concentrations became slightly higher following the Mid-Brunhes Transition at ~430 ka ago, even as permafrost persistence during these interglacials became more common, consistent with ice core evidence for minimal permafrost contribution to recent glacial-interglacial CH4 increases (33, 39). The persistence of deep permafrost through the past several glacial-interglacial cycles also implies long-term sequestration of a carbon pool preserved in permafrost that aggraded syngenetically with sediment accumulation over this time (3, 12, 4042). Although speculative, perhaps such sequestration of carbon not released back to the carbon cycle during warm interglacials contributed to the longer, deeper glaciations of the 100-ka world by slightly altering the thresholds and feedbacks needed to kick-start deglaciation.

These hypotheses highlight the complexity of the long-term interactions between climate, permafrost, and the carbon cycle (35, 36, 43) and the challenges of extrapolating interpretations of past permafrost dynamics to the future, given the different time scales of forcings, initial conditions, and local effects. Studies of contemporary permafrost dynamics and carbon cycling are important for evaluating potential contribution of the permafrost-carbon climate feedback to anticipated warming over the next several centuries in our high CO2 world (1, 30), but better understanding of past permafrost response to persistent warming—through field studies of ancient permafrost (3, 40, 44), development of new proxies for permafrost degradation (22), and numerical modeling of long-term permafrost-carbon feedbacks—is also needed to help resolve the carbon cycle conundrums identified here and thus clarify the role of permafrost and associated carbon stocks in past and future climate change.


Speleothem collection

Speleothems from the southern Canadian Rocky Mountains were collected by D. Ford, R. Harmon, P. Thompson, and S. Worthington in 1968–1974 and 1984–1985. Samples from the Nahanni area were collected by D. Ford, R. Harmon, and G. Brook in 1972–1973. B. Lauriol collected the northern Yukon speleothems between 1985 and 1995.

U-Th dating

U-Th samples were processed at Massachusetts Institute of Technology (MIT) in a clean laboratory equipped with HEPA-filtered air and ULPA-filtered laminar flow hoods. Samples were first dissolved in nitric acid along with a U-Th isotope tracer (229Th-233U-236U). After converting to hydrochloric acid, iron was added in solution, and the pH was raised by the addition of ammonium hydroxide to precipitate iron oxides that scavenge U and Th from solution. Calcium and other divalent cations, which primarily remain in solution at this step, were then removed by repeated centrifugation. Samples were then dissolved in nitric acid and then purified by ion exchange chromatography using custom-made columns with 500-μl bed volumes containing AG 1-X8 resin (100 to 200 mesh; Cl form) and a polyethylene frit. Following elution of iron and divalent cations with 8 M nitric acid, Th was eluted with 6 N hydrochloric acid, and U was eluted with water. Samples were then converted to 0.3 N nitric acid for analysis.

The U-Th isotope tracer was prepared by mixing a 229Th solution with the 3636a 233U-236U solution from the Institute for Reference Materials and Measurements. It was calibrated by mixing with U-Th gravimetric standards prepared from solid uranium (CRM 112a) and thorium (Ames Th bar) metals at MIT, followed by measurements of 238U/236U and 232Th/229Th in the combined solution. During the course of this study, four aliquots of the GC-1 secular equilibrium standard (21) were prepared and analyzed using the same methods used for samples in this study, giving a mean 230Th/238U activity ratio of 1.001 ± 0.002 and a mean 234U/238U activity ratio of 1.0000 ± 0.0008.

Samples were analyzed on a Nu Plasma II-ES multi-collector ICP-MS at MIT equipped with an Aridus II desolvating nebulizer and a Savillex C-Flow PFA micronebulizer (100 μl/min). For most samples, uranium isotopes were measured in a static analysis with 238U, 236U, 235U, and 233U measured on Faraday cups and 234U measured on a secondary electron multiplier (SEM). The Faraday cup measuring 233U was equipped with a 1012-ohm resistor, while the other cups had a 1011-ohm resistor. Samples were bracketed by measurements of unspiked CRM 112a measured with the same detectors. Measurements at masses 237, 236.5, 234.5, and 233.5 on the SEM were performed using a dynamic routine after each on-mass measurement. Signals at 237 and 236.5 were used to fit a power law to the tail from 238U; this relationship was then applied to estimate the down-mass tails from 236U and 235U and to remove the effects of tailing on all masses except 238U. Tailing at 234U was consistently 2 to 4 per mil of the total 234U intensity. Once per day, upmass tailing and hydride formation were measured by measurement at masses 238.5, 239, and 239.5. The combined upmass tail and hydride were consistently 1 × 10−6 to 2 × 10−6 relative to the 238U beam. This signal was multiplied by the 233U intensity and removed from the 234U beam for each sample. Ion counter yield was determined using 234U/238U ratios measured in bracketing CRM 112a standard measurements after correction for tailing as described above and correction for mass bias using measured 238U/235U ratios in bracketing standards. Backgrounds for both on-mass and tailing measurements were determined by measurements of blank acid using the same detector array before each analysis.

Thorium isotope analyses were performed in static mode with 232Th measured in a Faraday cup with a 1011-ohm resistor, 230Th measured in an SEM, and 229Th measured in a Faraday cup with a 1012-ohm resistor. Measurements were bracketed by measurements of an internal 232Th-230Th-229Th standard (MITh-1) calibrated by measurement of an admixed 233U-236U solution (IRMM 3636a). Standard analyses were used to determine mass bias and ion counter yield. Tailing from 232Th on 230Th was measured once per day using dynamic measurement of masses 230.5 and 229.5 on the SEM; as tailing from 232Th on 230Th was consistently <0.3 parts per million (ppm) and we excluded samples with 230Th/232Th ratios <100 ppm, variations in this tailing comprised a minor portion of uncertainty for our samples.

For a small number of samples, all isotopes of U and Th fractions were measured using Faraday cups in a static analysis to reduce uncertainties associated with ion counter yield variability. For uranium, 238U was collected in a Faraday cup equipped with a 1010-ohm resistor, and 234U was collected in a Faraday cup equipped with a 1012-ohm resistor. All other cups had 1011-ohm resistors. For thorium, 230Th was measured in a cup with a 1012-ohm resistor. Uranium and thorium analyses were bracketed by the same standards as above but with higher concentrations. Tailing of higher masses on 234U was estimated from 234U/238U ratios in bracketing standards after correction for mass bias.

Uncertainties reported for isotope ratios incorporate uncertainties associated with background, tailing, mass bias, ion counter yield variability, spike composition, and procedural blanks. Total procedural blanks averaged <5-pg 238U and <0.05-fg 230Th (as compared to typical U contents of 50 to 500 ng and 230Th contents of >1 pg per sample). U-Th ages were calculated using the 238U decay constant of Jaffey et al. (45) and the 234U and 230Th decay constants of Cheng et al. (21). U-Th ages were corrected for initial 230Th assuming an initial 230Th/232Th ratio of 4.4 (±2.2) × 10−6 (atomic); samples with 230Th/232Th < 100 × 10−6 were excluded. U-Th age uncertainties incorporate analytical, initial 230Th, and decay constant uncertainties and were calculated with a Monte Carlo routine using MATLAB software. Because of the SEM methods used for most samples in this study, we consider the upper limit of data in which we place high confidence to be 500 ka.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank T. Schuur, C. Yonge, S. Ewing, and G. Zazula for discussions and C. Riggs for ArcGIS support. D.C.F. and B.L. thank the people of Old Crow (Yukon) for support during field work. Thoughtful reviews by I. Orland and an anonymous referee strengthened the manuscript. Funding: Research was supported by NSF ARC-1607816 and 1607968, NSERC Discovery Grants, and the Polar Continental Shelf Program. Author contributions: J.D.S., D.M., and C.I.W. designed the study. D.M. oversaw laboratory work. N.B.-C. performed sample analyses with B.H. and I.T. A.V.R. contributed to data interpretation. D.C.F. and B.L. investigated caves and collected samples. All authors contributed to the preparation of the manuscript. 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.

Stay Connected to Science Advances

Navigate This Article