Decadal Changes in Soil and Atmosphere Temperature Differences Linked With Environment Shifts Over Northern Eurasia

The difference between soil and air temperatures (ΔT) in a specified time is dependent on meteorological conditions, properties of soil and land covers. Understanding ΔT is critical in assessing land–atmosphere thermal interactions in changing environment. However, systematic knowledge of interannual variations and responses of ΔT to the environmental changes (e.g., snow cover and soil moisture) at decadal scales remain limited. Here, variations of the mean annual air and soil temperatures, and ΔT were investigated at 217 sites in northern Eurasia during 1981–2015. It is found that changes in the mean annual air and soil temperatures were inconsistent as the average increase in soil temperature was generally less than that of air. The relationships between trends in soil and air temperatures were significant in the upper ground (0.2 and 0.8 m) over the seasonal frost region but insignificant in the permafrost regions and deeper ground. During the period, widespread changes in ΔT occurred and closely responded to the environmental changes, but the relationships varied with soil depth. Among the tested factors, trends in snow cover thickness dominantly control trends in ΔT, followed by trends in snow cover duration and solar radiation. Both linear and nonlinear analyses indicate enhanced relationships between trends in snow depth and ΔT as depth increases. This study provides the first view of decadal trends in ΔT and conjunctions with the environmental changes during 1981–2015 over northern Eurasia. The findings are relevant to quantify land–atmosphere thermal interactions given impacts of future environmental changes.

. Systematic interannual changes in environmental variables discord trends in the mean annual air and soil temperatures by introducing nonatmospheric trends to soil temperature records, and hence causes temporal variations of ΔT (Bulygina et al., 2011;García-García et al., 2016;González-Rouco et al., 2003;Grundstein et al., 2005;Peng et al., 2016;Romanovsky et al., 2007). However, long-term interannual variations of ΔT and its response to the changing environmental conditions are little known yet as most previous efforts have been devoted to ΔT and its controls at seasonal (or diurnal) scale (Aalto et al., 2018;Bartlett et al., 2004;Wang et al., 2017).
Examining trends in ΔT and the role of environmental changes at large spatiotemporal scales are challenging due to data limitation in terms of the consistency and spatiotemporal coverage of observations. Thus, it has constrained only partial characterizations from case studies in short periods or by excluding key elements as snow cover has received more attention while, e.g., vegetation and soil moisture are rarely addressed simultaneously (Bartlett et al., 2004;Hrbáček et al., 2016;Jungqvist et al., 2014;Yazaki et al., 2013;Zhang, 2005). In northern Eurasia, multidepth soil temperatures and snow cover depth have been documented for decades, enabling us to investigate the interannual variability of ΔT at a large spatiotemporal scale.
Northern Eurasia is a key part of the global earth system, and regional climate warming is greater than the global average accompanied by increased frequency and intensity of extreme weather events (Blunden & Arndt, 2019;Groisman et al., 2016;Papalexiou & Montanari, 2019). The weather extremes even have prolonged impacts and feedback to land-atmosphere thermal interactions after taking place (Knapp et al., 2008). Previous attention has been mostly paid to ΔT between air and topsoil temperatures (<0.2 m) (Aalto et al., 2018;Bartlett et al., 2004;Luo et al., 2018;Zhang et al., 2018). In contrast, the knowledge of ΔT between air and deeper ground (>0.2 m) remains limited and is also of great importance as over half of soil carbon out of 2,344 Gt C at depths of up to 3 m is stocked in 0.2-2.0 m layer (Jobbágy & Jackson, 2000). Moreover, in frost-related regions, soil carbon decomposition in the deeper ground may proceed more rapidly than the upper under a warming climate (Bernal et al., 2016;Koven et al., 2015). Thus, examining long-term trends in ΔT and, in conjunction with the environmental changes, are essential in assessing land-climate thermal interactions under a changing climate.
Here, we provide the first holistic view of observed trends in ΔT at multiple soil depths and environmental variables at 217 sites in northern Eurasia during 1981-2015. Further, we use statistical methods to investigate links between trends in ΔT and environmental variables, including snow cover, rainfall, solar radiation, vegetation, and soil moisture during the period.

