| Introduction

The presence of ice wedges in permafrost regions is influenced by climatic conditions, soil properties, topography, and landscape history. Ice wedges typically form in regions with continuous permafrost, where mean annual air temperatures remain sufficiently cold, generally below -4°C to -6°C, to allow sufficient ground-cracking conditions and near-surface permafrost stability [1][2][3]. Soil moisture is also crucial, as ice wedges primarily develop in organic and organic-rich, finegrained saturated soils that undergo pronounced wintertime temperature fluctuations with extended cold snaps [4][5][6]. Thermal contraction cracking occurs during the frozen period, followed by snowmelt in the spring and subsequent filling of the cracks with water that freezes. Over time, repeatedthough often irregular-thermal contraction cracking allows thin ice veins to accumulate and gradually enlarge into interconnected ice wedge networks, driving the development of ice wedge polygonal terrain [4,7,8]. Topographic factors, such as low-lying areas with poor drainage, create conditions conducive to ice wedge formation but ice wedges can also form on hillslopes [9,10]. Additionally, changes in snow cover and vegetation, which can insulate the ground and limit the depth of freezing, affect the thermal regime of permafrost, further influencing ice wedge distribution [11].

The distribution of active and inactive ice wedges in the Arctic and sub-Arctic is closely tied to climate gradients, ecosystem properties, landscape features, landscape history, and disturbance regimes [2,3,12]. Active ice wedges are primarily found in regions with continuous and cold permafrost, where conditions still allow for the regular freeze-thaw cycles necessary for their growth and preservation [13]. Numerous examples of active ice wedge growth exist across the northern circumpolar permafrost region: Alaska [14], Canada [6], Siberia [15], and Svalbard [16]. Observations of inactive ice wedges or ice wedges preserved outside of environmental conditions conducive to their growth and/ or complete degradation are less common [2,13]. The presence of ice wedges in warm permafrost regions offers a valuable opportunity to study the vulnerability of different types of ground ice to rising air temperatures and to better understand the resilience and sensitivity of permafrost to climate warming and surface disturbances [17][18][19].

Warm permafrost, where ground temperatures are close to 0°C, is particularly vulnerable to degradation, which has direct implications for preservation of ice wedges at the southern fringe of the permafrost region. In these regions, ice wedges are preserved primarily through a delicate balance between ground temperature, soil moisture, and surface insulation [20]. Cold winters and cool summers are essential for maintaining ice wedge integrity, but even small increases in active layer thickness can lead to permafrost degradation, causing the melting of ice wedges [21,22]. In some cases, permafrost and ice wedges can remain warm but stable due to thick organic layers that buffer sub-surface temperature changes, helping to preserve ice wedges despite warming air temperatures [13,23,24].

Peatlands, characterized by saturated, organic-rich soils, create favorable conditions for ice wedge formation by retaining moisture and exhibiting low tensile strength when frozen-factors that promote frost cracking and ice wedge growth [13,25]. In peat-rich areas, the thermal and mechanical properties of organic matter are critical: frozen, saturated peat in winter conducts heat efficiently, facilitating ground cooling, while thawed, dry peat in summer acts as an insulator, limiting heat input. This seasonal contrast helps preserve relict permafrost and supports syngenetic ice wedge development, where wedges grow alongside accumulating peat [26][27][28][29]. However, in warming regions, these ice wedges are increasingly vulnerable to degradation. The dynamic interaction between organic soils and permafrost processes presents both challenges and opportunities for understanding ice wedge stability in a changing climate [30].

The Yukon-Kuskokwim Delta (YKD), spanning ~75,000 km 2 in southwestern Alaska, is located at the southern margin of the discontinuous permafrost zone and within the relict northern circumpolar Yedoma permafrost region ([31]; Figure 1). In our recent studies focused on assessing the agricultural potential of permafrost terrain in Bethel, Alaska, we synthesized permafrost information from several decades' worth of geotechnical reports. One of the reports indicated the presence of massive ground ice, which prompted subsequent remote sensing-based analyses and fieldwork focused on better characterizing nearsurface permafrost conditions in Bethel. Occurrence of massive ground ice in or around several adjacent villages (e.g., Chevak, Eek, Quinhagak, and Upper Kalskag) was also mentioned in available geotechnical reports [32]. However, most publications neglect the occurrence of ice wedges in the YKD, or mention the lack of extensive ice-wedge polygons as a defining characteristic of this predominantly tundra region [33][34][35].

In this study, we present the first map of ice wedge networks for the Bethel area based on analysis of 20-cm resolution orthoimagery from 2018. We provide ground truth for this map based on three field campaigns conducted in 2023 and 2024 focused on drilling to confirm the presence and morphology of wedge ice as well as coring to describe permafrost cryostratigraphy and ground-ice characteristics of frozen soils in polygonal centers. In addition, we incorporated ground penetrating radar (GPR) surveys to add context for our remote sensing-based ice wedge networks map. We also present data on the stable oxygen and hydrogen isotopic composition of wedge ice, radiocarbon ages of organic material preserved in near-surface permafrost, the presence of charcoal, and a near-surface tephra to date the timing of ice wedge formation and the potential of our studied permafrost sequences as a paleoclimate archive. Scaling our observations using very high-resolution satellite imagery shows that ice wedges are more widespread in the permafrost areas of the YKD than previously thought. The discovery of ice wedge networks on the YKD has implications for the resilience of ice-rich permafrost to warming while also contributing to assessments of the agricultural potential of the Bethel area and the broader YKD region.

| Study Area

The climate in the YKD region is subarctic with long, cold winters and short, mild summers. Historical records show a warming trend in mean annual air temperature (MAAT). In Bethel, the long-term (1923-2023) MAAT has already warmed from -1.2°C to -0.3°C as of the most recent standard climatological period (1991-2020) [36]. This warming has sped up since the 1950s, increasing at a rate of 0.4°C per decade, and is predicted to nearly double to 0.8°C per decade in the 21st century [37]. The warming is more pronounced in winter, with the cold season warming 60% faster than the warm season, leading to shorter winters and longer summers [38]. Annual precipitation has increased at a rate of 14.8 mm per decade historically and is expected to rise slightly slower at 13.8 mm per decade in the future [39]. Summer precipitation patterns in the YKD have shifted over the past four decades, with a decline in large-scale precipitation due to weakening Bering Sea storms and an increase in convective precipitation driven by enhanced local hydroclimatology [40]. These trends, derived from both NOAA historical data and CMIP5 climate model predictions, underscore the dynamic climatology of the YKD and how these rapid climate changes might be affecting the region's permafrost and ecological conditions. The region's vegetation is predominantly tundra, composed of low shrubs, sedges, grasses, and mosses, typical of wetland environments in Arctic and subarctic regions [41]. Tundra vegetation in Bethel forms a dense, continuous ground cover. Uplands support sedge and shrub tundra with species such as dwarf arctic birch (Betula nana), Labrador tea (Ledum palustre), and cotton grass (Eriophorum vaginatum). Along drainageways and pond margins, wet sedge tundra dominates, featuring Carex aquatilis and sphagnum. Aquatic plants, including Carex rostrata and Eriophorum scheuchzeri, are found around thermokarst lakes. These ecosystems thrive in poorly to very poorly drained soils due to the presence of near-surface permafrost [42].

