Does tall vegetation warm or cool the ground surface? Constraining the ground thermal impacts of upright vegetation in northern environments

Increased upright vegetation growth (i.e. trees and shrubs) in northern environments can profoundly impact ground surface thermal conditions through winter warming (e.g. enhanced snow trapping) and summer cooling (e.g. increased shading). The debate over these opposite effects emphasizes the need to better constrain net temperature impacts of upright vegetation on soils in northern environments. We generate a series of simulations with a widely-used permafrost model to partition the absolute warming and cooling impacts of upright vegetation on ground surface temperatures for a variety of shading scenarios, climates and surficial materials types (i.e. bedrock, mineral and organic soils). These scenarios simulate annual temperature differences between the air and ground surface caused by upright vegetation to provide likely ranges for the net effects induced by vegetation. These simulations showed that ground surface temperature warming in the winter mostly overwhelmed ground surface cooling in the thawing season even when simulations included extreme shading effects. Constraining the simulations to current best estimates of the possible summer cooling impact of vegetation yielded a dominant winter warming signal for most snow depths and climate types. Differences in the magnitude of air-surface temperature offsets between sites underlain by bedrock, mineral and organic soil highlights the importance of considering differences in unfrozen moisture content in areas where the ground freezes and thaws seasonally. The results of this study suggest that the net ground surface temperature impacts of increased snow trapping by vegetation will far exceed cooling caused by enhanced shading following increases in tall vegetation in most northern environments.


