Local-scale heterogeneity of soil thermal dynamics and controlling factors in a discontinuous permafrost region

In permafrost regions, the strong spatial and temporal variability in soil temperature cannot be explained by the weather forcing only. Understanding the local heterogeneity of soil thermal dynamics and their controls is essential to understand how permafrost systems respond to climate change and to develop process-based models or remote sensing products for predicting soil temperature. In this study, we analyzed soil temperature dynamics and their controls in a discontinuous permafrost region on the Seward Peninsula, Alaska. We acquired one-year temperature time series at multiple depths (at 5 or 10 cm intervals up to 85 cm depth) at 45 discrete locations across a 2.3 km2 watershed. We observed a larger spatial variability in winter temperatures than that in summer temperatures at all depths, with the former controlling most of the spatial variability in mean annual temperatures. We also observed a strong correlation between mean annual ground temperature at a depth of 85 cm and mean annual or winter season ground surface temperature across the 45 locations. We demonstrate that soils classified as cold, intermediate, or warm using hierarchical clustering of full-year temperature data closely match their co-located vegetation (graminoid tundra, dwarf shrub tundra, and tall shrub tundra, respectively). We show that the spatial heterogeneity in soil temperature is primarily driven by spatial heterogeneity in snow cover, which induces variable winter insulation and soil thermal diffusivity. These effects further extend to the subsequent summer by causing variable latent heat exchanges. Finally, we discuss the challenges of predicting soil temperatures from snow depth and vegetation height alone by considering the complexity observed in the field data and reproduced in a model sensitivity analysis.


