Canopy cover and microtopography control precipitation-enhanced thaw of ecosystem-protected permafrost

Northern high-latitudes are projected to get warmer and wetter, which will affect rates of permafrost thaw and mechanisms by which thaw occurs. To better understand the impact of rain, as well as other factors such as snow depth, canopy cover, and microtopography, we instrumented a degrading permafrost plateau in south-central Alaska with high-resolution soil temperature sensors. The site contains ecosystem-protected permafrost, which persists in unfavorable climates due to favorable ecologic conditions. Our study (2020–2022) captured three of the snowiest years and three of the four wettest years since the site was first studied in 2015. Average thaw rates along an across-site transect increased nine-fold from 6 ± 5 cm yr−1 (2015–2020) to 56 ± 12 cm yr−1 (2020–2022). This thaw was not uniform. Hummock locations, residing on topographic high points with relatively dense canopy, experienced only 8 ± 9 cm yr−1 of thaw, on average. Hollows, topographic low points with low canopy cover, and transition locations, which had canopy cover and elevation between hummocks and hollows, thawed 44 ± 6 cm yr−1 and 39 ± 13 cm yr−1, respectively. Mechanisms of thaw differed between these locations. Hollows had high warm-season soil moisture, which increased thermal conductivity, and deep cold-season snow coverage, which insulated soil. Transition locations thawed primarily due to thermal energy transported through subsurface taliks during individual rain events. Most increases in depth to permafrost occurred below the ∼45 cm thickness seasonally frozen layer, and therefore, expanded existing site taliks. Results highlight the importance of canopy cover and microtopography in controlling soil thermal inputs, the ability of subsurface runoff from individual rain events to trigger warming and thaw, and the acceleration of thaw caused by consecutive wet and snowy years. As northern high-latitudes become warmer and wetter, and weather events become more extreme, the importance of these controls on soil warming and thaw is likely to increase.