Permafrost in the region is poorly studied compared to other regions in the Arctic and sub-Arctic. The Brown et al. [43] circumarctic permafrost region map considered the YKD to be wholly within the continuous permafrost region; however, the updated Jorgenson et al. [44] permafrost map for Alaska remapped this region as part of the discontinuous permafrost zone. In addition, Pastick et al. [45] have modeled a relatively low probability of finding permafrost in the upper 1 m of the land surface across the YKD. The majority of permafrost-based field studies have occurred on the outer portion of the YKD in association with migratory bird studies and ecosystem processbased research [46][47][48]. Inland on the YKD, permafrost-based field studies have focused on disturbances associated with recent tundra fires [33,49] and carbon cycling [50,51]; however only scant information was presented on near-surface permafrost characteristics. Exploration for groundwater near Bethel in the mid-20th century found that thickness of permafrost extended down ~100 to 120 m below ground level [52]. Nearsurface permafrost studies in Bethel during the 1980s found that mean active layer thickness ranged from 38 cm in silty moist sedge tundra, to 64 cm in wet sedge tundra, to 103 cm in sandy soils, emphasizing the impact of organic layer thickness and soil texture on thaw depth [42].

| Methods

| Geotechnical Report Synthesis

Geotechnical and permafrost data were extracted from 18 historical geotechnical drilling reports from the Bethel area, covering the years 1967 to 2004. These reports, available as digital scans, were processed using Adobe Acrobat Pro 2020, specifically the Scan & OCR tool, to convert them into searchable text. Key data-including permafrost presence, depth to frozen ground, depth to the base of frozen ground, talik presence, visible ice, massive ice, and borehole depth-were manually extracted and compiled into a spreadsheet [53]. For data filtering, boreholes with depths of 3 m or less were excluded from further analysis, as distinguishing between seasonal frost and permafrost at these shallow depths was unreliable due to the timing of most drilling efforts. Boreholes deeper than this threshold were analyzed for permafrost and talik presence. Classification was guided by regional permafrost conditions, with depth to frozen ground serving as a primary criterion. Additional recoding ensured consistency across geotechnical reports. Lastly, boreholes were classified based on ice content, differentiating between visible and massive ice. Borehole locations were determined from maps within the reports, which were georeferenced to a basemap image in ArcGIS Pro 3.4. ArcGIS was also used to map and visualize the extracted data attributes (Figure 2). To ensure accuracy, extracted coordinates were cross-referenced with the original scanned reports.

| Mapping Trough Networks

Trough networks were mapped manually in a 60 km 2 area, covering the greater Bethel area, using a 20 cm resolution colorinfrared orthoimage from 18 August 2018 provided by Kodiak Mapping LLC. The networks were manually mapped in ArcGIS Pro 3.4 at a scale of 1:1000 as polyline vector features. Trough networks were aggregated into individual clusters of interconnected lines and later classified as either being an ice wedge network or a sand wedge (troughs resulting from thermal contraction cracks that fill with sand) network based on drilling during our fieldwork studies and association with stabilized sand dunes.

Ice wedge networks were also manually mapped across the YKD using a systematic approach that combined transect-based and targeted point sampling methods. Similar manual approaches have been used to map thermokarst landforms and to develop training datasets in the Arctic [54][55][56][57]. Transects were spaced 35 km apart in a north-south direction, with ice wedge point locations identified along each of the 13, east-west transect lines.

Mapping was conducted at a scale of 1:2000, with a field of view extending 210 m on either side of the transect lines, covering a total mapped area of ~2000 km 2 , or ~3% of the YKD. Ice wedge networks were mapped as point vector features, representing individual ice wedge networks or clusters of networks where they were abundant. The mapping utilized very high-resolution satellite imagery available as basemaps in ArcGIS Pro 3.4. The source imagery acquisition years for the study area range from 2012 to 2024.

The ice wedge mapping results were overlain on a 5 m IfSAR elevation dataset and the 2016 U.S. Geological Survey National Land Cover Dataset (NLCD; [58]) to create the ice wedge zone map for the YKD. We overlaid the point data with 5 m resolution IfSAR-derived elevation data to identify elevation thresholds associated with mapped ice wedge networks. Based on the distribution of confirmed occurrences, we excluded areas below 4 m and above 80 m in elevation, as these primarily fall outside the observed range for ice wedge formation in the YKD. Next, we incorporated the 2016 NLCD dataset to refine land cover suitability where ice wedges were mapped. Only land cover classes where ice wedges were confirmed to occur-primarily dwarf shrub and sedge/herbaceous tundra-were retained. Active or recently active riverine corridors were excluded based on their land cover classification and morphology. Finally, we removed geomorphologically unsuitable terrain by masking coastal volcanic areas and rounded mountain features using digitized units from the Coonrad [59] geologic map. These areas lack the surface conditions necessary for polygonal ground development. The resulting product is a refined ice wedge zone map representing regions where conditions are potentially most favorable for ice wedge development and preservation across the YKD.

| Field-Based Studies

Our permafrost studies in Bethel, Alaska, incorporated three primary field campaigns between 2023 and 2024. In May 2023, five boreholes were drilled using a 7.5-cm-diameter SIPRE corer in polygon centers and troughs, alongside 20 additional boreholes with a Kovacs auger across ice wedge troughs and rims to measure ice-wedge geometries. During this field campaign, we also installed a ground temperature array in a polygon center using iButtons (±0.5°C resolution) placed within a 1.5-in. diameter PVC tube at depths of 15 cm and 100 cm. The June 2024 campaign involved one SIPRE-corer borehole in an ice wedge polygon center and 14 Kovacs boreholes in ice wedge troughs, while in September-October 2024, we drilled two SIPRE-corer and nine Kovacs-auger boreholes in polygon centers and ice wedge troughs, respectively. The boreholes reached depths of up to 5.5 m with the SIPRE corer and up to 7.5 m with the Kovacs auger, enabling sampling across various depths for ground-ice content, stable water isotopes of wedge ice, and radiocarbon dating. Boreholes were strategically placed across four profiles in three primary study areas, which allowed for measurement of ice-wedge dimensions, including widths and vertical extents, and coring of ice wedge polygon centers. Frozen cores were described and photographed in the field. Cryostratigraphic descriptions were based on classifications of massive ground ice and cryostructures (patterns formed by ice inclusions in the frozen soil) adapted from Russian and North American literature [7,60,61].