Introduction
The Arctic region is experiencing rapid warming in response to global climate change (Serreze andBarry 2011, Biskaborn et al 2019).Soil warming in permafrost regions enhances microbial decomposition of accumulated organic carbon, causing the substantial release of greenhouse gasses into the atmosphere (Zimov et al 2006, Tarnocai et al 2009, Schuur et al 2015).Warming-stimulated plant growth (e.g.shrub expansion), on the other hand, offsets soil carbon emission by plant carbon uptake and accumulation (Mekonnen et al 2021).Thus understanding soil thermal dynamics and their controls is critical to accurately predicting the biogeochemical responses of permafrost ecosystems to climate change (Jorgenson et al 2010, Loranty et al 2018).
One major challenge in quantifying current and future impacts of climate change on permafrost regions is the difficulty of understanding the controls of local heterogeneity in above-and belowground properties and processes, and representing them adequately in climate or ecosystem models (Riseborough et al 2008, Lara et al 2020).Local heterogeneity of snow, vegetation, and soil characteristics at the meter scale can lead to greater variability in soil temperatures than latitudinal air temperature gradients at the 100 kilometer scale (Shirley et al 2022).Impacts from multiple factors complicate the understanding and prediction of soil thermal dynamics.For example, heterogeneous snow cover impacts winter soil temperatures through various insulating effects (e.g.Zhang 2005, Romanovsky et al 2010a, Gisnås et al 2014).Vegetation influences the soil temperature by impacting snow distribution in the winter (e.g.Sturm et al 2001, Komarov andSturm 2023) and ground shading in the summer (e.g.Blok et al 2010).Topographic variables (e.g.elevation, aspect, and slope) impact soil temperatures by controlling snow accumulation and solar radiation (e.g.Gubler et al 2011).Beyond above-ground controls, subsurface heat processes (e.g.conduction, advection, and latent heat exchange) also affect soil temperatures but their mechanisms are more complex and dependent on multiple factors such as soil texture (e.While studies have investigated the interactions between soil temperature and various environmental factors, the large geographical extent and associated low spatial resolution of most studies have not been conducive to documenting how local scale heterogeneity influences these interactions (e.g.Romanovsky et al 2010b, Cable et al 2016, Kropp et al 2020).Most existing studies on local variability of soil temperatures only focus on ground surface temperatures (GSTs) (e.g.Gubler et al 2011, Grünberg et al 2020, Zhang et al 2021, Evans et al 2022), considering little on subsurface processes or deep soil temperatures that better represent long-term thermal states.In addition, the coupled impact of various processes on soil temperature varies across sites and studies, and thus indicates the need for observations across various environments.In discontinuous permafrost regions, where soil temperatures close to the 0 • C threshold represent a tipping point for several hydro-biogeochemical processes (e.g.Natali et al 2019), field-based investigations of interactions between above-and below-ground processes are still lacking (e.g.Uhlemann et al 2021).Improving our understanding of controls and indicators of soil thermal dynamics is crucial for developing processbased models that reproduce various trajectories in transitional landscapes, as well as remote sensing products that can be used to estimate soil thermal and physical states over large spatial extents.
In this study, we aim to answer the following questions related to local heterogeneity (within km-scale spatial extent) in soil temperature in a discontinuous permafrost region: 1.How is spatial variability in GST related to spatial variability in ground temperature (GT) at different times and time scales?2. How do above-and below-ground processes impact soil temperatures seasonally and annually?3. Can soil thermal states be estimated from the above-ground dynamics such as vegetation type, vegetation height, and snow depth?

Study site and datasets
This study was performed at a watershed (∼2.3 km 2 area) located about 40 km northwest of Nome, Alaska on the southern Seward Peninsula (figure 1(a)).The watershed's elevation ranges from 60 to 280 m and it is characterized as a discontinuous permafrost region (Jorgenson et al 2008) where near-surface permafrost and talik co-exist (Léger et al 2019, Uhlemann et al 2021, Bennett et al 2022).From 2013 to 2017, the mean annual air temperature is −1.02 • C and the yearly rain and snow precipitation are 450 mm, and 1705 mm, respectively (Léger et al 2019).The site has mostly south-east facing slopes and the general geology consists of a thin layer of organic materials overlying silts, with schist bedrock beneath (Uhlemann et al 2023).
Vertically-resolved profiles of soil temperature were acquired every 15 min at discrete locations across the watershed (figure 1(b)) using distributed temperature profiling (DTP) systems.Each DTP system included a data acquisition logger connected to a 120 cm long thin (1 cm diameter) probe containing 15 digital sensors recording temperature with factoryassured sensor accuracy of 0.1 • C (Dafflon et al 2022).The sensors along the probe were spaced 5 cm from the first (shallowest) to the 7th sensor and 10 cm from the 7th to the 15th sensor.At each location, a guide hole was drilled and the probe was inserted vertically into the ground with the second sensor from the top of the probe located slightly (1-2 cm) below the ground surface (defined as 0 cm depth) to capture the GST.This enables temperature measurements from 5 cm above ground surface down to 105 cm depth.In early October 2019, 93 DTP probes were deployed along multiple transects that spanned topography or vegetation cover gradients.Due to the pandemic, the probes were not accessed and maintained until June 2021, causing many probe failures.We therefore only focus on 45 locations (figure 1(b)) that produced one entire year time series from 4 October 2019 to 4 October 2020 at all depths after linear interpolation on consecutive missing data less than 2 h.The selected locations mainly vary by vegetation and snow accumulation and thus do not fully account for other factors such as geomorphological features.Long-term soil records of GST and GT at ∼1 m depth at six locations across the watershed show similar annual variation from 2016 to 2022 (Romanovsky et al 2022), making our analysis of one-year temperature time series a good representation of the general thermal dynamics (supplementary section S1).
We classified the vegetation types at each DTP location into four categories based on the circumpolar arctic vegetation map (CAVM) (Walker et al 2005): (1) graminoid tundra (CAVM unit G3: nontussock sedge, dwarf-shrub, moss tundra) (2) wetland complex (CAVM unit W2: sedge, moss, dwarfshrub wetland complex), (3) dwarf shrub tundra (CAVM unit S1: erect dwarf-shrub, moss tundra) and (4) tall shrub tundra (taller than 40 cm, corresponding to CAVM unit S2: low-shrub tundra) (figure 1(c)).A local weather station within the watershed (figure 1(b)) provided meteorological records, which included air temperature at 1.5 m above the ground, ultrasonic snow depth, and cumulative rainfall (Busey et al 2017).In order to evaluate the influence of snow depth, vegetation height, and soil moisture content on soil temperature variability, we use data collected for these environmental factors from a different year as proxies as the site was inaccessible during 2020.UAV images collected on 1 April 2022 and 19 July 2017 were used to derive snow depth and vegetation height, respectively, and volumetric soil moisture content of the top 20 cm was collected in June 2022.Additional information on these datasets and justification for their use as proxies can be found in supplementary section S2.

Data analysis
Temperature time series were pre-processed to adjust the sensor depth of the DTP systems that rose upward due to frost jacking (supplementary section S3).The DTP systems in 28 locations were impacted by frost jacking, leading to an average and maximum displacement of 9.5 cm and 20 cm respectively at the end of summer.We therefore focused our analysis on the temperature data to a maximum depth of 85 cm to ensure a complete one-year time series at all selected locations.Using the multi-depth temperature time series, we first calculated the mean temperature at different time scales including daily, annual, and seasonal.We converted the mean daily temperature time series to the daily freeze/thaw index time series, which used −1, 0, and 1 to represent the status of frozen (<−0.3 • C), freezing/thawing (between −0.3 • C and 0 • C), and thawed (>0 • C), respectively.We used −0.3 • C as the threshold to define the complete frozen status considering the depression point of water mixed with porous media as commonly used in the literature (e.g.Cable et al 2016, Farquharson et al 2022)).From the acquired freeze/thaw index time series, we determined the frozen/thawed depth and first snow-free date.We then assessed the decoupling between air and GST using the freezing (or thawing) n-factors calculated as ground surface freezing (or thawing) degree days divided by air freezing (or thawing) degree days.We also estimated the summer soil thermal diffusivity from pairs of temperature time series at 5 cm depth intervals up to 20 cm depth using an analytical method (Hinkel 1997) (supplementary section S4).
To understand the heterogeneity and controls of soil thermal dynamics across the measurement locations, we grouped locations with similar temperature responses using hierarchical clustering following Cable et al (2016).We used a one-year time series of mean daily temperature and daily freeze/thaw index at all depths as input to compute the pairwise Euclidian distance between each location.The agglomerative hierarchical clustering was performed using ward linkage and then visualized as a dendrogram, which displays the similarity of the input datasets.Note that we normalized the temperature across all locations at each depth to avoid large absolute temperature differences in the shallow soils dominating the clustering process.The addition of the freezethaw index time series helped differentiate frozen and thawed status for temperatures around 0 • C. Additional information on hierarchical clustering is presented in supplementary section S5.While other studies have grouped locations based on vegetation types (e.g.Grünberg et al 2020, Kropp et al 2020), we grouped the locations based on their temperature characteristics, as vegetation types are not the only factor determining soil thermal dynamics.
We compared relationships of soil temperature, snow depth, and vegetation height observed in our field data to those obtained from processed-based model simulations from a recently published sensitivity analysis study (Shirley et al 2022).Shirley et al (2022) applied a global sensitivity analysis to ecosys, a process-rich terrestrial ecosystem model to understand how variability in soil properties, landscape position, and weather forcing affect variability in soil thermal regimes, vegetation dynamics, and carbon cycling.They applied ecosys to simulate the response of a dynamic plant canopy and multi-layer soil column to hourly meteorological inputs from year 1836-2019.Forcing and parameterization for the model sensitivity analysis were configured to represent the conditions in the same watershed studied in this work, including variation in 24 factors related to soil properties, boundary conditions, and weather forcings (supplementary section S6).We examined the relationship between the mean annual GT, latewinter snow depth, and vegetation height in 2015-2018 from their 604 simulation runs and compared it to the relationship developed using our field data.We use this comparison to examine the extent to which GT can be predicted by the above-ground metrics.