Data Compilation
We extracted daily air and soil temperatures, snow depth, and rainfall observations collected during 1981-2015 at weather stations operated by the All-Russia Research Institute of Hydrometeorological Information -World Data Center (RIHMI-WDC). The observation network covers the majority of land in northern Eurasia at latitudes from 42 to 70°N and longitudes from 28 to 170°E (Figure 1). The air temperature was measured 2 meters above ground, and soil temperatures were measured under natural cover of grasses and snow at depths of 0.2, 0.8, and 1.6 m (Gilichinsky et al., 1998). The measurement accuracies of temperatures and snow depth are within 0.1°C and 0.1 cm, respectively. Both datasets were released to the public after extensive quality control checks, and only qualified data were used in the study (Text S1 and Table S1) (Bulygina & Razuvaev, 2012;Sherstiukov, 2012).
The daily measurements were then computed into annual values: mean annual soil temperature (MAST), mean annual air temperature (MAAT), the annual rainfall (Rainfall, liquid precipitation), the annual mean of daily snow depth (AMSD, an average of daily snow depth during a year), snow cover duration (SCD, a sum of the days that snow exists during a year). In computing annual values, we defined a qualified year as a calendar year, during which at least 300-day data are available. On average, there are over 360-day measurements available for computing annual values in the datasets (full summary in Table S2). Apart from in situ observations, we extracted monthly gridded (9 km spatial resolution) attributes from ERA5-Land products, which were then computed to annual values: surface net solar radiation (SolarRad, absorbed solar radiation), leaf area index (LAI), volumetric water content (VWC, including the ice content) (Hersbach et al., 2020). The VWC is available at three layers (VWC 0.2 , VWC 0.8 , and VWC 1.6 , corresponding to soil depth intervals of 0-0.2 m, 0-0.8 m, and 0-1.6 m, respectively).

Data Processing and Analyses
We investigated trends in air, soil temperatures, and other environmental variables during 1981-2015. Trends in individual variables were calculated as the slopes of linear least squares regressions at each location. The criteria in computing trends are the following; 1). A station has to have a minimum of 10-year CHEN ET AL.  qualified data. 2). At least 90% of the years are needed during computing period (for example, if observations at a station spanned across 10 years during 1981-1990, at least 9-year data are needed for computing trends). 3). Trends of all variables were calculated in a common set of years for consistency in analyses. We then defined the annual difference (ΔT, °C) between MAST and MAAT as ( Figure 1):

ΔT = MAST -MAAT
We examined trends in ΔT (ΔT 0.2 , ΔT 0.8 , and ΔT 1.6 ), which is derived from MAAT and MASTs at multiple soil depths (0.2, 0.8, and 1.6 m). In total, we assembled 184, 194, and 196 sites in calculating trends of ΔT 0.2 , ΔT 0.8 , and ΔT 1.6 , respectively, due to the different availabilities of soil temperature data with depth. The average numbers of years used for calculating trends range from 27.4 to 30.2 years at depths during 1981-2015 (details in Figure S1).
Then, trends in ΔT were linked to the environmental changes by generalized additive mixed model (GAMM). Unlike linear regression, GAMM captures the impact of changing environment on ΔT trends by smooth functions, which can be linear or nonlinear depending on the underlying patterns in the data (Text S2) (Hastie & Tibshirani, 1990;Wood, 2017). The number of years was included as a random term because of the different years used for computing trends at the sites ( Figure S1). The collinearity diagnostics by multiple indices including pairwise correlations, condition index (CI) (Belsley et al., 2005;Douglass et al., 2003), variance inflation factor (VIF) (Hair et al., 1995), and tolerance suggest that no highly correlated variables were included in the models (Table S3, Figure S2) (Dormann et al., 2013). The response curves of trends in ΔT and other variables were fitted and optimized by the restricted maximum likelihood method (REML) (Patterson & Thompson, 1971;Wood, 2017) and visualized with 95% confidence intervals. The statistical significance of each variable across models is reported by three levels: very significant: p ≤ 0.001; significant 0.001< p ≤ 0.05; insignificant; p > 0.05. The detailed summary statistics of models are available in Table S4. To characterize the contribution of each variable in the models, we use Pearson correlation coefficient between fitted values of the models and predictions where the variable under investigation has been randomly permutated (Breiman, 2001;Thuiller et al., 2009). Then, the importance of the variable is simply calculated by the following equation: Variable importance = 1 -|Pearson correlation (prediction original varaible , prediction one variable permutated )| Consequently, the magnitude of variable importance is constrained from 0 to 1, in which the higher, the more influential the variable is. Finally, the variable importance was averaged over 100 permutations with an error bar representing a 95% confidence interval over the permutations. All analyses were undertaken on R platform (v.3.6.3) (R core team, 2017) with packages of ggplot2 (Wickham, Chang, et al., 2020), ggpubr (Kassambara, 2020), dplyr (Wickham, François, et al., 2020), cowplot (Wilke, 2019), mgcv (Wood, 2019), mgcViz (Fasiolo et al., 2020), biomod2 (Thuiller et al., 2020).