Thaw depth measurements were conducted along two, 100 m long transects during early October 2024. Measurements were taken on October 1, 2024 at the BET-31/-32 transect and on October 3, 2024 at the BET-7 transect. A 1.5-m long graduated probe was used to measure thaw depth to the nearest centimeter at 1 m intervals along each 100 m long transect. For instances where thaw depth exceeded 1.5 m, measurements were recorded as greater than 1.5 m. At the BET-31/-32 transect, a deeper thaw feature existed along the transect so we used the Kovacs auger to drill through unfrozen sediments to the permafrost table, enabling the determination of thaw depth beyond 1.5 m in this case. Ground surface elevations come from a differential GPS survey at BET-7 and an uncrewed aerial vehicle (UAV) survey at BET-31/-32.

We utilized a Ground Explorer ground-penetrating radar (GPR; MALÅ, Brookvale, NSW, Australia) equipped with a 160 MHz shielded antenna to non-invasively measure thaw depths and identify subsurface ice wedges during early October 2024 fieldwork along the two 100 m transects. The shielded antenna was hand-towed along each transect using a common-offset configuration with a trace spacing of 0.5 m, triggered by an attached encoding wheel. The GPR system's built-in GPS recorded location data corresponding to each radar trace. Data processing in ReflexW (Sandmeier Geophysical Research, Karlsruhe, Germany) included distance corrections using GPS coordinates, airwave alignment, low-frequency noise filtering, and background removal. Manual thaw depth probe measurements collected on the same days as the GPR surveys were used to calibrate the radar velocity within the thawed active layer [62].

| Lab Analyses of Samples

To analyze the ground-ice content, a total of 44 soil samples were taken from the SIPRE and Kovacs boreholes. The samples were collected from the frozen part of the active layer (n = 4) and perennially frozen soils (n = 40) for gravimetric and volumetric moisture content and excess ground-ice content measurements; ice contents were determined according to Shur et al. [63]. These analyses provide insights into ground-ice distribution and overall ice content within various permafrost layers.