Heterogeneity of thermal states
The soil temperature across the studied watershed exhibits large temporal and spatial variability (figures 2(c) and (d)).The mean annual temperature of all monitoring locations has a mean and standard deviation of 1.18 ± 1.12 • C, 0.57 ± 1.28 • C, and 0.73 ± 1.21 • C at a depth of 0 cm, 20 cm, and 85 cm, respectively.The largest spatial variability is observed in the winter season, especially between 1 January and 15 March when the snow accumulates (figures 2(a) and (d)).A short period of elevated variability at shallow depths is found between May 1st and 1 July, which corresponds to the various timing of snowmelt completion (figures 2(b) and (d)).In summer, temperature spatial variability is low at 5 cm above and at the ground surface, but slightly increases with increasing depth (figure 2(d)).The smallest spatial variability is observed between 1 November and 1 January, as well as between 15 April and 20 June, which corresponds to the soil freezing and thawing periods with temperatures close to zero. Figure 2(e) presents how well the spatial variability in temperatures at various times, time scales, and depths correlates to the spatial variability in the mean annual ground temperature at a depth of 85 cm (MAGT85) that better represents the long-term thermal states.Mean annual and winter temperatures show strong correlations (r > 0.75) with MAGT85 across 45 locations, while summer temperatures at shallow depths show weak or negative correlations with MAGT85 (figure 2(e)).