Introduction
Satellite remote sensing shows that high-latitude landscapes are greening in response to climate warming and Arctic amplification processes , Fraser et al 2011, Ju and Masek 2016, Vickers et al 2016, Arndt et al 2019. Warmer summer air temperatures are changing northern vegetation (Post et al 2009, Olthof andPouliot 2010) by increasing the abundance and density of shrubs in mountain and tundra environments (Tape et al 2006) and increasing ecosystem net primary productivity (Euskirchen et al 2009). Fieldbased empirical studies have shown increased shrub abundance in tundra vegetation (Myers-Smith et al 2019), with plot-scale data demonstrating biomewide trends of increasing plant canopy height, abundance of evergreen, prostrate and upright shrubs, and a decrease in bare ground (Elmendorf et al 2012). The combination of increased shrub biomass, abundance and cover is colloquially termed shrubification (Myers-Smith et al 2011). Shrubification can influence biophysical environments by modifying snow depth and other associated hydrologic elements, surface albedo (Myers-Smith et al 2011), northern wildlife (Rickbeil et al 2018), the diversity and abundance of understory species, ground temperatures (Elmendorf et al 2012) and permafrost (Wilcox et al 2019). Land cover changes can markedly impact ground thermal conditions (Myers-Smith and Hik 2013) through various interrelated factors by: (a) changing surface organic cover (Shur and Jorgenson 2007, Jorgenson et al 2010, Smith and Riseborough 2010a; (b) altering snow cover and associated processes (Sturm et al 2001b, Lawrence and Swenson 2011, Wang et al 2019; and (c) changing maximum thaw depths (Fisher et al 2016, Wilcox et al 2019. Whether shrubification in mountain and tundra environments will result in a positive or negative ground temperature feedback remains somewhat contested in the interdisciplinary literature (Myers-Smith et al 2011, Loranty et al 2018. Empirical studies have largely shown that increased shrub cover results in warmer winter ground temperatures (Sturm et al 2001b, Myers-Smith and  shrub removal has been shown to advance ground thaw (Kade and Walker 2008, Blok et al 2010, Nauta et al 2015. The net impact of shrub encroachment will likely be heterogeneous with the combined influence of peat, organic layer thickness and vegetation on soil temperatures possibly modifying both climate warming effects (Yi et al 2007) and maximum summer thaw depths (Walker et al 2003). The opposing impacts of freezing season warming (e.g. due to snow) and thawing season cooling (e.g. due to shading) emphasize the need to better constrain the total potential ground surface temperature effect of increasing vegetation in mountainous, Subarctic and Arctic environments.
This study uses ground thermal modelling to investigate the theoretical range of ground surface temperature responses to differences in snow thickness and vegetation shading caused by increases in upright vegetation. Here, an increase in upright vegetation refers to changes in the height or biomass of shrubs and trees due to greening and shrubification processes. The process-based modelling we employ enables quantification of theoretical constraints on the absolute magnitude of air-surface temperature differences because of upright vegetation during the thawing (warm) and freezing (cold) seasons. This modelling is guided by outputs from one-dimensional ground thermal modelling and is parameterized by a meta-analysis of the permafrost and ecology literature which enables us to evaluate thermal impacts across a range of soil types and climates.

Methods
Approaches to modelling ground surface temperatures range in complexity from fully coupled surface energy balance solutions (process-based) to localscale statistical fits of field data (empirical modelling) . Process-based models are preferred in permafrost modelling studies because empirical models do not explicitly consider physical heat transfer theory and ground thermal dynamics (Etzelmüller 2013). Process-based ground thermal models use either analytical (intermediate complexity) or numerical (high complexity) frameworks with the latter being more computationally onerous due to fewer implicit assumptions and thus more complex inputs (Etzelmüller 2013, Stieglitz et al 2003, Riseborough et al 2012. Nevertheless, intermediate-complexity analytical models generally consider physical processes like numerical models while allowing efficient calculations and parameterization from field data . Such analytical models have become widespread for applications ranging from site specific modelling of the ground thermal regime (e.g. Wright et al 2003, Bevington and Lewkowicz 2015, Ferreira et al 2017 to global mapping of permafrost distribution (e.g. Obu et al 2019). In this study, we combine numerical modelling outputs with empirical data from field studies to partition theoretical ground surface temperature impacts of upright vegetation growth (figure 1).
We use the surface offset (SO) as presented in the temperature at the top of permafrost (TTOP) model (Romanovsky and Osterkamp 1995, Smith and Riseborough 1996, Riseborough 2004 to quantify annual temperature differences between the air and ground surface. Application of a TTOP framework enables partitioning of the SO into components summarizing air-surface temperature (t) differences during the thawing (t > 0 • C) and freezing (t < 0 • C) seasons, including implied warming impacts from snow and cooling impacts from canopy shading. Our estimate of the potential cooling impact of vegetation during the thawing season (e.g. shading and organic mat related) on ground surface temperatures employs simple heat transfer theory constrained by local climate, and literature-informed (vegetative) cooling scenarios. Simulations of air-surface temperature differences caused by snow during the freezing season are derived from one-dimensional energy balance simulations previously generated for a range of index soil and climate conditions found in Canada (Riseborough 2004).

Theoretical framework
The TTOP model is a widely used analytical framework that has been validated against heat conduction models (Romanovsky and Osterkamp 1995, Smith and Riseborough 1996, Riseborough 2004) and empirical field data (Wright et al 2003, Juliussen and  and the empirical model-based thawing season offset (TSO) estimations to elucidate the total annual air-surface temperature difference (i.e. the surface offset (SO)). Humlum 2007, Gisnås et al 2013, Bevington and Lewkowicz 2015, Obu et al 2019. TTOP provides an equilibrium estimate of the mean annual ground temperature at the top of permafrost or at the base of the freeze-thaw layer in seasonally frozen ground (Smith and Riseborough 2002). The typical equation for TTOP and the variant used in this study is expressed in appendix 1. TTOP uses transfer functions known as n-factors (Lunardini 1978) such as the thawing n-factor (i.e. ratio of thawing degree days at ground surface to thawing degree days in the air; TDDs/TDDa) to convert accumulated degree days ( • C) in the air to corresponding values at the ground surface. These transfer functions account for the cumulative impacts of snow, vegetation cover and structure, and moisture on the ground surface (Smith and Riseborough 2002). The mean ground temperature at TTOP can be calculated using frozen and thawed thermal conductivities for the ground, or via a thermal offset summarizing the annual temperature influence of soil characteristics in the active layer or freeze-thaw layer (Burn and Smith 1988, Romanovsky andOsterkamp 1995).
Following a TTOP framework, the mean annual ground surface temperature (MAGST) can be expressed as equations (1) or (2). Combining equations (1) and (2) gives equation (3) for the SO which determines the mean annual temperature difference between the air and the ground surface (Smith and Riseborough 2002).
where: MAGST = Mean annual ground surface temperature ( • C) TDDs = Cumulative thawing degree days recorded at the ground surface ( • C days) FDDs = Cumulative freezing degree days recorded at the ground surface ( • C days) P = Period (typically, 365 days) where: MAAT = Mean annual air temperature ( • C) or (TDDa − FDDa) ÷ P SO = Mean annual surface offset ( • C) TDDa = Cumulative thawing degree days recorded in the air ( • C days) FDDa = Cumulative freezing degree days recorded in the air ( • C days) The SO (equation (3)) can be further partitioned into two components (equation (4)) reflecting differences between the air and ground surface during the freezing (t < 0 • C) (FDDa, FDDs) and thawing (t > 0 • C) (TDDa, TDDs) seasons. Functionally, these two separate components correspond to the nival offset (NVO; • C) (equation (5)) discussed by Smith and Riseborough (2002) and the thawing season offset (TSO; • C) (equation (6)) described by To quantify air-surface temperature differences, many prior permafrost studies have used freezing (nf) and thawing (nt) n-factors and the original TTOP formulation (appendix 1) (e.g. Karunaratne and Burn 2003, Wright et al 2003, Juliussen and Humlum 2007, Bevington and Lewkowicz 2015, Ferreira et al 2017, Obu et al 2019 rather than the NVO/TSO formulation given above. However, numerical model experiments performed by Riseborough (2004) suggest that NVO may be more intuitively linked to air temperature than the freezing n-factor. The NVO/TSO formulation also provides a direct temperature effect estimate as opposed to n-factors which require climate inputs to quantify absolute temperature impacts. For these reasons, we use the NVO and TSO to elucidate freezing season (e.g. snow) and thawing season (e.g. shading & organic mat) air-surface temperature differences in modelling scenarios (figure 1).

Simulating air-surface temperature differences in the thawing season
Thawing season air-surface temperature differences (TSO) were modelled by adjusting cumulative thawing degree days ( • C days) in the air (TDDa) with five scenarios meant to represent the net impacts of different magnitudes of vegetative cooling (e.g. due to differing vegetation canopy or structural characteristics) at the ground surface (figure 2(a)). TDDa were derived using a sine curve fit to mean annual air temperature and annual temperature amplitude following Riseborough et al (2012) (supplemental S1 (available online at stacks.iop.org/ERL/16/054077/ mmedia)). Simulated TDDa used MAATs ranging from −16 • C to +6 • C (intervals of 0.5 • C) and assumed a temperature amplitude of 20 • C (Riseborough 2004). The minimum and maximum TDDa for a given MAAT, when used to calculate TSO (equation (6)) , provide theoretical limits for the net air-surface temperature differences that could arise from different vegetation characteristics in the thawing season ( figure 1(b)). Prior empirical estimates of nt for upright vegetation show that the values are usually between 0.5 and 1.0, with a best estimate around 0.75 (e.g. partial cooling) (table 1) (e.g. Kropp et al 2021). This partial cooling (nt = 0.75) scenario was included in addition to two more extreme cooling scenarios where only a quarter of TDDa (nt = 0.25; large cooling) and half the TDDa (nt = 0.5; half cooling) reach the ground surface (figures 2(a) and (b)). These two more unlikely (extreme) scenarios help us explore the limits of possible temperature effects of vegetation in the thawing season.

Simulating air-surface temperature differences in the freezing season
Snow cover impacts on ground surface temperatures depend largely on snow characteristics including depth, density and duration (Sturm et al 2001a, Johansson et al 2013, Gisnås et al 2014, local climate conditions and the soil layers' characteristics (Osterkamp and Romanovsky 1998, Romanovsky and Osterkamp 2000, Riseborough 2004, Throop et al 2012. TTOPbased studies typically have related freezing n-factors (nf) to snow depths or ecosystem types using field data (Ménard et al 1998, Juliussen and Humlum 2007, Karunaratne et al 2008, Davesne et al 2017, Gisnås et al 2016 or one-dimensional model simulations stratified by climate types (Henry and Smith 2001, Smith and Riseborough 2002, Way and Lewkowicz 2016, Obu et al 2019. A critical challenge to both approaches is the role of unfrozen moisture and latent heat in the active layer which modifies the thermal impact of snow for different surficial types (Romanovsky and Osterkamp 2000, Riseborough 2002, 2004, Zhang 2005, Throop et al 2012. In our model scenarios, we determined NVO for three standard surficial types (bedrock, mineral soil, and organic soil) meant to broadly represent a range of moisture/soil thermal characteristics in northern environments using experiments originally performed by Riseborough (2004) (figures 2(c) and Figure 2. (a) Simulated thawing degree days ( • C days; TDDa) in the air ( • C days; TDDa) and at the ground surface used for five thawing season cooling scenarios. TDDa are presented as thawing n-factor (nt) scenarios, where nt = 1.0 is equivalent to the thawing degree days in the air (no cooling). Thawing n-factors and thawing degree days at the ground surface are stratified by mean annual air temperature ( • C); (b) simulated thawing season offsets ( • C; TSO) calculated using equation (6), the thawing degree days in the air and ground surface for warm-season cooling scenarios presented in (a); interpolated nival offset ( • C; NVO) surfaces derived from Riseborough (2004) TONE simulations, with thermal and moisture properties for (c) bedrock and (d) mineral soil. Surfaces were generated using MAATs ranging from −18 • C to +6 • C and snow depths from 0 to 1.5 m. Contours for NVO use a 1 • C interval. Surface plots were generated using the fields package (Nychka et al 2020) in R.
(d)). Each surficial type corresponds to a different moisture regime ranging from high (i.e. organic soil), to moderate (i.e. mineral soil) to negligible (i.e. bedrock) moisture contents. These simulations were originally prepared with the one-dimensional finite element thermal conduction model TONE (Goodrich 1978(Goodrich , 1982 which has been thoroughly validated for use in northern environments with empirical data, soil analytical solutions and surface energy balance modelling (Zhang et al 1996, 2003, Zhang and Stamnes 1998, Oelke et al 2003, Riseborough 2004, Smith and Riseborough 2010b, O'Neill and Burn 2017. The original TONE simulations were performed using MAATs ranging from −18 • C to +6 • C and snow depths from 0 to 1.5 m (see supplemental S2 for model details). These simulations were not available in tabular form but were precisely digitized from Riseborough (2004) and then interpolated using a thin plate spline from the fields package (Nychka et al 2020) in R (figures 2(c) and (d); see supplemental S3 for full graphics).
Estimates of NVO were generated for bedrock, mineral, and organic soil for the same MAATs used in the TSO experiments (−16 • C to +6 • C, 0.5 • C intervals) and for late-winter snow depths ranging from 0 to 1.5 m (0.1 m intervals). Snow cover duration varied depending on the air temperature (e.g. snow cover duration increased with cooler MAATs) and numerical simulations used a mean snow density of 250 kg m −3 ; therefore, applying this method across different snow densities (e.g. ecosystem types or climate classes) would require normalization to equivalent snow depth (e.g. Riseborough 2004, Way andLewkowicz 2016). For our analysis, we focused on differences between bedrock and mineral soil. Organic soil simulations are only presented in the Online Supplement because of similarities between results for mineral and organic soils, though divergence does occur at MAATs colder than −8 • C (supplemental S4).

Model simulations
The results in this section are presented in the form of SOs (i.e. SO or NVO + TSO) which sum the model-based estimates for the warming effect of snow cover (NVO) and the scenario-based cooling effects attributed to vegetation (TSO) (e.g. figure 1). Simulated SOs show the mean annual air-surface temperature difference, with positive values when warming exceeds cooling (NVO > TSO) and negative Table 1. Summary of prior studies providing empirical estimates of the thawing n-factor (nt). Details are provided for study authors and locations, with the ranges of nt values measured and nt observations for upright vegetation presented. values when cooling exceeds warming (TSO > NVO).

Region of study
In total, we derived NVOs for 45 different MAATs and 16 different late-winter snow depths yielding a total of 720 NVOs for each of the three soil types (n = 2160 NVOs total). Adding the NVO estimates to TSO values from the five vegetation cooling scenarios (e.g. figure 2(b)) resulted in a total of 10 800 estimates of SO ( figure 3). Individual results for all cooling scenarios and surficial types are presented in the supplemental documentation (supplemental S5-S9).

SO simulations for bedrock, mineral and organic soil
SO simulations for bedrock, mineral and organic soil produced median estimates above 0 • C for each of the cooling scenarios meant to represent different ranges of vegetation coverage (figure 3). Mean SOs in all scenarios were also greater than 0 • C and in most scenarios were statistically different across surficial material types (supplemental S10). SOs were closer to 0 • C for bedrock under each cooling scenario when compared to mineral and organic soil. SOs were also always below +10 • C in bedrock simulations, whereas mineral and organic soils produced SOs exceeding +10 • C in many cooling scenarios. The large cooling (nt = 0.25), half-cooling (nt = 0.5), and partial cooling (nt = 0.75) scenarios produced nearly the same number of SO simulations above 10 • C as those below 0 • C for both mineral and organic soils. SOs for the maximum cooling scenario exhibited the largest range of any scenario, though bedrock (∼−10 to +8 • C) showed a smaller range in SOs than mineral (∼−10 to +13 • C) and organic soil (∼−10 to +12 • C). The no cooling scenario (nt = 1.0) resulted in the smallest range in SOs, with bedrock (∼0 to +10 • C) again having a smaller range than mineral and organic soils (∼0 to +14 • C). The partial cooling scenario (i.e. best-estimate of thawing season cooling effects of upright vegetation; nt = 0.75; table 1) produced only a limited number of simulations that caused an overall cooling impact (i.e. negative SOs) ( figure 3). Generally, decreasing the thawing season cooling effect (i.e. increasing nt) resulted in three main outcomes for all surficial material types: (a) an increased number of simulations with positive SOs; (b) a decreased total range of SOs, albeit less pronounced for bedrock; and (c) more simulations with positive SOs for mineral and organic soils than bedrock.

SOs, late-winter snow depth and mean annual air temperature
Simulated SOs for each scenario, as a function of snow depth and mean annual air temperature, are presented in supplemental S5-S9 with figures 4 and 5 depicting only the maximum (nt = 0.0) and partial (nt = 0.75) cooling scenarios for bedrock and mineral soil. As a simplification, MAAT is discussed here in the context of typical air temperature ranges used for permafrost zonation: (a) continuous permafrost associated with MAATs <−8 • C; (b) discontinuous permafrost represented by MAATs between −8 • C and 0 • C; and (c) mostly seasonal ground freezing for MAATs >0 • C. SOs for the complete cooling scenario (nt = 0.0; figure 4) differed between bedrock and mineral soil, with the mineral soil simulations producing larger SOs than bedrock (i.e. reaching +12 • C and +8 • C, respectively) at deeper snow depths and colder MAATs. The largest positive SOs for bedrock (+6 • C to +8 • C) were produced at MAATs typically found in the continuous permafrost zone at late-winter snow depths exceeding ∼0.8 m. Comparatively, for mineral soil the largest SOs (+6 • C to +12 • C) were simulated at snow depths exceeding ∼0.5 m and MAATs less than −6 • C (i.e. within climates characteristic of continuous and discontinuous permafrost zones). Negative SOs were more frequently produced in the seasonally frozen ground zone for both surficial types but also at shallower snow depths (∼0.1 m to ∼ 0.6 m) at MAATs typical of the discontinuous permafrost zone (e.g. 0 • C to −8 • C). Bedrock, however, typically produced negative SOs in deeper snow depths than mineral soil. For example, at ∼−5 • C, negative SOs occurred at snow depths as high as ∼0.4 m for bedrock, while 0.1 m was sufficiently thick for warming to exceed cooling in mineral soil. In the climatic zone typical of discontinuous permafrost zone, TSO had a stronger impact on bedrock at shallower snow depths (e.g. between ∼0.2 m and 0.6 m); however, deeper snow depths only produced negative SOs at MAATs between ∼−3 • C and 0 • C.
The partial cooling scenario (i.e. nt = 0.75; figure 5) showed greater similarity in the SO patterns for bedrock and mineral soil. NVO was greater than the TSO at most late-winter snow depths and MAATs, with negative SOs constrained to a narrow range of conditions for both surficial types. The largest positive SOs for bedrock and mineral soil (+8 • C and +14 • C, respectively) occurred at deeper snow depths and MAATs typical for the continuous permafrost zone. Negative SOs were produced for both surficial types in climatic conditions characteristic of discontinuous permafrost, at late-winter snow depths below 0.1 m. At MAATs usually in the zone of seasonal ground freezing, negative SOs only occurred where snow depth was less than 0.2 m and 0.4 m for mineral soil and bedrock, respectively. However, where MAATs exceeded +5 • C, negative SOs occurred at deeper snow depths. Overall, negative SOs were generally simulated in warmer climates where snow depths were shallow (i.e. ∼<0.1 m and ∼<0.2 m in mineral soil and bedrock, respectively).

Implications of vegetation growth impacts on ground surface temperatures
Several implications can be drawn from the simulation results presented in this study. First, at MAATs usually associated with continuous permafrost, it is exceedingly difficult, even under the largest vegetation cover cooling scenarios (e.g. nt = 0.00, nt = 0.25), to produce MAGSTs that are appreciably colder than MAATs when the late-winter snow depth exceeds 50 cm. The modelled warming effect of thicker snow in the freezing season typically overwhelms cooling effects during the thawing season in most scenarios, including most simulations with no thawing season temperatures reaching the ground surface (nt = 0.0). This result aligns with a recent northern hemisphere-wide observational study that showed that cold-season effects on ground thermal regimes are larger than warm-season effects in permafrost environments (Kropp et al 2021). Second, the situations where SOs can be negative are limited to those with sizable vegetation-induced cooling, shallow snow cover and warmer MAATs. This agrees with a variety of empirical and modelling studies finding that local ground temperature variability in northern environments is mostly controlled by snow characteristics (  Further, our modelling scenarios likely underestimated the magnitude of shrubification-and upright vegetation-induced ground surface temperature warming by not including the impact that enhanced upright vegetation growth would have on snow accumulation or the thermal characteristics of the snowpack. A myriad of studies have demonstrated that taller vegetation is associated with greater snow accumulation and extensive drifting around dense patches (Nicholson and Granberg 1973, Nicholson 1979, Sturm et al 2001a, Marsh et al 2010, Palmer et al 2012, Lantz et al 2013, Myers-Smith and Hik 2013, Jean and Payette 2014, Paradis et al 2016, Frost et al 2018, Pelletier et al 2019. Further, snow density in areas of upright vegetation is typically lower than for tundra snow (Sturm et al 2001a, Loranty et al 2018, thus a transition to less dense snow would favor less winter heat escape and enhanced air-surface temperature differences in the freezing season. Finally, the simulations that we present highlight the challenges in deciphering short-term ground surface temperature data in the absence of a robust understanding of local snow conditions and geomorphological context. Snow cover can influence the ground thermal regime, independent of the near surface air temperatures (Stieglitz et al 2003) through a variety of factors such as timing, duration, accumulation, density and melting, in addition to interactions between snow cover and other variables such as vegetation and geographic location (Zhang 2005). Trend assessment of ground surface temperature data under climate change scenarios may thus be susceptible to bias as climate driven reductions to snow cover could diminish or even mask potential ground temperature warming that may otherwise be occurring (e.g. Zhang et al 2008).
Overall, our modelling scenarios suggest that even a small change in late-winter snow depth (e.g. 10 cm) can be sufficient to overcome a large vegetative cooling effect or alternatively a summer warming effect induced by climate change. Therefore, climate driven reductions and increases to snow cover have the potential to either dampen or amplify ground surface warming over differing timescales, respectively (Zhang 2005, Stieglitz et al 2003. These results highlight the need for enhanced deployment of low-cost methods for monitoring snow depth and snow thermal characteristics at remote field sites (e.g. Danby and Hik 2007, Lewkowicz 2008, Staub and Delaloye 2017. Interactions between soils, vegetation species, climate and the chemical composition of litter can also affect the rate of soil organic carbon sequestration and the quantity of carbon stock in soils (Lal 2005). The soil micrometeorology results from our study may have implications for soil carbon cycling because changes in soil temperature can:

Limitations
There are several drawbacks to the overall approach used in this study warranting discussion. Most importantly, we characterize total effects with a process-based analytical model that, although guided by numerical model simulations, was not explicitly run using a fully-coupled atmospheric-land surface energy balance scheme. The use of late-winter snow depth as a single summary of snow cover effects on the ground surface will not include the potential impacts of short duration Goodrich 1982, 1978 winter warming events or early-onset fall snow accumulation (e.g. Riseborough 2004, O'Neill and Burn 2017, Jan and Painter 2020. However, late-winter snow thickness is generally found to be amongst the strongest predictors of ground temperature variability in field and model-based studies (Nicholson 1979, Ménard et al 1998, Zhang 2005 Application of our NVO estimation method to simulations using a surface energy balance approach  showed that our SO simulations had a mean absolute difference of 0.9 • C ± 0.6 • C across 595 model years (supplemental S11).
Another limitation is related to trade-offs between the thermal conductivity and moisture contents between bedrock, mineral and organic soils on air-surface temperature differences. In our study, the larger modelled negative SOs for bedrock are likely due to lower unfrozen moisture contents and minimal latent heat release during active layer freeze up. However, both thermal conductivity and moisture contents can influence the ground thermal regime (e.g. Jusak et al 2016) with thermal conductivity varying across soil textures, vegetation covers and organic layers (Gray et al 1988). A lower thermal conductivity of organic soils compared to mineral soils could result in reduced vertical heat fluxes within the soil column and greater insulation of the ground from warmer summer air temperatures (O'Donnell et al 2009). In turn, unfrozen moisture may mitigate ground surface cooling in frozen ground, increase heat fluxes originating from the ground and strengthen the insulating effects of snow cover (Romanovsky and Osterkamp 2000). For example, snow can have a greater influence on air and ground temperature relations when moisture contents in the active-layer are high (Throop et al 2012). While soil properties are a main influence on the ground thermal regime, our method did not allow to disaggregate between the specific impacts of soil moisture content and thermal conductivity and thus could not be used to investigate the trade-offs between the two soil properties. Disaggregating across different soil properties may be important for future work as variations in vegetation structure and composition are the result of interactions between various parameters including, climate, hydrology, thermal conditions, disturbance factors and microbial conditions (Walvoord and Kurylyk 2016). For example, changes in moisture content and/or surface water conditions associated with permafrost could change vegetation composition by shifting from black spruce to bogs (Baltzer et al 2014).
Additionally, the analytical model framework we use (TTOP) is derived from a model that assumes equilibrium rather than transient conditions. Application of this model to real-world conditions can be problematic in ice rich permafrost environments undergoing rapid environmental change (Riseborough 2007. In general, colder MAATs will result in colder permafrost temperatures (Throop et al 2012), but various factors can make the relationship less direct. Permafrost temperatures are influenced by climate as a first-order control while microclimate is a control on ground surface temperatures, with vegetation, moisture contents, snowpack characteristics, topography and earth materials impacting temperatures at a local scale (Goodrich 1982, Judge 1973, Throop et al 2012. We mostly reconciled this issue by examining ground surface temperatures rather than deeper surficial layers, but the effects of near-surface temperatures on the active layer remain under parameterized (Riseborough 2004).
Finally, the approach we use to characterize vegetation-related cooling does not disaggregate the vegetation-specific impacts on the local surface energy balance. Aggregation could thus average out opposing effects making it difficult to isolate specific processes, though the overall results would be unchanged. For example, shrubs can reduce active layer thickness (i.e. cool the ground) during the snowfree season by shading while also advancing snowmelt timing when their branches protrude through the snow surface and subsequently allow for earlier active layer thaw (i.e. ground warming) in the spring (Wilcox et al 2019). Here, the method will not disaggregate between both opposing processes. Rather, the result will reflect the net yearly impact on the ground surface temperature regime. It should be noted that the range of thawing season cooling simulations we use is not representative of the real-world ranges observed in field studies (refer to table 1). We conservatively decided to exclude scenarios where nt = >1 (common in empirical studies) and included those where nt = <0.5 (uncommon in empirical studies) (see table 1) which biases our results towards negative SOs (e.g. TSO > NVO) and overestimates the cooling effect of upright vegetation in the thawing season. Nevertheless, these extreme scenarios can help provide constraints on the maximum effects that could be (but are typically not) observed in the field.

Conclusion
In this study, we presented a range of simulations that explore air-surface temperature differences in thawing and freezing seasons for a range of theoretical surficial types and vegetation cooling scenarios. Significant contrasts (related to unfrozen moisture content) were observed between the estimated temperature effects of snow thickness for different surficial types highlighting the importance of geomorphological context for understanding thermal impacts of vegetation and snow cover at the ground surface. Most simulations were unable to produce a thawing season vegetation cooling effect large enough to overcome the warming impact of snow cover; however, scenarios with near-surface bedrock had smaller differences between snow and vegetation effects. Under most scenarios, the impact of snow cover was dominant at late-winter snow depths greater than 20 cm even with substantial thawing season cooling assumptions. Using best estimates from the literature, our results overwhelmingly show net warming due to snow cover that exceeds cooling effects induced by vegetation during the thawing season. Scenarios where ground surface temperatures were cooler than air temperatures often required unrealistic assumptions of a near total decoupling of air and ground surface temperatures during the thawing season.
We thus hypothesize that an annual cooling impact of enhanced vegetation growth is most likely to be observed where late-winter snow depths are shallow with extensive summer shading, coarser (drier) near-surface soils, and warmer mean annual air temperatures. The insights gained from this modelling experiment agree with prior empirical studies suggesting that the overall impact of taller vegetation will largely be long-term ground temperature warming. The influence of winter snow thickness on ground surface temperatures remains a considerable source of uncertainty in characterizing northern environmental change and may mask subtle changes in local temperature and moisture regimes. These results further highlight the necessity of enhanced snow and soil characteristic modelling and associated field data collection at northern ecological field sites.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.  where: kt = Thermal conductivity of thawed ground (Wm −1 K −1 ) kf = Thermal conductivity of frozen ground (Wm −1 K −1 ) FDDa = Cumulative freezing degree days recorded in the air ( • C days) FDDs = Cumulative freezing degree days recorded at the ground surface ( • C days) MAGST = Mean annual ground surface temperature ( • C) MAAT = Mean annual air temperature ( • C) nf = Freezing n -factor (unitless) calculated as FDDs ÷ FDDa nt = Thawing n -factor (unitless) calculated as TDDs÷ TDDa NVO = Annual nival offset ( • C) P = Period (days) calculated typically as 365 SO = Annual surface offset ( • C) TDDa = Cumulative thawing degree days recorded in the air ( • C days) TDDs = Cumulative thawing degree days recorded at the ground surface ( • C days) TO = Thermal offset ( • C) TSO = Annual thawing season offset ( • C).