Past climate and continentality inferred from ice wedges at Batagay megaslump in the Northern Hemisphere's most continental region, Yana Highlands, interior Yakutia

Ice wedges in the Yana Highlands of interior Yakutia – the most continental region of the Northern Hemisphere – were investigated to elucidate changes in winter climate and continentality that have taken place since the Middle Pleistocene. The Batagay megaslump exposes ice wedges and composite wedges that were sampled from three cryostratigraphic units: the lower ice complex of likely pre-Marine Isotope Stage (MIS) 6 age, the upper ice complex (Yedoma) and the upper sand unit (both MIS 3 to 2). A terrace of the nearby Adycha River provides a Late Holocene (MIS 1) ice wedge that serves as a modern reference for interpretation. The stable-isotope composition of ice wedges in the MIS 3 upper ice complex at Batagay is more depleted (mean δ¹⁸O about −35 ‰) than those from 17 other ice-wedge study sites across coastal and central Yakutia. This observation points to lower winter temperatures and therefore higher continentality in the Yana Highlands during MIS 3. Likewise, more depleted isotope values are found in Holocene wedge ice (mean δ¹⁸O about −29 ‰) compared to other sites in Yakutia. Ice-wedge isotopic signatures of the lower ice complex (mean δ¹⁸O about −33 ‰) and of the MIS 3–2 upper sand unit (mean δ¹⁸O from about −33 ‰ to −30 ‰) are less distinctive regionally. The latter unit preserves traces of fast formation in rapidly accumulating sand sheets and of post-depositional isotopic fractionation.

1 Introduction

Figure 1 (a) Location map of study region in northeast Siberia showing all ice-wedge study sites mentioned. (b) Study sites in the Yana Highlands (background image: Sentinel-2 mosaic; 7 and 21 August 2017).

Interior Yakutia is currently the most continental region of the Northern Hemisphere. In Verkhoyansk, the cold pole of the Northern Hemisphere, mean monthly air temperature ranges exceed 60 ^∘C and absolute temperature ranges can reach 100 ^∘C (Lydolph, 1985). Mean annual precipitation is only 155 mm (Lydolph, 1985). Scarce stable-isotope data indicate that modern winter precipitation in this region has the most depleted isotopic composition in Siberia (Kurita et al., 2004, 2005) according to the “continental effect” in stable-isotope fractionation, i.e. the farthest distance from the source region (Dansgaard, 1964). However, warm summers enrich stable isotopes in summer precipitation at interior sites compared to coastal sites in the north, which also reflects the strong continentality of interior Yakutia. During Late Pleistocene cold stages, when lower sea level exposed larger areas of continental shelves, interior Yakutia probably experienced even stronger continentality, with further decreased annual precipitation and increased annual temperature amplitudes. Such increased amplitudes must have resulted from either lower winter temperatures, higher summer temperatures, or a combination of both.

