Impacts of snow on soil temperature observed across the circumpolar north

Air temperatures (T_A) across northern high latitudes are increasing twice as fast as the global average (IPCC 2013), which could have significant impacts on permafrost (Chadburn et al 2017), infrastructure (Nelson et al 2001), and the large stock of organic carbon stored in frozen soils (Schuur et al 2015). These impacts are mainly driven by changes in soil temperature (T_S). Thus, understanding the response of T_S to T_A increase is critical to predict the magnitude of these impacts. Seasonal snow cover can cause significant differences between T_A and T_S (dT_SA = T_S−T_A), known as surface offset, through thermal insulation and reflection of solar radiation (Smith and Riseborough 2002, Zhang 2005). Since snow conditions are highly variable, observations and modelling studies show that annual dT_SA differs not only spatially (Zhang 2005, Morse et al 2012, Palmer et al 2012, Throop et al 2012), but also varies with time (Zhang et al 2001, Stieglitz et al 2003, Zhang et al 2005, Isard et al 2007, Osterkamp 2007, Romanovsky et al 2007, Woodbury et al 2009, Lawrence and Slater 2010, Qian et al 2011, Park et al 2014, Wang et al 2017). Previous observation-based investigations of dT_SA are limited to local or regional scales (Zhang et al 2001, Beltrami and Kellman 2003, Stieglitz et al 2003, Isard et al 2007, Osterkamp 2007, Romanovsky et al 2007, Sherstyukov 2008, Woodbury et al 2009, Qian et al 2011, Morse et al 2012, Park et al 2014, Streletskiy et al 2015, Wang et al 2017). Modelling studies can consider the entire circumpolar region but the results differ widely due to model structure and input data (Lawrence and Slater 2010, Wang et al 2016). This uncertainty limits our capacity to assess the impacts of climate change on permafrost and northern ecosystems and affects reconstruction of historical climate change from deep ground (>10 m) temperature profiles (Beltrami and Kellman 2003, Mann and Schmidt 2003, Bartlett et al 2005).

In this study, we compiled observations at 588 sites from climate stations and boreholes across the northern high latitudes. Based on this large number of site observations, we analyzed the spatial and temporal variations of near surface T_S (measured at about 0.2 m below surface), near surface T_A (measured at about 2 m above surface), dT_SA and the relationships with snow conditions. We use the data to explore the following questions. (1) What is the general range of mean dT_SA, and are there any broad scale spatial patterns evident across the circumpolar north? (2) How does annual dT_SA change with climate warming, and is there any evidence that changes in dT_SA affected the response of T_S to T_A? (3) How are the spatial and temporal variations in dT_SA related to snow conditions at the circumpolar scale?

The data used in this study include observations from climate stations and boreholes across the northern high latitudes. We only included sites with latitudes greater than 45 °N. Daily data from Russian meteorological stations were obtained from the All-Russian Research Institute of Hydrometeorological Information-World Data Centre (Sherstiukov 2012a). Daily data from Canadian climate stations were provided by Environment and Climate Change Canada. Borehole ground temperature and T_A data were obtained from Global Terrestrial Network for Permafrost (GTN-P) and the literature.

The Russian climate station dataset includes 458 stations for T_S, 619 stations for snow depth, 599 stations for T_A and precipitation. In combination, 264 stations have measurements of T_S, T_A, snow depth, and precipitation. A detailed description of the T_S dataset is provided in Sherstiukov (2012a). T_S was measured at 12 depths (0.02, 0.05, 0.1, 0.15, 0.2, 0.4, 0.6, 0.8, 1.2, 1.6, 2.4, and 3.2 m). The records are sparse at the top four depths. Therefore, we used observations at 0.2 m as the near surface T_S for the analysis. The observations were mainly from 1985–2011, and a third to a quarter of the stations also had observations from 1963–1976 and 1984. Sherstiukov (2012b) conducted a series of quality checks on the T_S dataset. Each record was given a quality flag but no corrections. In this study, we corrected the obvious errors based on the quality flags but strove to avoid over-correction. Obvious errors corrected included sudden sign changes and significant value changes. We performed the correction by comparing the flagged records with the records on the preceding and the following days. On average, about one to two records of T_S at 0.2 m depth were corrected for each station.

In Canada, daily T_S was measured at about 80 stations across the country by Environment and Climate Change Canada. Most of these stations are in the south of 60 °N. T_S was measured at depths of 0.05, 0.1, 0.2, 0.5, 1.0, and 1.5 m during various periods between 1958 and 2008, mainly 1964–1999. T_S was recorded twice daily, morning (0800 h) and afternoon (1500 h). We only used morning observations in the analysis as there were more missing data in the afternoon observations (Qian et al 2011). In total, 70 stations have at least one year of complete observations for both T_A and T_S. Among these, 69 stations also have at least one year of snow depth observations.