Classification of thermal dynamics and their link to vegetation type
Soil thermal dynamics observed at the 45 monitoring locations are classified via hierarchical clustering into three different thermal groups, referred to as cold, intermediate, and warm groups based on their mean annual temperature (figures 3(a) and (b)).The three thermal groups show a strong match with the vegetation types, with the cold, intermediate, and warm groups dominated by graminoid tundra, dwarf shrub tundra, and tall shrub tundra, respectively (figure 3(a)).Wetland complex locations show more complex interaction with microtopography or tall shrubs and they will not be discussed extensively due to the limited number of locations.The difference in temperature between the three thermal groups is larger in winter (December to March) than in summer (June to September) (figures 3(c)-(f)).The temperature in the warm group mostly stays above 0 • C at deeper depths (e.g.85 cm) during the winter (figure 3(f)).Summer temperature variability is low at 0 cm but increases considerably with depth (figures 3(c)-(f)).Frozen and thawed depths also show large differences between the three thermal groups (figure 3(g)).In early November, all locations freeze at similar timing and speed down to about 10 cm depth.Later, the frozen depths diverge with the median frozen depths of cold and intermediate groups reaching 85 cm on 7 January and 10 February, respectively.In contrast, the median frozen depth of the warm group only reaches 35 cm by 1 March.Early thawing occurs at similar timing and speed for all thermal groups, reaching 40 cm by mid-June.Because deep soils remain unfrozen throughout the winter in the warm group, thaw depth abruptly reaches 85 cm in mid-June.In contrast, cold and intermediate groups remain frozen at 85 cm for a month or longer.

Variability of thermal metrics and environmental factors
The controls on the three classified thermal groups are evaluated using various metrics and properties (figure 4).The difference between the MAGT85 and mean annual GST (MAGST) is similar for all groups, but the difference between the mean annual air temperature and MAGST increases from cold to warm groups (figure 4(a)).Freezing n-factors decrease from cold to warm groups and, particularly in deeper soils, the number of frozen days decreases (figures 4(b) and (c)).Correspondingly, the snow depth increases from cold to warm groups but with large overlaps, resulting from uncertainty linked to the use of the snow depth proxy and/or the impact of other environmental factors on the soil temperature (figure 4(g)).
The estimated snow-free date varies over 30 d across all locations with the warm group having significantly later snowmelt completion than the other two groups (figure 4(d)).In summer, the three groups show similar daily mean temperatures at shallow depth (figure 3(c)) and 5 cm above the ground surface (figure 4(e)), and the thawing n-factor does not vary significantly between the three groups (p > 0.05) (figure 4(b)).The warm group, which is covered by taller vegetation (25th-75th percentile: 0.30-1.64m) than the cold and intermediate groups (25th-75th percentile: 0.13-0.22m, 0.05-0.13m, respectively), is characterized by lower maximum and higher minimum daily temperatures than the other groups (figures 4(e) and (f)).
Soil moisture (top 20 cm) is significantly higher in the cold group than in intermediate and warm groups and shows a strong correlation with the vegetation type (e.g.graminoid tundra and wetland complex soils are typically wetter than tall and dwarf shrub soils).Summer soil thermal diffusivity increases from 0 to 20 cm depths, but does not differ significantly between the thermal groups at any depth.The median thermal diffusivity of all locations gradually increases from 0.05 mm 2 s −1 at 0-5 cm depth to 0.3 mm 2 s −1 at 15-20 cm depth, which falls within the thermal diffusivity range of peats and clay/silt materials, respectively.

