Northeast Siberian Permafrost Ice‐Wedge Stable Isotopes Depict Pronounced Last Glacial Maximum Winter Cooling

Stable isotopes (δ18O, δD) of wedge ice hold potential to reconstruct past winter climate conditions. Here, we present records of the marine isotope stages (MIS) 3 and 2 including the last Glacial maximum (LGM) from Bol’shoy Lyakhovsky Island (NE Siberia). MIS 3 wedge ice dated from 52 to 40 Kyr b2k varies between −32 and −29‰ in δ18O. Colder LGM conditions are implied by δ18O of −37‰ around 25 Kyr b2k. Similar Deuterium excess values indicate comparable moisture sources during MIS 3 and MIS 2. Regional LGM climate reconstructions depend on the seasonal resolution of the proxies and model simulations. Our wedge‐ice record reflects coldest winters during global minima in atmospheric CO2 and sea level. The extreme LGM winter cooling is not represented in model projections of global LGM climate where West Beringia shows noticeably little cooling or even warming in mean annual temperatures compared to the late Holocene.

temperatures around 21 Kyr BP in West Beringia based on pollen records and modeling efforts (Weitzel et al., 2021). Similar-to-modern LGM summer temperatures between 20 and 18 Kyr BP have been reported by V. D. Meyer et al. (2017) on Kamchatka Peninsula in far-east Russia based on biomarker proxy data from marine sediments.
LGM proxy records from West Beringia are generally sparse and have methodological limitations such as proxy uncertainties and seasonal biases which challenge the comprehensive understanding of LGM timing, magnitude, and spatial patterns on larger scales.
To reconstruct Beringian paleo-environments and climate such as those of the LGM, proxy data from permafrost deposits and ground ice are essential. Permafrost preserves both summer-specific fossils Kienast et al., 2011) and winter-specific wedge ice (e.g., Opel et al., 2018). The use of wedge ice stable oxygen (δ 18 O) and hydrogen (δD) isotopes as a winter temperature proxy arises from its distinct formation: Spring snowmelt fills polygonal winter thermal contraction cracks in the ground (i.e., ice-wedge polygons) and re-freezes, thus forming vertical ice veins (Leffingwell, 1915) that contain the isotopic composition of the winter time snowpack. By repeated frost cracking over long periods and with ongoing sedimentation, ice wedges grow syngenetically in width and height. Variations in stable water isotopes preserve WETTERICH ET AL.  (Dyke et al., 2003;Manley & Kaufman, 2002;Svendsen et al., 2004), exposed shelf areas during LGM sea level lowstand of −120 m, and the position of Bol'shoy Lyakhovsky Island (red circle); (b) map of the study site (Landsat 5, June 27, 2009). The white line indicates the transect shown in (b); (c) coastal relief profile (after Andreev et al. (2009) with ice-wedge (triangles) and sediment (squares) profiles. (d) photograph of MIS 2 Yedoma Ice Complex hosting LGM ice wedges (dotted red rectangle in (c). information on air temperature and moisture source for the cold season, that is the meteorological winter and spring (December to May, DJFMAM) as defined by H. Meyer et al. (2015).
To use ice wedges in paleoclimate research, chronologies are obtained from radiocarbon dating of fossil remains, CO 2 , or dissolved organic carbon (DOC) preserved inside the wedge ice (e.g., Lachniet et al., 2012). Host sediments provide further indication of the ice-wedge formation age (Opel, Wetterich, et al., 2017). Dating of ice wedges can provide up to centennial resolution as shown for the Younger Dryas (H. Meyer et al., 2010). Stacked records of dated Holocene stable isotope samples combined from multiple sites (Holland et al., 2020;H. Meyer et al., 2015) or from a single site (Opel, Laepple, et al., 2017) represent state-ofthe-art approaches to deduce reliable ice-wedge chronologies. However, wedge-ice chronologies from the late Pleistocene are less constrained (H. Opel, Wetterich, et al., 2017 and scarce (Porter & Opel, 2020), but crucial for a better regional and seasonal assessment of LGM climate.
Here, we present radiocarbon ages and stable isotope data of ice wedges from Yedoma Ice Complex of MIS 3-2 age exposed at the southern coast of Bol'shoy Lyakhovsky Island (New Siberian Archipelago) ( Figure 1). Specifically, we seek (1) to quantify isotopic changes in glacial winter climate conditions from interstadial MIS 3 to stadial MIS 2, (2) to constrain the timing of the coldest winter conditions during the LGM, and (3) to present new insights into LGM cooling and regional seasonality patterns in Beringia.

Ice-Wedge Sampling
Bol'shoy Lyakhovsky is the southernmost island of the New Siberian Archipelago at the Dmitry Laptev Strait and located in West Beringia (Figure 1a). The southern coast near the Zimov'e River mouth exposes late Quaternary permafrost and has been intensively studied in recent decades (e.g., Schirrmeister et al., 2011;Tumskoy et al., 2012).
The Bol'shoy Lyakhovsky permafrost sequences exposed in coastal bluffs discontinuously cover approximately the last 200 Kyr Wetterich et al., 2019). Four distinct Ice Complex generations of MIS 7, MIS 5, MIS 3, and MIS 2 age are exposed, intersected by floodplain deposits of MIS 6 and MIS 4 age as well as by MIS 5 thermokarst deposits, and covered by Holocene deposits Wetterich et al., 2009Wetterich et al., , 2016Wetterich et al., , 2019. Substantial parts of the exposed permafrost belong to the Yedoma Ice Complex of MIS 3-2 age (H. Meyer, Dereviagin, Siegert, Schirrmeister, et al., 2002;Wetterich et al., 2011Wetterich et al., , 2014, which is the focus of this study. Modern wedge-ice growth takes preferentially place in thermokarst basins and on the Zimov'e river floodplain. In total, seven ice wedge profiles from MIS 3-2 Yedoma Ice Complex are considered in our study (Figures 1b  and 1c). The profiles exposed over a total distance of about 2.5 km were sampled in 1999, 2007, and 2014 either by ice screw or by chainsaw. More details are given in the supporting information S1.

Radiocarbon Dating
The radiocarbon dating of organic material from host deposits and wedge-ice inclusions was performed at CologneAMS (University of Cologne, Germany; Rethemeyer et al., 2019), at Leibniz Laboratory for Radiometric Dating and Stable Isotope Research (Kiel, Germany; Grootes et al., 2004) and at MICADAS (Mini Carbon Dating System, AWI, Bremerhaven, Germany; Opel et al., 2019). In total, eight new radiocarbon dates from distinct floral or faunal macrofossils (>1 mm in size) and seven re-calibrated dates from H. Meyer, Dereviagin, Siegert, Schirrmeister, et al. (2002) and Wetterich et al. (2011) are considered (Table S1). Ages are given as median values of calibrated years before the year 2000 (yr b2k) using the IntCal20 calibration data set (Reimer et al., 2020).

Ice-Wedge Stable Isotopes
Melted ice samples were analyzed for oxygen (δ 18 O) and hydrogen (δD) stable isotopes using a Finnigan MAT Delta-S mass spectrometer (1σ < 0.1‰ for δ 18 O, 1σ < 0.8‰ for δD; H. Meyer et al., 2000). The values are reported as per mil (‰) difference from the VSMOW (Vienna Standard Mean Ocean Water). The deuterium excess (d) is calculated according to Dansgaard (1964) in Equation 1: In total, 357 samples from seven ice wedges were analyzed (Table S2; Figure 2). The study focuses on δ 18 O as it is more commonly used in Arctic ground ice paleoclimate studies. Respective δD values are given in Table S2.

Ice-Wedge Formation Ages
The formation of the MIS 3 Yedoma Ice Complex of Bol'shoy Lyakhovsky is constrained by numerous radiocarbon dates from frozen deposits. Two 15-m-thick vertical sequences reflect continuous Yedoma formation between >50 and 32.7 Kyr b2k (LYA-97; Andreev et al., 2009) (c) LGM (bluish symbols) ages. Directly dated wedge-ice R9 (black diamonds) and TZ-2-4 (black crosses; H. Meyer, Dereviagin, Siegert, Schirrmeister, et al., 2002) are shown for comparison in (a). L7-07 data shown in (c) are from Wetterich et al. (2011). Radiocarbon-dated ice wedge profiles of (d) δ 18 O and (e) deuterium excess (d) delineating the MIS 3 to MIS 2 transition toward LGM are shown. Orange asterisks in (d) indicate the dated samples. Profile means and 1σ rages in (e) are indicated by lines and dotted lines, respectively. Note that data from L14-13 and from L14-07 derive from ice wedges that were sampled only in their central parts. similar geomorphologic positions west and east of the Zimov'e River ( Figure 1c). Scarce radiocarbon ages from MIS 3 Yedoma wedge ice range from 51.7 to 39.5 Kyr b2k (H. Meyer, Dereviagin, Siegert, Schirrmeister, et al., 2002). Newly obtained dates from woody inclusions in the older part of ice wedge L6 fit into that range (43.3 Kyr b2k in L6-A and 43.9 Kyr b2k in L6-B) ( Table S1).
Evidence of continued Ice Complex accumulation during MIS 2 on Bol'shoy Lyakhovsky was only found in slope positions of the Zimov'e River valley, but not on top of nearby MIS 3 Ice Complex (Figure 1c). This change in accumulation area can be explained by a lowered erosion base during LGM sea level lowstand and corresponding changes in the hydrological system of areas with higher relief inclination leading to preferential Ice Complex formation in valley areas . On the Zimov'e River valley slope, previously published radiocarbon dates  and new data from in situ Rangifer tarandus bones confirm MIS 2 Yedoma Ice Complex formation from about 30 to 26 Kyr b2k. However, one age of 34 Kyr b2k (Equus sp.) bone is interpreted as re-deposited, possibly by downslope slumping from above-lying MIS 3 deposits (Table S1). New dates from MIS 2 Yedoma ice wedges show ages of 27.5 Kyr b2k (L14-13), 25.9 and 24 Kyr b2k (both L14-07; Table S1). Thus, the overall MIS 2 Ice Complex formation spans from 30 to at least 24 Kyr b2k.
Syngenetic wedge ice is younger than the host deposits at the same elevation level due to the ice-wedge formation in vertical frost cracks within older sediments (Porter & Opel, 2020). Such age offset is systematic and defined by the cracking depth and permafrost aggradation rates. In the MIS 2 Yedoma record of Bol'shoy Lyakhovsky Island, three in situ R. tarandus bones sampled at 7 masl range between 26.2 and 26 Kyr b2k, while the dates of ice wedge L14-07 at the same sampling height are slightly younger with 25.9 Kyr b2k and 24 Kyr b2k (Table S2). The corresponding ages of MIS 2 host deposits and wedge ice suggest reliable age information.

Stable Isotope Composition of Wedge Ice
Yedoma wedge ice of MIS 3 and MIS 2 age from Bol'shoy Lyakhovsky has distinct stable isotope compositions (Table S2) , we hypothesize a change in the prevailing frost cracking dynamics possibly due to a change in the polygonal pattern, and, hence, rather lateral growth for ice wedge L6, although the L6-C part has not been dated. We regard the L6-B age as derived from re-deposited organic material, and consequently omit it from interpretation. Irregular growth differing from the central ice-wedge growth model (Lachenbruch, 1962) has been previously reported (Lachniet et al., 2012;Mackay & Burn, 2002;Opel et al., 2011), and is thought to explain the directed isotopic trend from L6-A to L6-C by a presumably lateral growth direction. The addition of younger ice laterally (L6-C) is also supported by field observations ( Figure S3). The L6-C wedge ice is isotopically very close to the L14-13 wedge ice record (27.5 Kyr b2k; δ 18 O of −34.1 ± 0.2‰, d of 5.9 ± 0.4‰; Table S2) although the mean d value in L6-C (7.2 ± 0.9‰) is slightly higher (Figure 2e). Both ice wedges, L14-13 and L6-C were sampled in similar stratigraphic and morphologic positions west and east of the Zimov'e River mouth, respectively ( Figure 1c). Thus, given the very similar isotopic composition, we assume a formation time for both ice wedges L14-13 and L6-C during the early MIS 2 reflecting a winter cooling trend.
LGM wedge ice was encountered in ice wedge profile L14-07 (25.9 and 24 Kyr b2k) as well as in ice wedges L7-07-1 and L7-07-2 of the same isotopic composition (Figure 2c). These records represent LGM winter conditions with most depleted mean δ 18 O values of up to −37.4 ± 0.4‰ and mean d between 6.1 and 7.3‰ (Table S2).

Paleoclimatic Implications of Wedge-Ice δ 18 O and Deuterium Excess
The variations in ice wedge stable isotopes from Bol'shoy Lyakhovsky correspond to large-scale interstadial MIS 3 and stadial MIS 2 climate conditions, including the LGM (Figure 3) (Table S1) on Bol'shoy Lyakhovsky Island shown together with modeled cold-season surface temperatures of the Laptev Sea region (DJFMAM, Ganopolski & Calov, 2011). The LGM period is shown as gray-shaded area.
imply relatively higher winter temperatures and indicate more variable interstadial conditions between >51.7 and 39.5 Kyr b2k.
The deuterium excess in wedge ice is influenced by local conditions (Opel et al., 2018) as well as by isotopic smoothing with time, but is commonly applied to deduce changes in moisture generation and transport. Thus, ice wedge d might serve as a potential indicator of shifting moisture sources related to changing major atmospheric circulation patterns (H. Meyer et al., 2010) or sea ice cover. The mean d values in MIS 3-2 wedge ice overlap. The MIS 3 mean values of d vary only between 6.8 ± 0.7‰ and 9.3 ± 0.9‰ and those for MIS 2 even less between 5.9 ± 0.4 ‰ and 7.3 ± 0.4‰ (Table S2). We therefore infer that winter moisture sources for Bol'shoy Lyakhovsky Island remained relatively stable during MIS 3-2. Consequently, variations in MIS 3-2 wedge-ice δ 18 O are thought to reflect winter temperature changes. We assume the North Atlantic and potentially the North Pacific as main moisture sources with varying contributions to the precipitation eventually forming the ice wedges. To which extent both moisture sources finally contributed to winter precipitation of our study region during MIS 3-2 is yet unknown. Further, minor contributions of proximal primary and secondary sources are possible but unlikely due to the persistent sea ice cover and frozen open water bodies in the winter season.
The interstadial MIS 3 wedge ice of Bol'shoy Lyakhovsky varies over about 3‰ in δ 18 O, which is less than half the amplitude seen in the NGRIP ice core for the same time slice that is driven by prominent Dansgaard-Oeschger (D-O) events: rapid climate oscillations that originated in the North Atlantic region (Figure 3c;Dansgaard et al., 1993). Those are not identified but only hypothesized in Siberian climate archives yet (Baumer et al., 2021). The temporal resolution of ice wedge stable isotope records would not allow recognizing D-O events yet even though they might have impacted East Siberian climate dynamics. The different range in MIS 3 isotopes might also relate to the differences in seasonality represented in the two archives and spatial coverage. Bol'shoy Lyakhovsky wedge ice reflects regional winter temperatures while the Greenlandic ice rather represents the North Atlantic annual temperatures and does not specifically resolve the winter-season. The NGRIP LGM signal is not characterized by most depleted isotope values, but varies within its MIS 3 range. In contrast, Bol'shoy Lyakhovsky LGM wedge ice δ 18 O values are up to 9‰ lower than MIS 3 maxima (Figure 3d).
In general, distinct LGM wedge-ice records showing extremely cold winter conditions are rare on the circum-arctic scale (Porter & Opel, 2020). In eastern Beringia, Lachniet et al. (2012) found a 6‰ decrease in δ 18 O within one ice wedge that formed between 28 and 22 Kyr b2k, reflecting LGM cooling. More often, similar MIS 3 and MIS 2 isotopic signatures are found (Porter & Opel, 2020 et al., 2020). Hitherto, it is not sufficiently resolved yet, whether this is due to a less cold LGM climate in the region or poor LGM representation in the studied ice wedge profiles. Therefore, the well-dated Bol'shoy Lyakhovsky LGM wedge-ice record is exceptional and requires further explanation. Porter and Opel (2020) emphasize that rare LGM wedge-ice records might be explained by local winter conditions affecting ice-wedge formation when thick snow cover potentially prevented frost cracking. Alternatively, thin or no snow cover and deflation at exposed terrain surfaces by enhanced winter storm activity (Guthrie, 2001) might not have provided sufficient melt water for wedge-ice formation during the LGM. For our record, we consider that the slope position of MIS 2 Yedoma Ice Complex may have promoted accumulation of snow drift potentially leading to an overrepresentation of isotopically more depleted precipitation from the coldest period of the winter time that became preserved in the LGM wedge ice. If so, during MIS 3 and MIS 2 the local conditions of Ice Complex formation differed considerably and could therefore have superimposed the general interstadial-stadial winter climate pattern.
However, both the interstadial and the stadial climate patterns as reflected in our wedge-ice records are supported by variations in global atmospheric CO 2 concentrations (Lüthi et al., 2008) and in sea level reconstructions (Spratt & Lisiecki, 2016) (Figure 3). The Bol'shoy Lyakhovsky ice-wedge record is further in line with modeled cold-season air temperatures of the Laptev Sea region (Ganopolski & Calov, 2011). This model output is based on an Earth system model of intermediate complexity, mainly forced by Earth's orbital parameters and atmospheric concentration of major greenhouse gases (Ganopolski & Calov, 2011;Petoukhov et al., 2000). Although a calibrated relationship between wedge-ice δ 18 O and cold-season air temperature has not been developed yet for northern Eurasia (Opel et al., 2018), the isotopic depletion of 9‰ in δ 18 O between warmest interstadial MIS 3 and coldest stadial MIS 2 periods would correspond to modeled cold-season LGM temperatures up to 6°C lower than during MIS 3 ( Figure 3) if temperature alone controlled the wedge-ice δ 18 O.
Holocene ice-wedge records show increasing δ 18 O corresponding to the mid to late Holocene winter warming that is related to increasing cold-season insolation and partially by rising greenhouse gas concentrations (Holland et al., 2020;H. Meyer et al., 2015;Opel, Laepple, et al., 2017). In contrast, the LGM winter cooling appears not to be controlled by cold-season insolation, which increased over this period ( Figure 3). Instead, the overall LGM cooling was driven by a strong minimum in summer and annual insolation in mid to high northern latitudes due to orbital cycles (Milankovitch, 1941). Between about 53 and 22 Kyr b2k, the decrease in summer insolation of about 15 W m −2 dominated over the winter insolation increase of about 9 W m −2 (Figure 3) leading finally to overall full glacial conditions of the LGM. Less pronounced orbitally driven seasonality during MIS 2 would imply smaller temperature amplitudes between LGM summers and winters in contrast to higher seasonality with higher seasonal temperature amplitudes during MIS 3. However, this potential effect seems superimposed by increased continentality with a stronger Siberian High and stronger winter cooling for MIS 2 in the study region as a result of the northern hemispheric configuration of ice sheets, sea ice coverage on the Arctic Ocean and global sea level lowstand (Guthrie, 2001;Murton et al., 2015).
The substantial LGM winter cooling under highly continental conditions as reflected in the Bol'shoy Lyakhovsky ice-wedge stable isotope record complements the results of recent LGM climate modeling studies. Modeled annual LGM temperatures (Tierney et al., 2020) and summer-season LGM temperatures (Weitzel et al., 2021) as well as biomarker-based summer-season temperature reconstructions (V. D. Meyer et al., 2017) indicate only slightly colder or even warmer LGM conditions than today in East Siberia (West Beringia) for periods slightly post-dating the present winter-season record. By highlighting the seasonal biases of the different studies, the apparent differences imply the need for well-dated regionally constrained seasonal-scale temperature reconstructions and climate model output for holistically assessing the climate conditions of particular time periods such as the LGM.

Conclusions
Directly dated ice wedges are the only conclusive winter climate archives for the late Pleistocene in Beringia. However, winter climate variability between interstadial MIS 3 and stadial MIS 2, including the last Glacial maximum, has not yet been well constrained by radiocarbon-dated ice wedge stable isotope records. One exception from West Beringia is found on Bol'shoy Lyakhovsky Island. Here, interstadial MIS 3 winter climate variability is reflected by mean δ 18 O values between −32 and −29‰ (dated to 51.7 to 39.5 Kyr b2k) followed by early MIS 2 cooling reflected by mean δ 18 O values from −35 to −34‰ (dated to 27.5 Kyr b2k) and the full LGM with −37‰ in δ 18 O (dated to 25.9 to 24 Kyr b2k). However, as chronologically less constrained records further west in the Central Laptev Sea coastal region do not show this prominent LGM winter temperature minimum, and local conditions in snow distribution and consecutive wedge-ice formation might have affected the finally preserved isotopic signal in addition to climate.
Increased continentality in the study region was caused by the LGM sea level minimum. Together with the northern hemispheric configuration of ice sheets and a persistent Arctic Ocean sea ice cover in winter, a stronger Siberian High and thus stronger winter cooling are likely. The Bol'shoy Lyakhovsky LGM icewedge record reflects coldest winter conditions in NE Siberia during global minima in atmospheric CO 2 and sea level. Cool and dry LGM summers as deduced from biological proxy data and exceptionally cold LGM winters as seen in the present study correspond to low seasonality with smaller seasonal temperature amplitudes compared to MIS 3.