Trends in Air and Soil Temperatures From 1981 to 2015
The increases in soil and air temperatures during 1981-2015 were significant and varied spatially . For the region as a whole, warming rate of soil temperature at 0.2 m by 0.36 ± 0.03°C decade −1 (Mean ± standard error of the mean) was greater than the deeper ground (1.6 m, 0.27 ± 0.02°C decade −1 ). The increase in air temperature at 0.34 ± 0.01°C decade −1 was slightly less than soil warming at 0.2 m but greater than the warming in the deeper ground. Regionally, increases in soil temperatures at multiple depths were generally less than that of air temperature except at a depth of 0.2 m in the seasonal frost region. Of note, in the discontinuous permafrost region, the warming rate of soil temperature at 0.2 m was only about half of that of air; however, soil temperatures warmed at a similar rate to air temperature at greater depth. This difference can be attributed to the latent heat effects due to freeze-thaw cycles of water, which is more significant close to the surface than at greater depth. Therefore, change in MAAT cannot account for the warming of MASTs with depth and such inconsistency consequently caused interannual variations in ΔT during the period.
The linear relationships indicate that trends in air and soil temperatures at 0.2 and 0.8 m were significantly related over the seasonal frost region ( Figure 3 and Table 1). However, in the permafrost regions and a The total numbers of sites are reported in brackets. The extent of frozen ground (Continuous permafrost region (extent of permafrost ≥ 90%), Discontinuous permafrost region (0% < extent of permafrost < 90%)) is derived from permafrost map of the International Permafrost Association (Brown et al., 1997).  Table 1. deeper depth (1.6 m), there were no clear relationships between trends in air and soil temperatures. In the region as a whole, the relationships between trends in soil and air temperature were significant at depths of 0.2 and 0.8 m, and insignificant at 1.6 m.

Decadal Trends in ΔT during 1981-2015
Although the average trends in ΔT are slightly below the null (0°C decade −1 ) at multiple depths, they are spatially heterogeneous across 217 sites (Figure 4). For example, trends in ΔT 0.2 show the greatest ranging from −2.3 to over 1.0°C decade −1 while similar variability can be found in trends of ΔT 0.8 and ΔT 1.6. Across the entire region, trends in ΔT in the upper ground decrease from south to north as they significantly relate to the latitudes of locations (Figures 5c and 5d). As depth increases, the relationship between trends in ΔT and latitude diminishes and becomes less significant. There are no clear relationships between trends in ΔT and MASTs (p > 0.1), whereas trends in ΔT account for 0.7 of the trends in MAST at 0.2 m and account less in the deeper ground (0.8 and 1.6 m) (Figures 5a-5d).

