High Arctic Vegetation Communities With a Thick Moss Layer Slow Active Layer Thaw

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ High Arctic Vegetation Communities With a Thick Moss Layer Slow Active Layer Thaw © 2024. The Author(s) Published version Schuuring, Sil; Halvorsen, Rune; Bronken Eidesen, Pernille; Niittynen, Pekka; Kemppinen, Julia; Lang, Simone I. Schuuring, S., Halvorsen, R., Bronken Eidesen, P., Niittynen, P., Kemppinen, J., & Lang, S. I. (2024). High Arctic Vegetation Communities With a Thick Moss Layer Slow Active Layer Thaw. Journal of Geophysical Research: Biogeosciences, 129(8), Article e2023JG007880. https://doi.org/10.1029/2023jg007880 2024 High Arctic Vegetation Communities With a Thick Moss Layer Slow Active Layer Thaw Sil Schuuring 1,2 , Rune Halvorsen 2 , Pernille Bronken Eidesen 1,3 , Pekka Niittynen 4 , Julia Kemppinen 5 , and Simone I. Lang 1 1 Department of Arctic Biology, University Centre in Svalbard, Longyearbyen, Norway, 2 Natural History Museum, University of Oslo, Oslo, Norway, 3 Department of Biosciences, The Faculty of Mathematics and Natural Sciences, University of Oslo, Oslo, Norway, 4 Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland, 5 Geography Research Unit, University of Oulu, Oulu, Finland Abstract Svalbards permafrost is thawing as a direct consequence of climate change. In the Low Arctic, vegetation has been shown to slow down and reduce the active layer thaw, yet it is unknown whether this also applies to High Arctic regions like Svalbard where vegetation is smaller, sparser, and thus likely less able to insulate the soil. Therefore, it remains unknown which components of High Arctic vegetation impact active layer thaw and at which temporal scale this insulation could be effective. Such knowledge is necessary to predict and understand future changes in active layer in a changing Arctic. In this study we used frost tubes placed in study grids located in Svalbard with known vegetation composition, to monitor the progression of active layer thaw and analyze the relationship between vegetation composition, vegetation structure and snow conditions, and active layer thaw early in summer. We found that moss thickness, shrub and forb height, and vascular vegetation cover delayed soil thaw immediately after snow melt. These insulating effects attenuated as thaw progressed, until no effect on thaw depth was present after 8 weeks. High Arctic mosses are expected to decline due to climate change, which could lead to a loss in insulating capacity, potentially accelerating early summer active layer thaw. This may have important repercussions for a wide range of ecosystem functions such as plant phenology and decomposition processes. Plain Language Summary Temperatures are rising in the Arctic, causing increased thaw of the layer of soil located above the permanently frozen ground. In Low Arctic regions vegetation cools the soil, which reduces the thawing. So far, we do not know whether the small plants growing in the High Arctic may be able to slow or reduce thaw. We measured soil thaw throughout the summer in High Arctic Svalbard in locations where vegetation composition is known. We also measured thickness of the moss layer, height of plants and snow depth. We found that moss thickness was the strongest factor in insulating the soil. Also the cover of plants, height of shrubs and forbs, and height of grass‐like plants slowed soil thaw in the early summer. The insulating effects became less over time and no effects were found 8 weeks after onset of thaw. As climate change is causing changes in the Arctic vegetation, mosses and small shrubs are expected to decrease. As we found these to be the most important factors in insulating the soil, a future decrease in mosses and small shrubs may cause accelerated soil thaw at the start of summer. 1. Introduction Throughout the Arctic, climate change is leading to degradation of permafrost through deepening of the active layer and rise of permafrost temperature (Smith et al., 2022; Thoman et al., 2022). In western and central Siberia, the total active layer thickness (ALT) has increased with 1.3 cm/year, and similar rates have been found in Scandinavia (Smith et al., 2022). This thawing of the permafrost and the associated increase in ALT is expected to have a major effect both locally and globally through the release of greenhouse gases, damage to infrastructure and release of harmful chemicals (Hjort et al., 2022; Miner et al., 2021,2022). Global warming has also led to large‐scale changes in pan‐Arctic vegetation (Bjorkman et al., 2020). Shrubs and graminoids are generally expanding as a result of increasing temperatures, while the cover of mosses and lichens is decreasing (Elmendorf et al., 2012; Lang et al., 2012). Furthermore, northward expansion of shrubs and trees is being observed, along with an overall increase in NDVI (Myers‐Smith et al., 2020). RESEARCH ARTICLE 10.1029/2023JG007880 Key Points: •High Arctic vegetation slows active layer thaw in early summer after snow melt •Mosses show a stronger negative relation with thaw depth than vascular vegetation •Factors influencing active layer thaw change over time in early summer Supporting Information: Supporting Information may be found in the online version of this article. Correspondence to: S. Schuuring, [email protected] Citation: Schuuring, S., Halvorsen, R., Bronken Eidesen, P., Niittynen, P., Kemppinen, J., & Lang, S. I. (2024). High Arctic vegetation communities with a thick moss layer slow active layer thaw. Journal of Geophysical Research: Biogeosciences, 129, e2023JG007880. https://doi.org/10. 1029/2023JG007880 Received 3 NOV 2023 Accepted 13 JUL 2024 Author Contributions: Conceptualization: Rune Halvorsen, Pernille Bronken Eidesen, Simone I. Lang Data curation: Rune Halvorsen Formal analysis: Rune Halvorsen Funding acquisition: Simone I. Lang Investigation: Simone I. Lang Methodology: Rune Halvorsen, Pernille Bronken Eidesen, Simone I. Lang Supervision: Rune Halvorsen, Pernille Bronken Eidesen, Simone I. Lang Writing – review & editing: Rune Halvorsen, Pernille Bronken Eidesen, Pekka Niittynen, Julia Kemppinen, Simone I. Lang © 2024. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. SCHUURING ET AL. 1 of 12 There are indications that this climate‐driven increase in vegetation may mitigate permafrost degradation. Oehri et al. (2022) demonstrated that vegetation type is one of the main predictors of soil energy balance both in summer and on an annual basis. During summer, vegetation generally has a cooling effect on the soil by shading from incoming sunlight and creating an insulating layer (Blok et al., 2010; Jorgenson et al., 2010; Magnússon et al., 2020). While these effects are mainly demonstrated for shrubs, mosses may also play an important role in reducing thaw depth (TD). A thick mat of moss lowers thermal conductivity, reduces soil temperatures, and protects the permafrost against increasing mean air temperatures and temperature fluctuations (Blok et al., 2011; Gornall et al., 2007; Porada et al., 2016; Soudzilovskaia et al., 2013). These summer cooling effects are often balanced by winter warming. Higher and denser shrub vegetation traps more snow, insulating the soil against the colder air, and thus results in higher soil temperatures (Grünberg et al., 2020; Sturm et al., 2001). In spring, shrubs protrude above the snow and absorb solar radiation, which can lead to accelerated snow melt and earlier exposure of the soil (Wilcox et al., 2019). Abiotic factors are other drivers of active layer thaw. Snow has been demonstrated to influence active layer thaw in early spring. Late‐lying snow insulates the soil against higher air temperatures in spring and reflects most of the incoming radiation, which leads to delayed soil thaw (Callaghan et al., 2011). Soil moisture also plays an important role in permafrost dynamics. High soil moisture leads to slower refreezing in autumn, as more energy is required for the water to freeze (Romanovsky & Osterkamp, 2000). Therefore, increased soil moisture after extreme rainfall in summer delays autumn freezing (Magnússon et al., 2022). However, this does not result in slower thawing in spring as thaw occurs more abruptly, and a higher heat capacity and latent heat loss of the soil is counteracted by a higher soil thermal conductivity caused by infiltration of meltwater and rain (Hinkel & Outcalt, 1994; Loranty et al., 2018). The above‐mentioned studies of vegetation‐permafrost interactions have so far mostly focused on Low Arctic regions where the vegetation is largely dominated by shrubs and the ecosystem is considered to modify or protect the permafrost (Heijmans et al., 2022; Ran et al., 2021; Shur & Jorgenson, 2007). In High Arctic regions, where mosses dominate and vascular plants are smaller and occur more sparsely than in the Low Arctic, the relationship between vegetation and active layer thaw is less understood (Ran et al., 2021; Shur & Jorgenson, 2007). While Gornall et al. (2007) studied the relation between moss thickness and thaw depth in an experimental setup, observations over a natural gradient are still lacking. Accordingly, Heijmans et al. (2022) have emphasized the importance of effects of vegetation and the need for further studies in the High Arctic. In addition, most studies have concentrated on the effect on ALT at the end of the thawing season (e.g., Blok et al., 2010). Yet the onset and progression of soil thaw determines a wide range of ecosystem processes, from plant phenology to decomposition processes (Bring et al., 2016; Hobbie et al., 2000; Schaub et al., 2022). In the light of climate change and a changing vegetation, changes in timing of thaw may thus be crucial for ecosystem functioning, and shedding light on which factors are most important in insulating the ground during the early growing season, will help us understand future implications of climate change. As of now it remains unclear which factors are most important on a temporal scale from the onset of thaw and during the thawing season, and whether these factors change throughout this period. In this study we investigate factors influencing active layer thaw during early summer in High Arctic Svalbard, focusing on the role of vegetation and using fine‐scale study grids along topographical gradients to cover a wide range of natural variation in vegetation and abiotic properties. Based on existing knowledge, we expect negative relationships between the progression of active layer thaw and (a) the thickness of the moss layer; and (b) the cover and height of High Arctic vascular vegetation. We also hypothesize (c) that factors most strongly correlated with active layer thaw will change during the season. 2. Materials and Methods 2.1. Location The study was conducted in four study grids in Endalen and Adventdalen, located in the vicinity of Longyearbyen, Svalbard (78.20°N 15.63°E, Figure 1, Table S1 in Supporting Information S1). Mean annual temperature is −3.9°C, with an average July temperature of 7.0°C (1991–2020). Mean annual precipitation is 218 mm, of which 52 mm falls in summer (June, July and August) (1991–2020) (Norwegian Centre for Climate Services, 2023). Each grid had the dimension of 8 ×20 m, and consisted of 160 plots, with a plot size of 1 m ×1 m (Kemppinen, Niittynen, le Roux et al., 2021; Niittynen et al., 2020b). On a subset of 20 plots in each grid (80 plots in total), a Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 2 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License detailed vegetation survey was conducted (i.e., species level cover estimates of vascular plants, bryophytes, and lichens), and used in this study. The grid was divided into sectors based on vegetation types and environmental factors (e.g., mesotopography) and plots were randomly selected within all those sectors with the condition that plots were not adjacent to one another. Thus, the randomly selected plots would cover the grids maximally both in terms of space and vegetation types. A short description of each grid is given in Table S1 in Supporting Information S1. 2.2. Abiotic Measurements 2.2.1. Thaw Depth A large amount of belowground rocks made manual probing to measure active layer depth not possible. To overcome this, we installed frost tubes in each plot in the autumn of 2020. Installation of frost tubes also allowed for repeated measurements within the same plot without causing disturbance or destruction of the vegetation. Frost tubes are found to have a linear deviation when measuring thaw depth, but as we are comparing between plots this was deemed to be acceptable (Iwata et al., 2012). The frost tubes were constructed according to the frost tube protocol (GLOBE, 2019). Each frost tube consisted of an PVC tube, 2 m long with an outer diameter of 20 mm, a screw cap on the top and a watertight seal at the bottom (Figure 2). The tube contained an PVC hose filled with water, sealed at both ends. To install the tubes, holes with a 2 cm diameter were drilled to a soil depth of 100–140 cm, into which the tube was lowered. If any air gaps existed between the tube and the soil existed these were filled with local soil. Monitoring of TD was performed once per week in 2021, from the first snow melt (May 31) until the onset of freezing (October 3). Measurement of TD started when less than 25% of the surface of the plot was covered by snow. To measure TD, the distance from the top of the hose to the thaw front was measured in the field in cm, and thaw depth calculated by subtracting the tube extension above the ground from this measurement. Figure 1. Location of study grids 1 to 4 in the vicinity of Longyearbyen. Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 3 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2.2.2. Soil Moisture In all plots, surface level volumetric soil moisture (%) in the top 6 cm was recorded at the same time as TD using handheld sensors (delta‐T ML2x moisture probe, calibration setting mineral soil) as soon as the soil thaw progressed far enough to insert the probe. The average of three measurements in each plot was used for further analyses. In order to correlate the moisture values to the progressing TD, a cumulative soil moisture value was chosen to function as an index to indicate soil moisture conditions in the period prior to the measurement. This cumulative value was calculated by adding the mean of the three soil moisture measurements to the sum of means recorded earlier in the same plot. 2.2.3. Relative Elevation, Microtopography and Aspect The elevation of each plot relative to a fixed reference point outside each grid was measured in m using a Topcon GTS‐226 total station (Topcon Corporation, Japan). Relative plot elevation was obtained by subtracting the lowest value in each grid from plot values. Microtopography was described by an index of surface roughness, obtained as the standard deviation of distances (in cm) between 25 regularly spaced points on a frame, placed over each 1 m ×1 m plot, and the soil surface. Aspect was measured in degrees in 2018 (Kemppinen, Niittynen, le Roux et al., 2021; Niittynen et al., 2020b). 2.2.4. Snow Depth and Snow Melt Snow depth was measured in cm using an avalanche probe, once in February and March 2021 (26‐02‐2021 and 31‐03‐2021). Snow cover (%) was monitored twice per week starting 25‐05‐2021 until all snow had disappeared from the research plots. 2.2.5. Air Temperature Air temperature data (°C) were obtained from the Adventdalen weather station using a PT1000 thermometer located 2 m above ground level (UNIS, n.d.). The weather station was located between 2 and 3.5 km from the grids (Table S1 in Supporting Information S1). Measurements were collected hourly and used to calculate daily averages. Cumulative air temperature was calculated by summing up daily positive mean temperatures between grid sampling days. 2.3. Biotic Measurements 2.3.1. Vegetation In July 2018, data were collected on vegetation composition (in %; including species level coverage estimates on vascular plants, lichens, mosses and rock cover, see Niittynen et al. (2020b) for methodology and Niittynen et al. (2020a) for data). An overview of dominant species can be found in Table S1 in Supporting Information S1. Organic layer depth (cm) included the uppermost peaty organic layer and the underlying soil horizon that was clearly influenced by organic matter (Kemppinen, Niittynen, le Roux et al., 2021; Niittynen et al., 2020b). Vegetation height measurements (cm) were carried out in July 2020 by placing a 1 m ×1 m frame with 25 regularly spaced points over each plot. At each point, the tallest forb or dwarf shrub (excluding the infloresce), the tallest graminoid and moss layer thickness were measured. Graminoids were recorded separately from forbs and shrubs as they clearly protruded above average vegetation height. Moss thickness was measured at each point using a ruler to millimetre accuracy, measuring the height of the green part of the moss layer. 2.4. Data Analysis Kendall's rank correlation coefficients τ (Kendall, 1938) between predictor variables were calculated as an initial analytic step. Time of snowmelt and onset of soil thaw varied considerably among plots within and among grids Figure 2. Illustration of frost tube as installed in the field. Distances are indications, precise measurements varied among plots. The inner blue part of the frost tube denotes liquid water (T>0°C) while the white part denotes ice (T<0°C). Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 4 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License (28 days). In order to account for such variations, the timing of TD measurements was standardised to days since complete plot snow melt, and grouped into weekly intervals (Week of thaw: W OT ). that is, W OT 1 corresponds to the week in which a plot became snow‐free. As the active layer progressed deeper than the installation depth of the majority of the tubes, data on TD after the tubes completely thawed and total ALT was unavailable in these plots. Therefore, focus was put on the first eight weeks of thaw, when the data set was mostly complete (56 out of 80) and balanced. For each pair of predictor variables, we tested the hypothesis "Kendall's τ=0" against the two‐tailed alternative hypothesis. No correction for multiple testing was applied as the intention was to evaluate the specific paired relationships. The response variable (TD) and several predictor variables showed strongly skewed frequency distributions. In order to comply with assumptions of linear statistical models, a zero‐skewness transformation (Økland et al., 2001,2006) was carried out, separately for each subset. This transformation leveled out differences in magnitude between variables, thus allowing comparisons among different variables recorded in the same W OT . Results for the same variable in different W OT were visualized using back‐transformed predictions from the resulting Linear Mixed Models (LMM) in all graphs. Thus, a given variable could be compared over time. LMMs with TD as response variable and grid as random factor were obtained separately for each single predictor variable as well as for the set of all predictors. A Fisher's combining probabilities test (Fisher, 1970) was carried out for the results of the LMMs for each W OT in order to test for an overall relationship between predictors and TD. Variables were selected for the multi‐predictor model by forward selection using the Bayesian Index Criterion (BIC). All analyses were performed using R version 4.1.2 (R Core Team, 2021), LMMs using the package nlme (Pinheiro et al., 2021) and graphs were made with the ggplot2 package (Wickham, 2016). 3. Results 3.1. Effects of Snow Conditions The day of snowmelt varied considerably among plots, both between and within grids. The first plots became snow‐free on the 31st of May, while the last were uncovered on the 28th of June. Most predictors were significantly correlated (P<0.05) with snow depth and/or day of snowmelt, except for total vegetation cover, graminoid height, slope angle and aspect, while moss thickness only showed a significant correlation with maximum snow depth (Table S2). The strongest correlation among any pair of predictors was found between day of snowmelt and snow depth (τ=0.832, P<0.001). Strong correlations were also found between day of snowmelt and moss cover, vascular plant cover, lichen cover and rock cover (τ=0.326, P<0.001; τ= − 0.345, P<0.001; τ= − 0.211, P=0.014; τ= − 0.243, P=0.006 respectively). Both maximum snow depth and day of snowmelt were negatively correlated with vascular plant‐ and lichen cover and positively correlated with moss cover, while the correlation with total vegetation cover was neutral. Both day of snowmelt and maximum snow depth were positively correlated with forb and shrub height (τ=0.206, P=0.011 and τ=0.228, P=0.003, respectively). Snow depth was positively correlated with moss thickness (τ=0.162, P=0.035) while day of snowmelt showed no correlation (τ=0.121, P=0.135 for day of snowmelt). A strong negative correlation was also found with relative elevation (τ= − 0.337, P<0.001 for day of snowmelt; τ= − 0.330, P<0.001 for snow depth). Both snow‐condition variables were strongly positively correlated with air temperature and soil moisture over the different weeks of thaw (Table S3 in Supporting Information S1). 3.2. Relationship Between Variables and Progress of Soil Thaw Most variables were significantly correlated with TD at least once during W OT 1–7 (Table 1). Only moss cover, rock cover, cumulative soil moisture and microtopography showed no significant correlation. A Fisher's combination test carried out over all LMMs separately for each W OT showed that TD was affected by the measured variables during W OT 1–7, but not in W OT 8 (Table 2). Average moss thickness had a significant negative relation to TD in six out of eight weeks (Table 1). It showed the strongest slope and R 2 values W OT 3, 4, 5 and 7, and only the strongest slope value in W OT 1 and W OT 6. Highest R 2 values in those weeks of thaw were obtained for cumulative air temperature (W OT 1) and day of snowmelt (W OT 6), both positively relating to TD when significant. Forb and shrub height showed a similar trend as moss thickness, negatively relating to TD in six out of eight weeks with slightly lower slope and R 2 values. Overall, relation between cover values and TD showed lower slopes than the relation between TD and the various Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 5 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License vegetation heights. Moss cover had no significant relation to TD, and vascular plant cover only showed a significant negative relation during W OT 4–7. Organic layer depth also showed a significant negative relationship to TD in W OT 1–5 and 7 (and a Pvalue slightly above 0.05 in W OT 6) with slopes comparable to cumulative air temperature, while R 2 values were relatively low. Similarly, graminoid height was negatively related to TD in W OT 3, 4, 6 and 7, yet with a low R 2 . Table 1 Results of Linear Mixed Models (LMM) of Single Variables Throughout the Early Summer, With TD as a Response Variable Slope pMarginal R 2 Slope pMarginal R 2 Slope pMarginal R 2 Slope pMarginal R 2 W OT 1 W OT 2 W OT 3 W OT 4 Moss Cover −0.09 0.365 0.011 −0.1 0.393 0.008 −0.08 0.421 0.007 −0.08 0.381 0.008 Vascular plant cover −0.04 0.762 0.001 −0.16 0.096 0.035 −0.19 0.064 0.031 −0.29 0.002 0.085 Lichen cover −0.01 0.895 0 −0.22 0.011 0.068 −0.18 0.053 0.031 −0.14 0.099 0.023 Rock cover 0 0.086 0.039 0 0.515 0.004 0 0.336 0.01 0 0.085 0.030 Moss thickness −0.27 0.028 0.063 0.01 0.905 0 −0.49 <0.001 0.212 −0.46 <0.001 0.213 Forb and shrub height −0.2 0.104 0.035 −0.44 <0.001 0.166 −0.41 0.003 0.126 −0.33 0.006 0.097 Graminoid height −0.23 0.098 0.037 −0.06 0.627 0.002 −0.27 0.036 0.048 −0.3 0.010 0.068 Cumulative soil moisture 0.1 0.586 0.006 0.18 0.088 0.04 0.12 0.347 0.008 0.1 0.376 0.007 Relative elevation −0.13 0.209 0.023 −0.45 <0.001 0.29 −0.16 0.074 0.028 −0.14 0.092 0.024 Cumulative mean air temperature 0.25 0.042 0.089 0.25 0.014 0.126 0.39 0.001 0.182 0.28 0.001 0.141 Organic layer depth −0.24 0.028 0.064 −0.37 0.001 0.099 −0.24 0.027 0.055 −0.28 0.005 0.084 Microtopography −0.01 0.936 0 −0.13 0.204 0.019 −0.03 0.821 0 0.03 0.801 0.001 Day of snowmelt 0.16 0.082 0.057 −0.07 0.502 0.005 0.38 <0.001 0.203 0.33 <0.001 0.192 Max. snow depth 0.06 0.514 0.006 0.39 <0.001 0.357 0.33 0.002 0.113 0.3 0.002 0.115 Aspect 0.02 0.863 0 0.26 0.012 0.106 −0.12 0.497 0.02 −0.11 0.509 0.018 Slope angle −0.04 0.736 0.002 −0.35 0.025 0.073 −0.17 0.141 0.021 −0.12 0.256 0.012 W OT 5 W OT 6 W OT 7 W OT 8 Moss cover −0.08 0.344 0.011 0.01 0.942 0 −0.08 0.387 0.011 0.02 0.818 0.001 Vascular plant cover −0.25 0.004 0.085 −0.3 0.003 0.094 −0.3 0.004 0.104 −0.2 0.081 0.054 Lichen cover −0.08 0.310 0.010 −0.13 0.163 0.019 −0.04 0.729 0.001 −0.14 0.236 0.026 Rock cover 0 0.125 0.027 0 0.151 0.024 0 0.160 0.029 0 0.740 0.002 Moss thickness −0.41 <0.001 0.213 −0.34 0.010 0.113 −0.39 0.002 0.176 −0.14 0.236 0.027 Forb and shrub height −0.25 0.027 0.072 −0.32 0.017 0.088 −0.38 0.004 0.144 −0.18 0.196 0.037 Graminoid height −0.19 0.082 0.038 −0.33 0.009 0.083 −0.36 0.005 0.111 −0.24 0.069 0.061 Cumulative soil moisture 0.08 0.433 0.007 0.1 0.360 0.009 0.08 0.443 0.008 0.16 0.157 0.037 Relative elevation −0.02 0.832 0.001 0 0.965 0 0.1 0.294 0.019 0.13 0.144 0.038 Cumulative mean air temperature 0.18 0.020 0.087 0.22 0.014 0.097 0.14 0.137 0.042 0.08 0.269 0.022 Organic layer depth −0.23 0.013 0.075 −0.21 0.053 0.047 −0.22 0.039 0.066 −0.05 0.654 0.004 Microtopography 0.02 0.851 0 0.04 0.700 0.001 0.08 0.496 0.006 0.1 0.383 0.014 Day of snowmelt 0.24 0.003 0.145 0.28 0.002 0.153 0.18 0.059 0.069 0.12 0.177 0.033 Max. snow depth 0.23 0.011 0.086 0.29 0.003 0.115 0.2 0.055 0.060 0.13 0.241 0.025 Aspect −0.03 0.840 0.002 0.01 0.962 0 0.07 0.668 0.009 0.02 0.834 0.001 Slope angle −0.05 0.577 0.004 −0.11 0.315 0.012 −0.05 0.674 0.003 0.06 0.593 0.005 Note. Values in bold indicate significant tests (P<0.05). Note that slope values can only be compared within the same WOT, as data were zero‐skewness transformed separately (see Figure 3for comparisons of each variable over time). Values in italics represent the strongest slope and the highest R 2 value in each WOT. Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 6 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Cumulative air temperature consistently showed a positive relationship with TD until W OT 7, maximum snow depth and day of snowmelt were positively related to TD in W OT 2–6 and W OT 3–6, respectively. Day of snowmelt was significantly related to TD throughout this period, even though the data were corrected for timing of snow melt. Maximum snow depth had the highest R 2 value of all variables in W OT 2, and had a relatively steep slope, while day of snowmelt had the highest R 2 in W OT 6. Terrain variables showed less strong relationships with TD than vegetation or climate variables. Relative elevation, slope angle and aspect were significantly related to TD in W OT 2 only, while no significant relationship was found for rock cover and microtopography. Relative elevation was strongly related to TD in W OT 2, with a high R 2 and the highest slope value for TD. After back‐transformation of the LMM results, results showed a relatively flat slope of TD with respect to the predictor variables in W OT 1 and 2, followed by the steepest slopes throughout the research period, observed in W OT 3 and 4 (Figure 3and Figure S5 in Supporting Information S1). After that, the response attenuated again during W OT 5–7 and by W OT 8 all variables showed no significant relations (Table 1). Thus, the strongest relation of TD to the predictor variables was observed 3–4 weeks after onset of soil thaw. Forward selection of predictor variables in LMMs showed that average moss thickness and day of snow melt were the generally most important drivers of TD (Table 3). These predictors were selected in W OT 3–6 (together with relative elevation in W OT 6 and 7). Moss thickness showed a negative relationship to TD, whereas day of snow melt was positively related to TD. The negative slope of moss thickness was stronger than the positive slope of timing of snow melt except for W OT 6, where the order was reversed, and in W OT 7, where the slopes were almost equal in magnitude although in different directions. Furthermore, day of snow melt and average moss thickness each had higher slopes when selected together compared to being the single predictor (Table 1). Moss thickness was also selected in W OT 1, as the only significant variable and W OT 7 together with maximum snow depth and Figure 3. Predictions from selected linear mixed models, based upon results of biotic predictors in Table 1. Lines show the model, while scatterplots portray the data on which the models were based.The figures show predictions back‐transformed to the original scales used for recording the variables. Variables were included if at least two models were significant, and maximum R 2 was >0.1. Pand R 2 values are the same as in Table 1. Table 2 Results of Fisher's Combining Probabilities Test, Carried Out for the Weekly LMM Results Presented in Table 1 W OT 1 W OT 2 W OT 3 W OT 4 W OT 5 W OT 6 W OT 7 W OT 8 χ 2 48.065 114.266 111.570 125.562 88.799 87.256 79.517 39.7428 df 32 32 32 32 32 32 32 32 P0.034 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 0.163 Note. Bold values indigate significant tests (P<0.05). Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 7 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License relative elevation. In W OT 2, the selected variables were forb and shrub height and maximum snow depth. W OT 1 and 8 had lower R 2 values, compared to W OT 2–7. Marginal R 2 values reached as high as 0.580. In W OT 8, no significant variables were selected, corresponding to the lack of significance found in the Fisher's combining probabilities test. 4. Discussion 4.1. Both Mosses and Vascular Vegetation Slow Thaw Depth Progression The results confirm our hypothesis that both High Arctic mosses and High Arctic vascular vegetation negatively affect TD at the start of the summer. Overall, moss thickness explained the variation in TD best of the two, which were negatively correlated. While vascular plant height and cover also were significantly correlated with TD, the relationship was overall less strong compared to moss thickness. This corresponds with the findings of Jaroszynska et al. (2023), who found that bryophytes significantly reduce soil temperatures compared to shrubs and graminoids, especially in colder, alpine climates and on warm, sunny days. The high insulating effect of mosses is hypothesized to be due to their higher moisture content and water holding capacity (Soudzilovskaia et al., 2013), which increases their heat capacity and promotes evaporative cooling, therefore acting as a buffer during periods Table 3 Results of Forward Selection of Predictors in Linear Mixed Models With Grid as a Random Factor W OT Selected variables BIC Marginal R 2 Conditional R 2 1Predictor Moss thickness 16.071 0.063 0.063 Slope −0.27 P<0.001 2Predictor Forb and shrub height Max. snow depth −26.982 0.580 0.583 Slope −0.469 0.367 P<0.001 <0.001 3Predictor Day of snowmelt Moss thickness −24.161 0.405 0.623 Slope 0.405 −0.536 P<0.001 <0.001 4Predictor Day of snowmelt Moss thickness −39.488 0.448 0.590 Slope 0.361 −0.54 P<0.001 <0.001 5Predictor Moss thickness Day of snowmelt −45.516 0.434 0.482 Slope −0.528 0.284 P<0.001 <0.001 6Predictor Day of snowmelt Moss thickness Relative Elevation −21.843 0.510 0.639 Slope 0.558 −0.438 0.305 P<0.001 <0.001 0.004 7Predictor Moss thickness Max. snow depth Relative Elevation −20.900 0.435 0.435 Slope −0.457 0.469 0.249 P<0.001 <0.001 0.003 8Predictor Graminoid height −5.704 0.061 0.094 Slope −0.224 P0.07 Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007880 SCHUURING ET AL. 8 of 12 21698961, 2024, 8, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JG007880 by University Of Jyväskylä Library, Wiley Online Library on [07/08/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License