The Global Terrestrial Network for Permafrost provides a comprehensive database for permafrost monitoring parameters, including ground temperatures at various depths, T_A and active-layer thickness. We downloaded (performed in March 2017) data from sites with at least one year of complete observations of T_A and/or ground temperatures at about 0.2 m depth. Some boreholes are near climate stations, so we estimated T_A for these boreholes from the nearby climate station data if T_A was not measured at the site. In total we compiled 145 boreholes with at least one year of complete observations of T_A and ground temperatures at about 0.2 m depth (101 boreholes in Russia, 29 boreholes in Alaska, 15 boreholes in Canada). We also included 109 sites from recently published papers: 53 sites across the Mackenzie Valley Corridor in northwestern Canada (Duchesne et al 2014, Wolfe et al 2010); one site in Peel Plateau in northwestern Canada (averaged for two spruce forest observation sites) (O'Neill et al 2015); two sites for alluvial and uplands in the Kendall Island Bird Sanctuary in the outer Mackenzie Delta (averaged for seven and eight observation sites, respectively) (Morse et al 2012); three sites for black spruce forest, white spruce forest and open black spruce peatland, respectively, in the Great Slave Region (each site was averaged from three nearby observation sites) (Morse et al 2016); three sites for bog, fen and palsa, respectively, in the southern Hudson Bay Lowlands (each site was averaged from two nearby observation sites) (Ou et al 2016); 24 sites in Labrador and Quebec in eastern Canada (measurements at barren and rock sites were excluded in boreal regions and the sites very close to one another were averaged as one site) (Way and Lewkowicz 2018); 11 sites across permafrost regions in Canada (Throop et al 2012); one site in the Lena River Delta, Russia (Boike et al 2013); and 11 sites in southern Norway (Farbrot et al 2011). Together with observations at climate stations, we compiled 588 sites with at least one year of complete observations of both T_A and T_S. The longest dataset available is 38 years. There are 442 sites with at least five years of complete observations although the years may not be continuous. Figure 1(a) shows the distribution of the sites and the mean T_A.

Zoom In Zoom Out Reset image size

Using daily snow depth observations at climate stations, we calculated annual snow cover duration, winter mean snow depth (last December, January and February) and annual mean snow depth (sum of daily snow depths divided by the total number of days in a year). In the text, annual T_A or annual mean T_A (or other variables) is the average of a year at a site, and mean T_A or mean annual T_A (or other variables) is the average of all the years of observations at a site. Annual averages are based on calendar year except for the correlation analysis of annual snow conditions described in the last two paragraphs in section 3.4, where a year is defined from October 1st to September 30th in the following year. Correlation analysis is based on Pearson momentum correlation coefficients (R). Simple linear regressions were calculated based on least squares to estimate temporal trends and linear relations between two variables. We also calculated mean absolute deviation (MAD) to represent inter-annual variations for annual T_A and T_S. These statistics were calculated for sites with at least five years of complete observations. We used a marginal significance level (one-sided p < 0.1) in trends and correlation analysis to maximize the number of sites for spatial coverage. The general patterns are similar when we limited to sites with longer observations and higher statistical significance levels.

3.1. Spatial distribution of mean dT_SA

Mean T_A decreases from south to north (figure 1(a)). While mean T_S shares a similar pattern, it exhibits a greater variation than that of T_A due to local variations in snow, soil and vegetation conditions. Figure 1(b) shows the spatial distribution of mean dT_SA. Mean dT_SA ranges from about 0.0 °C–12.9 °C, with an average of 4.4 °C. Mean dT_SA shows a non-linear relationship with mean T_A: It increases from less than 2 °C to up to 12.9 °C (with an average of 7.3 °C) until mean T_A reaches about −10 °C, then it decreases to about 2 °C until mean T_A is about 7 °C, after which it has no obvious trend (figure 2(a)). Mean dT_SA is greatest in the low arctic and northern boreal regions. It decreases both northward and southward (figures 1(b) and 2). Because of these different patterns, the spatial gradients of mean T_S and T_A vary among the different biomes. The north-south gradient of mean T_S is about one and a half times of the gradient of mean T_A in the arctic, but only about two thirds of that of mean T_A in boreal regions on average. In northern temperate regions, the spatial gradients of mean T_S and T_A are similar, especially in southern areas.

Zoom In Zoom Out Reset image size

3.2. Temporal variations of annual T_S and dT_SA with T_A

To investigate the response of T_S to changes in T_A, we calculated simple linear regressions and Pearson correlation coefficients between annual T_S and annual T_A, and between annual dT_SA and annual T_A based on the time series data for each site (table 1). Since T_A is the primary driver of the variations in T_S, annual T_S is positively and significantly (p < 0.1) correlated with annual T_A for 83% of the sites. However, the relationship between annual T_S and T_A varies with biomes. The R-values tend to be lower in boreal regions than in the arctic and the northern temperate regions (figure 3(a) and table 1). Annual dT_SA is significantly (p < 0.1) but negatively correlated with T_A for 83% of the sites, and the correlation tends to be stronger (closer to −1) in boreal regions than in the arctic and the northern temperate regions (figure 3(b) and table 1). This result indicates a weaker coupling between annual T_S and T_A associated with a stronger correlation between annual dT_SA and T_A in boreal regions than in the arctic and the northern temperate regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Slope coefficients (K) of linear regressions between annual T_S and T_A are less than 1.0 for almost all the sites (only 11 sites are exceptions). The values of K in boreal regions are smaller than in the arctic and the northern temperate regions (figure 3(c), table 1). Conversely, the values of K between annual dT_SA and T_A are negative, and they are smaller (closer to −1) in boreal regions than in the arctic and the northern temperate regions (figure 3(d), table 1). This analysis suggests that annual T_S is less responsive to changes in annual T_A in boreal regions than in the arctic and the northern temperate regions. However, annual dT_SA is more sensitive to changes in T_A in boreal regions than in other regions.

3.3. Impacts of dT_SA on the inter-annual variations and long-term trends of annual T_S

3.4. The causes of the spatial and temporal differences in dT_SA

Across the northern high latitudes, mean annual dT_SA is closely correlated with mean winter dT_SA (R = 0.894, N = 554) but the correlation with mean summer dT_SA is poor and negative (R = −0.265, N = 554 sites), indicating that mean annual dT_SA is mainly determined by winter conditions. Mean winter dT_SA is negatively correlated with mean winter T_A (R = −0.740, N = 554) while the correlation between mean summer dT_SA and mean summer T_A is poor (R = 0.358, N = 554).

Observations of T_A, T_S and snow depth at climate stations in Russia and Canada indicate that mean winter dT_SA is correlated with mean winter snow depth (R = 0.661, N = 333 stations, figure 5(a)). Mean dT_SA is also closely correlated with mean annual snow depth (R = 0.678, N = 329 stations, figure 5(b)). For stations with three to eight months of snow cover in a year, mean dT_SA increases linearly with mean snow cover duration (R = 0.797, N = 295 stations) (figure 5(c)). Mean annual snow depth tends to be low when mean T_A is very low or very high (figure 6(a)), and snow cover duration decreases with the increase in mean T_A (figure 6(b)). Thus, the large dT_SA in the low arctic and northern boreal regions is mainly due to the long duration of relatively thick snow cover in this region. In the high arctic, thin snow cover due to low snowfall, wind redistribution (Morse et al 2012, Palmer et al 2012) and development of wind slab (Derksen et al 2014) reduce the insulation effects on the ground. In the northern temperate regions, short duration of snow cover reduces annual dT_SA.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The year-to-year variation of annual dT_SA is also positively correlated with snow depth and snow cover duration for most of the stations (table 3), indicating that snow condition has significant impacts on the temporal variation of annual dT_SA. Annual mean snow depth and annual snow cover duration (from October 1st to September 30th in the following year) are also positively correlated to each other. Snow effects on the response of T_S to changes in T_A depend on the changes in snow conditions with T_A (due to snowfall and snowmelt) and the consequent changes in dT_SA. Annual snow depths and snow cover durations generally have negative correlations with annual T_A, especially for snow cover durations (table 3) (except for some stations in Siberia and one station in northeast Canada where the correlations are positive). However, the strength of the correlations vary with biomes. In the arctic, the correlations between snow conditions (depth and duration) and annual T_A are generally not statistically significant. The correlations are stronger (more negative) in the northern temperate than in boreal regions. However, annual dT_SA is more sensitive (higher K-values) to changes in snow depth in boreal regions than in the northern temperate regions (table 3). This snow depth effect on dT_SA probably is the main reason for the reduced response of annual T_S to changes in T_A in boreal regions.

Since the periods of observation are short, the temporal trends of dT_SA and snow conditions are not significant for most of the stations. For stations with significant trends (p < 0.1) for both annual dT_SA and snow conditions, the trends of annual dT_SA are positively correlated with the trends of snow conditions (with winter mean snow depth: R = 0.865, N = 40 stations; with annual mean snow depth: R = 0.585, N = 46 stations; with annual snow cover duration: R = 0.669, N = 42 stations).

The effects of snow on ground temperature (the surface offset) have been recognized in previous studies (Smith and Riseborough 2002, Zhang 2005, Palmer et al 2012). Several observation-based studies have assessed the differences between T_S and T_A and their temporal dynamics but only for limited sites and areas. Modelling studies can cover large areas, but their results vary widely due to model structure, resolutions and inputs (Lawrence and Slater 2010, Wang et al 2016). This study synthesizes a large number of observations across the northern high latitudes. The data show that there is a large offset between T_S and T_A in the low arctic and northern boreal regions across the circumpolar, but small offsets in the high arctic, southern boreal and northern temperate regions. In boreal regions, a thick and relatively long snow cover duration and its strong impacts on annual dT_SA variation decouple T_S and T_A and dampen the response of T_S to climate warming (increase in T_A).

The results of this study agree qualitatively with previous observational analysis and modelling studies, which have shown that snow is the major factor differentiating mean T_S from T_A in permafrost regions (e.g. Smith and Riseborough 2002). The insulation effect of snow on annual T_S is at least twice that of the shading and cooling effects of vegetation in summer (Matyshak et al 2015). When snow cover is thin or absent in winter, the annual mean dT_SA can be very small (Lacelle et al 2016). Based on observations at ten sites across Canadian permafrost regions, Throop et al (2012) indicate that dT_SA is the smallest in the arctic sites, and the largest in the boreal of the Mackenzie Valley. A modelling study by Lawrence and Slater (2010) shows that the warming of T_S is less than that of T_A for most of the northern high latitudes during 1950–2100. The ratio of the warming between T_S and T_A is smaller in the boreal than in the arctic and the northern temperate regions (figure 1(c) in Lawrence and Slater 2010).

Climate stations account for about half of the sites reported in this study. They are usually located in flat areas with short vegetation, which can be very different from natural conditions, especially in forest regions. In boreal regions, vegetation is denser and higher than in the arctic, organic layers and mosses are usually deep and wet, and ground ice is common in near surface permafrost. Thus, T_S in natural boreal areas may be even lower and less responsive to changes in T_A than measured at climate stations. On the other hand, arctic regions have low and sparse vegetation, and organic and mineral deposits are usually thin. Thus, the response of T_S to changes in T_A is likely to be direct and rapid (Throop et al 2012).

The results of this study are useful for understanding the current distribution of permafrost and its response to climate warming, for assessing the impacts of climate change on ecosystems and soil organic carbon, for improving models about snow effects on T_S, and for reconstruction of historical climate based on ground temperature profiles at the northern high latitudes. For example, figure 2 and table 2 indicate that the increase of annual T_S with time is smaller than that of T_A in boreal regions. Thus, the warming and degradation of permafrost in boreal regions will be less than the estimates directly from T_A. In contrast, near surface T_S in arctic regions is highly responsive to changes in T_A, therefore these regions can be very sensitive to climate warming. On the other hand, reconstructions of historical climate from deep ground temperature profiles would likely underestimate the changes in T_A in boreal regions. For instance, the estimated warming of near-surface T_S during 1950–2000 in Canada (Beltrami et al 2003) is smaller than that of T_A observed at climate stations, especially in boreal regions (Vincent et al 2015, Zhang et al 2000).

Our results show the insulative effects of snow on T_S at the circumpolar scale. The data also show strong site variations due to local and regional ground and climate conditions. Significant changes in snowfall could directly affect snow depth, thus altering dT_SA and the response of T_S to T_A. Changes in soil and vegetation conditions due to fire disturbances and gradual responses to climate warming could also affect T_S and its response to T_A.

This study compiled T_S and T_A observations from 588 climate stations and boreholes across the northern high latitudes. This broad observational data coverage analyzed in this paper shows a large offset between mean T_S and T_A in the low arctic and northern boreal regions across the circumpolar north. The offset decreases both northward and southward. Thus, the north-south gradient of mean T_S is about one and a half times of that of T_A in the arctic, but only about two-thirds of T_A in boreal regions. Further south, the spatial gradients of mean T_S and T_A are similar. This pattern is closely related to snow conditions. Thick and persistent snow cover and its strong impacts on annual dT_SA in boreal regions weaken the coupling between annual T_S and T_A, and the response of annual T_S to variation in annual T_A is smaller in the boreal than in the arctic and the northern temperate regions. Consequently, the inter-annual variation and the trends of T_S are smaller than that of T_A in boreal regions. This systematic and significant differences in the relationship between T_S and T_A across the circumpolar north is important for understanding and assessing the impacts of climate change and for reconstruction of historical climate based on ground temperature profiles for the northern high latitudes.

The authors thank Steve Wolfe and Peter Morse for stimulating discussion and their careful review of the earlier versions of the paper. This study was funded by Polar Knowledge Canada Science and Technology Program (project 186). This study also contributes to a project affiliated to the Arctic Boreal Vulnerability Experiment (ABoVE), a NASA Terrestrial Ecology program.