Linear Relationships Between Trends in ΔT and Environmental Variables
The environmental conditions in the area have changed over the past few decades. Snow characteristics have changed, including a minor increase in AMSD (0.06 ± 0.08 cm decade −1 , Mean ± SEM) but decrease in SCD (−2.39 ± 0.23 days decade −1 ) (Figures 6a, S3a and S3b). The reversed changes between the two snow parameters demonstrate previous studies on increased snowfall extremes and maximum snow depth in the area (Bulygina et al., 2011;Donat et al., 2016). SolarRad has widely increased at over 95% sites with an average of 0.06 w m −2 decade −1 (Figures 6a and Figure S3c). Annual rainfall also significantly increased on an average of 5.20 ± 1.20 mm decade −1 (Figures 6a and S3d). Nevertheless, soil moisture loss is significant at multiple depths over the region, which is likely attributable to the alterations in the evaporation rate and infiltration of soil moisture to the deeper ground (Figures 6a, S3e-S3g, and S4) (Andresen et al., 2020).
Linear regressions demonstrate that trends in AMSD and SCD have significant relationships with trends in ΔT (Table 2). At a depth of 0.2 m, trends in SCD and AMSD explain the variance of ΔT trends at a similar level, while trends in SCD has a greater agreement with ΔT trends at 0.8 m than AMSD. The opposite can be found at 1.6 m, where trends in AMSD explain 30.4 % variance of ΔT trends, which is more than three times the variance attributed to SCD trends (10.0%). Trends in ΔT at 0.2 m significantly relate to trends in SolarRad, but the relationships are not significant at the deeper ground. There are no clear relationships between trends in LAI or Rainfall and ΔT. However, responses of ΔT to the changes in VWC are significant at depths of 0.2 and 0.8, but not at 1.6 m.

Trends in ΔT Linked With the Environmental Changes
To link decadal trends in ΔT and the changing environmental variables, we ran general additive mixed regression models with trends in environmental variables (AMSD, SCD, SolarRad, Rainfall, VWC, and LAI) as fixed smooth terms and the number of years used for computing trends as a random-effect term. It is noted that relationships between observed trends in variables and ΔT are nonlinear, suggesting that results from GAMMs are potentially more robust than the results from linear regressions that cannot capture nonlinear relationships due to the nature of the method (Figures 6b-6d). For example, GAMMs indicate significant relationships between trends in SolarRad and ΔT 0.8 /ΔT 1.6 , which, however, are not reflected by linear regressions (Figures 6b-6d and Table 2). Similarly, trend in ΔT 1.6 significantly responds to VWC trend from GAMM but not as significant in linear regressions.
Trends in ΔT at multiple depths significantly relate to trends in environmental variables by explained variances ranging from 0.38 to 0.56 ( Figure 6 and Table S4). Trends in ΔT are primarily attributed to the snow cover changes. Trends in AMSD are more influential at the deeper ground with over twofold variable importance of ΔT 1.6 trends (0.79 ± 0.006, Mean ± 0.95 CI) than that in ΔT 0.2 (0.31 ± 0.003), which agrees with linear analyses (Table 2 and Figure 6e). Of note, the negative trends in ΔT (meaning soil-atmosphere temperature difference became smaller) associate with the decrease in AMSD consistently and a turnover (negative to positive) of trends in ΔT occurs along with a shift of trend from decreasing to increasing AMSD. Compared with trends in AMSD, the relationships between trends in ΔT and SCD are flatter but still significant.
Trends in SolarRad significantly relate to trends in ΔT at multiple depths, and the influences decrease with depth (Figures 6b-6d). Trend in rainfall is the second influential factor for ΔT 0.2 trend, while the relationships at other soil depths are not as clear. Notably, trends in SolarRad and rainfall poorly explain positive ΔT trends as curves are mostly below the null (Figures 6b-6d). In contrast, trends in soil VWC have significant nonlinear relationships with ΔT trends at all depths, but the response curves are mainly above the null. Moreover, trends in ΔT and LAI are less related with large uncertainties, but the relationship at a depth of 1.6 m is significant and negative, suggesting greening of vegetation may be a contributor to the decrease in ΔT 1.6 by cooling ground even under a warming climate. In contrast, such relationships are less significant in the upper ground at decadal scales.

