Impacts of snow on soil temperature observed across the circumpolar north

Climate warming has significant impacts on permafrost, infrastructure and soil organic carbon at the northern high latitudes. These impacts are mainly driven by changes in soil temperature (TS). Snow insulation can cause significant differences between TS and air temperature (TA), and our understanding about this effect through space and time is currently limited. In this study, we compiled soil and air temperature observations (measured at about 0.2 m depth and 2 m height, respectively) at 588 sites from climate stations and boreholes across the northern high latitudes. Analysis of this circumpolar dataset demonstrates the large offset between mean TS and TA in the low arctic and northern boreal regions. The offset decreases both northward and southward due to changes in snow conditions. Correlation analysis shows that the coupling between annual TS and TA is weaker, and the response of annual TS to changes in TA is smaller in boreal regions than in the arctic and the northern temperate regions. Consequently, the inter-annual variation and the increasing trends of annual TS are smaller than that of TA in boreal regions. The systematic and significant differences in the relationship between TS and TA 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.


Introduction
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 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, 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 andSlater 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?

Data and processing
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-1976and 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 .
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.

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 nonlinear 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.

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.
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.

Impacts of dT SA on the inter-annual variations and long-term trends of annual T S
A negative correlation between annual dT SA and T A or a smaller response of annual T S to the changes in annual T A has two implications. First, the interannual fluctuation of T S will be smaller than that of T A . Second, if there is a long-term increasing trend in annual T A , the increase of annual T S will be less than that of T A . To confirm these patterns, we calculated MAD and linear temporal trends for annual T S and T A based on the time series data for each site with at least five years of observations. The ratio between MAD of annual T S and MAD of annual T A increases linearly with the K between annual T S and T A and with the R between annual dT SA and T A (figures 4(a) and (b)). The ratio between the trend of annual T S and the trend of annual T A also increases linearly with the K between annual T S and T A and with the R between annual dT SA and T A (figures 4(c) and (d)). This result indicates that the inter-annual variation and the trend of T S are closely related to the response of T S to changes in annual T A and the correlation between dT SA and T A . When the response of annual T S to Figure 3. The spatial distributions of R and K calculated based on time series data for each site with at least five years of observations. (a) R between annual T S and annual T A ; (b) R between annual dT SA and annual T A ; (c) K between annual dT S and annual T A when the correlation is significant at p < 0.1; and (d) K between annual dT SA and annual T A when the correlation is significant at p < 0.1. For a site, the K between annual T S and annual T A equals the K between annual dT SA and annual T A minus 1. The data are divided into 15 classes using natural breaks provided by ArcGIS. The grey area is for land or glaciers, and the light green area is for boreal regions. The blue and black curves are the southern boundaries of continuous and isolated patches of permafrost, respectively. change in T A is small (small K-values), the inter-annual variation and the trend of T S are smaller than that of T A . Similarly, when annual dT SA and T A are strongly correlated, the inter-annual variation and the trend of T S are smaller than that of T A .
For different biomes, the ratio between MAD of annual T S and MAD of annual T A tends to be smaller in boreal regions than in the arctic and the northern temperate regions (table 2). Similarly, the ratio between the trend of annual T S and the trend of annual T A tends to be smaller in boreal regions than in the arctic and the northern temperate regions (table 2). These results are consistent with the weaker coupling between annual T S and T A (stronger correlation between annual dT SA and T A ) and smaller K-values between annual T S and T A in the boreal than in the arctic and the northern temperate regions (table 1).

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)).

Figure 4. The distributions of the ratios between the mean absolute deviation (MAD) of annual T S and MAD of T A with (a) the K between annual T S and T A and (b) the R between annual dT S and T A .
And the distributions of the ratios between the trend of annual T S and the trend of annual T A with (c) the K between annual T S and T A and (d) the R between annual dT S and T A . Each circle represents one site with significant correlations (p < 0.1) for the K between annual T S and T A and for the trends of annual T S and T A . They are calculated based on the time series of data for each site with at least five years of observations. Table 1. The correlation between annual T S and annual T A and between annual dT SA and annual T A for sites with at least five years of observations.

Biomes
The    Correlation between annual mean snow depth and annual mean T A Average ( 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 ( 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).

Discussion
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 andSlater 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 Figure 5. Scatter graphs (a) between mean winter dT SA and mean winter snow depth, (b) between mean annual dT SA and mean annual snow depth, and (c) between mean annual dT SA and mean annual snow cover duration. Each circle represents one site (averaged for all the years with data available).
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 Figure 6. The distribution of mean T A with (a) mean annual snow depth, and (b) mean annual snow cover duration. Each circle represents one site (averaged for all the years with data available). The data are from 329 climate stations with T S observations in Russia and Canada. 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  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 .

Conclusions
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.