Introduction
Permafrost covers approximately 24% of the terrestrial landscape of the Northern Hemisphere (Brown et al 1997).Atmospheric temperatures in the Northern Hemisphere have warmed faster than the global average, and this amplified warming is expected to strengthen in the future (Rantanen et al 2022).Trends in permafrost temperature are consistent with trends in atmospheric temperature, with widespread warming observed since 1980 (Smith et al 2022).However, other environmental factors, such as vegetation, soil properties, hydrology, topography, local climate, and snow, also regulate permafrost temperature and thaw progression (Engstrom et  Rainfall is one such environmental factor that contributes to permafrost dynamics.Rainfall can accelerate permafrost thaw through advective heat transport (Iijima et al 2010, Connon et al 2018, Devoie et al 2019, 2021, Douglas et al 2020, Mekonnen et al 2021, Magnússon et al 2022) and by increasing soil thermal conductivity (Farouki 1981, Hinkel and Nelson 2003, Subin et al 2013), but rain can also slow thaw since wetter soils require more energy to warm due to greater heat capacity (Farouki 1981, Hinkel et al 2001, Subin et al 2013, Clayton et al 2021) and experience increased evaporative cooling (Hamm and Frampton 2021).Northern high-latitude precipitation is expected to increase (Bintanja andSelten 2014, Bintanja et al 2020), which may enhance the role of rain in permafrost progression.
We hypothesize that ecosystem-protected permafrost, the warmest and most vulnerable of the permafrost classes, will be most sensitive to these changes.Ecosystem-protected permafrost is located in regions with a mean annual air temperature (MAAT) above 0 • C where permafrost persists in unfavorable climates due to favorable ecologic conditions (Shur and Jorgenson 2007).Understanding how ecosystemprotected permafrost degrades is crucial because ecosystem-protected permafrost acts as a forecaster and predictor of future thaw in more northern latitudes (Beilman et al 2001) which are becoming more climatically similar to existing ecosystem-protected areas (Cohen et al 2014, Kusunoki et al 2015, Liu et al 2021, Wang et al 2021, Rantanen et al 2022, Thackeray et al 2022).Degradational features such as taliks are common in these warm permafrost systems (Harris 1988, Connon et al 2018, Devoie et al 2019).Taliks are layers of year-round unfrozen ground that have been shown to enhance permafrost thaw.One potential mechanism is the enhancement of advective heat transport with rain by providing subsurface lateral flow paths for runoff (Connon et al 2018, Devoie et al 2019, 2021).
To better understand the impact of rain, as well as other environmental factors affecting permafrost in the ecosystem-protected zone, we instrumented a rapidly degrading permafrost plateau in southcentral Alaska with high-resolution soil temperature sensors and measured environmental variables including canopy cover, snow, and microtopography.Our study took place from 2019 to 2022 and contained the wettest and snowiest years on record since 2015.

Site description
The site is located in the western Kenai Peninsula lowlands (60.47 • N, −150.73 • W) (figures 1(A) and (B)) within Kenai National Wildlife Refuge (permit: 2018-Res-RBNeumann-6276).Fieldwork took place May 2019 to September 2022.The plateau is surrounded by and contains wetland features (figure 1(C)) and experiences a semi-continental climate due to Cook Inlet to the northwest, Kachemak Bay to the south, and its location within the rain shadow of the Kenai mountains to the east.Plateau and wetland soil is peat-dominated to 3 m depth (Jones et al 2016).In 2015, the average seasonally frozen ground thickness was 45 cm, which overlayed continuous taliks (Jones et al 2016).Average depth to permafrost measured in 2015 was 1.37 m with an average permafrost thickness of 5.61 m (Jones et al 2016).Black spruce trees (Picea mariana) grow on the most elevated sections (hummocks) of the permafrost plateau (figure 1(D)).Plateau lowland areas (hollows) contain no live trees.

Soil temperature measurements
Soil temperatures were measured using 160 cm long Distributed Temperature Profilers (DTPs) with sensors every 10 cm (Dafflon et al 2022, Wielandt et al 2023) installed in August 2020.DTPs have an accuracy of ±0.05 • C and a 15 min measurement interval.Data were depth-corrected annually for frost heave and subsidence and linearly interpolated for locationto-location consistency.

Geospatial measurements and relative-elevations
A digital elevation model was created in MATLAB (R2022b) using 1620 real-time kinematic points recorded using Emlid Reach RS2 units in November 2020 (figure 1(E)).Relative-elevations were determined by comparing DTP elevations to the average elevation of the surrounding 5 m radius area.Relative-elevations were used to define locations as a hummock, transition, or hollow (figures 1(C) and (E)).Hummock locations were >0.2 m above, transition locations were between 0.2 m above and 0.1 m below, and hollow locations were >0.1 m below their surroundings (figure 1(E)).

Meteorological measurements and data
Precipitation was measured with an annually calibrated HOBO Rain Gauge installed April 2019 in a 5 m clearing.Air temperature was recorded with a HOBO Temperature/RH Logger installed April 2019 using dense black spruce as a radiation shield (Lundquist and Huggett 2008).Data were gap-filled using two National Weather Service stations, COCOAKKP7 and COOPKEYA2, located within 6 km (NOAA 2023).Air temperature and precipitation data from 2015 to 2019 were sourced from interagency weather station SDFA2 located 16 km west.Trends in SDFA2 data were consistent with in situ meteorological measurements recorded during the study period.Snow depth was measured with an avalanche probe in March 2022.Snow water equivalent data were sourced from Moose Pen SNOTEL site 32 km northeast (USDA 2023).Cumulative shortwave radiation was calculated using HPEval, a MATLAB-based software tool that uses high-resolution hemispherical photos to estimate shortwave radiation below forest canopy (Jonas et al 2020).

Hydrologic measurements
In-Situ pressure transducers and barometer were deployed in fully-screened PVC wells annually from May to September (2019)(2020)(2021)(2022).Sensor data were verified with manual groundwater measurements taken 2-3 times annually between May and September with a Solinst Water Level Meter.Soil moisture was measured continuously at one hummock (90 cm length sensor) and one transition (120 cm length sensor) location (figure 1(C)) using GroPoint Multi-Depth Soil Moisture Sensors attached to GroPoint Data Loggers installed in August 2020.Profilers measured a composite value of percent water every 15 cm along the length of the sensor at a 30 min interval.Outputs were calibrated for peat soils (Yoshikawa et al 2004) and a manufacturer-provided dielectric correction (1.019) was applied for ⩽0 • C soils.

Thermal conductivity modeling
Soil thermal conductivity was estimated using a physically-based model (Zhao et al 2019) that used degree of saturation, dry thermal conductivity, and wet thermal conductivity as inputs.Degree of saturation was calculated by normalizing measured water content by the highest observed water content, 96.7%.Dry and wet thermal conductivities were estimated through weighted averages of the thermal conductivities of air (0.025 W m -1 K -1 ), water (0.579 W m −1 K -1 at 10 • C), and organic matter (∼0.074W m −1 K -1 ).Dry thermal conductivity (0.027 W m −1 K -1 ) assumed 96.7% air and 3.3% organic matter.Wet thermal conductivity (0.562 W m −1 K -1 ) assumed 96.7% water and 3.3% organic matter.For ⩽0 • C soils, we calculated thermal conductivities assuming no ice content and maximum possible ice content.Maximum ice content was estimated as the difference between the current and highest observed water content of the measurement location.The thermal conductivity of 0

Canopy measurements
Canopy cover percentage was determined through image analysis of upward-facing smartphone photos taken with a Moment 14 mm fisheye lens 1 m above the soil surface in July 2022.Photos were analyzed in ImageJ by cropping the photo into a circle, splitting the photo into color panes, thresholding the blue pane into sky and vegetation pixels, and calculating percent cover from the relative proportion of pixels (Beckschäfer 2015).

Depth to frost table measurements
Permafrost table depth was estimated annually through manual probing.Probe length was 2.2 m from 2019-2021 and extended to 3.1 m in 2022.Measurements took place annually in mid-September.DTP location probing included four replicates within a 30 cm radius.Across-site transect (figure 1(C)) probing included three replicates, mimicking the methods and locations of Jones et al (2016).

Statistical methods
Environmental variables and thermal response metrics were compared for hummock, transition, and hollow locations using one-way ANOVA followed by Tukey HSD post-hoc using α = 0.05.Data were log 10 transformed to increase normality and reduce skew before ANOVA analysis.Data presented in figures are shown in non-transformed space for visual clarity.Correlations were determined by fitting data with a linear regression model and using two-sided hypothesis testing with α = 0.05.

Site meteorological conditions
From 2015 to 2022, MAAT was +3.3 • C with a mean annual precipitation of 480 mm (NOAA 2023).From 2020 to 2022, the site experienced the three snowiest years (figure 2(A)) and the second, third, and fourth wettest years (figure 2(B)) since 2015.Average air temperature did not vary significantly (SI figure S1).Three large warm-season rain events, defined as events with >50 mm of rain over two weeks, took place during our study in July 2020 (58 mm), July 2021 (54 mm), and July 2022 (77 mm) (figure 2(Bi,ii).The July 2020 event took place before the installation of site DTPs and soil moisture profilers and was therefore not included in further analyses.

Permafrost and subsurface temperature regimes within microtopographic groups
Hummock locations had significantly shallower depths to the permafrost table (figure 5(A), p < 0.001) and significantly slower rates of annual permafrost thaw (figure 5(B), p < 0.05) than transition and hollow locations.In 2022, average depth to permafrost at hummock locations was 97 ± 39 cm compared to 252 ± 53 cm and 288 ± 15 cm at transition and hollow locations, respectively.Average annual thaw rates between 2020 and 2022 at hummock locations were 8 ± 9 cm yr −1 compared to 39 ± 13 cm yr −1 and 44 ± 6 cm yr −1 at transition and hollow locations, respectively.
At a depth of 10 cm, cold-season soil temperatures of hummock locations were significantly colder (−3.5 ± 1.7 • C) than hollow locations (−0.6 ± 0.4 • C), while transition locations (−1.7 ± 1.1 • C) were statistically indistinguishable from either group (figure 5(C), p < 0.05).A comparison depth of 10 cm was chosen to maximize sensitivity to surface-atmosphere interactions while avoiding diel temperature oscillations.In hummock and transition locations, freezing temperatures (−0.05 • C) reached significantly deeper depths of 28 ± 11 cm and 20 ± 11 cm, respectively, than in hollow locations where −0.05 • C reached a depth of 3 ± 6 cm (figure 5(E), p < 0.001).A threshold of −0.05 • C was used to track extent of freezing since DTPs have a confidence interval of ±0.05 • C.
Hummock locations were significantly cooler (5.2 ± 1.7 • C) than hollows (9.1 ± 0.7 • C) during the warm-season at 10 cm depth, while temperatures at transition locations (7.9 ± 0.4 • C) were statistically indistinguishable from either group (figure 5(D),  p < 0.05).In hummock locations, warm temperatures (+3 • C) only reached a depth of 20 ± 17 cm, while in transition and hollow locations, warm temperatures reached depths of 69 ± 33 cm and 101 ± 9 cm, respectively.The depth of warming in hummock locations was significantly shallower than that of transition and hollow locations (figure 5(F), p < 0.05).A threshold of +3 • C was the warmest temperature that did not penetrate beyond the 160 cm sensor length of any location.

Thermal and hydrologic responses to large rain events
Large rain events (>50 mm) occurred in July 2021 and 2022 (figure 2(B)).During these events, a greater thickness of the soil column at transition locations warmed from ⩽0 • C to above-freezing temperatures than at hummock or hollow locations in both 2021 (figure 6(A), p < 0.01) and 2022 (figure 6(B), p < 0.01).During the 2021 rain event, an average of 107 ± 20 cm of soil within transition locations warmed from ⩽0 • C to >0 • C compared to 31 ± 42 cm and 6 ± 4 cm of soil in hollow and hummock locations, respectively (figure 6(A)).A similar thickness of transition soil warmed from ⩽0 • C to >0 • C during the 2022 rain event, 106 ± 14 cm, compared to <1 cm and 10 ± 1 cm of soil in hollow and hummock locations, respectively (figure 6(B)).The observed thickness of transition-location warming and thaw that occurred during both rain events (106-107 cm) included over 40% of total abovepermafrost soil during both years (252 ± 53 cm) (figures 5(A) and 6).Smaller warm-season rain events did not cause notable warming.Transition location #2 had a greater proportion of surrounding higher elevation area than transition location #1.In the context of hydrologic flow, the proportion of surrounding higher elevation area can be represented as the ratio of area contributing recharge to lower elevations (Area C ,) to area receiving recharge (Area R ).Transition location #1 had 2 m 2 of contributing area per 1 m 2 of receiving area leading to an Area C :Area R ratio of 2:1, while transition location #2 had 6 m 2 of contributing area per 1 m 2 of receiving area leading to an Area C :Area R ratio of 6:1.Observed responses of all example locations were similar in 2022 (SI figure S3).Soil temperature figures, thermal data, and environmental data for all site locations are available in SI sections S5-S19.
Before the rain events, the hummock location had an estimated thermal conductivity of 0.265 ± 0.044 W m −1 K -1 in 2021 and 0.244 ± 0.053 W m −1 K -1 in 2022 (table 1).Estimated pre-event thermal conductivity for the transition location ranged from 0.263 to 0.627 W m −1 K -1 and 0.287-0.556W m −1 K -1 in 2021 and 2022, respectively (table 1).This range was due to uncertainty associated with ice content (see methods 2.6).Soil moisture was not logged at hollows, so thermal conductivity could not be calculated.Due to presumably high soil moisture, hollow locations likely had greater near-surface thermal conductivity than transition or hummock locations.Hummock thermal conductivity remained unchanged during the 2021 rain event and increased to 0.297 ± 0.069 W m −1 K -1 (22% change) during the 2022 event (table 1).Thermal conductivity of the transition location converged to within 1% of the hummock location after the 2021 event (0.267 ± 0.071 W m −1 K -1 ) and within 5% after the event in 2022 (0.307 ± 0.061 W m −1 K -1 ).Possible changes in transition-location thermal conductivity due to the rain events ranged from a decrease of 45%-57% (high pre-event ice content) to an increase of 1%-7% (ice-free).Table 1.Average depth to water table (WT), near-surface water content (θ), and near-surface thermal conductivity (λ) in response to July 2021 and 2022 rain events.Change in depth to water table is averaged for two hollow locations (see figure 8

Accelerated thaw due to consecutive rainy and snowy years
As noted in the results, our study captured notably wet and snowy years (figure 2) accompanied by a nine-fold increase in thaw rates along the across-site transect (figure 3) which expanded existing taliks.These taliks begin isolated, progress to connected taliks, and mature into flowthrough taliks (Devoie et al 2021).Flow-through taliks that are hydrologically connected to wetlands often experience enhanced thaw that is nearly irreversible (Sjöberg et al 2016, Devoie et al 2019).This process could partially explain rapid thaw observed at the wetland-plateau interface (figure 3).Hummock locations thawed slowly (figure 5(B)) during the study period while all other locations, including the across-site transect which contained no hummock locations, thawed rapidly (figures 3 and 5(B)).This outcome emphasizes the stability of hummock locations and susceptibility of transition and hollow locations to rain and snow (figure 5  Increased thermal conductivity due to increased soil moisture does not appear to be the dominant mechanism.Soil-atmospheric thermal interactions are controlled by near-surface thermal conductivity (El-Din 1999).Transition-location thermal conductivity likely decreased during the rain events.Depending on pre-event ice content, the instrumented transition location experienced between a ∼57% decrease and ∼7% increase in near-surface thermal conductivity (table 1).Regardless, there was no evidence of increased top-down warming or thaw during the rain events (figures 7(C) and 8(B)).During the 2021 event, soils at 20 to 30 cm depth of transition location #1 remained ⩽0 • C while soils from 40 to 140 cm warmed (figures 7(C) and 8(B)).Thus, our data indicates a mechanism that enhances deep soil warming, as opposed to enhancing surface soil-atmospheric energy transport, is responsible for the observed warming and thaw at transition locations.

Transition locations thawed rapidly due to advective heat transport by rain
The pattern of warming and thaw indicates advective heat transport as the primary mechanism.Rain temperature, to a first approximation, is the same as air temperature (Byers et al 1949).When rain temperature is greater than soil temperature, rain events can add thermal energy to soil through thermal conduction, rapidly warming soils (Neumann et al 2019).During the large July rain events, air temperature was substantially warmer than soil temperature at hummock and transition locations, and slightly warmer than hollowlocation soil (figure 7).Transition locations presumably received runoff from surrounding higherelevation hummocks and released water into lowerelevation hollows.Subsurface runoff increases energy available for warming and thaw, explaining why transition locations warmed during rain events, but hummocks, which did not receive runoff, did not experience enhanced warming.Once flow-limiting soil ice melted within transition locations, subsurface runoff flowed into taliks, warming deeper soil layers.Hollow locations received runoff, but had warm soils before rain events, which minimized the soil-runoff thermal gradient.
Energy contributed by advective heat transport depends on runoff water temperature and volume, while energy needed to warm and thaw soils depends on pre-event soil temperature, water content, and ice content.A key relationship that connects energy contributed to energy needed in this context is the ratio between the land area contributing runoff (Area C ) to the land area receiving runoff (Area R ).The contributing area influences the volume of runoff generated while the receiving area influences the volume of water and ice that must be warmed and thawed.Transition location #1, which had an Area R :Area C ratio of 1:2, had a smaller magnitude thermal response to rain (figure 7(C)) than transition location #2 (figure 7(D)), which had a ratio of 1:6.The larger Area R :Area C ratio of transition location #2 presumably corresponded with relatively greater runoff volume, and therefore, more energy for warming and thaw.
Figure 9 explores the ability of runoff from rain to cause the observed transition location warming and thaw.During these events, ∼100 cm of soil with 86% water content (figure 8(B)) and an unknown amount of ground ice warmed from from 0 • C to 1 • C (SI section 1, figures 7(C) and (D)).Since runoff temperature and ice content were unknown quantities, ranges of these values were used as figure axes (figure 9).Minimum and maximum possible runoff temperatures range from the minimum observed event soil temperature, 0 • C, to the average air temperature during the event, 15 • C. At the site, Area R :Area C varied from 1:2 to 1:6, depending on Figure 9. Energy analysis of rain-induced warming and thaw.Dashed black lines mark scenarios where the energy stored in runoff is equal to the energy needed to warm and thaw 100 cm of transition-location soil for ratios of soil area receiving runoff (AreaR) to soil area contributing runoff (AreaC) of 1:2, 1:4, and 1:6.Below each line, the observed amount of warming is not possible with energy solely from runoff.Above each line, there is a surplus of energy.In the orange region, all scenarios have excess runoff energy.In the green section, all scenarios had insufficient energy to cause the oberved warming and thaw.In the blue region, some scenarios have surplus energy while others have insufficient energy based on AreaR:AreaC ratio.The x-axis is the percent ice contained in the layer of transition-location soil.The y-axis is the temperature of runoff as it reaches the transition location.This analysis assumes a homogeneous soil profile.The full analysis is shown in SI Section 1. the transition location.Three different Area R :Area C scenarios are presented in the figure, 1:2, 1:4, and 1:6 (figure 9).Along each scenario line, energy in subsurface runoff equals energy required for the observed warming and thaw.Above each line, there is a surplus of runoff energy.Below each line, there is insufficient energy to cause the observed warming and thaw.In the orange figure section, all Area R :Area C scenarios had excess energy.In the blue section, some scenarios had excess energy while others had insufficient energy.In the green section, all scenarios had insufficient energy.See SI section 1 for a full calculation description.
Figure 9 calculations confirm that energy in raininduced subsurface runoff was sufficient to cause the observed warming and thaw if average ground ice content was low.Using the average Area R :Area C ratio of 1:4 and an estimated runoff temperature of 7.5 • C, our analysis suggests sufficient energy to warm and thaw 100 cm of soil containing an average of 1.5% ground ice from 0 • C to 1 • C. Runoff temperature was estimated as 7.5 ± 5 • C by mass averaging known masses of 1 • C hummock soil water and 15 • C rainwater.Since all above-permafrost soil except the ∼45 cm seasonally frozen layer is ice-free talik (Jones et al 2016), ice content is not homogeneously distributed, but concentrated in the top 45 cm of the soil profile.With ice constrained to the seasonally frozen layer, energy from the same 1:4 Area R :Area C ratio and 7.5 • C runoff was sufficient to warm and thaw 100 cm of soil containing 3.5% of average ground ice in the seasonally frozen layer.The highest observed Area R :Area C ratio of 1:6, coupled with a 7.5 • C runoff temperature, had sufficient energy to warm and thaw 100 cm of soil with up to 6% of average ground ice in the seasonally frozen layer.Since rain-enhanced thaw depends on runoff, small rain events are not expected to cause notable warming and thaw.To create largescale advective heat transport, rain events must be sufficiently large to offset losses from evapotranspiration and add sufficient moisture to reach the point of free drainage.
Multiple field and modeling studies have highlighted the ability of rain to enhance warming and thaw (Iijima et al 2010, Connon et al 2018, Devoie et al 2019, 2021, Douglas et al 2020, Mekonnen et al 2021, Magnússon et al 2022), while others have reported a minimal (or cooling) impact (Hamm andFrampton 2021, Gao andCoon 2022).The ability of a single rain event to cause soil warming through advective heat transport has been shown in a wetland environment (Neumann et al 2019).Our study highlights the impact of a single rain event on warming and thaw of a warm permafrost plateau containing taliks.Transition locations, which our results indicate are vulnerable to the thermal inputs from rain, are highlighted in dark red and labeled as 'Rain Enhanced Thaw' in figure 10, which presents a conceptual overview of factors affecting permafrost thaw at our site.

Hummock and hollow location thaw was controlled by canopy cover and microtopography
During our study, hummock locations thawed gradually (8 ± 9 cm yr −1 ) while hollow locations thawed rapidly (44 ± 6 cm yr −1 ) (figure 5(B)).Observed hummock soil temperatures were always colder than hollow temperatures (figures 5(C) and (D)), with cold temperatures reaching deeper depths at hummock locations during the cold-season (figure 5(E)), and warm temperatures unable to deeply penetrate during the warm-season (figure 5(F)).These temperature regimes were expected since hollow locations had deeper snow (figure 4  melt-inducing longwave radiation to snowpack (Woo and Giesbrecht 2000, Lundquist et al 2013, Webster et al 2016).Microtopography can affect snow depth through accumulation in microtopographic depressions during wind-driven snow transport (Clark et al 2011).These relationships are reflected in our data establishing correlations between both canopy and relative-elevation with snow depth (figures 4(D) and (E)).Black spruce can reduce soil moisture through transpiration and direct interception of precipitation (Perron et al 2023), and higher relative-elevation locations have reduced soil moisture due to exiting runoff (Engstrom et al 2005, Graham et al 2022).In this way, dense canopy and high relative-elevation keep hummock locations drier (table 1), and therefore, cooler during the warm-season (figures 5(D) and (F)).These cooling effects from high canopy coverage and high relative-elevation are forms of ecosystem protection and are referred to as 'Canopy Protection' and 'Topographic Protection' in figure 10.
Without topographic variation, canopy cover would theoretically be the primary control on both soil moisture and snow depth.In this scenario, areas with dense canopy would thaw slowly due to lower soil moisture and shallower snow depth, while areas with minimal canopy cover would thaw more rapidly.Low canopy areas would, in turn, subside over time, becoming localized low points that accumulate moisture and snow.These wet depressions would continue to thaw and become wetter until a lack of oxygen due to soil saturation would cause back spruce death (Krause and Lemay 2022), leading to less canopy at topographic low points, and perpetuating the feedback loop (Dearborn et al 2021, Haynes et al 2021).The eventual result of this theoretical model for landscape progression is a positive correlation between relative-elevation and canopy cover like our site (SI figure S1).Without any black spruce, this scenario would advance more rapidly and uniformly.

Implications of a warmer and wetter future
Observed warming and thaw at this ecosystemprotected site provides a preview of expected future permafrost progression regimes in more northern latitudes.Rain-induced energy transport to transition locations provided sufficient energy to cause warming and thaw of low ice-content soil, due in part to hydrologically connected taliks (figures 9 and 10).In the absence of taliks, this mechanism could, over time, cause enhanced thaw and the development of taliks.As larger rain events become more common in northern high-latitudes (Cohen et al 2014, Kusunoki et al 2015, Liu et al 2021, Wang et al 2021, Thackeray et al 2022), and air temperatures continue to rise (Rantanen et al 2022), the energy contribution from this mechanism is expected to increase.This increase, however, is not expected across the entire permafrost zone.Precipitation changes are highly variable across northern high latitudes (Bintanja et al 2020).In addition, permafrost sites located in continental climates, such as the semi-continental climate of our site, have been shown to experience the greatest rain-induced warming through a combination of enhanced summer precipitation and warmer summer air temperatures, while maritime climates tend to experience a slight cooling effect from heavy rainfall (Hamm et al 2023).
Our analyses indicate the nine-fold increase in thaw rates observed between 2020 and 2022 (figure 3) was facilitated by low ice content and the presence of taliks within the permafrost plateau.Decades of exposure to a MAAT above 0 • C presumably caused permafrost soils to warm, taliks to form, and soil to lose ice content.This prior slow warming allowed energy inputs from three snowy and wet years, and two large warm-season rain events, to overcome the low remaining resistance to thaw.As temperatures in northern high-latitudes continue to warm more quickly than the global average (Rantanen et al 2022), causing more of the permafrost zone to experience mean annual temperatures above 0 • C, more permafrost will likely experience reduced ice content and become vulnerable to similar rapid thaw events.For example, in response to increased winter precipitation and snow at a colder and more stable site in northern Sweden, Webster et al (2016) observed permafrost soils warm while the late-season thaw depth remained relatively stable.This site is now more vulnerable to the type of future thaw observed at our site.

Conclusion
Based on high-resolution soil temperature data, repeat permafrost surveys, and measurements of environmental variables, we determined that the thaw of our ecosystem-protected permafrost site, and expansion of taliks, accelerated nine-fold during three consecutive wet and snowy years.Neither thaw rates, nor mechanisms of thaw, were uniform across the site.Hummock locations experienced only 8 ± 9 cm yr −1 of thaw, while hollow and transition locations thawed 44 ± 6 cm yr −1 and 39 ± 13 cm yr −1 , respectively (figure 5(B)).
Mechanisms of thaw between hollow and transition locations differed.Hollow locations had high thermal conductivity due to high soil moisture (table 1) and ample cold-season insulation from deep snow (figure 3(C)), while transition locations experienced substantial warming and thaw from thermal energy transported through subsurface taliks by runoff during individual rain events (figures 6 and 7(C), (D) and 9).As northern high-latitudes continue to become warmer and wetter, and rain events become more extreme, these observed mechanisms of warming and thaw are expected to increase.

Figure 1 .
Figure 1.Site overview.(A) Alaska map with an inset of the Kenai Peninsula shown in panel (B).(B) Kenai Peninsula map with an arrow indicating an inset of the site location shown in panel (C).(C) A Google Earth image of the site with an inset of our core instrumentation area shown in panel (E).Instrumented high-elevation (hummock) locations are green, mid-elevation (transition) locations are orange, and low-elevation (hollow) locations are light blue.Wetland and permafrost areas along with rain gauge, temperature/RH sensor, and two soil moisture sensor locations are labeled.The across-site permafrost transect is shown as a white line.(D) A field photo of plateau topography and snow cover heterogeneity.(E) A digital elevation model of our core instrumentation area.Shapes denote the instrumented locations shown in the inset of panel (C).Yellow topographic contour lines label every 0.5 m of elevation change.Names and locations of individual DTPs are shown in SI figure S4.

Figure 2 .
Figure 2. Meteorological forcing from 2015 to 2022.Line graphs show stacked annual time series of (A) snow water equivalent and (B) cumulative precipitation.Data for 2020 is blue, 2021 is orange, and 2022 is green.Data from 2015-2019 are gray or black with varied line styles displayed in the legend.The snow-to-rain transition period is shown as a gray box in Panel (B).Black-bordered orange (i) and green (ii) boxes in Panel (B) show periods of intense rain during July 2021 and 2022, respectively.

Figure 3 .
Figure 3. Permafrost thaw along the across-site transect shown in figure 1(C).Depths to the permafrost table measured in 2015 by Jones et al (2016) are shown as black circles.Our measurements from September 2020 and September 2022 are shown as blue diamonds and green squares, respectively.Surface elevation is marked with a thin solid black line based on data from Jones et al (2016).Maximum probe depth is marked with a dashed black line.Gray vertical bars label locations that had permafrost deeper than our 3.1 m probe in 2022.Wetland locations (positions 23-63) are labeled.

Figure 4 .
Figure 4. Environmental variables across the site.(A)-(C) Box and whisker plots showing: (A) relative-elevation of sensor locations, (B) canopy percentage, and (C) snow depth.On each box plot, the central line indicates the median, 'x' marks the mean, edges of the box indicate the 25th and 75th percentiles, whiskers indicate the most extreme non-outlier data points, and a red plus marker indicates an outlier greater than 2.7 standard deviations from the mean.Letters on top of boxes indicate data sets with statistically significant differences in mean as determined through one-way ANOVA.(D), (E) Scatter plots showing: (D) snow depth versus relative-elevation and (E) snow depth versus canopy cover.

Figure 7
Figure 7 shows the soil thermal impact of the 2021 rain event on example locations from each microtopographic group.Hummock and hollow locations experienced short-term perturbations in soil temperature where temperatures briefly warmed and then returned to near-linear rates of seasonal warming (figures 7(B) and (E)).At transition location #1, soil temperatures between 40-140 cm increased from 0 • C to <1 • C positive temperatures while soils at 20-30 cm depths remained ⩽0 • C (figure 7(C)).Transition location #2 had a greater response than transition location #1.Soils between 20-120 cm depth at transition location #2 warmed from ⩽0 • C to 1 • C-4 • C (figure 7(D)).Transition location #2 had a greater proportion of surrounding higher elevation area than transition location #1.In the context of hydrologic flow, the proportion of surrounding higher elevation area can be represented as the ratio of area contributing recharge to lower elevations (Area C ,) to area receiving recharge

Figure 5 .
Figure 5. Permafrost progression and thermal metrics among relative-elevation groups.Box and whisker plots showing: (A) Depth to the permafrost table measured on 9 September 2022, (B) average annual deepening of the permafrost table between 2020 and 2022, (C) minimum and (D) maximum temperature reached at a depth of 10 cm between 2020 to 2022, (E), (F) maximum penetration of (E) −0.5 • C and (F) +3 • C. Boxplot statistics described in figure 4.

Figure
Figure Thermal response to rain among microtopographic groups.Box and whisker plots showing: (A), (B) soil profile thickness of hollow, transition, and hummock locations that warmed from ⩽0 • C to above-freezing temperatures from one week before the largest rain event of the year to one week after the conclusion of the rain event in (A) 2021 and (B) 2022.Boxplot statistics are described in figure4.

Figure 7 .
Figure 7. Soil temperature responses to the rain event in July 2021.(A) Air temperature (red, left axis) and cumulative precipitation (blue, right axis).(B)-(E) Soil temperature at depths from 10 to 140 cm below the soil surface (depths are labeled on the figure and colored according to the legend) for (B) hummock location, (C) transition location #1 with zoomed inset (i), (D) transition location #2, and (E) hollow location.The hummock location and transition location #1 were instrumented with soil moisture sensors, see figure 8. Y-axis scales vary by location.

Figure 8 .
Figure 8. Hydrologic response to rain events in July 2021 (left column) and July 2022 (right column).Colormap plots show soil moisture as percent water for: (A) hummock location and (B) transition location, which are the hummock and transition #1 locations in figure 7(B), (C).A transparent white mask shows areas with a soil temperature of ⩽0 • C. (C) Water table depth of two hollow locations (green and light blue, left axis) and cumulative precipitation (dark blue, right axis).Y-axis scales vary by location.
(C)).Soil moisture and thermal conductivity values are average values for soil between 0 and 10 cm from one transition and one hummock location before (3 July) and after (2 August) the rain events in 2021 and 2022 (see figures 8(A) and (B)).Thermal conductivity values of soils at ⩽0 • C are given as a range corresponding to maximum and minimum potential ice content (see methods 2.6).Standard deviations are given for average values.
(B)).Increased thickness of taliks, which accompanied permafrost thaw at most site locations, has biogeochemical consequences such as increased methane emissions (Heslop et al 2015, O'Neill et al 2020), thermal consequences such as greater thaw rates (Connon et al 2018, Devoie et al 2019), and hydrologic consequences including increased connectivity (Connon et al 2018, Devoie et al 2019, 2021).It is known that rain and snow accelerate thaw.During three rainy and snowy years like our site experienced, Iijima et al (2010) observed a transition from permafrost stability to thaw (∼5 cm yr −1 ) at a colder (−8 • C MAAT) discontinuous permafrost site.Especially rainy years accelerated observed thaw rates in Interior, AK (Douglas et al 2020), and rain from a single wet summer increased thaw for multiple years at a Siberian tundra site (Magnússon et al 2022).Seasonally frozen layer thickness has been shown to increase with snow depth in both field (Johansson et al 2013, O'Neill and Burn 2017) and modeling (Zhou et al 2013, Park et al 2015) studies.At our site, locations with greater snow depth (figure 4(C)) thawed more rapidly (figure 5(B)), supporting the notion that snowy conditions between 2020 and 2022 exacerbated permafrost degradation.Composite effects of consecutive wet and snowy years likely contributed to the observed increase in thaw.
Over 100 cm of above-permafrost transition soil warmed from ⩽0 • C to >0 • C during individual rain events in 2021 and 2022 (figures 6, 7(C) and (D)).Two potential mechanisms of rain-induced energy transport could have caused this warming: increased soil thermal conductivity due to increased soil moisture (Farouki 1981, Hinkel and Nelson 2003, Subin et al 2013) and advective heat transport (Iijima et al 2010, Neumann et al 2019, Douglas et al 2020, Mekonnen et al 2021, Magnússon et al 2022).

Figure 10 .
Figure10.Conceptual overview.The top section shows a permafrost plateau cross-section.Light red represents areas expected to thaw rapidly with or without large rain events, and dark red represents areas only expected to thaw rapidly in the presence of large rain events.Light blue geometric shapes represent soil ice content.A blue triangle marks the top of the water table.Thick black arrows show subsurface flow paths of rain-induced runoff that can warm soil through advective heat transport.A thin black arrow shows the thaw that took place between 2020 and 2022.A dashed black line represents the depth to permafrost before 2020.The seasonally frozen layer is shown in light tan, the talik layer is shown in dark tan, and permafrost is shown in gray.The bottom section shows the relative magnitude of canopy protection, topographic protection, rain-enhanced thaw, and annual thaw rates among hummock, transition, and hollow locations.Canopy protection and topographic protection refer to the slowing of permafrost thaw caused by dense canopy cover and high relative-elevation, both of which reduce snow depth and soil moisture, which are known to accelerate permafrost thaw (Discussion 4.3).