Winter soil temperature varies with canopy cover in Siberian larch forests

In the Arctic, winter soil temperatures exert strong control over mean annual soil temperature and winter CO2 emissions. In tundra ecosystems there is evidence that plant canopy influences on snow accumulation alter winter soil temperatures. By comparison, there has been relatively little research examining the impacts of heterogeneity in boreal forest cover on soil temperatures. Using seven years of data from six sites in northeastern Siberia that vary in stem density we show that snow-depth and forest canopy cover exert equally strong control on cumulative soil freezing degrees days (FDDsoil). Together snow depth and canopy cover explain approximately 75% of the variance in linear models of FDDsoil and freezing n-factors (n f; calculated as the quotient of FDDsoil and FDDair), across sites and years. Including variables related to air temperature, or antecedent soil temperatures does not substantially improve models. The observed increase in FDDsoil with canopy cover suggests that canopy interception of snow or thermal conduction through trees may be important for winter soil temperature dynamics in forested ecosystems underlain by continuous permafrost. Our results imply that changes in Siberian larch forest cover that arise from climate warming or fire regime changes may have important impacts on winter soil temperature dynamics.


Introduction
Amplified warming in the Arctic (Rantanen et al 2022) is driving widespread changes across the Arctic (Post et al 2019).Key changes include altered vegetation productivity (Frost and Epstein 2014, Myers-Smith et al 2020, Berner and Goetz 2022), increased fire frequency and extent (Descals et al 2022), and permafrost thaw leading to increased carbon emissions, altered hydrology, and rapid geomorphic change associated with ground ice degradation (Liljedahl et al 2016, Turetsky et al 2020, Schuur et al 2022).These changes impact wildlife, traditional livelihoods, and human infrastructure (Post et al 2019).Precipitation change also impacts these processes in important ways; however, the magnitude and direction of precipitation change with warming are heterogeneous, dependent on the seasonality of temperature change, and generally less well understood (Box et al 2019).
Snow cover change is particularly important due to its influence on a series of biogeophysical, biogeochemical, and hydrological processes.Snow modulates permafrost temperatures by insulating soils from frigid winter air temperatures (Stieglitz et al 2003, Zhang 2005, Zhang et al 2018), an effect demonstrated by numerous field experiments that use snow fences to warm soils and thaw permafrost (Walker et al 1999, Hinkel and Hurd 2006, Natali et al 2011, O'Neill and Burn 2017).In the Arctic, winter temperatures typically drive mean annual soil temperature (Fedorov-Davydov et al 2018, Kropp et al 2020), and winter carbon dioxide (CO 2 ) emissions constitute a substantial component of the regional carbon budget (Zimov et al 1996, 1999, Natali et al 2019).In addition to enhanced carbon emissions, warming associated with increased snow cover can initiate irreversible permafrost thaw and thermokarst formation, as well as alter seasonal hydrology through enhanced runoff and soil moisture inputs (Turetsky et al 2020).
Northeastern Siberia is underlain by ice-rich, carbon-dense permafrost (Strauss et al 2017), meaning that thaw is typically abrupt, irreversible, and accompanied by high carbon emissions (Nitzbon et al 2020, Turetsky et al 2020).Snowfall has increased in northeastern Siberia over the past four decades (Pulliainen et al 2020), and these increases are projected to continue, though future predictions are poorly constrained (Bintanja and Selten 2014).Indeed, recent years with anomalously high snowfall in the region have been associated with increases in the depth of the seasonally thawed active layer and increased atmospheric CO 2 and methane (CH 4 ) concentrations (Anisimov and Zimov 2020).However, the impacts of increased snowfall on permafrost varied with a range of factors including vegetation cover, highlighting the importance of vegetationsnow interactions.Research investigating the influence of vegetation on snow accumulation typically focuses on differences between key vegetation types (Sturm et al 2001, Lundquist et al 2013, Grünberg et al 2020, Douglas and Zhang 2021).Northeastern Siberia is dominated by sparse larch forests, and while research has examined variability in snow accumulation across gradients in temperate forest cover (Sun et al 2022), there is a relative paucity of information on how heterogeneity in boreal tree composition and density influences snowpack dynamics (Smith et al 2022).
Interactions between vegetation and snow cover affect soil temperature, though in the Arctic, most studies documenting vegetation-snow-soil feedbacks focus on soil warming due to snow-trapping by shrubs (Myers-Smith andHik 2013, Frost et al 2018).While this interaction typically occurs at the patch scale (i.e.tens of meters), similar effects have been observed at the ecosystem or landscape scale (Paradis et al 2016, Lafleur and Humphreys 2018, Grünberg et al 2020).In line with these observations are several data synthesis and modeling studies indicating positive associations between winter soil temperature and vegetation height across Arctic and boreal regions (Zhang et al 2018, Kropp et al 2020, Way and Lapalme 2021).However, conventional knowledge suggests that interception by increased forest cover will reduce snow accumulation and also soil temperature.This holds true in studies where forest patches are sufficiently large to preclude snow-trapping effects (Jean and Payette 2014, Roy-Léveillée et al 2014).These inconsistencies are undoubtedly attributable in part to differences in the scale of observation and analysis; patterns revealed by large synthesis data sets are driven largely by climatic gradients, whereas site-level analyses reflect interactions between process drivers and local ecological and topographic heterogeneity.
Despite the clear importance of snow-vegetationsoil process interactions in permafrost ecosystems, research examining winter soil temperature variability across forest cover gradients is limited (Smith et al 2022).In this study, we used data collected from six study sites situated across a gradient of forest cover in northeastern Siberia to quantify variation in winter soil temperature dynamics in relation to canopy cover.Canopy cover at the sites ranges from 13% to >70%; permafrost active layer depth and understory tall shrub cover decrease with increasing canopy cover, while understory moss cover increases with canopy cover (Paulson et al 2021, Bendavid et al 2023).Data were collected over seven years, during which there was considerable variation in snow accumulation, allowing us to investigate the influence of snow cover variation on winter soil temperature as well.

Study site
The study site lies within the perimeter of a wildfire that occurred ca.1940 near the settlement of Cherskiy, Sakha Republic, Russia (68.75 • N, 161.48 • E; figure 1).Fire is common in the region (Talucci et al 2022), and variation in canopy cover is characteristic of early-to mid-successional forests (Alexander et al 2012).Cajander larch (Larix cajanderi Mayr.), a deciduous conifer, is monodominant in the region, and mean stand age at the study site is ca.60 years.Within the fire perimeter, canopy cover and stem density ranged from low-(<25% and 83 stems ha −1 , respectively) to high-density (>75% and 42 000 stems ha −1 , respectively) (Paulson et al 2021).In low density stands, the understory was dominated by erect deciduous shrubs (e.g.Betula nana subsp.Exilis and Salix spp.), while mosses (e.g.Sanionia spp.) dominated high-density stands (Paulson et al 2021).The region is underlain by continuous permafrost, and soils are characterized by a well-defined soil organic layer (SOL) composed of undecomposed moss, plant litter, and fine roots, that sits atop a silt-clay mineral soil layer.For the period 1980-2014 mean January, July, and annual temperatures were −32.7 • C, 13.2 • C, and −10.

Data collection
In July 2014, understory micrometeorological instrumentation was installed at six sites within the fire perimeter that spanned a range of tree density and canopy cover.At each site, a pyranometer (mounted 1 m aboveground) was used to measure photosynthetically active radiation (QSO-S, Decagon Devices, Pullman Washington), air temperature and relative humidity (either VP-3 or VP-4, Decagon Devices, Pullman Washington).Two integrated soil moisture, temperature, and conductivity sensors (GS-3, Decagon Devices, Pullman Washington) were installed to monitor soil microclimate conditions.At each site, one sensor was installed at the interface between the organic and mineral soil horizons (8-18 cm depth), and the second sensor was installed 5 cm below this.Initial data analysis revealed close agreement between both sensors at each site, so in the summer of 2015 the second sensor was repositioned to 10 cm below the organic-mineral soil interface.Data was logged hourly to EM-50 data loggers (Decagon Devices, Pullman Washington).Annual data downloads occurred during summer field campaigns, and there are occasional gaps due to power failures and wildlife disruption.Data are archived at the Arctic Data Center (Loranty and Alexander 2023).For this study, we analyzed data from the upper sensor located at the organic-mineral horizon interface.To help corroborate patterns observed from these sensors, we also include soil temperature data at a depth of 50 cm from the two endmember sites with the highest and lowest canopy cover collected using the same types of sensors (GS-3) and data loggers (EM-50) as part of a separate study at the sites (Kropp et al 2019).
To quantify interannual variability in snow cover, we used snow depth data recorded at the weather station in Cherskiy (WMO station ID 25123, downloaded from https://rp5.ru/Weather_archive_in_Cherskiy).Note: this provided a single data set that we assumed to be representative of conditions at all study sites, which were located ∼5 km from the station.Micrometeorological data were quality screened before submission to the Arctic Data Center, which is documented in the code archived with the data.Snow depth data were manually quality screened, and corrections are documented in the code cited below.
Vegetation and soil characteristics at 25 stands within the fire perimeter, including those instrumented for micrometeorology, were measured during summer field seasons from 2010 to 2017.At each stand, a series of measurements were taken along three parallel 30 m long belt transects spaced ∼30 m apart.Data were averaged to the stand-level, and variables used in this study include canopy cover measured using a hemispherical densiometer and SOL depth.

Analyses
Using data from the six micrometeorological stations, we calculated a series of metrics to characterize variation in soil temperature conditions and associated microclimatic drivers between sites and across years.In order to focus on the winter period, we analyzed data for each water year (e.g. 1 October 1-September 31).To characterize winter soil conditions, we calculated freezing degree days (FDD), or the sum of mean daily temperatures below zero for each water year, for soil (FDD soil ) and air temperatures (FDD air ).Similarly, in order to explore whether soil or air temperature conditions during the previous growing season impacted winter soil temperatures, we calculated thawing degree days (TDD) for the soil (TDD soil ) and air (TDD air ) as the sum mean daily temperatures above freezing during the water year.We used the mean of the top 5th percentile of daily snow depth measurements (SD 5 ) for each water year to quantify interannual variability in snow cover.We chose this as a metric of maximum snow depth to avoid spurious effects of a single maximum value, and also because data gaps prevented the use of seasonally integrated metrics.Lastly, we calculated freezing n-factors (n f ) as the quotient of FDD soil and FDD air in order to characterize differences in the winter energy balance between sites and across years.We also note that, unlike FDD soil , n f values are directly comparable across sites and can be compared with values from other studies.Our dataset included full records for the 2015-2021 water years, or 42 site-years of data, however data gaps required removal of six additional site years.We removed water years that were missing more than 20 d of data, and note that in some cases this may have been limited to specific sensors (e.g.only air or soil temperature data).Summaries of all site years of data can be found in table S1.
To determine the drivers of winter soil temperatures, we constructed a series of ordinary least squares linear models in R (R Core Team 2023).We began by comparing FDD soil with a series of independent variables that included FDD air , TDD soil and TDD air from the previous growing season (i.e.previous water year), SD 5 , and canopy cover.After this, we constructed multivariate linear models, retaining only those variables that were significantly related to FDD soil .In addition, we developed models with n f that excluded FDD air because it is used to calculate n f .We used Akike's Information Criterion corrected for small sample size (AICc) to compare models.All code used in analyses and data that are not archived elsewhere can be accessed at https://github.com/mloranty/dg_winter.

Results
We observed a high degree of variability in microclimate between sites and across the study period (table 1, figure 2).Winter soil temperatures varied substantially between sites, with mean FDD soil across the study period ranging from 419 to 1537 degree days (table 1).At our two endmember sites, FDD soil at 50 cm depth for the 2017 water year was 159 and 955.Summer soil temperatures also exhibited considerable variability, with TDD soil ranging from 190 to 569 degree days (table 1).By contrast air temperatures were relatively less variable, with site averages of FDD air and TDD air over the study period ranging from 4408-4799 and 1106-1148 degree days, respectively (table 1).SD 5 (mean of the top 5th percentile of seasonal daily snow depth values) ranged from a low of 460 mm in 2021 to a high of 851 mm in 2018 (figure 2).Mean Annual Soil Temperature at the sites was strongly correlated with FDD soil (figure 3; p < 0.01; Adj.R 2 = 0.89).
Variation in FDD soil across sites and years was related to a combination of climatological and ecological factors (figure 4, table 2).Conditions during the previous growing season were least important: previous year TDD soil was not significantly related to FDD soil (p > 0.05), while previous year TDD air was inversely related to FDD soil but explained relatively little variance (p < 0.05; Adj R 2 = 0.13).Winter air temperature (FDD air ) was strongly positively associated with FDD soil (p < 0.001; Adj.R 2 = 0.44).Canopy cover was also significantly positively related to FDD soil (p < 0.001; Adj.R 2 = 0.42), while SD 5 was inversely related to FDD soil (p < 0.01; Adj.R 2 = 0.25).
Multiple regression models showed that Canopy Cover and SD 5 (snow depth) explain the majority of variability in FDD soil and n f , though including FDD air and previous year TDD air yielded slight improvements in model performance (table 3).Including combinations of significant variables from individual regression models explained from 72% to 84% of the variance in FDD soil .The best model included all four significant variables; FDD air , previous year TDD air , SD 5 , and Canopy Cover (Adj R 2 = 0.84, AICc = 362.3).Removing previous year TDD air only slightly reduced the explanatory power (Adj R 2 = 0.83, AICc = 472.4),while the simplest model that included only SD 5 and Canopy Cover still explained nearly three quarters of the variance in FDD soil (Adj R 2 = 0.72, AICc = 507.6).Variation in n f is shown in figure 5 and associated models are described in table 4. SD 5 and Canopy Cover explained a similar amount of variability in n f (Adj R 2 = 0.74, AICc = −77.8),and including previous year TDD air in the model provided a relatively small increase in model performance (Adj R 2 = 0.79, AICc = −91.7).

Discussion
We observed substantial variation in winter soil temperature, quantified using FDD soil and n f , across six e n f = freezing n-factor (FDD soil /FDD air ).
f TDD = Thawing degree days.g This site had FDD soil = 159 at 50 cm depth for the 2017 water year.h This site had FDD soil = 955 at 50 cm depth for the 2017 water year.larch stands in northeastern Siberia with canopy cover ranging from 13% to 71%.The majority of variation in winter soil temperature was explained by canopy cover and interannual variation in snow depth.Given the long winters at this site, mean annual soil temperature is closely related to FDD soil , indicating that annual soil temperatures are driven largely by winter temperatures.Differences in permafrost thaw depth previously reported for some of these sites (Kropp et al 2019) align with the soil temperature differences we describe here; dense stands with colder soils have shallower permafrost thaw.In the absence of detailed site-level snow depth measurements the process drivers of these temperature differences are unclear.Nonetheless these results have several important implications for understanding permafrost change in boreal forests.

Why are soils colder under dense larch canopies in winter?
One potential explanation is that denser canopies intercept more snow, leading to shallower less insulative snowpacks beneath the canopy.This is typically inferred via differences in measured snow depth between forests and adjacent open areas, attributed to canopy interception and subsequent sublimation or redistribution by wind (Hedstrom and Pomeroy 1998).The magnitude of forest induced snow depth reductions depends on a number of biological, physical, and topographic factors.Canopy     Additionally, Siberian larch branches are covered in epiphytic lichen which may enhance interception, and timelapse photography from our study sites illustrates appreciable canopy interception (figure 6).
The extent to which canopy interception reduces snow accumulation depends on a variety of factors (Pomeroy et al 1998).During large precipitation events the amount of intercepted snow may exceed canopy holding capacity, and unloading may occur.Therefore, canopy induced snow depth reductions are typically greater with multiple smaller precipitation events, where the canopy loading capacity is not exceeded and intercepted snow can be sublimated or redistributed by wind (Moeser et al 2016).Related, absence of wind, due to either weather or topographic position, can also limit the magnitude of snow depth reduction.Interannual weather variability and landscape position are therefore important mediators of canopy interception impacts on snow accumulation.
Meteorological and topographic factors also interact with vegetation distribution at the landscape scale to control patterns of snow redistribution.In the Arctic, the snow-trapping effects of shrub canopies are a prominent and well documented example of this process.Patches of shrubs in graminoid dominated tundra ecosystems have canopies taller than surrounding grasses that act as windbreaks, and snow typically accumulates up to the height of the shrub canopy (Sturm et al 2001, Roy-Léveillée et al 2014, Lafleur and Humphreys 2018).These patch-scale effects have been observed at the ecosystem-scale as well, where tall shrub tundra ecosystems exhibit greater snow holding capacity than adjacent low shrub and graminoid tundra ecosystems (Lafleur and Humphreys 2018).In the latter case, snow depth typically reached shrub canopy height (e.g.∼30-60 cm), as snow was not limiting.The effects of forest cover on snow distribution in boreal forests are less clear.For example, a study in northern Canada observed greater snow depths in small forest patches surrounded by shrub tundra, while larger forest patches limited accumulation of snow, and in shrub dominated areas with no trees, snow depth was positively correlated with shrub canopy height (Roy-Léveillée et al 2014).A similar study in northern Quebec found deeper snow on forested palsas compared to unforested palsas (peat mounds with frozen permanently frozen center, upwards of 150 m in diameter); however, canopy gaps on forested palsas had deeper snow than adjacent forested areas (Jean and Payette 2014).These complex patterns illustrate that the spatial arrangement of vegetation on the landscape is also an important determinant of snow accumulation.Landscape scale heterogeneity in Siberian boreal forests is strongly influenced by wildfire influences on post-fire regrowth dynamics (Alexander et al 2018, Talucci et al 2020); thus, changes in the fire regime that drive forest structural dynamics may influence snow depth distribution patterns.
Snow density is also affected by vegetation cover, and is related to the thermal properties of snow that can affect soil temperature.Denser snow conducts heat more effectively (Zhang 2005), and redistribution and drifting caused by wind can lead to formation of high-density wind slab, which is common in tundra and high Arctic ecosystems (Derksen et al 2014).The study by Roy-Léveillée et al that observed variation in snow depth related to tree and shrub height also found that snow density was inversely associated with canopy height (Roy-Léveillée et al 2014), likely due to increased presence of wind slab (Zhang 2005).Alternatively, depth hoar, snow characterized by large crystals that metamorphose due to steep temperature gradients within the snowpack, has low thermal conductivity and is also affected by vegetation.In particular, the presence of erect shrubs is known to enhance depth hoar formation, especially relative to areas dominated by graminoids and/or cryptogams (Paradis et al 2016).A recent study detailing a transect in eastern Canada from boreal to high-Arctic latitudes found that shrubs were a primary determinant of snow thermal conductivity via impacts on depth hoar formation (Royer et al 2021).One possible explanation for the winter soil temperature patterns we observed is that greater shrub abundance in low density stands at our study site (Paulson et al 2021, Bendavid et al 2023) could facilitate depth hoar formation.
An alternative explanation unrelated to snow cover is that trees may act as thermal bridges that conduct heat from the soil to the atmosphere during the winter.At our research sites, patterns in larch stem density, basal area, and aboveground biomass were positively correlated with canopy cover (Loranty et al 2018, Paulson et al 2021).A recent study documented that thermal bridging by shrub stems reduces winter soil temperatures in tundra ecosystems in northern Canada (Domine et al 2022).Furthermore, this study found that the pattern reverses with the onset of the growing season; once air temperatures exceed soil temperatures, shrubs conduct heat into the soils, leading to a warming effect.At our sites, this might explain why soils thaw at the same time regardless of canopy cover (e.g.figure 2).However, to our knowledge there have not been any studies documenting thermal bridging by trees in boreal forests.It is also possible the patterns in winter soil temperature we observed are caused by a combination of thermal bridging and canopy impacts on snow cover.

Implications for ongoing change in the Arctic
We have demonstrated a link between canopy cover, a vegetation characteristic that can be measured remotely across large areas, and soil temperature metrics that are useful for quantifying important climate feedbacks.This may provide new opportunities for understanding landscape scale variability in winter CO 2 emissions and vulnerability to permafrost thaw.Our results imply that permafrost soils in Siberian larch forests with denser canopies may be more resilient to ongoing winter warming and/or increases in snow depth than those with sparse canopies.With colder soils, denser forests may also have lower winter CO 2 emissions.However, we do note that while winter fluxes decline with soil temperature, they also increase with canopy cover and leaf area index (Natali et al 2019).
Our results may also be useful for anticipating how ongoing vegetation change may mediate climatic influences on soil and permafrost temperatures.For example, if spectral greening observed in the region (Berner and Goetz 2022) is associated with forest infilling or increases in canopy cover, then the rate of soil warming with climate may slow as the cooling effect of forest canopies becomes more prevalent.
Similarly, recent large fire years in the region (Talucci et al 2022) coupled with increased larch recruitment with fire severity (Alexander et al 2018) may result in denser forest regrowth across the region, which in turn, could accelerate post-fire soil temperature recovery relative to areas with low-density regrowth.

Conclusions
We investigated soil temperature dynamics across canopy cover gradient at six larch stands in northeastern Siberia over seven years with substantial interannual variability in snow depth.We observed a roughly three-fold increase in soil FDD and freezing n-factors and our results demonstrate that nearly three quarters of this variability is explained by forest canopy cover and interannual variation in snow depth.The processes by which canopy cover variability influence soil temperature remain unclear, highlighting a need for additional research.Our results indicate that spatial heterogeneity in forest cover may serve as a useful proxy for winter and annual soil temperature variability in this region.Ongoing changes in forest cover, associated with climate warming and fire regime changes, as we all as changes in snow cover will interact to control the coupling between air and soil temperature (i.e.n-factor) in the future.

Figure 1 .
Figure 1.Map of the study site showing stand locations with percent canopy cover denoted (A).Representative photographs showing high (B), medium (C), and low (D) density stands.Oblique aerial photographs of high (E) and low (F) density stands.Photos M. Loranty.
Additional variables are described in detail in several previous publications (Loranty et al 2018, Paulson et al 2021, Hewitt et al 2022, Bendavid et al 2023), and the full dataset is archived at the Arctic Data Center.

Figure 2 .
Figure 2. Timeseries of soil temperature at the organic-mineral interface, and snow cover for six larch forests of varying density located within a single fire perimeter in northeastern Siberia.Soil temperature line colors correspond to stand canopy cover, and snow depth is shown in blue.

Figure 3 .
Figure 3. Relationship between mean annual soil temperature (MAST) and soil freezing degree days (FDD soil ) across a forest density gradient in northeastern Siberia.Each point represents a single water year for an individual site, and color represents percent canopy cover.The intercept and slope values are 1.47 and −0.003, respectively (p < 0.001, Adj.R 2 = 0.89).

Figure 4 .
Figure 4. Individual linear regressions between FDD soil and previous growing season TDD soil (a), previous growing season TDD air (b), FDD air (c), snow depth (d), and canopy cover (e).Regression lines are shown for significant relationships (p < 0.05), and the associated regression statistics are found in table 2. Color of individual points corresponds to canopy cover for each site.

Figure 5 .
Figure 5. Scatterplots showing relationship between n f and canopy cover (a), and snow depth (b).

Figure 6 .
Figure 6.Snow intercepted by larch canopies in northeastern Siberia.The top photo shows a high-density stand located near our study sites on 10 December 2016, and the bottom shows a low-density stand that was included in this study on 17 November 2016.Power failures and lens obstructions (e.g.condensation and snow) prevented inclusion of concurrent photos.Photo credit H Kropp.

Table 1 .
Summary of canopy cover, soil organic layer thickness, and micrometeorological variability for the six forested larch sites spanning a density gradient in northeastern Siberia.
a Soil organic layer.b Mean annual air temperature.c Mean annual soil temperature.d FDD = Freezing degree days.

Table 2 .
Individual regression summary for soil freezing degree days (FDD soil ).

Table 3 .
Multiple regression summary for soil freezing degree days (FDD soil ).

Table 4 .
Model fit summaries for multiple regression models for freezing n-factor (n f ).ModelDF Adjusted R 2 AICc n f ∼ TDD air + snow + CC