Comparison of relationships observed in field data and in ecosystem model sensitivity analysis
We compared relationships observed using field data to those obtained from processed-based model simulations from a recently published sensitivity analysis study using ecosys (Shirley et al 2022).Field data show that deep snowpacks are associated with tall vegetation in the watershed, but in short vegetation snow depth is highly variable (ranging from less than 10 cm to more than 1 m) (figure 5(c)).The relationship between snow depth and soil temperature observed in the field data is similar to that produced by the model sensitivity analysis, although with less scatter since environmental factors were varied over large ranges in the sensitivity analysis (figure 5(a)).A binary relationship between vegetation height and soil temperature is observed in the field data, with soil temperature exhibiting strong variability in vegetation less than 50 cm and reaching a plateau of around 2 • C at vegetation higher than 50 cm (figure 5(b)).The process model does not represent canopy interception of wind-blown snow, and the relationship between vegetation height and soil temperature is weaker (figure 5(b)).

Impacts of above-ground processes on soil thermal dynamics
Our datasets demonstrate that the spatial variability in soil thermal states at an annual scale is primarily influenced by the above-ground decoupling  S4) and therefore we use the combined data from multiple years to present a general trend of these relationships.Note that the scattering observed in the relationship is not explained by varied air temperature inputs (figure S4).
between mean annual air temperature and MAGST rather than the below-ground decoupling between MAGST and MAGT85 (figure 4(a)).Besides, the spatial variability in soil mean annual temperature is mainly determined by variability in winter temperature rather than summer temperature, highlighting the strong control of winter processes (figures 2(e) and 3).The thermal insulation effect of snow is demonstrated by the sharp difference of freezing nfactor between the three classified thermal groups (figure 4(b)) (Kade et al 2006).The warm group is covered by taller vegetation and snow depth than the other two clustered groups (figures 4(f) and (g)).This observation is consistent with previous studies that have linked snow accumulation by tall vegetation with soil warming in arctic tundra and boreal ecosystems (e.g.Grünberg et al 2020, Kropp et al 2020).The cold group has significantly higher freezing nfactors than the intermediate group, but vegetation height and snow depth proxies are similar between the two groups (figures 4(b), (f), and (g)).Besides the uncertainty linked to using snow depth from another year as proxy, this may also be due to reduced snow density inside the dense dwarf shrub canopy that leads to enhanced insulation near the ground surface.Additionally, the UAV survey may fail to sense the ground surface due to the dense canopy, leading to the underestimation of vegetation height and snow depth (Harder et al 2020).Despite the strong association between vegetation type and soil temperature, the presence of multiple vegetation types within each classified thermal group highlights the role of other factors determining soil temperatures (sites with symbols in figure 3(a)).For example, we observe that two locations on local high microtopography (therefore with reduced snow depth) are classified into the cold group.We also observed low-stature vegetation locations surrounded by or close to groups of tall shrubs show increased soil temperature, likely due to the increased snow depth by shrub fencing (Sturm et al 2001).
Spatial variability in the timing and amplitude of summer processes shows some impact on shortterm variation in soil temperature, but very limited impact on the mean annual temperature.In the thawing season, snowmelt timing variability causes only a temporary increase in spatial variability of shallow soil temperatures (figure 2(d)).The earlier first snow-free day of cold and intermediate groups leads to higher ground heat fluxes and slightly warmer shallow soil temperatures than the warm group from late May to early July (figure 3(d)).Considering the entire warm season, less variable thawing n-factors between the three groups demonstrate small variability in the decoupling between summer air temperature and GST driven by vegetation shading and surface soil and organic materials, which is consistent with previous studies (e.g.Paradis et al 2016, Grünberg et al 2020, Zhang et al 2021).Although the temperature at 5 cm above the ground surface during daytime reveals a canopy shading effect of tall shrubs, this effect is partly balanced by warmer night temperatures possibly due to the attenuated wind (figure 4(e)).Consequently, summer GSTs do not show significant differences between the three groups (figure 3(c)).

