Source: https://www.nature.com/articles/s41467-019-09622-y?error=cookies_not_supported&code=16ecc016-b37c-471d-89fa-daa3a32575ec
Timestamp: 2019-04-20 08:36:58+00:00

Document:
Eastern Beringia is one of the few Western Arctic regions where full Holocene climate reconstructions are possible. However, most full Holocene reconstructions in Eastern Beringia are based either on pollen or midges, which show conflicting early Holocene summer temperature histories. This discrepancy precludes understanding the factors that drove past (and potentially future) climate change and calls for independent proxies to advance the debate. We present a ~13.6 ka summer temperature reconstruction in central Yukon, part of Eastern Beringia, using precipitation isotopes in syngenetic permafrost. The reconstruction shows that early Holocene summers were consistently warmer than the Holocene mean, as supported by midges, and a thermal maximum at ~7.6–6.6 ka BP. This maximum was followed by a ~6 ka cooling, and later abruptly reversed by industrial-era warming leading to a modern climate that is unprecedented in the Holocene context and exceeds the Holocene thermal maximum by +1.7 ± 0.7 °C.
The Arctic has warmed faster than any other region this past century, owing to climate-cryosphere feedbacks1. This warming is driving a range of biotic and abiotic changes that directly impact northern communities and ecosystems2,3, as well as enhanced C fluxes to the atmosphere from thawing permafrost4 and ice loss from glaciers and ice caps5, which have implications for global climate and sea level, respectively. The broad-reaching effects of Arctic change highlight the need to advance our knowledge of the Arctic system.
Current knowledge of the Arctic system is informed by a sparse network of relatively short (~50–100 years) instrumental records. However, the Arctic is transitioning to a warmer state that appears to have no analogue in the historical record. Natural proxy archives (e.g. ice cores and lake sediments) offer an alternative source of data on past climate and environmental change extending back thousands of years, with analogues from the distant past that may be especially relevant for anticipating the form of a future, warmer Arctic6. Arctic proxy reconstructions for the Common Era (CE, last 2 ka)7 and Holocene (last ~11.7 ka)8,9 broadly agree that the late Holocene was characterised by a long cooling trend that ended with abrupt warming during the industrial era. Many proxy reconstructions agree that industrial-era warming was exceptional in overall rate, but not magnitude in the Holocene context10,11, while others suggest that industrial-era warming has been exceptional in both regards6,12,13. Differences can be explained by spatial heterogeneity in the climate field or proxy-related noise (e.g. chronological error, secondary climate sensitivities and differences in seasonality). Noise is inherent to all proxies, but recovery of the climate signal can be optimised through replication and a diverse, multi-proxy approach.
Eastern Beringia, the contiguous area of Yukon and Alaska that was ice free during the last glaciation, is one Arctic borderland where the Holocene climate history is still debated. A review by Kaufman et al.9 showed that nearly all quantitative, summer-specific temperature reconstructions from Eastern Beringia with early Holocene coverage are based on two proxy types—pollen and midges (Chironomidae) in lake sediments. Regionally, pollen and midges agree there was a climate optimum at ~7–5 ka BP, followed by a long-term cooling. However, these proxies show large differences before 8 ka BP with pollen indicating that early Holocene summers were persistently colder than the Holocene mean (ΔT range = –3.3 to –0.5 °C; mean ΔT = –1.6 °C), and midges showing that early Holocene summers were similar to the Holocene mean (ΔT range = –0.1 to + 0.4 °C; mean ΔT = + 0.1 °C). This discrepancy confounds our understanding of the drivers of early Holocene climate change and highlights the need for independent proxies to resolve this debate.
Precipitation isotope ratios—2H/1H (δDprecip) and 18O/16O (δ18Oprecip)—in ice core records have a long history of use as a temperature proxy14,15,16 and have played a crucial role in shaping our knowledge of Quaternary climates in glaciated locales. More recently, this concept has been applied to Pleistocene and Holocene ice wedges (snowmelt derived) in permafrost environments in parts of Siberia12,17 and Eastern Beringia18,19,20 to reconstruct winter δprecip (δDprecip or δ18Oprecip) and temperature. Relict pore ice in syngenetic (aggrading) permafrost also shows potential as a δprecip proxy19, although there are few pore ice studies with a paleoclimate focus, and only one example from Siberia21 where a Holocene δpore ice series was developed and used for qualitative (warmer vs. cooler) temperature inferences.
In this study, we present the first full Holocene, quantitative temperature reconstruction from a sequence of relict pore ice in syngenetic permafrost from a soligenous peatland called DHP174 (Dempster Highway Peatland, near 174 km marker) in central Yukon (Fig. 1a), part of Eastern Beringia. Surface peat accumulation is the main driver of permafrost aggradation at this site, which simultaneously archives a warm-season blend of precipitation (integrated as pore water) that reaches the base of the seasonally thawed active layer in the summer. This reconstruction is enabled by a significant positive association between local air temperatures and precipitation isotope ratios (Fig. 1b; Supplementary Fig. 1). Our pore ice record is conceptually similar to an ice core record in that it preserves a stratigraphic profile of ancient precipitation, but distinct because of its seasonal bias toward summer precipitation and a fractionation that occurs during freezing and causes relict pore ice to be enriched in δD and δ18O relative to the initial pore waters22. Based on this record, we develop a full Holocene δDprecip (pore ice)-based summer ΔT reconstruction, and an abiotic perspective to advance our understanding of early Holocene climate in Eastern Beringia.
The DHP174 peat-sediment composite record (5.32 m composite depth) is the product of two core profiles that overlap stratigraphically and includes three lithostratigraphic units (Fig. 1c) that were deposited over the last ~15.6 ka (see Methods). Unit 1 (5.32–4.78 m) is silty loess and was deposited between 15.6 and 14.3 ka cal BP. Unit 2 (4.78–4.38 m) is a graminoid-rich silty unit that aggraded relatively slowly (~0.1 cm decade−1 average) from 14.3 to 11.2 ka cal BP. Unit 2 overlaps with the Younger Dryas interval (~12.9–11.7 ka cal BP), when the climate of northern Alaska was much colder than today18,23. Finally, Unit 3 (4.38 m to the surface, 10.5 ka cal BP to present) is primarily Sphagnum moss, reflecting an aggrading peatland phase. The peat sequence accumulated at a relatively stable rate of ~0.4 cm decade−1, with the exception of higher rates at the time of peatland inception (4.38–4.13 m, ~2 cm decade−1) and a brief slowdown near the core top (0.21–0.15 m, ~0.1 cm decade−1) during the period ~1000–200 cal yr BP.
Peat fibres were exceptionally well preserved throughout the profile (Supplementary Fig. 2a–c), which suggests a relatively uninterrupted aggradation of permafrost since inception, and a syngenetic origin of relict pore ice in this sequence. Downcore ice content is primarily structureless pore ice (Supplementary Fig. 2) representing pore waters that froze in situ24. Rare ice lenses (typically <0.5 cm thick; Supplementary Fig. 2b) were encountered in Unit 1 (n = 7) and Unit 3 (n = 2) but were not sampled. Ice lenses represent active layer waters drawn to the freezing front by cryosuction24, a process that is common in fine-grained sediments and may explain their higher number in Unit 1. Pore ice and lens ice derive from potentially different active layer water and cannot be directly compared. In this study, we focus exclusively on the dominant form of ice—relict pore ice.
In contrast to the ~15.6-ka-long peat-sediment chronology described above, the relict pore ice hosted in the same peat-sediment sequence spans the last ~13.6 ka cal BP. The shorter relict pore ice chronology relates to the maximum depth at which active layer pore ice forms and becomes ‘relict’, which equals the maximum active layer thickness (ALTmax). ALTmax is defined as the thickest active layer likely to occur in a given climatic period and is ~58 cm for the study site today (Supplementary Note 1). Because relict pore ice forms below the surface it will always be younger than the host sediments. For example, ca. modern relict ice (–68 cal yr BP; 58 cm below ground) is hosted in fossil peats that grew at the surface ~2065 cal yr BP (i.e. the depositional age of the peats). Relict pore ice in the basal sediments is likewise ~2 ka younger than the depositional age of the hosting sediments (see Methods for full details on the relict pore ice chronology).
Relict pore ice δD and δ18O series from the upper and lower core profiles are compared and combined in a composite series in this section. For the upper core, relict pore ice includes all pore ice from depths greater than ALTmax or ~58 cm (i.e. excluding all active layer pore ice). For the lower core profile, we consider all pore ice at or below a composite depth of 290 cm (~90 cm below the lower coring surface) to be relict pore ice. The deeper cut-off depth to relict pore ice in the lower core avoids pore ice that has been overprinted by modern waters following disturbance of the lower surface (Supplementary Fig. 3).
The upper and lower relict pore ice records collectively span the last ~13.6 ka cal BP. The two records show small systematic differences in their mean δD and δ18O values during the period of overlap from 8.5 to 6.4 ka cal BP (composite depths 380–290 cm), the upper core being relatively depleted compared to the lower core in δD and δ18O (Fig. 2a, b). There is no offset in dexcess (dexcess = δD – 8 × δ18O) (Fig. 2c) The δD and δ18O differences are small in absolute terms (the mean offsets are 1.35‰ and 0.22‰ for δD and δ18O, respectively) and are not unexpected for different sites across a landscape given the same source water. As we discuss below, relict pore ice is enriched in heavy isotopologues relative to the initial pore waters due a fractionation that occurs during the freezing process22. The overall enrichment (εice-water) depends on the rates of freezing and molecular diffusivity. For our site, the lower active layer is saturated in the late summer and, thus, the rate of freezing will be buffered by latent heat effects25. However, water table depth and total latent heat content are unlikely to be spatially homogenous due to variations in microtopography and slope, and therefore the rate of freezing and εice-water are also unlikely to be homogenous. Because the two core locations (~50 m apart) receive the same precipitation and runoff, spatial variability in εice-water is a probable explanation for the small offsets in δpore ice.
Regardless of origin, the offsets are systematic and the cores show coherent variations at 102–103 year timescales (Fig. 2a, b), which suggests that both cores are in fact archiving the same source water. To enable a direct comparison of the two series in a composite record, we adjusted the upper core δDpore ice and δ18Opore ice series by adding a constant equal to half of the mean δD and δ18O offsets, and the lower core δDpore ice and δ18Opore ice by subtracting half of the mean δD and δ18O offsets (Fig. 2a, b). This aligns both series to a common mean value during their period of overlap without modifying the individual trends, and, thus, the composite δDpore ice and δ18Opore ice records (Fig. 2d, e) also accurately depict the raw δpore ice trends. These adjustments are small in absolute terms and do not influence trends in dexcess (Fig. 2c–f) nor conclusions on seasonality of relict ice, which we based on raw δpore ice (next section).
The composite record covers the last ~13.6 ka with an average time-step of ~158 years between pore ice samples since 10.5 ka cal BP; however, from 13.6–10.5 ka cal BP, the average time-step is relatively coarse (~1020 years between samples) due to low accumulation rates just prior to the peatland initiation.
Co-isotope (δD and δ18O) measurements of relict pore ice plot tightly on the Local Meteoric Water Line (LMWL; Fig. 3a), indicating that the pore ice is derived from meteoric waters that froze in situ rather than evaporatively enriched soil waters, which would plot below the LMWL along a shallow evaporation line26. However, relict pore ice older than 5 ka BP plots above the modern LMWL (Fig. 3b), which may reflect changing Holocene boundary conditions as we discuss in the next section. The remainder of this section focuses on seasonality of ca. modern (top-of-permafrost) δpore ice.
Relict pore ice is ‘first ice’ (the first ice to form in a bottom-up or top-down freezeback sequence) that once formed at the bottom of a paleo-active layer. First ice is enriched in heavy isotopes relative to the initial pore waters22. Based on an analysis of DHP174-13 active layer ice, the enrichment is ~5.8‰ for δD and ~0.8‰ for δ18O, and is not expected to vary significantly through time due to complacent freezeback conditions in a saturated active layer and unchanging peat lithology for most of the record (Supplementary Note 2). To establish the δprecip seasonality of ca. modern δpore ice, we account for these enrichments to estimate the initial δpore water, which can be compared to the δprecip climatology of the nearest GNIP site at Mayo, Yukon.
Mean δDpore water and δ18Opore water estimates for the uppermost 3 cm of DHP174-13 relict pore ice (58–61 cm; representing the last ~3 decades of syngenetic permafrost) are −161.5 and −20.2‰, respectively, with a mean dexcess value of +0.1‰. Compared to local meteoric waters from Mayo, our δDpore water estimate is enriched relative to mean annual δDprecip and within range of thaw season (May–September) δDprecip values (Fig. 4a). The pore water dexcess estimate is also biased toward warm-season precipitation values (i.e. depleted relative to annual precipitation dexcess and enriched compared to the warmest months; Fig. 4b). Both comparisons suggest relict pore ice derives mainly from warm-season precipitation. However, a caveat of these comparisons is that the temporal integrations of the Mayo δDprecip record (1985–1989, variable replication depending on month) and ca. modern pore ice (last ~3 decades) are unequal, which may impact how closely they agree. Mayo δDprecip data generally follow expectations of more depleted values during cold winter months and enriched values during warm summer months (as in Fig. 1b), but discrepancies do occur (e.g., peak δDprecip in September, which is not the warmest month), which are likely attributable to the limited replication in the Mayo record. The proposed seasonality of relict pore ice can also be corroborated by comparing δpore ice-based temperature estimates (Methods) with the local climatology as constrained by climate data from the Ogilvie River station (Supplementary Note 3), ~17 km north of the DHP174 site. The average temperature estimate is 6.3 ± 3.4 °C (1σ) based on δDpore ice (Fig. 4c), and 7.0 ± 3.3 °C (1σ) based on δ18Opore ice (Supplementary Fig. 4c), which are the same within error. Both estimates are bracketed by mean May–September temperatures in the study area and provide further confirmation that relict pore ice derives primarily from warm-season precipitation.
A warm-season seasonality is also expected for climatological and hydrological reasons. Roughly 65% of the annual precipitation budget is delivered during the thaw season, the bulk of which is rain (Fig. 4d). Secondly, the typical active layer hydrology style of raised northern peat plateaus favours a warm-season seasonality. At the time of snowmelt, the active layer is largely undeveloped (e.g. 0–10 cm), which forces meltwaters to drain near the surface where saturated hydraulic conductivity (Ks) is high (10–1000 m day−1, 0–10 cm depth)27. Conversely, as the active layer thaws in the summer, precipitation percolates to greater depths and will drain slowly due to lower Ks in the deeper active layer (0.5–5 m day−1 below 20 cm)27. This hydrological style favours greater representation of warm-season precipitation in the lower active layer, where pore waters eventually freeze and become relict following permafrost aggradation.
Finally, we note that the composite record of pore ice (Fig. 2d, e) is characterised by relatively small inter-sample δpore ice variability (on average 0.4‰ for δD and 0.1‰ for δ18O) in comparison to the large inter-annual (weather-related) variability in δprecip that is evident in local GNIP records for most months (see δDprecip ± 1σ; Fig. 4a). The absence of extreme variability in δpore ice suggests that each pore ice sample represents a multi-year seasonal blend of precipitation, whereby weather-related fluctuations in δprecip are attenuated through the mixing of precipitation from several warm-season months, and multiple years facilitated by repeated freeze–thaw cycles. This is logical given the slow peatland aggradation rate (~0.4 cm decade−1), which dictates that up to ~20–30 years of mostly warm-season precipitation can mix in the lower active layer before the permafrost table aggrades 1 cm, thereby trapping pore ice situated at the former permafrost table as constrained by ALTmax. While ALTmax will not be reached every thaw season (Supplementary Note 1), it is the infrequent, warmer-than-normal summers that can reanimate and mix with pore ice that froze at intermediate thaw depths (<ALTmax) during prior, cooler summers. Thus, natural averaging of multi-year pore ice in the active layer ensures a climatic timescale integration (depending on the aggradation rate) of δprecip and representation of cool and warm summers in any given sample of δpore ice, which explains the lack of extreme inter-sample δpore ice variability noted above.
Over Holocene and longer timescales, the isotopic composition of seawater is variable. The composite δpore ice records were corrected for changes in δseawater (Fig. 2d–f; see Methods) in order to isolate δpore ice variability related to changes in climate. The δDpore ice (seawater corrected) record (Fig. 2d) shows a rising trend through the Pleistocene–Holocene transition reaching some of the most enriched values of the early Holocene at ~10 ka cal BP, then a lesser but still positive trend between 10 and 6.6 ka cal BP, followed by an ~6-ka-long negative trend leading to the most depleted values on record in the early CE (~2 ka BP), and last an abrupt trend reversal that began ~150 yr BP and has culminated in the most enriched δDpore ice values at present. It is possible this reversal started before 150 yr BP, but the 1000–150 yr BP portion of our chronology is not well resolved due to the slow peat accumulation during this period, as discussed above. The reason for the slow accumulation is not known but could reflect a local climate limitation to peat growth (e.g. more arid) during the late pre-industrial era.
The δ18Opore ice (seawater corrected) record (Fig. 2e) documents similar low- and high-frequency trends, but a more pronounced positive trend during the early Holocene between ~9.2 and 5 ka cal BP. This difference is reflected in the dexcess record (Fig. 2f), which transitions from high mean values (~4 ‰) at 9.2 ka cal BP to low mean values (~0‰) by 5 ka cal BP. dexcess did not change much since the mid-Holocene. Under today’s δprecip boundary conditions (i.e. as defined by the present-day LMWL), dexcess values >4‰ are common during the cold-season months (October–April), and so the dexcess change in our pore ice record could be interpreted as a change in the precipitation seasonality (e.g. ratio of warm-to-cold-season precipitation) integrated in pore ice; however, this interpretation would also require a constant LMWL, which the co-isotope results (Fig. 3b) show is an invalid assumption. A precipitation seasonality mechanism is also unlikely given the active layer hydrology of northern peat plateaus, as discussed above, which precludes snowmelt from the deep active layer where relict pore ice is created. The Holocene dexcess trend is attributed to a transition from a paleo-LMWL (with a more positive intercept compared to today) defined by early Holocene relict pore ice (>~9.2 ka cal BP) to the contemporary LWML defined by modern GNIP data, but which aptly describes late Holocene relict pore ice younger than 5 ka BP (Fig. 3b). Continental boundary conditions such as local orographic effects and mean weather trajectories define the slope of the LWML22, whereas the intercept of the LMWL (which points to the dexcess of the first evaporate) is influenced by the mean evaporative boundary conditions at the moisture source28. Given the similar slopes but different intercepts of the paleo vs. modern LMWLs, we conclude that a change in evaporative boundary conditions at the moisture source is the most likely explanation for the Holocene dexcess transition.
The dexcess parameter reflects non-equilibrium partitioning of HDO relative to the slower moving H218O isotopologue in the evaporate leaving the ocean. dexcess is sensitive to factors that influence the evaporation rate including sea surface temperature (SST) and relative humidity (RH)28, which has made it a useful proxy in ice core studies for reconstructing changes that are far removed from the ice core site14. The North Pacific is the primary moisture source for central Yukon today19. Assuming a stable moisture source through the Holocene, the observed dexcess trend could reflect a cooling in North Pacific SST and/or rise in RH in the moisture source region. SST proxies encompassing a large area of the North Pacific broadly agree that SSTs peaked around ~10 ka and cooled by ~1–2 K by the mid-Holocene29,30. However, spatial coherence diminishes in the late Holocene with some SST reconstructions, indicating warming and others cooling. The balance of paleo-SST evidence is inconclusive, but hints at a possible and spatially complex relation between dexcess and SST.
Alternatively, the dexcess trend could indicate a geographic shift in the moisture source region with potentially different mean SST or RH at the different localities. The initial timing of the dexcess trend is coincident with the final collapse of the Laurentide ice sheet (LIS) from ~11–7 ka BP31. The LIS strongly impacted atmospheric circulation and climate in western Eurasia32. Likewise, climate modelling studies suggest that the collapse of the LIS had a warming effect on E. Beringia33, although the dynamical response of summertime atmospheric circulation patterns in the North Pacific sector to LIS collapse is unknown.
The co-isotope results indicate a transition from one set of δprecip boundary conditions to another. Although the physical mechanisms behind this transition are not well known, further interrogation of possible changes—such as modified SST or RH at a constant moisture source, or a geographic change in moisture source (e.g. due to altered atmospheric circulation pattern in the summer months)—will benefit from greater replication of the paleo-SST field and isotope-enabled GCM (General Circulation Model) experiments.
Under modern δprecip boundary conditions, δDprecip and δ18Oprecip are both temperature-sensitive proxies (Fig. 1b; Supplementary Fig. 1). However, if the modern δprecip-temperature transfer functions (both of which are linear) were applied to the δDpore ice and δ18Opore ice records, they would yield different reconstructions given their contrasting trends (Fig. 2d, e)—and clearly both cannot be correct. As evidenced by the early to late Holocene LMWL transition (as discussed above), δprecip boundary conditions did change and were likely driven by changes in evaporative conditions at the moisture source. The slow-moving H218O isotopologue is especially sensitive to non-equilibrium effects driven by changes in evaporation rate28. By comparison, HDO is relatively insensitive to this effect. For that reason, we consider δDpore ice to be a potentially more reliable local temperature proxy and use it as the basis for our ΔT anomaly reconstruction (see Methods).
The ΔT reconstruction (expressed as anomalies w.r.t. Holocene mean; Fig. 5d) shows a strong warming trend during the Pleistocene–Holocene transition, and variable ΔT anomalies ranging from –0.5 to +0.5 °C during the early Holocene (10.4–9.6 ka cal BP) with a mean ΔT of +0.1 °C. Local summer ΔT reconstructions from pollen and midges show major differences during this interval, with pollen indicating cool ΔT values from –2.2 to –1.2 °C, and midges indicating warm values from +0.1 to +0.4 °C (Fig. 5c). This difference may be due to limitations of the modern analogue technique (MAT), which is the quantitative method used to estimate summer temperatures from the E. Beringian pollen records. As discussed by Kaufman et al.9, it is notable that qualitative interpretations based on pollen in this region have reached the opposite conclusion that the early Holocene was warm. A warm climate is also indicated by a diverse suite of qualitative proxy evidence in E. Beringia, including high frequencies of thaw lake and peatland initiation9; a multi-proxy ‘all-temperature’ composite record that depicts the dominant mode of variability of all-temperature-sensitive marine and terrestrial records (Fig. 5b); and warmer SSTs as discussed above. However, not all of these supporting proxies are necessarily sensitive to summer temperatures exclusively. Last, a brief thermal maximum at ~10 ka cal BP is also evident in a multi-proxy temperature record representing the northern extratropics (30–90°N) (Fig. 5a), suggesting this phase of early Holocene warmth was widespread. The balance of evidence favours a warm early Holocene. It remains unclear why MAT-based interpretations of pollen in this region, on average, suggest colder summers. Further studies are needed to understand this discrepancy, especially given the popular use of pollen data in Quaternary studies.
From 9.6 to 6.6 ka cal BP, our reconstruction shows a sustained warming with ΔT anomalies ranging from –0.1 to +0.4 °C, and an average ΔT of 0.2 °C (Fig. 5d). Local maxima in ΔT occur at 9.6 ka cal BP (+0.4 °C), 8.8 ka cal BP (+0.4 °C), 7.6 ka cal BP (+0.4 °C) and 6.8–6.6 ka cal BP (+0.3 °C). Individual ΔT anomalies should not be over-interpreted; however, it is notable that the +0.4 °C anomaly at 8.8 ka cal BP coincides with a deep thaw unconformity preserved in permafrost near Mayo, which is dated to 8.87 ± 0.2 ka BP and has been linked to a brief period of exceptionally warm conditions34. We also note that the later warm interval between 7.6 and 6.6 ka cal BP coincides with a period of Arctic treeline expansion in the Mackenzie Delta region (~400 km N.E. of DHP174) that began ~8 ka BP35. The 7.6–6.6 ka cal BP warm interval also matches a thermal maximum in the extratropical composite (Fig. 5a). The E. Beringia ‘all-temperature’ composite shows a slightly later thermal maximum at ~6.5–6.0 ka BP (Fig. 5b). The regional pollen composite shows warm summers (+0.47 °C) from 7.5 to 7.0 ka BP (comparable to the pore ice ΔT anomaly at ~7.6 ka BP), but slightly warmer summers (+0.55 °C) from 6.0 to 5.5 ka BP (Fig. 5c). Similarly, the midge composite shows warm summers (+0.19 °C) from 7.5 to 6.5 ka BP, but even warmer summers (+0.34 °C) from 6.0 to 5.5 ka BP. Differences in magnitude and timing of early-middle Holocene peak warmth between the summer-specific proxies (i.e. pollen, midges and δDpore ice) are relatively small, and to some extent are likely explained by the unequal spatial representation of the different records. Additional pore ice records across the region may improve inter-proxy coherence.
From 6.6 to 1.2 ka cal BP, our reconstruction shows a long-term cooling trend of –0.16 °C ka−1 (p = 6.82 × 10–5) (Fig. 5d). Prominent cooling trends over the same interval are documented in the E. Beringia multi-proxy composite and extratropical multi-proxy composite (Fig. 5a, b). The extratropical composite shows a cooling trend of the same magnitude (–0.17 °C ka−1) between 6.6 and 1.2 ka BP. Regional pollen and midge composites also show cooling trends through this interval, but at lesser rates of change (–0.04 and –0.05 °C ka−1, respectively; Fig. 5c). Differences in cooling rates between the records may be explained by seasonality or spatial (e.g. continentality) inequities.
From ~1200 to 150 cal yr BP our reconstruction shows a gradual warming trend. A similar warming trend over the last ~1–2 ka is observed in the E. Beringia all-temperature composite (Fig. 5b). However, a warming trend during this interval is relatively uncommon to CE temperature reconstructions globally36, and is not supported by the multi-proxy composites from 30 to 90°N (Fig. 5a) and 60 to 90°N (Fig. 5e), which show a cooling trend during this time. Due to the slow peat accumulation rates from ~1000 to 150 cal yr BP (~0.1 cm ka−1), the late pre-industrial era is not well resolved in our pore ice record; therefore, our reconstruction should not be used for detailed estimates of temperature change during this interval.
Finally, our reconstruction shows accelerated warming since ~150 cal yr BP, culminating in the warmest summer climate at present, which (based on the last ~1 decade of relict pore ice, or the uppermost 1 cm of relict pore ice) we estimate is 2.1 ± 0.7 °C warmer than the Holocene average or 1.7 ± 0.7 °C warmer than the peak ΔT anomalies during the early-middle Holocene thermal maximum (Fig. 5d). The average ∆T anomaly for the last ~30 years of relict pore ice (uppermost 5 cm of relict pore ice) is +1.9 ± 0.7 °C warmer than the Holocene mean. The last ~100 years of our record is especially well resolved (e.g. sub-decadal for the top 5 cm of relict pore ice) due to the finer pore ice sampling resolution through this interval, and the historic cryptotephra (1912 CE Novarupta-Katmai at 10 cm depth) and post-bomb 14C dating (−35 cal yr BP at 4–5 cm depth) that constrain this part of the chronology (Tables S1, S2). Our reconstruction shows ~2 °C of warming over the last 100 years, which matches the overall increase in observed May–September temperature anomalies from this region37 (Fig. 5d). Similar warming rates are evident in the extratropical and Arctic composites (Fig. 5a–e), and temperature estimates from the Agassiz ice core ‘melt’ record, Ellesmere Island10. However, the magnitude of warming in the industrial era vs. early Holocene is still debated. The extratropical multi-proxy composite suggests the modern climate has not yet exceeded early Holocene peak warmth11, while some land-only (North America and Europe) pollen-based temperature reconstructions suggest the last decade was ~1.2 °C warmer than the Holocene mean13. Retreating glaciers on Baffin Island in the Canadian Arctic have exposed rooted plant fossils with non-finite 14C ages, which indicates that recent climate warming at that locality is unprecedented since the last interglacial6. While these select examples do not resolve the debate over whether today’s climate is exceptional everywhere in the Holocene context, they do support the idea that the magnitude of early Holocene warming has varied spatially. At our site in easternmost Beringia, the present-day summer climate does appear to be exceptional in the last 13.6 ka cal BP, and our reconstruction suggests that summer temperatures in recent decades have consistently exceeded the Holocene Thermal Maximum since about the mid-20th century.
Variability in Holocene summer climate in this region owes to multiple factors. The centroid of the July–August insolation curves and subsequent decline is closely aligned with the low-frequency trend in our reconstruction (Fig. 5d–f), and implies that summer insolation played an important role in pacing summer temperatures. The importance of insolation in driving Holocene climate in this region is corroborated by numerical climate model simulations33. However, other factors likely had a moderating effect. For example, the flooding of the Bering Strait (as late as 11 ka BP38) and development of thaw lakes in Yukon and Alaska in the early–middle Holocene probably had a regional cooling effect on summer climate33, while the LIS collapse and increases in shrub and boreal taxa likely had a warming effect33,39. The integration of all relevant forcings and feedbacks likely explains some of the more complex variabilities in early Holocene climate that diverge from a simple insolation-based interpretation. Finally, the industrial-era warming in our record (~2 °C warming over the last century) is coeval with the rapid increase in global pCO2 (Fig. 5g), which is an important driver of 20th century warming at the global scale40. In addition to the pCO2 forcing, Arctic warming is amplified above the global average due to reductions in sea ice concentrations, especially in the Beaufort and Chukchi Seas (north of Alaska)41, which likely accounts for at least part of the 20th century warming trend in our record.
The rate and magnitude of recent climate warming has implications for the stability of ice-rich permafrost in E. Beringia, which is both widespread and sensitive to climate warming3. Cryostratigraphic evidence for regional failure of ice-rich permafrost in E. Beringia is linked to past hypsithermals, including thermokarst deposits in Alaska and Yukon dating to Marine Isotope Stage 5e42, and thaw unconformities and thaw lake initiation dating to the early Holocene in parts of Yukon and the Mackenzie Delta region34,43,44,45. However, the warming threshold for such a widespread permafrost response is not well known. Our finding that today’s summer climate (last ~30 years) exceeds the early–middle Holocene thermal maximum by +1.7 ± 0.7 °C implies that the threshold for ice-rich permafrost stability in this region has already been crossed. This finding is consistent with recent evidence of climate-driven acceleration of thaw slump activity3,46 and ice wedge network degradation47 in the Mackenzie Delta region and Western Arctic more broadly, which have been linked to climate warming. Continued warming and permafrost degradation in this region threatens to expose vast pools of frozen organic matter in near-surface permafrost to microbial activity and, if released to the global atmosphere as CO2 or CH4, has the potential to amplify future warming trends4.
Ice-rich peat cores were collected from a peatland we refer to as DHP174 (65.21°N, 138.32°W; Fig. 1a–f), situated in the continuous permafrost zone. The peatland is dominated by Sphagnum spp. moss, and includes sparse forb and lichen groundcover, and an open canopy of Picea mariana trees. DHP174 is a soligenous peatland, sustained by precipitation and seasonal runoff from a south-flanking hillslope. There is ~15–20 m of relief across the ~100–150 m peatland surface, measured from the bottom of the hillslope to the northern edge of the peatland48. The slope ensures continuous drainage of surface and active layer runoff throughout the thaw season, and well-drained conditions in the uppermost 20–30 cm of peat. Saturated conditions during the summer occur below 20–30 cm depth and are maintained by elevated summer precipitation totals (Fig. 4d) and low saturated hydraulic conductivities (e.g. ~ 1 m day−1) that are typical below 30 cm depth in northern peat plateaus27. The DHP174 peatland sits on an alluvial terrace in a valley that drains into Engineer Creek (Fig. 1) and is underlain by permafrost, which limits runoff to the seasonally thawed active layer.
For the water isotope analysis presented in this study, we collected two 10 cm diameter core profiles using a portable permafrost drilling system49. The first core profile (DHP174-13, 3.83 m long) is from the ‘undisturbed’ upper surface, and the second core profile (DHP174-13L, 3.32 m long) was collected from an ‘excavated’ lower surface (Fig. 1c); the cores overlap stratigraphically, and the lower core includes the basal peats and underlying strata. The lower surface was excavated by the highway maintenance authority before 2007 (Google Earth time stamp) to thermally erode the underlying permafrost in order to expand the roadside quarry shown in the satellite images (Fig. 1d, e). At the time of coring (June 2013), frozen ground was encountered 16 and 24 cm below ground at the upper and lower coring sites, respectively, which are also the starting depths of the cores.
The frozen cores were split along their long axis with a bandsaw. The two core profiles were sub-sampled into 1 cm increments. For the DHP174-13 core, every tenth centimetre increment of pore ice was analysed for water isotopes starting from the core top (16 cm below the surface) to core bottom (382 cm composite depth); the only exception being DHP174-13 core depths 16–76 cm, which were analysed at near 1 cm resolution to allow for a more detailed analysis of pore ice in the active layer and across the active layer–permafrost interface. For the DHP174-13L core, every tenth centimetre increment of pore ice was analysed for water isotopes starting at a depth of 30 cm below the lower coring surface (230 cm composite depth; see Methods—chronology) down to the basal sediments at 530 cm composite depth). In total, 104 pore ice samples were analysed for water isotope ratios, 43 of which are active layer or transient layer pore ice, and 71 of which are relict pore ice.
The frozen samples were thawed at ~4 °C in plastic bags and the liquid was filtered into 2 mL autosampler vials. δD and δ18O ratios of the water were measured at the University of Alberta using a Picarro L2130-i analyser. Unknowns and standards were injected eight times, and mean raw values were calculated from the last four injections. Raw isotope ratios were corrected for inter-sample memory effects that are common to cavity ring-down spectroscopy systems50. We used certified water reference standards (USGS-45 and USGS-46) to normalise raw isotope ratios to the Vienna Standard Mean Ocean Water-Standard Light Antarctic Precipitation (VSMOW-SLAP) scale. Analytical error is 0.6‰ for δD and 0.1‰ for δ18O based on routine analysis of an internal deionised water standard.
A composite Bayesian age-depth model for the peat-sediment sequence was produced using radiocarbon (14C) dates and cryptotephra isochrons (see Davies51 for full details of the composite chronology). Depths are reported as core composite depths (ccd). The Bayesian age-depth model (see Supplementary Fig. 5) was constructed using the OxCal Poisson process model52.
Seventeen macrofossil samples were picked for 14C dating from the upper monolith and core profiles (Supplementary Fig. 1). This includes three dates from a previously collected core that was not used for pore ice analysis (DHP174-12, cored within 1.5 m of DHP174-13), but is included in the age model as it has precise modern cryptotephra and post-bomb 14C dates for its surface peat. Dates from all three cores were included in a composite profile where cryptotephra tie-points were used to correlate between the individual core profiles.
Where possible, Sphagnum remains were preferred as they can provide reliable 14C dates53. In the mineragenic, non-peat core sections at the base of the DHP174-13L profile, terrestrial plant macrofossils (twigs and leaves) were used for 14C dating. Samples were pre-treated at the University of Alberta following standard procedures54, and analysed for 14C content at the W. M. Keck Carbon Cycle AMS Laboratory (University of California, Irvine, CA, USA). Two secondary standards were pre-treated concurrently (a last interglacial non-finite-age wood, AVR-PAL-07; and a middle Holocene wood, FIRI-F) and the 14C results were within expected ranges. The unknown 14C dates were calibrated as appropriate using Bomb13NH155 and IntCal1356 calibration curves with Oxcal v4.3.257.
Glass shard concentration profiles for cryptotephra were produced at 1 cm resolution using standard methods58. Three samples geochemically correlated with known eruptions using major element EPMA glass analyses (Supplementary Fig. 2; see Davies51) were included in the composite age-depth model to improve model accuracy. Katmai-Novarupta (10 ccd) is a historic eruption with an observed occurrence date of 1912 CE. The Aniakchak CFE II (116 ccd) and Hayes set H unit F2 (133 ccd) eruptions are widely reported in Alaska and beyond, and have well-constrained Bayesian age estimates59.
Meteoric waters that reach the maximum thaw depth of ~58 cm (see maximum active layer thickness or ALTmax in Supplementary Note 1) will eventually freeze in situ as pore ice during the fall freezeback and become ‘relict’ as permafrost aggradation occurs. Since permafrost aggradation is paced by the accumulation of materials at the surface (peat or sediment), the age of relict pore ice at any given depth below the permafrost table can be referenced from the peat-sediment age-depth model by accounting for the 58 cm offset between relict pore ice and its contemporaneous surface. For example, ca. modern relict pore ice at 58 cm composite depth is age-equivalent to the modern surface (i.e. –63 cal yr BP in the age-depth model); likewise, relict pore ice at 530 cm composite depth is age-equivalent to a paleosurface at a composite depth of 472 cm (13,585 cal yr BP in our age-depth model). We assume a static ALTmax of 58 cm through time. Although ALTmax potentially did vary throughout the Holocene due to changes in summer climate, it likely remained within 3–6 cm of the modern ALTmax value and does not significantly influence our pore ice chronology or conclusions (Supplementary Note 1).
The δprecip isoscape of northern Canada is spatially complex due to regional differences in mean weather system trajectories, continentality, local orographic effects and varying degrees of continental moisture recycling22,60. Because each study area has its own potentially unique set of boundary conditions, δprecip signals preserved in ground ice are best interpreted in the context of local δprecip systematics22. We demonstrate the meteoric origins and seasonality of DHP174 pore ice by comparison with monthly δD and δ18O data61 from the nearest GNIP station at Mayo, Yukon, ~210 km south of DHP174 (Fig. 1a). Mayo and DHP174 are situated at similar elevations (600 and 720 m, respectively) in central Yukon in the western foothills of the Ogilvie Mountains. The Mayo data were used to constrain the LMWL (Fig. 3a), and the temperature-δprecip regression lines (Fig. 1b; Supplementary Fig. 1). Three Mayo records that plotted far below the LMWL (see Fig. 3a) were excluded from the regressions because of suspected evaporative enrichment. Extreme leverage points (>3σ Cook’s distance) were given a weight of zero in the calibration to avoid undue bias on the regression lines. By comparison with monthly GNIP data from Whitehorse (Fig. 1a), a site that is ~325 km south of Mayo and has a similar LMWL22, the Mayo temperature-δprecip regressions are found to be skilful in predicting δprecip at other continental sites in Yukon (Supplementary Fig. 1), which supports their application to more proximal sites such as DHP174.
We assume the δ18Oseawater reconstruction by Rohling et al.62 and that δDseawater changes were ~8× larger than δ18Oseawater, as supported by data from Late Pleistocene marine pore waters63 and as assumed by others14.
Two types of temperature estimates were made based on the pore ice record—mean temperatures for the modern era (Fig. 4c) and a full Holocene reconstruction of temperature anomalies (ΔT) (Fig. 5d). The first was calculated based on δDpore ice and δ18Opore ice of ca. modern relict pore ice (58–61 cm, uppermost 3 cm of permafrost from the DHP174-13 core), and was used as a line of evidence to validate the proposed summer seasonality of pore ice. Mean temperature estimates depend on the slope (m) and intercept (b) of the δprecip-temperature regression line, and δpore ice values corrected for the first ice enrichment effect (εice-water; 5.8‰ and 0.8‰ for δDpore ice and δ18Opore ice, respectively; Supplementary Note 2). Uncertainties in the mean temperature estimates take into account the 1σ error in δDpore ice and δ18Opore ice measurements, and 1σ uncertainties in the regression slope and intercept. The general formula used to calculate mean temperatures based on either isotope is: T = ([δpore ice ± 1σ] – εice-water – [b ± 1σ])/[m ± 1σ].
Temperature anomalies (ΔT) were reconstructed for the full Holocene based on δDpore ice anomalies (ΔδDpore ice, relative to Holocene mean) and the δDprecip–temperature slope (m). Since εice-water is assumed to be a constant (Supplementary Note 2), ΔδDpore ice and ΔδDprecip are equal. The formula used for estimating ΔT is: ΔT = [ΔδDpore ice ± 1σ]/[m ± 1σ].
The DHP174 isotope data and reconstruction are provided in the online supplement.
Journal peer review information: Nature Communications thanks Thomas Opel and the other anonymous reviewer for their contribution to the peer review of this work.
Serreze, M. C., Barrett, A. P., Stroeve, J. C., Kindig, D. N. & Holland, M. M. The emergence of surface-based Arctic amplification. Cryosphere 3, 11–19 (2009).
Lantz, T. C., Marsh, P. & Kokelj, S. V. Recent shrub proliferation in the Mackenzie Delta Uplands and microclimatic implications. Ecosystems 16, 47–59 (2013).
Kokelj, S. V., Lantz, T. C., Tunnicliffe, J., Segal, R. & Lacelle, D. Climate-driven thaw of permafrost preserved glacial landscapes, northwestern Canada. Geology 45, 371–374 (2017).
Millan, R., Mouginot, J. & Rignot, E. Mass budget of the glaciers and ice caps of the Queen Elizabeth Islands, Canada, from 1991 to 2015. Environ. Res. Lett. 12, 024016 (2017).
Miller, G. H., Lehman, S. J., Refsnider, K. A., Southon, J. R. & Zhong, Y. Unprecedented recent summer warmth in Arctic Canada. Geophys. Res. Lett. 40, 5745–5751 (2013).
McKay, N. P. & Kaufman, D. S. An extended Arctic proxy temperature database for the past 2,000 years. Sci. Data 1, 140026 (2014).
Briner, J. P. et al. Holocene climate change in Arctic Canada and Greenland. Quat. Sci. Rev. 147, 340–364 (2016).
Kaufman, D. S. et al. Holocene climate changes in eastern Beringia (NW North America)—a systematic review of multi-proxy evidence. Quat. Sci. Rev. 147, 312–339 (2016).
Lecavalier, B. S. et al. High Arctic Holocene temperature record from the Agassiz ice cap and Greenland ice sheet evolution. Proc. Natl Acad. Sci. USA 114, 5952–5957 (2017).
Marcott, S. A., Shakun, J. D., Clark, P. U. & Mix, A. C. A reconstruction of regional and global temperature for the past 11,300 years. Science (80-.). 339, 1198–1201 (2013).
Meyer, H. et al. Long-term winter warming trend in the Siberian Arctic during the mid- to late Holocene. Nat. Geosci. 8, 122–125 (2015).
Marsicek, J., Shuman, B. N., Bartlein, P. J., Shafer, S. L. & Brewer, S. Reconciling divergent trends and millennial variations in Holocene temperatures. Nature 554, 92–96 (2018).
Stenni, B. et al. The deuterium excess records of EPICA Dome C and Dronning Maud Land ice cores (East Antarctica). Quat. Sci. Rev. 29, 146–159 (2010).
Johnsen, S. J. et al. Oxygen isotope and palaeotemperature records from six Greenland ice-core stations: Camp Century, Dye-3, GRIP, GISP2, Renland and NorthGRIP. J. Quat. Sci. 16, 299–307 (2001).
Dansgaard, W., Johnsen, S. J., Møller, J. & Langway, C. C. One thousand centuries of climatic record from Camp Century on the Greenland ice sheet. Science (80-.). 166, 377–380 (1969).
Opel, T. et al. Ground-ice stable isotopes and cryostratigraphy reflect late Quaternary palaeoclimate in the Northeast Siberian Arctic (Oyogos Yar coast, Dmitry Laptev Strait). Clim. Past 13, 587–611 (2017).
Meyer, H. et al. Permafrost evidence for severe winter cooling during the Younger Dryas in northern Alaska. Geophys. Res. Lett. 37, 1–5 (2010).
Porter, T. J. et al. Multiple water isotope proxy reconstruction of extremely low last glacial temperatures in Eastern Beringia (Western Arctic). Quat. Sci. Rev. 137, 113–125 (2016).
Lachniet, M. S., Lawson, D. E. & Sloat, A. R. Revised 14C dating of ice wedge growth in interior Alaska (USA) to MIS 2 reveals cold paleoclimate and carbon recycling in ancient permafrost terrain. Quat. Res. 78, 217–225 (2012).
Schwamborn, G., Meyer, H., Fedorov, G., Schirrmeister, L. & Hubberten, H.-W. Ground ice and slope sediments archiving late Quaternary paleoenvironment and paleoclimate signals at the margins of El’gygytgyn Impact Crater, NE Siberia. Quat. Res. 66, 259–272 (2006).
Lacelle, D. On the δ18O, δD and d-excess relations in meteoric precipitation and during equilibrium freezing: theoretical approach and field examples. Permafr. Periglac. Process. 22, 13–25 (2011).
Gaglioti, B. V. et al. Younger-Dryas cooling and sea-ice feedbacks were prominent features of the Pleistocene–Holocene transition in Arctic Alaska. Quat. Sci. Rev. 169, 330–343 (2017).
French, H. & Shur, Y. The principles of cryostratigraphy. Earth Sci. Rev. 101, 190–206 (2010).
Outcalt, S. I., Nelson, F. E. & Hinkel, K. M. The zero-curtain effect: Heat and mass transfer across an isothermal region in freezing soil. Water Resour. Res. 26, 1509–1516 (1990).
Dansgaard, W. Stable isotopes in precipitation. Tellus 16, 436–468 (1964).
Quinton, W. L., Hayashi, M. & Carey, S. K. Peat hydraulic conductivity in cold regions and its relation to pore size and geometry. Hydrol. Process. 22, 2829–2837 (2008).
Pfahl, S. & Sodemann, H. What controls deuterium excess in global precipitation? Clim. Past 10, 771–781 (2014).
Max, L. et al. Sea surface temperature variability and sea-ice extent in the subarctic northwest Pacific during the past 15,000 years. Paleoceanography 27, https://doi.org/10.1029/2012PA002292 (2012).
Ruan, J. et al. Holocene variability in sea surface temperature and sea ice extent in the northern Bering Sea: a multiple biomarker study. Org. Geochem. 113, 1–9 (2017).
Ullman, D. J. et al. Final Laurentide ice-sheet deglaciation and Holocene climate-sea level change. Quat. Sci. Rev. 152, 49–59 (2016).
Baker, J. L., Lachniet, M. S., Chervyatsova, O., Asmerom, Y. & Polyak, V. J. Holocene warming in western continental Eurasia driven by glacial retreat and greenhouse forcing. Nat. Geosci. 10, 430–436 (2017).
Bartlein, P. J. et al. Early-Holocene warming in Beringia and its mediation by sea-level and vegetation changes. Clim. Past 11, 1197–1222 (2015).
Burn, C. R., Michel, F. A. & Smith, M. W. Stratigraphic, isotopic, and mineralogical evidence for an early Holocene thaw unconformity at Mayo, Yukon Territory. Can. J. Earth. Sci. 23, 794–803 (1986).
Ritchie, J. C. & Hare, F. K. Late Quaternary vegetation and climate near the Arctic treeline of northwestern North America. Quat. Res. 1, 331–342 (1971).
PAGES2k Consortium. A global multiproxy database for temperature reconstructions of the Common Era. Sci. Data 4, 170088 (2017).
Hansen, J., Ruedy, R., Sato, M. & Lo, K. Global surface temperature change. Rev. Geophys. 48, https://doi.org/10.1029/2010RG000345 (2010).
Jakobsson, M. et al. Post-glacial flooding of the Bering Land Bridge dated to 11 cal ka BP based on new geophysical and sediment records. Clim. Past 13, 991–1005 (2017).
Foley, J. A., Kutzbach, J. E., Coe, M. T. & Levis, S. Feedbacks between climate and the boreal forests during the Holocene epoch. Nature 371, 52–54 (1994).
IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge Univ Press, Cambridge, 2013).
Wendler, G., Moore, B. & Galloway, K. Strong temperature increase and shrinking sea ice in Arctic Alaska. Open Atmos. Sci. J. 8, 7–15 (2014).
Reyes, A. V., Froese, D. G. & Jensen, B. J. L. Permafrost response to last interglacial warming: field evidence from non-glaciated Yukon and Alaska. Quat. Sci. Rev. 29, 3256–3274 (2010).
Burn, C. R. Cryostratigraphy, paleogeography, and climate change during the early Holocene warm interval, western Arctic coast, Canada. Can. J. Earth. Sci. 34, 912–925 (1997).
Murton, J. B. & French, H. Cryostructures in permafrost. Tuktoyaktuk coastlands, western Arctic Canada. Can. J. Earth Sci. 31, 737–747 (1994).
Harry, D. G., French, H. & Pollard, W. H. Massive ground ice ad ice-cored terrain near Sabine Point, Yukon Coastal Plain. Can. J. Earth Sci. 25, 1846–1856 (1988).
Segal, R., Lantz, T. & Kokelj, S. V. Acceleration of thaw slump activity in glaciated landscapes of the Western Canadian Arctic. Environ. Res. Lett. 11, 034025 (2016).
Fraser, R. et al. Climate sensitivity of high arctic permafrost terrain demonstrated by widespread ice-wedge thermokarst on Banks Island. Remote Sens. 10, 954 (2018).
Porter, C. et al. ArcticDEM. Harvard Database, V1. https://doi.org/10.7910/DVN/OHHUKH. Accessed 1 Sept 2018.
Calmels, F., Gagnon, O. & Allard, M. A portable Earth-drill system for permafrost studies. Permafr. Periglac. Process. 16, 311–315 (2005).
Schauer, A. J., Schoenemann, S. W. & Steig, E. J. Routine high-precision analysis of triple water-isotope ratios using cavity ring-down spectroscopy. Rapid Commun. Mass Spectrom. 30, 2059–2069 (2016).
Davies, L. J. The Development of a Holocene Cryptotephra Framework in Northwestern North America (University of Alberta, Edmonton, 2018).
Ramsey, C. B. Deposition models for chronological records. Quat. Sci. Rev. 27, 42–60 (2008).
Nilsson, M., Klarqvist, M., Bohlin, E. & Possnert, G. Variation in 14C age of macrofossils and different fractions of minute peat samples dated by AMS. The Holocene 11, 579–586 (2001).
Reyes, A. V. et al. A late–Middle Pleistocene (Marine Isotope Stage 6) vegetated surface buried by Old Crow tephra at the Palisades, interior Alaska. Quat. Sci. Rev. 29, 801–811 (2010).
Hua, Q. & Barbetti, M. Atmospheric radiocarbon for the period 1950–2010. Radiocarbon 55, 2059–2072 (2013).
Reimer, P. et al. IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP. Radiocarbon 55, 1869–1887 (2013).
Blockley, S. P. E. et al. A new and less destructive laboratory procedure for the physical separation of distal glass tephra shards from sediments. Quat. Sci. Rev. 24, 1952–1960 (2005).
Davies, L. J., Jensen, B. J. L., Froese, D. G. & Wallace, K. L. Late Pleistocene and Holocene tephrostratigraphy of interior Alaska and Yukon: Key beds and chronologies over the past 30,000 years. Quat. Sci. Rev. 146, 28–53 (2016).
Aemisegger, F. et al. Deuterium excess as a proxy for continental moisture recycling and plant transpiration. Atmos. Chem. Phys. 14, 4029–4054 (2014).
IAEA/WMO. Global Network of Isotopes in Precipitation. The GNIP Database https://nucleus.iaea.org/wiser/. (2017).
Rohling, E. J. et al. Sea-level and deep-sea-temperature variability over the past 5.3 million years. Nature 508, 477–482 (2014).
Schrag, D. P. et al. The oxygen isotopic composition of seawater during the Last Glacial Maximum. Quat. Sci. Rev. 21, 331–342 (2002).
Laskar, J., Fienga, A., Gastineau, M. & Manche, H. La2010: a new orbital solution for the long-term motion of the Earth. Astron. Astrophys. 532, A89 (2011).
Keeling, C. D. et al. in A History of Atmospheric CO 2 and its effects on Plants, Animals, and Ecosystems (eds Ehleringer, J. R., Cerling, T. E. & Dearing, M. D.) 83–113 (Springer, Berlin, 2005).
MacFarling Meure, C. et al. Law Dome CO2, CH4 and N2O ice core records extended to 2000 years BP. Geophys. Res. Lett. 33, L14810 (2006).
This work was supported by an I.W. Killam Memorial Postdoctoral Fellowship to T.J.P., and Natural Sciences and Engineering Council of Canada Discovery Grants to T.J.P. and D.F. We thank Fabrice Calmels, Sydney Clackett, Gerard Otiniano and Aleesha Bakkelund for assisting with the collection of active layer measurements.
T.J.P., D.G.F. and L.J.D. co-conceived the study. T.J.P., D.G.F., L.J.D. and S.B. participated in the fieldwork. L.J.D. sub-sampled the core profiles and developed the age-depth model. L.J.D., S.B., S.W.S. and E.J.S. performed water isotope analysis. T.J.P. and S.W.S. analysed the water isotope data and developed the climate interpretation and the reconstruction. T.J.P. wrote the majority of the manuscript, and D.G.F., L.J.D., S.B., S.W.S. and E.J.S. co-wrote the manuscript.
Correspondence to Trevor J. Porter.

References: in fine
 V. 
 V. 
 V. 
 V. 
 V.