We radiocarbon dated organic material by individually picking samples from various depths to establish an age framework for organic deposits within the permafrost. Dates were sampled by targeting stratigraphic boundaries associated with the buried peat layers and underlying mineral deposits. Samples for AMS radiocarbon analysis were washed in distilled reverse osmosis water, photographed, and oven-dried (50°C) prior to shipment to the NOSAMS Woods Hole radiocarbon lab. The resulting dates were calibrated with OxCal Online Calib version 4.4 (https:// c14. arch. ox. ac. uk/ oxcal/ OxCal. html), using the IntCal20 calibration dataset [64]. This dating provided constraining ages for organic layer development, charcoal horizons, tephra deposition, ice wedge development, and ice-rich late-Pleistocene permafrost.

Stable oxygen and hydrogen isotope analyses (δ 18 O and δ 2 H) were conducted on ice wedge samples at the Alaska Stable Isotope Facility (ASIF) at the University of Alaska Fairbanks (UAF) [65]. Samples were analyzed using Continuous Flow Isotope Ratio Mass Spectrometry (CF-IRMS) with a Thermo Delta V + mass spectrometer, connected to appropriate peripherals for isotope measurements. Isotope ratios were measured using cavity ringdown spectroscopy (CRDS), and results are reported in standard delta notation (‰) relative to Vienna Standard Mean Ocean Water (VSMOW). In addition, the local meteoric water line (LMWL) for Bethel was calculated using monthly δ 2 H and δ 18 O precipitation estimates from the Online Isotopes in Precipitation Calculator (OIPC), based on unweighted monthly averages for the site coordinates (60.7921°N, 161.84°W, 4 m elevation). Analytical accuracy was maintained by analyzing internal standards calibrated to VSMOW-SLAP scales before and after each batch of samples. The precision of replicate measurements was typically ±0.1 ‰ for δ 18 O and ±0.5‰ for δ 2 H. This isotopic data helps to constrain the seasonal and climatic origins of the water that formed the ice wedges. In particular, they provide insight into the maritime influence on Holocene ice wedge formation in the Bethel region. The tephra layer was sampled from the BET-21 core (68-79 cm depth) and analyzed at the U.S. Geological Survey Alaska Tephra Laboratory in Anchorage, Alaska. The sample was wet sieved to remove grains < 63 μm. Backscatter electron images were acquired using a TESCAN MIRA-4 fieldemission SEM in low vacuum (~30 Pa) mode at 15 kV accelerating voltage, 10 nA beam current, and a 15 mm working distance. Combined, these lab-based assessments enabled us to generate a multi-faceted profile of the permafrost characteristics in Bethel, Alaska, contributing to a more comprehensive understanding of ground-ice distribution, permafrost formation processes, and ecosystem evolution in this region.

| Results

| Synthesis From Geotechnical Reports

A total of 803 geotechnical boreholes were analyzed from sites within Bethel, Alaska, with an average depth of 5.9 m, ranging from 1.2 to 36.0 m. Of these, 762 boreholes were successfully geolocated (Figure 2). Due to insufficient depth and the seasonal differences in drilling campaigns, 79 boreholes (less than 3 m deep) were excluded from further analysis, as they lacked adequate information to assess near-surface permafrost conditions. Permafrost was identified in 576 boreholes, while 107 boreholes indicated no near-surface permafrost but were classified as containing a talik, given the regional permafrost thickness of approximately 100 m. Among the boreholes with permafrost, 44% exhibited near-surface permafrost with no sub-aerial talik, 42% showed both near-surface permafrost and the presence of a subaerial talik, and 14% showed no near-surface permafrost and the presence of a sub-aerial talik. Notably, 160 boreholes contained visible ice within the permafrost, and one borehole revealed the presence of massive ice.

| Remotely Sensing Trough Networks-Bethel, AK

The trough network mapping results reveal a wide occurrence of ice wedges and sand wedges within the 60 km 2 study area (Figure 3). A total of 466 ice wedge networks were identified compared to 123 sand wedge networks. Ice wedges exhibit a mean length of 110 m, with the longest feature reaching 1182 m, summing to a total length of 51 km. In contrast, sand wedges have a mean length of 87 m, with the longest extending 332 m and a cumulative length of 11 km. Elevation analysis reveals differences in landscape position, with ice wedges typically occurring at lower elevations, averaging 18 m asl and ranging from 6 to 52 m asl, while sand wedges are found at higher elevations, with a mean of 32 m asl and a range of 9 to 50 m asl. These distinctions in distribution, size, and elevation underscore the geomorphic and environmental variability influencing past wedge network formation and presence near Bethel, Alaska. This mapping effort also provided the information that guided our field studies in 2023 and 2024 (Figures 4 and 5).

| Near-Surface Permafrost Characteristics

The ground temperature data collected from June 7, 2023 to June 2, 2024 in an ice wedge polygon center in Bethel reveals the presence of warm, near-surface permafrost (Figure 6). The mean annual ground temperature at 15 cm depth within the active layer was +0.7°C, while at 100 cm depth it was -0.8°C, indicating very warm near-surface permafrost. The MAAT during this period was +0.15°C.

Thaw depth measurements conducted along the 100 m transects at BET-7 and BET-31/-32 reveal distinct patterns (Figure 7). At the BET-7 transect, the mean thaw depth was 59.7 cm, with observed values ranging from 41 to 99 cm. In contrast, the BET-31 transect exhibited greater variability, with a mean thaw depth of 82.0 cm, a minimum of 42 cm, and a maximum of 276 cm. Locations with deep thaw along BET-31/-32 correspond to icewedge troughs (Figure 7a), which indicates recent or ongoing ice-wedge degradation, while ice wedges at BET-7 are generally stable. These results indicate significant spatial differences in thaw depth, reflecting variability in local environmental and soil conditions across the transects and are important for assessing the vulnerability of near-surface wedge ice. To compare our probe-based thaw depth measurements in early October 2024 that correspond to the permafrost table, with the GPR measurements, we used the calibrated radar velocity for the thawed layer (0.046 m/ns) to convert the two-way travel time to depth (Figure 8). The GPR data shows that there are twice as many ice wedges along the 100 m transects than are apparent in the remote sensing imagery (Figure 8). This indicates that our mapping of ice wedges with the remote sensing data in the Bethel area should be viewed as a conservative estimate.

The depth to the tops of wedges that were studied in the field (Figure 9) were assessed in the center of 14 troughs (n = 14) and along transects across four troughs (n = 15). The depth to the tops of ice wedges, based on these measurements, had a mean value of 90.4 cm, ranging from a minimum of 76.0 cm to a maximum of 166.0 cm. The standard deviation was 17.3 cm, indicating moderate variability in the depths across the surveyed points. For the four transects where ice wedge geometry was determined, the average ice wedge widths varied from 1.0 to 1.4 m (1.2 m average) and the average vertical extents varied from 2.2 to 2.4 m (2.3 m average).

The cryostratigraphy of BET-7 reveals a stratified permafrost profile with distinct variations in ground ice content (Table 1; Supplemental Material). The uppermost 15 cm consists of an unfrozen active layer (ALU) of moss and lichens. Below this, from 15 to 58 cm, a frozen active layer (ALF-indicating the active layer was still frozen at the time of field sampling) of peat transitions into a frozen active layer and potentially transient layer (TL) from 58 to 78 cm of very fine sandy silty soils with peat. Gravimetric moisture content (GMC) in this layer (at 58-68 cm) is 170.5%, volumetric moisture content (VMC) is 79.1%, and no excess-ice content (EIC) was detected. From 78 to 157 cm, a quasi-syngenetically frozen intermediate layer (IL) dominates, composed of very fine sandy silty soils with peat inclusions (especially at 92-108 cm); with GMC ranging from 185.3% to 202.0%, VMC from 80.5% to 83.8%, and EIC from 35.2% to 44.4%. The tephra layer was detected at 78-92 cm. From 157 to 216 cm, a syngentically frozen peat layer was observed with GMC of 937.5%, VMC of 94.0%, and EIC of 17.0%. The deepest layer, from 216 to 265 cm, consists of very fine to fine sandy silty soils with lower GMC (52.8%-57.7%) and VMC (59.9%-62.1%), and lower but still relatively high EIC (21.5%-27.1%). We identify this layer as quasi-syngenetic permafrost (QSP)-a buried IL that developed after peat accumulation due to a decrease in the active-layer thickness.

The cryostratigraphy of BET-21 reveals a complex profile with significant variability in ground ice content (Table 1; Figure 10; Supplemental Material). The uppermost 12 cm consists of an ALU of moss and peat, transitioning into an icepoor ALF with peat down to 40 cm and silty sand down to 57 cm. Below this, an ice-poor TL and relatively ice-rich IL extends to 79 cm, containing the tephra (from ~68 cm). From 79 to 214 cm, SP composed of peat dominates, with GMC ranging from 270.9% to 612.4%, VMC from 81.9% to 91.1%, and EIC from 1.8% to 27.8%. Deeper layers from 214 to 386 cm are presumably QSP, consisting of the ice-rich sandy silt and silt, with GMC of 125.9% to 267.1%, VMC of 76.3% to 87.2%, and EIC of 39.3% to 61.5%. In the Kovacs borehole BET-21d, the same ice-rich unit was encountered below the ice wedge from 285 to 520 cm (Figure 10; Table 1; Supplemental Material) with GMC of 53.4% to 837.6%, VMC of 60.2% to 94.9%, and EIC of 17.2% to 81.0%. From 530 to 740 cm, ice-and organic-poor sandy silt and silty sand unit (presumably refrozen talik) was encountered with GMC of 28.5% to 58.2%, VMC of 44.6% to 59.8%, and EIC of 0.0% to 2.9%.

The cryostratigraphy of BET-31 and BET-32 reveals complex, layered permafrost profiles with alternating organic-rich and mineral units (Table 1; Supplemental Material). At BET-31, the ALU extends from 0 to 65 cm, comprising moss, peat, and from 7 cm very fine silty sand, transitioning into a TL from 72 to 98 cm with peat and silty sand, and a thin IL from 98 to 102 cm. Below this, SP (peat) dominates from 102 to 456 cm. Below the peat layer (456 to 552 cm), very fine to fine sand (presumably QSP) was encountered. Similarly, at BET-32, the ALU spans 0-44 cm with moss, peat, and sandy silty peat from 15 cm, followed by a TL with very fine silty sand from 44 to 51 cm, and an IL from 51 to 108 cm (same soil), which includes the tephra from at least 99 to 108 cm. Below 108 cm, SP consists of peat, which at 296-305 cm is underlain by buried IL (QSP) that consists of very fine to fine silty sand. Ice contents for the BET-31/BET-32 area were determined from the Kovacs borehole BET-31e for the sandy silt and silty sand unit below the peat layer (from 380 to 750 cm), with GMC of 29.7% to 254.3%, VMC of 45.7% to 85.0%, and EIC of 0.0% to 25.0% with ice contents generally decreasing with depth.

Our estimations of excess-ice contents of the upper permafrost in the study area revealed rather high values, especially in the upper 5-6 m (Figure 11). Most of the 40 permafrost samples were obtained from organic-rich mineral soils, and the average excess-ice content was 27.3% ± 20.5% vol (n = 40). Ten samples (25% of all samples) had excess-ice content of more than 40%, while only eight samples had excess-ice content < 10%, and most of these were obtained from depths > 6 m. Based on this data, we expect that the degradation of the upper permafrost in our study area, could result in thaw settlement exceeding 1.5 m based solely on pore and segregated ice contents. Factoring in the presence of ice wedges and measurements of their geometry we can expect even more localized subsidence along the trough networks given the average width (1.2 m) and average vertical extent (2.3 m) measured by drilling in the field.

| Radiocarbon Ages of Permafrost-Preserved Organics

Samples for radiocarbon dating were collected at depths that represented cryostratigraphic and environmental transitions, including changes in ground-ice characteristics, charcoal layers, peat composition transitions, tephra accumulation events, and the development of basal peat deposits (Table 2). At BET-7, samples from 19 and 57 cm dated to 169 and 3128 calendar years before present (cal yrs BP), respectively, and were within two charcoal layers, which represent past tundra fire events. A sample at 160 cm (6725 cal yrs BP) marks a transition from buried woody peat to organic-rich mineral soil. The basal peat contact is identified at 216 cm, and dated to 8090 cal yrs BP. In BET-21, near-surface peat buried by tephra at 78 cm is dated to 4832 cal yrs BP, while the basal peat contact at 210 cm is 11,461 cal yrs BP. The deepest sample, at 370 cm, dates to 18,739 cal yrs BP, representing detrital plant fragments in icerich late-Pleistocene permafrost (likely reworked Yedoma) below the basal peat. In BET-31, the near-surface peat transitions into a mixture of organic-rich sands at approximately 100 cm depth, with an associated age of 4652 cal yrs BP. The lower peat sequence transitions from woody to herbaceous peat at 334 cm, dated to 6874 cal yrs BP. This is followed by the basal peat contact at 455 cm, which dates to 9443 cal yrs BP. At 544 cm, organic-rich silty sand below the basal peat is dated to 10,345 cal yrs BP. At BET-32, near-surface peat buried by tephra at 107 cm is dated to 3644 cal yrs BP, while the basal peat contact at 292 cm is dated to 8408 cal yrs BP. The deepest sample, at 303 cm, represents organic-rich silty sand below the basal peat, dated to 8860 cal yrs BP. These dates provide a sequential timeline of environmental and stratigraphic changes, including fire events, vegetation transitions, tephra deposition, and basal peat development, underscoring the dynamic history of the region.

| Stable Oxygen and Hydrogen Isotope Analyses of Wedge Ice

The δ 18 O and δ 2 H values of ice wedges from 16 boreholes exhibit a tight clustering reflecting consistent isotopic signatures and ground ice source waters across the study area (Table 3). The δ 18 O values average -15.7‰ with a standard deviation of 0.9‰, while δ 2 H values average -113.2‰ with a standard deviation of 4.6‰. These low standard deviations indicate minimal variation in isotopic composition (i.e., close to analytical precision), suggesting relatively homogeneous formation and source water conditions. This isotopic consistency highlights stable winter precipitation inputs and uniform freezing processes across the sampled locations.

| Initial Assessment of the Permafrost Preserved Tephra

The tephra deposit preserved at three sites in this study is appreciably thick, averaging 10-15 cm as measured at three sites. The pure tephra is fine grained with grains as large as 500 μm and composed of abundant white pumices and clear glass shards. The glass shards and pumices are predominantly clear with rare brown glass. The glass shards represent highly fragmental forms including thin bubble wall fragments (Figure 12). Together the thickness, grainsize and fragmental forms suggest it was formed during a significant eruption. Characteristics of this tephra layer including its age are consistent with the 3600 cal yrs BP eruption of Aniakchak volcano located on the Alaska Peninsula [66]. Aniakchak volcano has had at least two caldera-forming eruptions (CFE) in postglacial time, Aniakchak I (CFE I) and Aniachak II (CFE II) [67]. The modern 10 km diameter Aniakchak caldera was formed during a climactic explosive eruption ca. ~3600 cal yrs BP (CFE II) that produced widespread tephra deposits [67,68]. The ash cloud from this eruption was directed northward and fall deposits have been recognized in geologic outcrops along the Alaska Peninsula and mainland Alaska, as far as the northern Seward Peninsula [69-72], Siberia's Chukotka Peninsula [73], and in marine sediment in the Arctic Ocean [74].

| Remotely Sensing ice Wedge Networks Across the YKD

Mapping ice wedge networks across the YKD revealed a total of 761 ice wedge network clusters (Figure 13). Elevation analysis indicated a mean elevation of 18.4 m asl, with a median elevation of 13.5 m asl, and a minimum elevation of 3.0 m asl. Notably, 95% of the mapped ice wedges were found at elevations ranging between 4 and 80 m asl, and more than half of the ice wedges occur between 4 and 20 m asl. Landcover analysis, based on the 2016 USGS National Land Cover Dataset (NLCD), showed that 97% of the mapped ice wedges intersected four primary landcover classes: dwarf shrub, shrub/scrub, grassland herbaceous, and sedge herbaceous. We presume that these areas are underlain by similar peat deposits. By integrating these thresholds for elevation and landcover types, an ice wedge zone map was created using a minimum mapping unit of 10 ha. This zone, defined by thresholds from both datasets, encompasses 32% of the YKD region, providing a refined representation of the areas where ice wedges may occur.

| Discussion

| Permafrost and Ground ice Conditions on the YKD

Our permafrost studies in Bethel provide the first direct observations of ice wedge networks on the YKD, quantify groundice characteristics and excess ice contents, and confirm the presence of ice-rich Pleistocene permafrost, which was likely reworked prior to the onset of the Holocene, at depths of ~2 to > 5 m below the ground surface. Permafrost conditions on the YKD have long been a subject of debate and uncertainty. Michaelides et al. [35] used the YKD region as a case study region for evaluating a SAR-based remote sensing algorithm focused on mapping active layer variability because they deemed that thermokarst subsidence was minimal in the YKD. They based their assertion on the assumption that soils on the YKD contain little excess ground ice in the form of ice wedges or layers [35]. Additionally, Frost et al. [33] have interpreted that long-term permafrost thaw effects due to fire on the YKD TABLE 1 | Generalized permafrost borehole data for BET-7, -21, -31, and -32 showing the elevation (m asl) of each site, the depth of the borehole, inferred active layer depth (permafrost table), depth to the base of the peat layer, and ground-ice content information (GMC-gravimetric moisture content, % wt; VMC-volumetric moisture content, % vol; EIC-excess-ice content, % vol) for perennially-frozen peat and underlying organic-rich mineral soils (mainly very fine silty sand with peat inclusions). ND stands for not determined as these cores have been archived for future whole core analysis. For detailed borehole information, please see the Supplemental Material.

Borehole ID

Date

Elevation (m asl) Borehole depth (cm) Permafrost table (cm) Basal peat (cm) GMC (%) VMC (%) EIC (%) Peat Mineral Peat Mineral Peat Mineral BET-7 5/3/2023 15.4 265 58 216 937.5 52.8-202.0 94.0 59.9-83.8 17.0 21.5-44.4 BET-21 6/4/2024 14.5 386 53 214 270.9-612.4 125.9-267.7 81.9-91.1 76.3-87.2 1.8-27.8 28.9-61.5 BET-31 10/1/2024 19.9 552 72 456 ND ND ND ND ND ND BET-32 10/2/2024 18.3 305 51 296 ND ND ND ND ND ND 10991530, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ppp.70004 by Benjamin Jones -University Of Alaska Fairbanks , Wiley Online Library on [07/07/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License were muted since the region lacks extensive ice-wedge polygons reducing the vulnerability of YKD ecosystems to secondary impacts due to post-fire thermokarst. However, both Frost et al. [33] and Ludwig et al. [50] did note the presence of discontinuous, ice-rich permafrost in the region. In addition, Strauss et al. [31] included large portions of the YKD region in their Yedoma map of the northern hemisphere; however, they lacked field evidence for this delineation (Figure 1). Our findings contribute to a clearer and more accurate understanding of permafrost conditions on the YKD (Figure 13), highlighting the presence of Holocene ice wedge networks and possible relict Yedoma preserved in warm permafrost. These observations call into question the near-surface permafrost thaw susceptibility, raising concerns about how the region will respond to future climate and fire-related disturbances that may compromise the ecosystem protection of this permafrost [75].

Understanding the presence of ice wedges at the southern fringe of permafrost is essential due to the heightened vulnerability of these transitional areas to climate change [19,76]. Our preliminary ice wedge zone map indicates that roughly 32% of the YKD occurs in areas mapped as potentially having ice wedges. We found that these networks predominantly align with a ~4 m asl topographic feature, interpreted by Dupré [77] as a ~3000 cal yrs BP shoreline demarcating the line of late Holocene delta progradation by the Yukon and Kuskokwim River systems. While ice wedges were detected inland of this boundary up to 80 m asl, they were mostly absent seaward of it. This map does not imply that ice wedges are uniformly distributed within the ice wedge zone but rather that Holocene and possibly late-Pleistocene [31,77] ice wedges may be present. While only one or two deep boreholes were collected per transect due to logistical constraints, we supplemented these with multiple shallow auger holes and ground-penetrating radar surveys to capture local variability in ice wedge geometry and near-surface stratigraphy. The consistency in peat composition, wedge ice geometry, stable isotopic signatures, and radiocarbon age-depth relationships across sites supports the interpretation that these cores are broadly representative of the dominant polygon types in the study area. However, we acknowledge that additional cores across a wider range of polygon morphologies would further strengthen our ability to generalize cryostratigraphic trends across the YKD. Future work could target more diverse microtopographic settings to refine interpretations of lateral and vertical variability in ice wedge development and peat accumulation.

As the boundary between stable permafrost and thawing permafrost shifts, degrading ice wedges can serve as early indicators of thawing thresholds being crossed [21]. In our Bethel study area, ice wedge troughs were commonly underlain by wedges averaging 1.2 m in width and 2.3 m in vertical extent, while excess-ice content in organic-rich permafrost ranged from 0% to over 80%, with a mean of 27.3 ± 20.5% (n = 40). The combined presence of wedge ice and pore ice suggests that localized thaw settlement could exceed 1.5 m in trough and polygon center areas, particularly in low-lying terrain. These conditions are widespread across the ~32% of the YKD mapped as the ice wedge zone. The melting of ground-ice often triggers a cascade of environmental effects, such as ground subsidence, which disrupts hydrological systems and alters plant and animal habitats [18,22,78]. This is particularly important for Arctic and subarctic ecosystems, where even small shifts in soil stability can lead to the creation of thermokarst features, impacting local water drainage patterns and vegetation [17,79]. Changes in remotely sensed NDVI from 1982 to 2021 for the YKD reveals decadalscale variability rather than a consistent trend, with an initial increase until the late 1990s, a decline from 1998 to 2010, and a partial rebound post-2010 [34]. These fluctuations indicate that temperature alone does not govern tundra vegetation productivity as inferred from NDVI trends. Instead, ice wedge melting, thermokarst formation, and altered hydrology might play a critical but overlooked role in shaping NDVI patterns on the YKD. Thawing of permafrost can also have implications for the release of methane and carbon dioxide, which are important greenhouse gases that provide a positive feedback to global warming [80]. In addition, ice wedge degradation can have profound implications for infrastructure and communities. Roads, buildings, and utilities are vulnerable to ground instability as ice wedges melt, which results in differential thaw settlement and uneven terrain and structural damage [81]. As climate models predict continued warming, understanding how these processes unfold in vulnerable areas is crucial for developing long-term adaptation strategies for infrastructure and planning. By documenting the distribution and morphologies of ice wedges and the historical climate envelopes of these ice-rich permafrost areas that are close to the 0°C isotherm, we can better predict how these areas will change in the coming decades. This is especially useful for areas like the Bethel region, which is considering policy decisions and permafrost thaw mitigation efforts at local and global scales [37].

| Ice Wedge Temperature Dependence

Ice wedge formation is commonly referenced as being limited to MAATs colder than -4°C to -6°C [3,12]. However, examples TABLE 2

| Radiocarbon dates for organic material sampled from BET-7, -21, -31, and -32. We targeted transitions in lithology, tephra deposition, and charcoal layers in our sample preparation. Samples were calibrated using the IntCal20 calibration curve [64] and the calibrated age is reported as the 2-sigma midpoint in calendar years BP. Analyses were conducted at the National Ocean Sciences Accelerator Mass Spectrometry Facility (NOSAMS) at the Woods Hole Oceanographic Institution (WHOI). For detailed borehole information, please see the Supplemental Material.

Borehole ID

Depth, cm

Lab ID

Radiocarbon

age, yrs BP Fraction modern Calibrated age (2-sigma midpoint), calendar yrs BP Soil unit Dated material Sample notes BET-7 19 OS-180799 215 ± 15 0.9738 169 Peat Charcoal A past tundra fire BET-7 57 OS-180796 2970 ± 15 0.691 3128 Peat Charcoal A past tundra fire BET-7 160 OS-180797 5900 ± 25 0.4795 6725 Peat Wood Transition from woody to herbaceous peat BET-7 215 OS-180798 7260 ± 30 0.4048 8090 Peat Twig Basal peat contact BET-21 78 OS-180802 4230 ± 25 0.5908 4832 Peat Twig Near-surface peat buried by tephra BET-21 210 OS-180957 10,000 ± 35 0.2879 11,461 Peat Twig Basal peat contact BET-21 370 OS-181026 15,400 ± 65 0.1474 18,739 Organic Silt Detrital plant fragments Ice-rich syngenetic permafrost-yedomabelow basal peat BET-31 105 OS-183544 4140 ± 20 0.5973 4652 Peat Twig Near-surface peat buried by eolian material BET-31 334 OS-182616 6040 ± 25 0.4716 6874 Peat Wood Transition from woody to herbaceous peat BET-31 455 OS-182617 8390 ± 30 0.3519 9443 Peat Twig Basal peat contact BET-31 544 OS-182618 9210 ± 35 0.3178 10,345 Organic silty sand Twig Organic-rich silty sand below basal peat BET-32 107 OS-182619 3410 ± 25 0.6541 3644 Peat Twig Near-surface peat buried by tephra BET-32 292 OS-182620 7610 ± 30 0.3880 8408 Peat Twig Basal peat contact BET-32 303 OS-182621 7990 ± 30 0.3700 8860 Organic silty sand Twig Organic-rich silty sand below basal peat

in the literature indicate that ice wedges can be present and active in MAATs as warm as -2°C to -3°C in the presence of peat [13,23]. Ice wedges can be found in warm permafrost peat deposits due to unique thermal and physical characteristics of organic soils compared to mineral soils. The high volumetric water content of frozen peat, combined with its relatively low bulk density and high porosity, allows it to hold significantly more water than fine-grained mineral soils [13]. This high water content, along with a greater thermal contraction coefficient for ice, leads to higher thermal stresses in frozen, saturated organic soils during winter temperature drops, increasing the likelihood of thermal contraction cracking, essential for ice wedge formation (Al Moussawi, 1988;Andersland & Ladanyi, 2004). Furthermore, the insulating properties of thick snow cover in peatlands inhibit the freeze back of the active layer, allowing the peat to stay close to 0°C throughout winter [13]. This thermal stability makes the peat more susceptible to thermal contraction and the formation of ice wedges, as opposed to relatively warm mineral soils, which refreeze more completely and quickly [13]. Our findings add a new observation to the literature in Alaska with ice wedges persisting in a region with a long-term MAAT of -1.2°C and a more recent climate normal period where air temperatures averaged -0.3°C.

Figure 6a shows air temperature data at BET-7 from June 2023 to June 2024, the rapidly cooling winter air temperatures represented are potentially conducive for ice wedge development. Future work will focus on the potential activity (i.e., cracking, growth) of ice wedges in Bethel and on the YKD.

| A Hidden Paleo-Archive on the YKD

Our ice wedge network studies in Bethel have unveiled a previously unreported paleo-environmental archive in the scientific literature. We discovered the oldest basal peat sample ages in the YKD in the set of ice wedge polygon centers that we cored between 2023 and 2024. From oldest to youngest, our four basal peat ages are 11,461 cal yrs BP (BET-21, 210 cm depth), 9443 cal yrs BP (BET-31, 455 cm depth), 8408 cal yrs BP (BET-32, 292 cm depth), and 8090 cal yrs BP (BET-7, 215 cm depth) (Table 2; Figure 14). Previously, the oldest basal peat dates reported in the YKD region were < 2000 years old [82,83], which is much younger than the early Holocene peak in new peatland development across Alaska [82,84,85]. Basal peat, which represents the oldest layers of organic material at the base of a peat deposit, indicates a shift in hydrology and climate that lead to conditions that slow decomposition and the preservation of organic matter. Basal peat dates provide important clues about past vegetation trends, hydrology, and permafrost stability, helping to reconstruct the broader landscape evolution of the region [86]. Because peats accumulate stratigraphically with time, they are a key source for reconstructing past ecological and climatic conditions since peat initiation. Further, this new early Holocene peat initiation date (e.g., BET-21) reveals greater variability (spatial and temporal) in local climate and landscape evolution (e.g., peat formation conditions) within the regional paleo record than previously realized.

The discovery of these old basal peats offers insights into the potential role of the YKD as a long-term carbon sink if similar thick, buried Holocene permafrost-preserved peat deposits are more widely distributed in the YKD region than previously thought. Throughout the Holocene, these ice wedge polygon peat ecosystems could have sequestered vast amounts of carbon [87]. By analyzing the timing of peat initiation and accumulation rates, that have been corrected for the presence of ground-ice, across a broader region, we can begin to estimate carbon sequestration trends and assess the contribution of these warm permafrost peatlands to atmospheric carbon storage over millennia [88][89][90]. This is particularly relevant in the context of current climate change, as warming temperatures and permafrost thaw threaten to destabilize these ice-and carbon-rich environments, potentially releasing this stored YKD peatland carbon back into the atmosphere [80,85,91]. Thus, basal peat radiocarbon dating in Western Alaska provides a valuable framework for understanding both past, current, and future ecological processes in Arctic and sub-Arctic regions.

The discovery of the Aniakchak CFE II tephra in our permafrost cores from the YKD is also of importance (Figure 12). First, it caps and preserves the buried Holocene peat and ice wedges at depths ranging from 80 to 100 cm below the ground surface. The presence of the tephra coupled with our radiocarbon dates provides a compelling story for identifying likely periods of past ice wedge activity. The transition from a woodrich peat to a more herbaceous peat at BET-31 (~6800 cal yrs BP) and at BET-7 (~6700 cal yrs BP) likely indicates a shift from early Holocene warmth to the onset of neoglacial cooling [92,93], and the formation of Holocene ice wedges post

TABLE 3 | Bethel ice wedge stable water isotope (δ 18 O and δ 2 H) values. Analyses were conducted at the Alaska Stable Isotope Facility, Institute of Northern Engineering, University of Alaska Fairbanks. For detailed borehole information, please see the Supplemental Material.

Borehole ID Depth, cm δ 18 O (‰ V-SMOW) δ 2 H (‰ V-SMOW)

thermal maximum. The ice wedge and ice-rich peat at BET-31 are above this horizon. Given our radiocarbon dates from peat deposits in various landscape positions near Bethel, and the lack of ice wedges in the delta plain that formed after ~3000 cal yrs BP [77] on the YKD (Figure 12), we can tentatively bracket Holocene ice wedge formation to between ~3000 and 6500 cal yrs BP. This timing also fits well with a reported cool period with more ice in the Bering Sea from 3200 to 5500 cal yrs BP [84].

Most previous near-surface permafrost studies and observations on the inner portions of the YKD have been restricted to the upper 1 m of the land surface [42,94]. This could be a result of encountering the 10 to 15 cm thick tephra deposit without knowing it. Our ability to core through this layer has revealed a potentially useful Late Pleistocene and early Holocene sedimentary archive which includes the presence of ice-rich, Late Pleistocene permafrost and the buried early Holocene peat that is not represented in global databases [95,96] (Figure 14). The presence of Pleistocene-aged, ice-rich permafrost preserved below the old basal peats on the YK Delta is also important for understanding the past, present, and potential impact of its degradation on the carbon cycle and past and future landscape evolution. The YK Delta region likely reflects a highly degraded Yedoma landscape with numerous thermokarst lakes and drained lake basins indicating that much of this carbonrich Yedoma had previously degraded during the Holocene. Peat accumulation at BET-21 (lowest landscape position) began ~11,500 cal yrs BP, probably preserving the underlying reworked Yedoma. If Yedoma was present at BET-7, -31, and -32 it would have completely degraded during the early Holocene warm period, delaying the onset of peat accumulation to between ~8000 to 9500 cal yrs BP. Our findings show that the cryptic Holocene wedges formed in peat accumulated above this legacy Yedoma landscape. Future studies of YKD permafrost can potentially elucidate the amount of carbon that had been previously mobilized during the Holocene Thermal Maximum (HTM) [92], particularly as the climate shifted from a continental to a more maritime climate as sea level rose post-Last Glacial Maximum (LGM). Recent studies have documented an increase in tundra fire activity on the YKD, correlating it with climate warming and vegetation changes. The YKD has experienced at least 250 fires since circa 1940 [97] and the region has a fire return interval of ~150 years over the last 1000 years [98]. The Tungak Lake (also on the YKD) fire-history record, the longest in Alaska, reveals minimal fire activity over the past ~35,000 years, with slightly increased burning between 25.5 and 14.0 kcal BP due to drier conditions from lower sea levels, followed by a decline in fires during the Holocene (~11.7 kcal BP onward), with only one fire recorded, as moisture increased despite greater tundra biomass [99]. Besides these two studies, very little is known about the paleo tundra fire regime on the YKD. Our findings at BET-7 show charcoal layers preserved in peat that date to 169 and 3128 cal yrs BP (Figure 14; Supplemental Material). The presence of charcoal in the peat reveals the location of pre-historic tundra fires on the YKD for the first time and provides a point fire return interval measurement of ~1500 years over the course of the late Holocene. The interplay between climate change, vegetation shifts, and increased fire activity in recent years on the YKD underscores the need for ongoing research to better understand longer-term fire regimes and permafrost interactions.

| Stable Isotope Ice Wedge Outliers in Alaska

The typical water isotopic values of ice wedges are generally measured using the ratios of stable oxygen (δ 18 O) and hydrogen (δ 2 H) isotopes. Ice wedges tend to have highly depleted values of δ 18 O and δ 2 H, reflecting cold climatic conditions of permafrost regions, and winter precipitation sources. However, stable water isotopes show that ice wedges in the Bethel area are enriched relative to other ice wedges in Alaska (Figure 15). The mean δ 18 O and δ 2 H values for the 16 ice wedges measured at depths from ~1 to 3 m below the ground surface are -15.7‰ and -113.1‰, respectively. Water isotopes in ice wedges can vary depending on the region and local climate and continentality of a site, but typical δ 18 O ice wedge values reported in Alaska, Canada, and northern Siberia range from -25‰ to -35‰, while δ 2 H values range from around -200‰ to -280‰ [100][101][102][103][104][105][106].

The more enriched ice wedge isotope values in Bethel likely indicate a more proximal precipitation source during their formation. Our results align more closely with the Holocene ice wedge data from Lorino, Chukotka, Russia than they do with other sites in Alaska. In Lorino, the average δ 18 O values of -18‰ to -15‰ were obtained from ice wedges that formed between 6700 to 7700 cal yrs BP (Vasil'chuk et al., 2024). The overlap in isotopic values between Lorino and Bethel likely reflect similar regional climatic trends that are characterized by maritime conditions with severe winters that promoted ice wedge formation during the Holocene. We hypothesize that the climate conditions favorable for ice wedge development in Bethel (3600 to 6000 cal yrs BP) likely occurred several thousand years after those observed in Lorino (6700 to 7700 cal yrs BP). These findings suggest that Bethel experienced comparable climatic influences, likely driven by shared atmospheric and oceanic circulation patterns, but possibly at different periods in the Holocene, driven by differences occurring on the east and west side of the Bering Sea. We interpret the Bethel ice wedges-dated between ~3600 and 6000 cal yrs BP-as predominantly syngenetic or quasisyngenetic in origin, based on their occurrence within thick, accumulating peat sequences and their stratigraphic position beneath the Aniakchak CFE II tephra. This suggests that wedge formation occurred concurrently with sediment accumulation, supporting the reliability of radiocarbon ages from surrounding organic material. The presence of ice wedges in Bethel and across the Yukon-Kuskokwim Delta provides valuable context for interpreting relict polygonal networks elsewhere in Alaska. These results expand the climatic envelope for ice wedge formation by documenting a relatively warm, low-elevation end member for Holocene permafrost environments [1,12,107,108].

| Conclusion

This study underscores the potential significance of ice wedge networks in Bethel, Alaska, and the broader Yukon-Kuskokwim Delta (YKD) as key indicators of permafrost dynamics and valuable paleoclimate archives. Our findings provide direct evidence of ice-rich permafrost in a region where near-surface permafrost characteristics have been historically understudied. Through detailed mapping, field investigations, and laboratory analyses, we sampled and described ice wedges in the Bethel area for the first time. Additionally, GPR transects revealed ice wedges beyond those detectable in remote sensing imagery, suggesting that their abundance may be underestimated. The isotopic composition of wedge ice, combined with radiocarbon dating, offers new insights into the region's paleoecological and paleoclimatic history. Furthermore, our research documents buried Holocene peat, the persistence of late-Pleistocene ice-rich permafrost, charcoal evidence of past tundra fires, and the presence of the Aniakchak CFE II tephra, shedding light on the cryptic nature of Holocene ice wedge networks in the region. Combined, our findings highlight the complex interplay between climate, vegetation, and permafrost evolution on the YKD, reinforcing the importance of ice wedge networks as both climate proxies and indicators of permafrost resilience in a warming Arctic.

The implications of our findings are several-fold. First, the presence of ice wedge networks at the southern fringe of permafrost highlights their resilience to Holocene conditions. Second, the potential presence of ice wedges should be considered when studying the tundra landscapes of the YKD region. Third, the identification of buried, Holocene-aged peat layers and the reworked Yedoma permafrost beneath ice wedge polygons reveals the YKD as a cryptic but possibly significant carbon-stock reservoir. These deposits potentially offer a valuable paleoenvironmental archive that helps reconstruct past climate dynamics and informs future projections of carbon release under warming scenarios in the discontinuous permafrost region. This research highlights the critical need for continued observation and exploration of ice wedge systems in transitional permafrost regions, as they serve as key benchmarks for permafrost in a rapidly changing Arctic and Subarctic environments.