Impacts of below-ground processes on soil thermal dynamics
The impact of above-ground processes on the soil temperature is further modulated by below-ground processes, which cause variability in the seasonal soil thermal states.In summer, the warm group is significantly warmer than the other two groups below 20 cm depth, even though their GSTs are similar to each other.(figure 3(c)).This cannot be explained by heat conduction alone given the comparable thermal diffusivity of shallow soils across three groups (figure 4(i)).Instead, latent heat is a more determinant factor as evidenced by shorter frozen periods of the warm group (figure 4(c)), which makes summer energy penetrate into the deep soil faster without being extensively used for thawing the frozen soils.In addition, the unfrozen soil can allow faster and deeper infiltration of water and associated heat (Uhlemann et al 2021).
Though the thermal diffusivity of thawed soil estimated during the summer does not vary significantly across the site, large spatial variability in the thermal diffusivity of winter soils increases soil temperature heterogeneity.Soils with high water content have a much higher thermal diffusivity when frozen, since the thermal diffusivity of ice (1.02 mm 2 s −1 ) is much larger than water (0.14 mm 2 s −1 ) (Mcgaw et al 1978, Farouki 1981) (figure 4(i)).The cold group is characterized by wet soils that freeze quickly since snow depths are shallow, leading to high thermal diffusivity that enhances heat transfer from the soil to the atmosphere.The intermediate and warm groups, in contrast, are characterized by drier soils that freeze more slowly or do not freeze due to more snow insulation, so thermal diffusivity remains lower and heat transfer is slower.The results described above demonstrate how positive feedbacks created by subsurface processes augment the impact of snow depth on winter soil temperatures.

Implications for predicting soil thermal states from GST and above-ground properties
Our analysis demonstrates that the relative spatial variability in MAGST strongly indicates the spatial variability in MAGT85, mainly due to their strong correlation in winter (figure 2(e)).This indicates that GST, nowadays frequently monitored with commercial temperature loggers (e.g.Gubler et al 2011, Grünberg et al 2020), offers valuable information on deeper soil thermal states when they cover winter or one-year period.However, summer GST is highly decoupled from MAGT85, meaning neither groundnor airborne-based summer GST measurements infer the spatial variability in deeper soil thermal states (figure 2(e)).The strong association between vegetation types and clustered thermal groups (figure 3(a)) indicates that remote sensing-based products associated with vegetation structure or type may provide valuable information on subsurface thermal regimes, consistent with previous findings (e.g.Cable et al 2016, Kropp et al 2020).However, challenges remain when predicting soil temperature from snow depth or vegetation height alone, as evident by their relatively weak correlations observed in both the field data and model sensitivity analysis (figure 5).While the field data-based relationship in this study may suffer from uncertainties linked to the use of snow depth proxy from another year, the relationship from sensitivity analysis shows the same level of complexity, which results from the varying influences of multiple hydrological processes/properties (soil texture, water dynamics, etc).

Conclusion
We analyzed the local-scale spatial heterogeneity of soil thermal dynamics and their controls in a discontinuous permafrost region.The large spatial variability in mean annual temperatures is primarily influenced by winter temperatures, in contrast to summer temperatures, which exhibit less spatial variability.The mean annual or winter season GSTs are well correlated to the MAGT85.We attribute the large spatial variability in soil thermal dynamics to the variable thermal decoupling from dynamic snowpacks, which modulates ground insulation, winter heat conduction, and the subsequent summer latent heat exchanges.In contrast, the summer season surface thermal decoupling from vegetation and organic layers does not vary significantly across the study locations.We illustrate that vegetation types are good indicators of soil relative thermal states.However, quantitative prediction of soil temperature from snow depth or vegetation height is challenging, as illustrated by both the field observations and model sensitivity analysis.High-resolution monitoring and surveys on subsurface characteristics remain necessary to further validate the model simulation and unravel their complicated impacts on soil temperatures.Besides, advances in monitoring snow temporal dynamics and snowpack properties may improve the prediction of soil thermal dynamics and facilitate the quantification of the impact of other environmental factors by disentangling the snow effects.
g. Johnson et al 2013), soil moisture (e.g.Subin et al 2013, Clayton et al 2021), and rainfall (e.g.Neumann et al 2019).In addition, distribution of water bodies (e.g.Walvoord and Kurylyk 2016), human interventions (e.g.Ehlers et al 2022), and wildfire (e.g.Holloway et al 2020) can further modify the soil thermal states in a permafrost region.