While biological proxies such as pollen or plant macrofossil remains can be used to reconstruct summer climate, winter conditions in interior Yakutia may be inferred only from ice wedges, which can be found in Middle and Late Pleistocene as well as Holocene permafrost deposits. Ice wedges provide winter temperature information due to their specific seasonality, with frost cracking in winter and crack infilling by snowmelt in spring (Opel et al., 2018). They integrate the stable-isotope composition of cold-season precipitation (Meyer et al., 2015), which in high northern latitudes is mainly a function of temperature (Dansgaard, 1964). However, a knowledge gap exists about the winter paleoclimate of interior Yakutia, because ice-wedge studies in the Yana Highlands – between the Verkhoyansky and Chersky ranges – began only recently (Vasil'chuk et al., 2017). In contrast, other sites in east Siberia (west Beringia) have been extensively studied over recent decades, including the coastal lowlands to the north (e.g. Meyer et al., 2002a, b; Opel et al., 2017b), the Kolyma region to the east (e.g. Vasil'chuk et al., 2001; Vasil'chuk, 2013; Vasil'chuk and Vasil'chuk, 2014) as well as central Yakutia to the south (e.g. Popp et al., 2006) (Fig. 1). Ice-wedge studies have been carried out also in east Beringia, e.g. in Alaska, Yukon, and the Northwest Territories (Kotler and Burn, 2000; Meyer et al., 2010b; Fritz et al., 2012; Lachniet et al., 2012; Porter et al., 2016).

Here, we report cryostratigraphic field observations of ice wedges and composite wedges made during a reconnaissance expedition from 27 July to 5 August 2017 together with stable-isotope data from the Batagay megaslump and Adycha River floodplain (Fig. 1). New radiocarbon ages from organic remains in ice wedges and host sediments help constrain the late Quaternary chronology. The observations and ages elucidate the history of wedge development, past winter temperatures, and continentality starting from the Middle Pleistocene.

2 Study site characteristics, stratigraphy, and chronology

Neither study site has been glaciated during at least the last 50 000 years, and the Batagay megaslump provides access to permafrost formations starting from the Middle Pleistocene (Ashastina et al., 2017; Murton et al., 2017). The present mean annual air temperature at Batagay is −15.4 ^∘C (mean of the coldest month: −47.7 ^∘C; mean of the warmest month: 15.4 ^∘C) and the mean annual precipitation is 170 to 200 mm (Ivanova, 2003). Prevailing wind directions are southwest in winter (January) and north or northeast in summer (July) (Murton et al., 2017). The mean annual temperature of the permafrost is −8.0 to −5.5 ^∘C and the active-layer depth between 0.2 and 0.4 m beneath forest and moss covers and 0.4–1.2 m beneath open sites (Ivanova, 2003).

The Batagay megaslump (67.58^∘ N, 134.77^∘ E), close to the city of Batagay on the east bank of the Yana River, is the world's largest known thaw slump (Fig. 2) (Ashastina et al., 2017; Murton et al., 2017). It is located on a hillslope in the taiga landscape (328 m above sea level for the highest part of the headwall) and exposes a sequence of permafrost deposits at least 60 m thick, with underlying bedrock cropping out in places (Kunitsky et al., 2013). The slump has formed during recent decades (Savvinov et al., 2018), with increasing headwall retreat rates of up to 30 m per year. In 2016, it reached a width of 840 m and a total area of >70 ha (Günther et al., 2016).

Figure 2 Overview photograph of the Batagay megaslump showing the study sites of Ashastina et al. (2017) and Murton et al. (2017) and this study. Photograph taken by Alexander Gabyshev in 2015.

The slump exposes Pleistocene and Holocene permafrost formations ranging in age from MIS 6 (or older) to MIS 1 (Ashastina et al., 2017, 2018; Murton et al., 2017). Above slate bedrock and a basal diamicton, four major cryostratigraphic units contain ice wedges and/or composite wedges (Fig. 3, Table 1).

Table 1 Cryostratigraphic units of the Batagay megaslump as presented by Ashastina et al. (2017) and Murton et al. (2017) and the reconciled stratigraphic framework used in this study and proposed for future studies.

Figure 3 (a) Panoramic photograph of the Batagay megaslump with main stratigraphic units and approximate sampling locations of ice wedges (blue arrows). Red circles and numbers indicate the approximate sampling locations of ice wedges used by Vasil'chuk et al. (2017). (b) Schematic stratigraphic section of the exposed units along the headwall of the Batagay megaslump with approximate positions of the studied sections.

The lowest unit is a pebbly dark sand, 3–7 m thick, with ice wedges at least 2–3 m high and 1 m wide, truncated by a thaw unconformity. This lower ice complex has neither been dated nor sampled previously for wedge ice.

The lower sand unit above reaches about 20 m in thickness and is composed of yellowish pore-ice cemented (fine) sand with grey horizontal bands. It contains tall, narrow syngenetic ice wedges up to 0.5 m wide. The middle part of the unit was dated to 142.8±25.3 and >123.2 kyr by optically stimulated luminescence (OSL) and to 210.0±23.0 kyr by infrared stimulated luminescence (IRSL) (Ashastina et al., 2017), pointing to its deposition during MIS 6. This age is supported by palaeo-ecological data (Ashastina et al., 2018). An organic-rich unit of silty sand with abundant wood and plant remains (in lenses up to 3 m thick) overlies a distinct erosional surface near the top of the lower sand unit. A wood fragment yielded a conventional radiocarbon age of 49.32±3.15 ¹⁴C kyr BP, outside of the range of calibration (Murton et al., 2017), though palaeo-ecological analysis points to a last interglacial age for the organic material (Ashastina et al., 2018).

The overlying upper ice complex is 20–25 m thick and dominated by huge syngenetic ice wedges up to few metres wide and at least several metres high within silty and sandy deposits. Finite radiocarbon ages from 49.0±2.0 to 12.66±0.05 ¹⁴C kyr BP (Ashastina et al., 2017) indicate sedimentation gaps or erosive events and reveal a MIS 3 to MIS 2 age for the upper part of the upper ice complex. This unit has numerous chronostratigraphic analogies in the coastal lowlands of northern Yakutia, i.e. the Yedoma Ice Complex (Schirrmeister et al., 2011b; Murton et al., 2015). The upper ice complex forms the highest Pleistocene stratigraphic unit in the upper central part (upslope) of the slump headwall, whereas downslope towards the slump mouth it is overlain and partially grades into the upper sand unit (Fig. 3).

The upper sand is up to about 20 m thick and consists of pore-ice cemented brown to grey sand with narrow (≤0.5 m wide) syngenetic ice wedges and composite wedges. Radiocarbon ages between 36.30±0.70 and 26.18±0.22 ¹⁴C kyr BP (Ashastina et al., 2017; Murton et al., 2017) indicate deposition during MIS 3 and MIS 2, similar to the upper ice complex, but in a more complex spatiotemporal pattern, likely due to varying palaeotopography, sediment supply, and accumulation rates.

A near-surface layer 1–1.5 m thick of brown sand and modern soil covers the sequence. This layer was dated to 295±30 ¹⁴C yr BP (Ashastina et al., 2017). Distinct (thermo-)erosional contacts can be found between the lower ice complex and the lower sand, the lower sand and the upper ice complex, and below the near-surface layer.

According to Murton et al. (2017), exposed floodplains of proximal rivers such as the Batagay and Yana, within 2 and 10 km, respectively, of the slump are the assumed major source of the sediments exposed in the headwall, which implies upslope directed transport by wind. Periglacial and nival processes on nearby hillslopes may also have contributed to sediment supply (Ashastina et al., 2017).

For a Late Holocene reference, we studied an ice wedge at the actively eroding Holocene to recent bank of a small cut-off channel of the Adycha River (67.66^∘ N, 135.69^∘ E; 138 m above sea level), about 40 km east of the Batagay megaslump (Fig. 1). The Adycha site is located about 3 km upstream of the 60–65 m high and 1.2 km long river-bluff stratigraphic section at Ulakhan Sullar (Kaplina et al., 1983; Sher et al., 2011; Germonpre et al., 2017). Organic-rich silty sands of alluvial origin are exposed in the Adycha riverbank section.

3 Material and methods

3.1 Fieldwork

Cryostratigraphic observations of ice wedges, composite wedges, unit contacts, and sediments at the Batagay megaslump and a terrace of the Adycha River provided a framework for sample selection. A total of six ice wedges and composite wedges from the lower ice complex, the upper ice complex, and the upper sand were described in detail and sampled at two sections of the slump (Figs. 2, 3). Section 1 was studied in a gully in the icy badlands of the slump (Fig. 4), and section 2 was studied in a headwall slope segment (Fig. 5). Composite wedges of the lower sand unit were inaccessible due to dangerous outcrop conditions, and no wedges were observed in the near-surface layer. At the Adycha River, we sampled a single ice wedge.

Figure 4 Schematic cryostratigraphic sketch of section 1, Batagay megaslump.

Figure 5 Schematic cryostratigraphic sketch of section 2, Batagay megaslump.

Ice samples about 4 cm wide, mostly in horizontal profiles, were obtained by chainsaw. Details about the studied ice and composite wedges are given in Table S1 in the Supplement. The ice samples were melted on site in freshly opened standard Whirl-Pak sample bags, and the meltwater filled up 30 mL PE bottles that were then tightly closed and stored cool until stable-isotope analysis.

In addition, water from small streams draining the slump was sampled in several parts of the slump floor and the main outflow stream (n=7) as well as from summer precipitation (n=4). Supernatant water from selected Batagay megaslump host sediment samples (n=9) was taken for stable-isotope analysis of intrasedimental pore and segregated ice.

3.2 Stable-isotope analysis

The stable oxygen (δ¹⁸O) and hydrogen (δD) isotope ratios of all samples were determined at the Stable Isotope Laboratory of the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research in Potsdam, using a Finnigan MAT Delta S mass spectrometer with an analytical precision of better than ±0.1 ‰ for δ¹⁸O and ±0.8 ‰ for δD (Meyer et al., 2000). The isotopic composition is expressed in δ per mille values (‰) relative to the Vienna Standard Mean Ocean Water (VSMOW) standard. The deuterium excess d was calculated as $d=\mathit{\delta }\mathrm{D}-\mathrm{8}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (Dansgaard, 1964). For three of seven studied wedges, samples were excluded from further analysis (five in total; Table S2) on the basis that they were lateral samples with high sediment content and/or had isotopic values distinctly different from the wedge centres and similar to host sediments, probably due to exchange processes between wedge ice and host sediments (Meyer et al., 2002a, 2010a).

3.3 Radiocarbon dating

We picked organic material from our ice-wedge and host sediment samples for radiocarbon dating at the MICADAS ¹⁴C dating facility of the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research in Bremerhaven. The organic material was subjected to chemical leaching and cleaning following the ABA (acid–base–acid) procedure. Material was first submerged in 1 M HCl (for 30 min at 60 ^∘C) to remove carbonate contamination. Following rinsing with Milli-Q water, samples were submerged in 1 M NaOH (for 30 min at 60 ^∘C) to leach out humic acids. The base leaching was repeated for a minimum of seven times or until no colouring of the solution was evident. Afterwards, the material was subjected to a final acid treatment (1 M HCl; 30 min at 60 ^∘C), rinsed to neutral with Milli-Q, and dried for a minimum of 12 h at 60 ^∘C. Milli-Q used for cleaning and acid/base solutions was pre-cleaned by liquid–liquid extraction with dichloromethane to remove remaining organic contamination.

Samples were packed in tin capsules and combusted individually using an Elementar vario ISOTOPE EA (elemental analyser). If samples were smaller than approximately 200 µgC, the CO₂ produced by combustion with the EA was directly injected into the hybrid ion source of the Ionplus MICADAS using the gas interface system (GIS) (Fahrni et al., 2013). For larger samples, the CO₂ produced was graphitized using the Ionplus AGE3 system (automated graphitization system; Wacker et al., 2010c). The radiocarbon content of the samples was determined alongside size-matched reference standards (oxalic acid; NIST 4990c) and blanks (phthalic anhydride; Sigma-Aldrich 320065) using the Ionplus MICADAS dating system (Synal et al., 2007; Wacker et al., 2010a), and blank correction and standard normalization were performed using the BATS software (Wacker et al., 2010b). Results are reported as F¹⁴C (Reimer et al., 2004) and conventional radiocarbon (¹⁴C) ages in years BP. Several analyses yielded lower ¹⁴C intensities than the blanks included in the respective sequence. In these cases, results are reported as F¹⁴C_sample<F¹⁴C_blanks. Conventional ¹⁴C ages were calibrated using Oxcal 4.3 (Bronk Ramsey, 2009) based on the IntCAL13 dataset (Reimer et al., 2013).

4 Results

4.1 Batagay megaslump

4.1.1 Field observations of ice wedges

The present study confirms cryostratigraphic observations about ice wedges and composite wedges from earlier studies (Ashastina et al., 2017; Murton et al., 2017) and provides new observations about ice-wedge relationships across contacts between the cryostratigraphic units, as outlined below.

In the lower ice complex, ice wedges are truncated by a thaw unconformity, though the toes of some narrow syngenetic composite wedges in the overlying lower sand unit extend down across the contact (Fig. S1 in the Supplement). A single example of clear ice sharply overlying the shoulder of an ice wedge in the lower ice complex (B17-IW1, Fig. 4) contained brown plant remains and round, few-millimetre diameter organic bodies identified as hare droppings. This ice lacked the foliation characteristic of wedge ice and is interpreted as pool ice. The basal contact of the lower ice complex was not observed.

In the lower and upper sand units, narrow syngenetic wedges vary in apparent width from a few centimetres to about 0.5 m. Width commonly varies irregularly with height along individual wedges, sometimes gradually, sometimes abruptly. Wedge height varies from a few metres to at least 12 m. Wedges tend to be oriented approximately at right angle to colour bands in the sand, such that wedges in the upper sand unit – whose bands dip downslope, parallel to the ground surface – are characteristically subvertical, with their tops inclined downslope (Fig. S2). The wedge infills grade between end members of ice-wedge ice and icy sand wedges, with slightly sandy ice-wedge ice (i.e. composite wedges) as the most common type.

Near the base of the upper ice complex, syngenetic ice wedges, a few metres wide, tend to narrow downwards in the lower several metres of the unit and terminate with irregular or flattish or U- to V-shaped toes over a vertical distance of about 1–3 m (Fig. S2). At the northeast part of the slump headwall, however, the base of the upper ice complex is a sharp contact that cuts down into and truncates bands in the underlying lower sand unit over a vertical distance of up to at least several metres (Fig. 3). As a result, the upper ice complex thickens substantially downslope.

The top of the upper ice complex displays a variety of wedge relationships with the overlying upper sand unit (Fig. S3). Some wide syngenetic wedges terminate along planar to gently undulating contacts that are horizontal to gently dipping and have apparent widths of about 1–3 m; some taper upward into narrow wedges characteristic of the upper sand unit; some taper irregularly upward, marked by shoulders up to about 0.5 m wide; and some end upwards with a narrow offshoot. These changes occur over a vertical distance of about 1–3 m.

4.1.2 Lower ice complex (pre-MIS 6)

A small truncated ice wedge (B17-IW1; 0.5 m wide, 1.1 m high) from the lower ice complex was sampled in section 1 (Fig. 4) near the slump bottom, about 50 m below the ground surface and about 5 m below the altitude level in the lower sand unit dated by luminescence to 142.8±25.3 and >123.2 kyr (OSL) and to 210.0±23.0 kyr (IRSL) (Ashastina et al., 2017). It shows a mean stable-isotope composition of −33.1 ‰ for δ¹⁸O, −256.7 ‰ for δD, and 8.0 ‰ for d (n=6), with very little variability (Table 2, Fig. 6). The regression line in a δ¹⁸O–δD bi-plot is $\mathit{\delta }\mathrm{D}=\mathrm{12.78}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{166.23}$ (R²=0.80).

Figure 6 Co-isotopic δ¹⁸O–δD diagram of ice wedges and intrasedimental ice sampled at the Batagay megaslump and the Adycha River. The data from the Batagay megaslump outflow samples record an integrated signal of the entire Batagay permafrost sequence and modern precipitation.

Table 2 Stable-isotope (δ¹⁸O, δD, and d) minimum, mean, and maximum values, and standard deviations (SDs) as well as slopes and intercepts in the δ¹⁸O–δD diagram for all ice and composite wedges, intrasedimental ice, outflow water, and rainwater, respectively.

n/a – not applicable

The isotopic composition of the intrasedimental ice (n=2) shows similar mean δ values (−32.5 ‰ for δ¹⁸O, −247.1 ‰ for δD) but higher d values (12.8 ‰). The co-isotopic regression line is $\mathit{\delta }\mathrm{D}=\mathrm{4.09}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{114.27}$.

Radiocarbon dating of twigs, Cyperaceae stems, and roots from a sediment sample 0.5 m above the truncated ice wedge yielded nonfinite ages of > 53 400 and > 37 500 ¹⁴C yr BP (Table 3). A hare dropping from the ice wedge was dated to 15 792 ± 358 ¹⁴C yr BP, which can be only explained by contamination with relocated material.

Table 3 Radiocarbon ages of organic remains in ice wedges (sample ID includes IW) and host sediments of the Batagay megaslump and at the Adycha River. The samples were radiocarbon dated as gas targets (lab ID ends with 1.1) and graphite targets (lab ID ends with 2.1).

4.1.3 Lower sand unit (MIS 6)

The isotopic composition of intrasedimental ice (n=2) from the lower sand unit shows very little variability (mean −32.3 ‰ for δ¹⁸O, −241.6 ‰ for δD, 16.6 ‰ for d (Table 2, Fig. 6). The co-isotopic regression line is $\mathit{\delta }\mathrm{D}=\mathrm{1.33}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{198.54}$.

4.1.4 Upper ice complex (MIS 3–2)

Two large syngenetic ice wedges from the overlying ice-rich upper ice complex in Sect. 2 (Fig. 5) in the southern part of the slump show more depleted isotope values as compared to the lower ice complex (Table 2, Fig. 6). Ice wedge B17-IW6 (about 0.5 m wide, sampled about 26 m below surface, b.s.) shows a mean isotopic composition of −35.1 ‰ for δ¹⁸O, −269.4 ‰ for δD, and 11.0 ‰ for d (n=4). The co-isotopic regression is $\mathit{\delta }\mathrm{D}=\mathrm{12.25}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{159.79}$ (R²=0.99). The mean stable-isotope values of ice wedge B17-IW5 (about 1.6 m wide and sampled about 20 m b.s.) are −34.9 ‰ for δ¹⁸O, −271.0 ‰ for δD, and 8.1 ‰ for d (n=12). The regression line in a δ¹⁸O–δD bi-plot is $\mathit{\delta }\mathrm{D}=\mathrm{7.41}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{12.36}$ (R²=0.99). The intrasedimental ice data (n=2) show distinctly enriched mean values (−26.2 ‰ for δ¹⁸O, −195.0 ‰ for δD, and 14.5 ‰ for d).

Unidentified plant (bract fragments and roots) and insect (complete pieces and fragments of elytron) remains from ice wedge B17-IW6 as well as Cyperaceae remains and roots from host sediment at the lowest studied level in section 2, about 26 m b.s., both revealed nonfinite ages (> 37 500 ¹⁴C yr BP; Table 3). At the level of about 20 m b.s., rootlets from the host sediment yielded ages of > 37 500 and 47 550 ± 677 ¹⁴C yr BP, whereas unidentified plant remains (twigs, roots, and florets) from within ice wedge B17-IW5 were dated to 24 858 ± 536 ¹⁴C yr BP.

4.1.5 Upper sand unit (MIS 3–2)

Two narrow sand–ice wedges (composite wedges) and one ice wedge were sampled from the rather ice-poor upper sand unit in section 2 (Fig. 5) at the southern part of the slump. Composite wedge B17-IW2 (0.25 m wide) was sampled 2.6 and 2.8 m b.s., and composite wedge B17-IW3 (0.2 m wide) was sampled 2.0 m b.s. The stable-isotope composition of the composite wedges from the upper sand unit differs from those of the lower ice complex and the upper ice complex (Table 2, Fig. 6). The mean values are −32.8 ‰ and −31.9 ‰ for δ¹⁸O, −247.5 ‰ and −240.5 ‰ for δD, and 15.1 ‰ and 14.7 ‰ for d (n=4 and n=2), respectively, and plot well above the global meteoric water line (GMWL, $\mathit{\delta }\mathrm{D}=\mathrm{8}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{10}$) (Craig, 1961). The corresponding regressions are $\mathit{\delta }\mathrm{D}=\mathrm{12.51}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{163.29}$ (R²=0.92) and $\mathit{\delta }\mathrm{D}=\mathrm{14.47}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{220.91}$, respectively.

Ice wedge B17-IW4 (about 0.4 m wide and containing a few thin sand veins) was sampled in the upper sand 2.0 m b.s. and exhibits a more enriched stable-isotope composition (mean values −29.9 ‰ for δ¹⁸O, −227.8 ‰ for δD, and 11.6 ‰ for d (n=6). The co-isotopic regression is $\mathit{\delta }\mathrm{D}=\mathrm{15.03}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{222.01}$ (R²=0.91) (Table 2, Fig. 6).

The intrasedimental ice (n=3) shows similar mean δ values (−30.6 ‰ for δ¹⁸O, −220.0 ‰ for δD) but higher d (24.4 ‰) values compared to the composite wedges and the ice wedges of the upper sand. The regression line in a δ¹⁸O–δD bi-plot is $\mathit{\delta }\mathrm{D}=\mathrm{10.61}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{104.19}$ (R²=0.99).

Two radiocarbon ages (analysed as gas) of small samples of Empetrum nigrum (crowberry) leaf, twigs, and Cyperaceae remains from host sediment about 11 m b.s. are nonfinite, whereas a graphite target yielded an age of 38 348 ± 236 ¹⁴C yr BP.

4.1.6 Slump outflow

The samples from the outflow stream draining the slump show average stable-isotope values of −31.7 ‰ for δ¹⁸O, −238.6 ‰ for δD, and 15.2 ‰ for d. They plot on a co-isotopic regression line of $\mathit{\delta }\mathrm{D}=\mathrm{8.84}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{41.83}$ (R²=0.97, n=7), well above the GMWL (Table 2, Fig. 6).

4.1.7 Summer precipitation

The samples of summer precipitation at the Batagay megaslump show much more enriched stable-isotope compositions than the wedges, as expected (Table 2). Mean values are −14.3 ‰ for δ¹⁸O, −112.0 ‰ for δD, and 2.6 ‰ for d (n=4). The samples plot on a co-isotopic regression line of $\mathit{\delta }\mathrm{D}=\mathrm{7.93}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{1.65}$ (R²=1).

4.2 Adycha River

Ice wedge A17-IW3 was about 0.3 m wide and sampled about 1.4 m above river level and about 2 m below the ground surface. Its mean stable-isotope values are −29.0 ‰ for δ¹⁸O, −226.3 ‰ for δD, and 5.5 ‰ for d (n=10, Table 2, Fig. 6). The regression line in a δ¹⁸O–δD bi-plot is $\mathit{\delta }\mathrm{D}=\mathrm{7.01}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{23.25}$ (R²=0.99).

Rootlets and larch needles from the host sediment were dated to 213±109 ¹⁴C yr BP and confirm the Late Holocene age of the alluvial sediments and the ice wedge.

5 Discussion

5.1 Chronostratigraphy of the Batagay megaslump

To date, the exceptionally long record of permafrost history under highly continental climate conditions as exposed in the Batagay megaslump is still not entirely established in terms of its geochronology. Previous studies by Ashastina et al. (2017) and Murton et al. (2017) provide the general differentiation into four main units with ice and/or composite wedges overlain by Late Holocene cover deposits, which are the lower ice complex (pre-MIS 6), the lower sand unit (MIS 6), the upper ice complex (MIS 3–2), and the upper sand unit (MIS 3–2). Thaw unconformities as observed in sampled sections of this study (Figs. 4, 5), but also in previous studies (Ashastina et al., 2017; Murton et al., 2017), question the chronological continuity of the archive. Such (thermo-)erosional features are often observed in late Quaternary permafrost chronologies (e.g. Wetterich et al., 2009) and are most likely caused by intense permafrost degradation during warm stages such as the last interglacial with widespread thermokarst (Reyes et al., 2010; e.g. Kienast et al., 2011) or exceptionally warm interstadials but might also relate to palaeotopography and respective accumulation and erosion areas.

The Batagay megaslump, on account of its huge dimensions and rapid erosion, holds potential to identify additional records and units in future studies. With the present study and additional radiocarbon dating, further confidence of the existing stratigraphy is reached, although the wedge-ice-derived age information requires careful interpretation. Due to the formation mechanism and downward directed growth of wedge ice, organic remains from inside the ice are commonly younger than those of the host deposits at the same altitude (i.e. sampling depth). The risk of dating redeposited organic material, that entered into or later froze onto the surface of ice recently exposed in a slump, is common, in particular in erosional features such as slump-floor gullies. It can be seen in the truncated lower ice complex (Fig. 4) where an age of about 16 ¹⁴C kyr BP was obtained from the ice wedge, while radiocarbon dates of the host deposit are infinite and luminescence ages of the overlying lower sand unit reach 142.8±25.3 and >123.2 kyr by OSL as well as 210.0±23.0 kyr by IRSL (Ashastina et al., 2017). We therefore excluded this anomalously young radiocarbon age from further interpretation.

A distinct unit of wood and plant remains up to 3 m thick and traceable along the exposure is situated near the top of the lower sand unit and below the upper ice complex. Radiocarbon dating of this layer revealed an infinite radiocarbon age of >44 ¹⁴C kyr BP (Ashastina et al., 2017) as well as an out-of-calibration age of 49.32±3.15 ¹⁴C kyr BP (Murton et al., 2017). The palaeo-botanic proxy data indicate a warm climate. The plant macro-fossil record is characterized by considerable species richness (34 species), mainly representing northern taiga forests, including larch and birch trees, shrub alder and raspberry, and further representatives of the taiga's understory. Hence, the unit is aligned to the last interglacial (Ashastina et al., 2018), which caused particular thaw of the underlying lower sand unit and formed a disconformity.

The expected difference between ice-hosted and sediment-hosted ages of the same sample depth is clearly seen in the upper ice complex, where organic remains from the ice wedge are dated to about 25 ¹⁴C kyr BP and the host deposits yield nonfinite to about 48 ¹⁴C kyr BP ages (Fig. 5). However, even though the age offset is very large and a contamination of the dated sample in the field cannot be fully excluded, such information may be reliable and lies within the previously dated overall accumulation period of the upper ice complex between about 49 and 13 ¹⁴C kyr BP (Ashastina et al., 2017; Murton et al., 2017).

The new radiocarbon ages from the upper sand unit largely confirm the assumed accumulation period of this unit from about 36 to 26 ¹⁴C kyr BP (Ashastina et al., 2017; Murton et al., 2017) but predate it by about 2000 years, at least with a date of about 38 ¹⁴C kyr BP and nonfinite ages of >37.5 ¹⁴C kyr BP from the host deposits. We note that the upper sand unit has been studied at two different sites of the slump in different altitude levels (north–northeast of the central headwall (Ashastina et al., 2017; Murton et al., 2017) and together with the upper ice complex more downslope in the southern part of the slump (Fig. 2), which complicates the correlation of dating results from different sampling sites. Furthermore, close to our section 2, the upper sand could not clearly be differentiated from the upper ice complex during fieldwork in 2014 (Ashastina et al., 2017). Hence, the more or less simultaneous accumulation of both the upper ice complex and the upper sand unit during MIS 3–2 but their differing spatial distribution in the slump highlights the importance of palaeo-topography in a hillslope setting and varying sediment sources and velocities of deposition and requires future sampling and analysis at higher vertical resolution.

Interestingly, no indications for an upper sand unit have been found at the top of the modern headwall on the upslope part of the slump. Additionally, there are no significant Holocene deposits exposed at all in the upper part of the slump. However, carcasses of Equus lenensis and juvenile Bison priscus (dated to 4.45±0.35 and $\mathrm{8.22}+\mathrm{0.45}/-\mathrm{0.40}$ ¹⁴C kyr BP, respectively) were found at the exit of the slump (Murton et al., 2017). These observations may be related to (a) the palaeo-topographic situation and transport and deposition regimes of windblown material during MIS 3–1, likely decreasing upslope, (b) widespread permafrost degradation since the Late Glacial–Early Holocene transition, and/or (c) changed hydrological and vegetation patterns in the Holocene preventing transport and/or deposition of windblown sediment.

The lower ice complex is not dated so far. Given the dating results of the overlying lower sand unit, it might correspond to the Yukagir Ice Complex from Bol'shoy Lyakhovsky Island, which has been dated by radioisotope disequilibria (²³⁰Th∕U) in peat to 178±14, 221±27 (Wetterich et al., 2019) and 201±3 kyr (Schirrmeister et al., 2002). Hence, this unit most likely has survived both the last and current interglacials. Dating such an old ice complex is challenging, as the dating approaches are limited and partly still in development. Independently dated tephra layers have been proven useful to date old permafrost in east Beringia, i.e. in Alaska (Schirrmeister et al., 2016) and in Yukon (Froese et al., 2008). However, no tephra layers have been identified yet in the study region due to the long distance from active volcanos. But with its clearly defined stratigraphy and the access to very old permafrost and ice wedges, the Batagay megaslump may be a suitable site to evaluate and further develop the ice-wedge dating approaches using ${}^{\mathrm{36}}\mathrm{Cl}/{\mathrm{Cl}}^{-}$ (Blinov et al., 2009) as well as uranium isotope methods (Ewing et al., 2015). To better constrain the chronostratigraphy of the exposed units, a systematic dating approach is needed, which should additionally include independent dating of host sediments, e.g. by luminescence.

5.2 Cryostratigraphy at the Batagay megaslump – implications for ice-wedge growth regime

The shape and composition of the wedges in the Batagay megaslump correspond to the assumed genesis of the respective units and may also provide paleoclimate information (Table 4). Both the lower and upper ice complex units are characterized by ice wedges up to several metres wide, whereas both the lower and upper sand units contain tall narrow wedges. The former indicates more stable surface conditions with lower rates of sand aggradation and sufficient meltwater supply during the formation of both ice complexes that allowed more horizontal growth of ice wedges. The tall narrow wedges of both sand units, in contrast, point to rapid upslope-directed aeolian deposition of sand and therefore a predominant vertical ice-wedge growth in persistent polygonal patterns. The downslope inclination of narrow syngenetic wedges in the upper sand unit indicates upward growth subvertically at a right angle to the aggrading depositional (hillslope) surface. The tall and narrow wedges in the sand units are similar in size and shape to chimney-like sand wedges in thick aeolian sand sheets in the Tuktoyaktuk Coastlands of western Arctic Canada, attributed to rapid aggradation of aeolian sand and rapid wedge growth (Murton and Bateman, 2007). The composite wedges in the lower and upper sand units imply dry conditions with only little meltwater supply. However, as most ice wedges exhibit a significant sediment content, a steady supply of windblown sediment and a rather thin snow cover during ice-wedge formation is likely.

Table 4 Summary of past climate in the Yana Highlands. Palaeo-ecological interpretations are from Ashastina et al. (2018) and based on fieldwork in 2014.

At the top of the upper ice complex, ice wedges tend to narrow and partly transform into narrow composite wedges of the upper sand unit (Fig. S3). This indicates a rather gradual change of the deposition regime towards higher sedimentation rates and possibly lower meltwater supply due to drier conditions. Interestingly, the polygonal pattern has not been affected by this change, as indicated by the consistent frost-cracking positions. The upward transition of wedges from the upper ice complex to upper sand unit, however, was interrupted episodically by thaw, producing a number of thaw unconformities at different depths. In contrast, the erosional event truncating the lower ice complex seems to have changed the polygonal pattern in which the lower sand unit has been deposited.

A major erosional surface attributed to gullying by water flowing down a palaeo-hillslope is inferred from the large concave-up lower contact of the upper ice complex in the northeast part of the headwall (Fig. 3). Numerous gullies on the present hillslope near the megaslump are indicated on satellite images between 1968 and 2010 (Kunitsky et al., 2013), which suggests that gullies are characteristic landforms of such terrain under present environmental conditions. Water is supplied to such gullies by snowmelt in spring and rainfall and melt of ground ice in summer. Erosion of the underlying lower sand unit provided substantial accommodation space for development of an unusually thick upper ice complex with wide ice wedges above it. Stratigraphically, this erosional surface is at a similar level as the upper woody debris layer that is thought to be of last interglacial age (Ashastina et al., 2017, 2018), which in turn suggests that the erosion also probably took place during warm (i.e. interglacial) conditions.

5.3 Yana Highlands' ice-wedge stable isotopes – regional paleoclimate implications

Unfortunately, no comprehensive modern precipitation stable-isotope data are available for the Yana Highlands. The nearest sites with available data are Zhigansk, about 500 km to the west (Kurita, 2011); Tiksi, about 500 km to the northwest (Kloss, 2008); and Yakutsk, about 650 km to the south (Kloss, 2008; Papina et al., 2017). These sites show similar isotopic characteristics for the cold season (October to March). Mean δ¹⁸O values decrease along the north–south transect from about −29 ‰ in Tiksi to −31 ‰ in Yakutsk, and d varies strongly without a clear geographical pattern and shows values between about 8 ‰ and 15 ‰. In general, the mean winter δ values plot around the GMWL. However, these sites are characterized by less cold winters compared to the Yana Highlands. Furthermore, they are located east and south, respectively, of the Verkhoyansky range (Fig. 1), which constitutes an orographic barrier for atmospheric moisture transport by the westerlies and thereby contributes to the extreme cold and dry winters of the Yana Highlands. Hence, more depleted δ¹⁸O values can be expected for winter precipitation.

The Yana Highlands' ice-wedge stable-isotope data form two major clusters (Fig. 6, Table 2). The regression of the first cluster (cluster 1) plots below but mainly parallel to the GMWL ($\mathit{\delta }\mathrm{D}=\mathrm{7.45}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{10.46}$, R²=1) and comprises ice wedges B17-IW1 (lower ice complex), B17-IW5, B17-IW6 (both upper ice complex), and A17-IW3 (recent Adycha riverbank). The regression of the second cluster (cluster 2) plots clearly above the GMWL ($\mathit{\delta }\mathrm{D}=\mathrm{6.88}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{21.72}$, R²=0.97) and comprises all data from the upper sand unit (composite wedges B17-IW2, B17-IW3, ice wedge B17-IW4). Furthermore, all except one of the intrasedimental ice samples as well as all slump outflow samples – interpreted as integrated signal of the entire exposed Batagay sequence – plot well above the GMWL. In summary, the higher the sediment content of our samples, the higher the d values and hence the offset above the GMWL. The composite and ice wedges from the upper sand unit, likely formed in a generally dry environment (Sect. 5.2), exhibit high d values (cluster 2). In contrast, the ice wedges from cluster 1, likely formed under moister conditions with more abundant meltwater supply, show lower d values. This suggests that the divergence from cluster 1 (or the GMWL), i.e. higher d, might be an indicator of drier conditions.

Generally, the high d values indicate a higher kinetic fractionation at the moisture source, e.g. due to initial evaporation under low relative humidity and/or high sea-surface temperature conditions (Pfahl and Sodemann, 2014). As we are dealing with cold-stage MIS 3–2 ground ice, a changed moisture source region might be more likely than generally higher sea-surface temperatures compared to the Late Holocene as reflected by the ice wedge A17-IW3. As the cold season in east Siberia is characterized by higher d values than the warm season, one might conclude that cold-season precipitation does not only feed the ice and composite wedges but also contributes significantly to the formation of intrasedimental ice. However, secondary fractionation during snow cover evolution under extremely continental conditions with a generally lower snow cover, and snowmelt containing only a fraction of the initial snow might also affect d preserved in wedge ice (Opel et al., 2018; Grinter et al., 2019). For intrasedimental ice, too, multiple freeze–thaw cycles in the active layer have to be taken into account.

The regional palaeoclimatic implications for the Yana Highlands drawn from our ice-wedge δ¹⁸O data are summarized in Table 4 together with palaeo-ecological results from Ashastina et al. (2018). The lowest winter temperatures of our dataset are indicated by δ¹⁸O values around −35 ‰ for ice wedges from the upper ice complex, which formed during the MIS 3 interstadial (Kargin interstadial in Siberian stratigraphy). Similar values including d have been reported by Vasil'chuk et al. (2017) for three undated ice wedges from the upper ice complex sampled close to the headwall of the Batagay megaslump (Fig. 3). Slightly higher temperatures during the formation of the lower ice complex before MIS 6 are indicated by δ¹⁸O values around −33 ‰. However, given the small size of the corresponding ice wedge, the old age and rather similar isotope values of intrasedimental ice in the host deposits some post-depositional alteration of the originally recorded signal (Meyer et al., 2002a, 2010a) cannot be excluded. Hence, these values should be treated with caution. Comparable and slightly higher winter temperatures can be inferred for composite wedges (δ¹⁸O values around −33 ‰ and −32 ‰) and an ice wedge (δ¹⁸O values about −30 ‰) from the upper sand unit, which is attributed to the transition from MIS 3 to 2. Again, similar values have been reported by Vasil'chuk et al. (2017) for an undated narrow ice wedge from the upper sand unit sampled in another part of the Batagay megaslump (Fig. 3). However, the palaeoclimatic significance of the rather small composite wedges with few samples has to be treated with caution. Ice wedge B17-IW4 might even be an epigenetic ice wedge of Holocene age, given the quite similar δ¹⁸O values (about −29 ‰) of the definitively Holocene ice wedge from the Adycha riverbank (A17-IW3), which indicate higher winter temperatures in the last centuries to millennia. These values are close to those of modern winter precipitation at Tiksi, Zhigansk, and Yakutsk. As modern ice wedges are enriched by several per mille in δ¹⁸O compared to modern winter snow cover (Grinter et al., 2019; Opel et al., 2018), we suggest that modern winter precipitation in the Yana Highlands is isotopically distinctly more depleted than that at Tiksi, Zhigansk, and Yakutsk.

No palaeoclimatic information based on ice-wedge isotopes is available yet from Batagay for the MIS 4 (Zyryan) and MIS 2 (Sartan) stadials. On Bol'shoy Lyakhovsky Island of the New Siberian archipelago, ice wedges from both periods show δ¹⁸O values about 6 ‰ lower than those of MIS 3 Ice Complex ice wedges (Wetterich et al., 2011, 2019; Opel et al., 2017b). Hence, even more depleted δ¹⁸O values might be expected in ice wedges of Last Glacial Maximum age at the Batagay megaslump. However, we note that distinctly depleted δ¹⁸O values for MIS 4 and 2 have been found only at Bol'shoy Lyakhovsky Island.

When interpreting ice-wedge stable isotopes in terms of absolute palaeotemperatures, one has to keep in mind that the isotopic composition of the oceanic moisture source has changed over glacial–interglacial timescales, with enrichment in δ¹⁸O in the ocean water during past cold stages. Correcting the ice-wedge stable isotopes for this effect requires detailed knowledge about the moisture source region and in particular a better constrained chronology of ice-wedge isotopes than is available from our dataset.

In general, the inferred lower winter temperatures meet the expectations of an increased continentality during the Late Pleistocene cold stages as well as partly drier conditions indicated by the formation of composite wedges. Palaeo-ecological analysis (plant macrofossil remains, pollen, and beetles) of Batagay megaslump deposits (Ashastina et al., 2018) indicates generally much drier conditions compared to recent times (Table 4), in particular during the cold stages. Tree and beetle indicators imply warm summers for most units, whereas macrofossil remains indicate cold to very cold winters for the upper ice complex. These patterns are generally in line with our ice-wedge interpretation and support the hypothesized higher continentality.

5.4 Large-scale implications for paleoclimate and past continentality

5.4.1 MIS 3 interstadial wedge ice records

To assess the climate and continentality across much of west Beringia during the MIS 3 interstadial, we compare our Yana Highlands ice-wedge stable-isotope data to a dataset of 17 other ice-wedge sites (Figs. 1, 7). Most of the wedges are from the Yedoma Ice Complex (Schirrmeister et al., 2011b, 2013), which is widely distributed in east and central Siberia. We selected all available ice wedges per study site that provide both δ¹⁸O and d data (Table S3). Additionally, all except two of the considered wedges or their host sediments have been directly dated by radiocarbon methods to between about 50 and 30 kyr ago (MIS 3 interstadial). We appreciate that an ice wedge may not contain a full record of this time period and that the climate was not uniform throughout this period. It is, however, not clear yet whether millennial-scale climate oscillations known from the North Atlantic region (e.g. Dansgaard et al., 1993) also characterized the MIS 3 in east Siberia. Following Murton et al. (2017), there are hints that the MIS 3 in northern areas of northeast Siberia was a time of general climate stability even though the interstadial climate was not monolithic. In particular, a warmer period around 40 kyr BP is well known from several palaeo-ecological proxies for some of the study sites (Wetterich et al., 2014) and from ice wedges of the well-dated Bykovsky Peninsula section (Meyer et al., 2002a). The Batagay upper ice complex ice wedge B17-IW5 shows δ¹⁸O values of around −35 ‰, which is 6 ‰ to 2 ‰ lower than MIS 3 ice-wedge δ¹⁸O data from northern Yakutia, the Kolyma region, and central Yakutia (Fig. 8). As hypothesized, the relatively depleted δ¹⁸O values from Batagay clearly point to lower winter temperatures during MIS 3 in the Yana Highlands than elsewhere, probably due to the more continental location framed by the partially glaciated Verkhoyansky and Chersky ranges reaching elevations of more than 2000 m above sea level. Even the ice-wedge δ¹⁸O data from highly continental sites in central Yakutia do not reach such low δ values as observed in the Batagay ice wedges. Interestingly, the lowest MIS 3 δ¹⁸O values beside Batagay originate from the northeastern-most site, Novaya Sibir' island. If an Atlantic moisture source for MIS 3 wedge ice is assumed, Novaya Sibir' would represent the farthest location from the source region and therefore show strongly depleted δ¹⁸O values, almost reaching those of Batagay.

Figure 7 Comparison of stable-isotope data from Siberian ice wedges attributed to the MIS 3 interstadial (about 50 to 30 kyr), comprising data from the western Laptev Sea (Magens, 2005), the Lena River Delta (Schirrmeister et al., 2003, 2011a; Wetterich et al., 2008; Thomas Opel, unpublished data), the central Laptev Sea (Meyer et al., 2002a; Hanno Meyer and Thomas Opel, unpublished data; Schirrmeister et al., 2017), the New Siberian Islands and the Laptev Strait (Lutz Schirrmeister and Hanno Meyer, unpublished data; Wetterich et al., 2014; Opel et al., 2017b), the Kolyma Lowland (Strauss, 2010), the Yana Highlands (this study; Vasil'chuk et al., 2017), and central Yakutia (Lutz Schirrmeister, unpublished data; Popp et al., 2006). The new data from Batagay are marked in red. Further information is given in Table S3.

The peculiarity of the Yana Highlands is also apparent in the d data of the MIS 3 ice wedges. The Batagay ice wedges exhibit the highest d values, between 8 ‰ and 11 ‰, in the dataset. Similar values have been found in ice wedges from central Yakutia, whereas today's coastal lowlands of northern Yakutia, the New Siberian Islands, and the Kolyma region sites show values between 2 ‰ and 7 ‰ (Fig. 9). Again, from all northern sites, Novaya Sibir' island shows the value closest to those of the Yana Highlands and central Yakutia.

Figure 8 Map of mean δ¹⁸O data from Siberian ice wedges attributed to the MIS 3 interstadial (about 50 to 30 kyr). Site IDs are given in Fig. 7 and further details in Table S3.

Figure 9 Map of mean d data from Siberian ice wedges attributed to the MIS 3 interstadial (about 50 to 30 kyr). Site IDs are given in Fig. 7 and further details in Table S3.

Little is known about d variation in northeast Siberian ice wedges in general, and particularly within specific time periods. Previous studies have interpreted changes in the d values between the Late Pleistocene (low values) and Holocene (higher values) as changes in the moisture sources and transport pathways following the retreat of the Eurasian ice sheets (Meyer et al., 2002a, b). Higher d values in central Yakutia and the Verkhoyansk mountain forelands have also been found in modern winter snow along a west-to-east transect over east Siberia from Yakutsk to Magadan (Kurita et al., 2005). This spatial pattern has been attributed to the influence of different moisture sources, with greater kinetic effects during initial evaporation, as resulting from a lower relative humidity (Kurita et al., 2005). Analogously, during MIS 3, the Yana Highlands and central Yakutia may have received precipitation from other moisture sources and/or transport pathways compared to more northern sites. Due to their more southern location, the Yana Highlands might have received a contribution of moisture from a more southerly source, e.g. the Mediterranean. Such moisture usually has higher d values compared to North Atlantic moisture (Gat et al., 2003) and may have reached this region via a different southern moisture transport pathway.

5.4.2 Holocene

The mean δ¹⁸O values of the Holocene ice wedge sampled at the Adycha River also point to lower winter temperatures – due to greater continentality – than other Holocene ice wedges in the dataset. The Adycha samples (mean δ¹⁸O −29 ‰) are more depleted than Holocene wedges from northeast Siberian coastal sites, which show typically mean δ¹⁸O values between −27 ‰ and −23 ‰ (e.g. Meyer et al., 2002b, 2015; Opel et al., 2017a; Wetterich et al., 2018). Holocene ice wedges from the Bykovsky Peninsula (Meyer et al., 2002a) show exceptional values (mean δ¹⁸O of −28 ‰ to −26 ‰), which do not fit into the regional pattern and may be related to local climate or environmental effects. In contrast, Holocene ice wedges from central Yakutia exhibit mean δ¹⁸O values of −28 to −26 ‰ (Popp et al., 2006), which are closer to those from the Yana Highlands and also reflect the greater continentality and lower winter temperatures than the Holocene coastal sites.

The d values from Holocene ice wedges in Siberia are highly variable. Mean d values at the coastal sites are higher (compared to MIS 3), between 4 ‰ and 12 ‰ (e.g. Opel et al., 2011; Meyer et al., 2002b; Wetterich et al., 2008, 2018). In contrast, in the Yana Highlands, the Late Holocene mean d value of 6 ‰ is lower compared to MIS 3, whereas Late Holocene d values in central Yakutia (between 6 ‰ and 11 ‰) are rather similar to MIS 3 values (Popp et al., 2006). The exceptional ice wedges from the Bykovsky Peninsula (Meyer et al., 2002a) even reach d values of 16 ‰. In summary, Holocene ice-wedge d values do not show a clear spatial pattern that is readily interpretable in terms of moisture sources or continentality.

6 Conclusions and outlook

Our stable-isotope data clearly show that winter temperatures in the Yana Highlands during MIS 3, represented by the upper ice complex (Yedoma) at Batagay, and the Holocene were distinctly lower than at other ice-wedge study sites in coastal and central Yakutia. This indicates the persistence of enhanced continentality of the Yana Highlands region during at least part of the Late Pleistocene in west Beringia and during the Holocene. The stable-isotope data from narrow composite wedges of the upper sand unit (MIS 3–2) and an old ice wedge from the lower ice complex (pre-MIS 6) are less indicative and require additional studies.

High-resolution systematic sampling and dating now need to be carried out for all cryostratigraphic units of the Batagay permafrost sequence to validate our findings of increased continentality during MIS 3 and the Holocene, to improve the temporal resolution of the Batagay ice-wedge record, and to elucidate the palaeoclimatic history from other time slices. Of particular interest are the lower and uppermost parts of the upper ice complex, likely representing MIS 4 and 2 stadials, as well as more detailed studies of the yet-undated lower ice complex. To establish reliable ice-wedge chronologies for the upper ice complex and upper sand unit, radiocarbon dating should include macro remains, dissolved organic carbon, and CO₂ from gas bubbles (Lachniet et al., 2012; Kim et al., 2019). Additionally, previous luminescence dating results (Ashastina et al., 2017) need to be validated.