Discussion
The increase in the mean annual air and soil temperatures occurred at 217 sites in northern Eurasia during 1981-2015. The average increase in mean annual soil temperature is slightly lower than the increase in mean annual air temperature except at a depth of 0.2 m (0.44°C decade −1 ), which is greater than the mean annual air temperature increase (0.36°C decade −1 ) in the seasonal frost region (Figure 2). Such discordance CHEN ET AL.  Violin plots display the densities of trends over the sites. Solid red and black lines depict trends on average and percentiles (5th and 95th), respectively. Averaged trends (Mean ± SEM) are reported at the top. The significance of averaged trends against null (representing no change: 0°C decade −1 ) was tested by the Wilcoxon method and marked with Asterisks (p < 0.05) at the top. between trends in air and soil temperatures demonstrates that air temperature change alone cannot explain soil warming over the region, consequently causing interannual variations of ΔT (Figure 2). It is found that trends in ΔT are spatially heterogeneous and significantly respond to trends in environmental variables, especially snow cover, solar radiation, and soil moisture.
The role of snow cover in decoupling soil and air temperatures has been widely examined (Chen et al., 2021;Romanovsky et al., 2007;Wang et al., 2017). Snow cover generally insulates the ground from coldness due to low thermal conductivity, whereas it does not always warm ground since long-lasting snow cover till next spring (or even summer) could have an overall cooling effect because of the albedo and energy required to melt its pack (Zhang, 2005). Nevertheless, the impacts of the different snow cover parameters (snow cover thickness and duration) on ground temperature are relatively poorly quantified yet (Isaksen et al., 2011). A modeling study by Bartlett et al. (2004) suggests that snow layer thickness is less critical than snow cover duration in modulating the mean annual surface ground temperature, while Romanovsky et al. (2007) highlights the influence of snow cover thickness on the mean annual soil temperature (at a depth of 1.6 m). We highlight that trends in ΔT have closer linear relationships with trends in AMSD than SCD except at depths of 0.2 and 0.8 m in the discontinuous and continuous permafrost regions, respectively (Table 2). In the whole area, trends in AMSD (R 2 = 0.304) account for three times CHEN ET AL. greater variance of ΔT 1.6 trends than trends in SCD (R 2 = 0.100), while a similar difference in variance of ΔT trends is observed at shallower depths. Meanwhile, nonlinear analyses by GAMM also suggest that ΔT trends more closely respond to trends in AMSD than SCD at three depths in the area. Moreover, both linear and nonlinear analyses indicate that changes in AMSD tend to be more influential on trends of ΔT in the deeper ground than the upper over northern Eurasia during 1981-2015 (Table 2 and Figure 6e).
Essentially, all variables affecting hydrothermal fluxes near the land-atmosphere interface and within the ground influence soil temperature. At middle and high latitudes, ΔT during the snow-free season is mainly governed by the amount of absorbed solar radiation, usually leading to warmer air temperature than soil temperature (Aalto et al., 2018;Bartlett et al., 2006;Putnam & Chapman, 1996). At an interannual time scale, our correlative analyses show that trends in ΔT are significantly related to the increase in SolarRad over the period in northern Eurasia (Figures 6a-6d). At a depth of 0.2 m, trend in SolarRad is the third important variable for explaining ΔT 0.2 trend, which is at a similar level of trends in AMSD and rainfall, but not as important in the deeper ground (Figure 6e). The response curves for ΔT 0.8 and ΔT 1.6 are relatively flatter, whereas the decrease in SolarRad can explain the decrease in ΔT 0.2 (Figures 6b-6d), suggesting that a decrease in SolarRad contributes to less warming of topsoil (0.2 m) temperature than air, but such influence dwindles at the deeper ground.
Annual rainfall has generally increased except at a small group of sites in the southern and southwestern parts of the area, which is in line with the previous studies (Bintanja & Andry, 2017;Groisman et al., 2016). Trends in rainfall have significant and nonlinear relationships with ΔT 0.2 trends but not in the deeper ground (Figures 6b-6d). It is also noticeable that the response curve for trends of ΔT 0.2 is largely below the null (trends in ΔT = 0°C decade −1 ), indicating that changes in rainfall barely explain the increase in ΔT (Figure 6b). The impacts of altered rainfall on soil temperature are complicated in terms of timing, frequency, and intensity. Peng et al. (2016) concluded that changes in rainfall might both negatively or positively affect trends in soil temperature across multiple models. Although annual rainfall has a slight effect on CHEN ET AL.   annual soil temperature (Karjalainen et al., 2019;Knight et al., 2018), Zhang et al. (2001) suggested that the changes in rainfall during summer season may play a more significant role in controlling summer soil temperature at a depth of 40 cm than the changes in air temperature at Irkutsk, Russia. Summer precipitation is likely to have a cooling effect on soil temperature by increasing surface wetness and soil moisture, thus leads to more energy consumption for evaporation, while rainfall in prewinter may warm soil temperature through the effects of the greater freezing latent heat (Iijima et al., 2010;Park et al., 2013;Westermann et al., 2011;Zhang et al., 2001). In the warmer permafrost region (mean annual ground temperature > −0.5°C at the depth of zero annual amplitude), intensive precipitation in summer results in abrupt rising of soil temperature at shallow depth and an earlier thawing of the active layer, but reduces seasonal variations of soil temperatures and decreases the mean annual soil temperature (Luo et al., 2020). However, in the permafrost regions where the active layer is thin and always at or near saturation, an increase in rainfall during summer may have little impact (Zhang et al., 2001). Nevertheless, rainfall is a key factor controlling active layer thickness (Cao et al., 2017;James et al., 2019;Karjalainen et al., 2019), implying that the active layer may deepen in the future under a warming climate. Consequently, the influences of rainfall on soil temperatures are likely to be strengthened through more significant soil moisture feedback and runoff generation-outflow mechanisms (Eltahir, 1998;Gao et al., 2016;Luo et al., 2003).
Soil moisture has decreased in the seasonal frost region; however, it increased at a large group of sites in the permafrost regions ( Figures S3f-S3h). In the frost-related region where soil moisture remains relatively high, it has a greater influence on land-atmosphere heat exchange through freeze-thaw cycles (Cheng & Wu, 2007;Liljedahl et al., 2016;Martin et al., 2019). A study has shown that soil moisture in topsoil layers decreases during summer season (JJA) in the permafrost region (Teufel & Sushama, 2019); thus, an increase in soil moisture at the sites located in the permafrost region can be partially attributable to the increases of rainfall in prewinter season (SON) and active layer thickening (Andresen et al., 2020). Trends in soil moisture significantly and nonlinearly relate to ΔT trends at multiple depths, but thick confidence bands in response curves indicate a large number of uncertainties (Figures 6b-6d). An increase in vegetation leads to a decrease in soil temperature during summer when canopies are present, while the net effect of vegetation on soil temperature remains unclear (Loranty et al., 2018). At decadal scales, it is found that the changes in vegetation have little impact on trends in ΔT 0.2 and ΔT 0.8 , but significantly affect ΔT 1.6 trends (Figures 6b-6d). This further implies that greening of vegetation has a cooling effect on the mean annual soil temperature at 1.6 m and mean annual air temperature over the area, which is in line with a recent study in North America ( The coefficient of determination (R2) is reported with a statistical significance level (p) in the bracket. The significant relationships (p < 0.05) are in bold.

Conclusions
This study provides the first holistic pictures of long-term trends in soil-atmosphere temperature differences in response to the changing environment at a large spatiotemporal scale. By quantifying trends in ΔT and the role of the environmental alterations from 1981 to 2015 in northern Eurasia, we conclude: 1. The warming rate of soil temperature is generally less than that of air temperature at depths of 0.8 and 1.6 m except for soil temperature at 0.2 m in the seasonal frost region that increases more rapidly than air temperature 2. Inconsistent warming rates of the mean annual air and soil temperatures lead to interannual variations of ΔT, which are spatially heterogeneous across northern Eurasia 3. Changes in snow cover and soil moisture explain a significant portion of the positive trends in ΔT, suggesting that increases in ΔT are mainly associated with alterations in snow cover and soil moisture, especially by snow cover thickness changes. Meanwhile, decreases in ΔT relate to more factors, including the changes in snow cover duration, solar radiation, rainfall, and vegetation A sparse observation network in such a large area brings uncertainties and may limit implementation of the results. Thus, it is recommended to develop the current database in multilayer soil temperatures, particularly in the continuous permafrost region, and enhance observational networks on the dynamics of soil properties (e.g., soil carbon stock) (which are not visible in our study) to systematically improve our knowledge of mitigated soil-atmosphere temperature difference in a changing world. The findings can improve our understanding of how the changing environment would influence land-climate thermal interactions at decadal scales under future climate change. This is also of relevance and interest in multiple disciplines as temperature regimes at the land-atmosphere interface are fundamental to understanding of subsoil biogeochemical cycles, roots development, and infrastructure stability.

Data Availability Statement
The soil and air temperatures, and snow depth data can be downloaded from the following URL (http:// meteo.ru/english/climate/cl_data.php). The gridded data can be downloaded from ERA5-Land product (https://cds.climate.copernicus.eu/cdsapp#!/home).