Figure 1 .
Figure 1.Overview of the study site and settings.(a) Location of the study site on the southern Seward Peninsula in Alaska (Google Earth imagery).(b) Locations of the weather station and distributed temperature profiling (DTP) systems across the watershed (UTM Zone 3 N, ESRI imagery).Contour lines depict elevation intervals of 20 m.(c) Photos of the DTP system (logger contained in the 7.5 cm long jar) deployed in four different vegetation environments.

Figure 2 .
Figure 2. Weather forcing and spatial variability in soil temperature across 45 locations in the studied watershed from October 2019 to October 2020.(a) Air temperature at 1.5 m, snow depth, and cumulative rainfall measured at the local weather station (Busey et al 2017); (b) Number of locations without snow cover in the snowmelt season inferred from the first snow-free date; (c) and (d) Mean and standard deviation of the mean daily, annual, winter (December to March), and summer (June to September) temperatures across the 45 locations; (e) Pearson correlation coefficient between the mean daily, annual, winter, and summer temperatures at each depth versus MAGT85 across the 45 locations.Scatter plots between mean annual, winter, or summer GST versus MAGT85 are shown in supplementary figure S5 as examples of strong and weak correlations.

Figure 3 .
Figure 3. Hierarchical clustering of the soil temperature data.(a) Dendrogram of identified soil thermal groups and comparison to collocated vegetation types.Solid symbols indicate special environmental settings.Local microtopographic high and low represent bumps and depressions occurring within a few meter scale, respectively.Locations A5, D9, and D6 have identified shallow permafrost based on temperatures at 85 cm depth.(b)-(g) Mean annual (b), mean winter and mean summer (c) temperature at different depths, time series of daily mean temperature at depths of 0 cm (d), 20 cm (e), and 85 cm (f), and time series of frozen/thawed depths (g) visualized as the median (solid line) and 25th to 75th percentile (shaded area) of each classified thermal group.Each time series is from October 2019 to October 2020.

Figure 4 .
Figure 4. Selected thermal metrics and environmental properties of three classified thermal groups with vegetation types overlaid.Letters (a to c) denote significant differences, i.e. two groups without the same letter are significantly different (p < 0.05).No letter means no significant difference between the three thermal groups.(a) Difference between mean annual ground surface temperature (MAGST) and mean annual air temperature (MAAT) and difference between mean annual ground temperature at a depth of 85 cm (MAGT85) and MAGST.(b) Freezing and thawing n-factors.(c) Number of frozen days in a year at 0, 20, and 85 cm depths.(d) First snow-free date plotted as the number of days relative to 1 May, 2020.(e) Mean, maximum, and minimum daily temperature at 5 cm above ground averaged over the summer (June-September).(f) Vegetation height.(g) Snow depth derived from the UAV image surveyed in the 2022 winter as a proxy.(h) Soil moisture.(i) Thermal diffusivity of thawed soils estimated from temperature time series in summer.Shading areas from bottom to top indicate the thermal diffusivity range of peats/organic materials, clay/silt, and water-saturated sand/gravels, respectively (Márquez et al 2016).The arrow shows an example of soil thermal diffusivity increase from unfrozen to frozen status (Mcgaw et al 1978).

Figure 5 .
Figure 5. Relationships between soil temperature, snow depth, and vegetation height from field data and from sensitivity analysis using model simulations (Shirley et al 2022).The sensitivity analysis includes 604 runs with different inputs involving 24 factors (supplementary table S1).The relationships from sensitivity analysis are established based on the mean annual ground temperature at a depth of 85 cm (MAGT85) of years 2015 (October 2015-October 2016), 2016, 2017, and 2018, and the averaged snow depth between 15 March and 15 April of corresponding winters.The relationships from different years show similar responses (supplementary figureS4) and therefore we use the combined data from multiple years to present a general trend of these relationships.Note that the scattering observed in the relationship is not explained by varied air temperature inputs (figureS4).