Observed Decrease in Soil and Atmosphere Temperature Coupling in Recent Decades Over Northern Eurasia

Coupling of soil‐air temperature is fundamentally relevant to various biophysical and biogeochemical functions near the land surface. However, coupling at an interannual scale is less addressed than at finer timescales, such as diurnal and seasonal scales. Here, we show that interannual variability of soil‐air temperature coupling decreased significantly during 1984–2013 across 196 stations in northern Eurasia. Coupling at a depth of 0.2 m remained significant in only 53% of stations during 1999–2013, while it was significant in 80% of stations during 1984–1998. Decreased coupling is mainly attributable to air temperature change in the snow‐free season and snow cover retreat. Moreover, change in the extreme snow cover frequency also contributes to the decline, while other climate extremes are less influential. The findings can help us to understand long‐term land‐climate thermal interactions and to evaluate soil temperature profiles simulated by land surface models.

. On a seasonal basis, snow cover retreat can enhance β due to a reduced insulation effect in winter, whereas rising air temperature during the snow-free season may dampen β because abrupt air temperature change is likely to be filtered as depth increases in soil (Cheruy et al., 2017;Tingjun, 2005;Zhan et al., 2019). Furthermore, with significant increases in the number of hot days and nights since the 1950s and alterations in other climate extremes (IPCC, 2013), the status of the annual soil-air temperature tracking ability remains uncertain. Hence, it seems reasonable to suppose that compound environmental and climatic changes should mediate β, whereas there are so far only limited examples of the dynamics of β at an interannual timescale. Nevertheless, evidence from local and regional studies has hinted that the response of soil temperature to air temperature changes might not be constant over the long term (Fang et al., 2019;Romanovsky et al., 2007;Zhan et al., 2019;Way et al., 2018). For instance, Beltrami et al. (2003) and Chen et al. (2021a) both highlighted that long-term trends in air and soil temperature were inconsistent, leading to an interannual variation of air and soil temperature difference, which suggests that the coupling between air and soil temperature may not be stable.
Accessing variations of the tracking ability of soil-air temperature is critical for the parameterization of land surface models (LSMs). The soil temperature simulated by LSMs has been evaluated extensively and behaves differently from diurnal to decadal timescales. In particular, previous studies have pointed to the limitation in capturing the interannual variability of simulated soil temperature by various LSMs2018 (Decharme et al., 2016;Orth et al., 2017;Yang et al., 2020;Yang & Zhang, 2018). In this regard, assessing the interannual variability of β by using in situ, independently observed soil and air temperature data could provide new insights into the evaluation of simulated soil temperature over a long-term period. It is also imperative to test the stability of β to better understand land-atmosphere thermal interactions and to further assess related ecosystem functions such as soil carbon stock given other environmental changes in the future Fernández-Martínez et al., 2019;Loranty et al., 2018).
Here, we investigate the interannual variability of soil-air temperature coupling during 1984-2013 over northern Eurasia by analyzing paired observations of climatic variables. We began by showing that interannual variability of β in the region was not stable by running 15-year moving windows over the period. We then examined the potential mechanisms of decreased β.

Data Compilation
The study region encompasses most land in northern Eurasia at latitudes from 42° to 70°N and longitudes from 28° to 170°E (Figure 1a). Daily meteorological (including soil and air temperature, snow depth, and precipitation) data collected during 1984-2013 were extracted from the All-Russia Research Institute of Hydrometeorological Information-World Data Center. Extensive quality control checks were undertaken on all data sets before publication (Bulygina & Razuvaev, 2012;Sherstiukov, 2012), and only qualified data were used in our analyses. Furthermore, daily measurements were computed into annual attributes: MAAT and MAST (°C); freezing and thawing degree days (FDD and TDD, °C, a sum of daily degrees below and above 0°C for a year); annual rainfall (mm), annual mean of daily snow depth (AMSD, cm), and snow cover duration (SCD, a sum of the days on which snow cover exists for a year). A minimum of 300 days of data was required for computing annual attributes and, on average, more than 362 days of data in a year were available across data sets (Table S1). Apart from the in situ observations, we also extracted monthly gridded data (9 km resolution) from the ERA5-Land product (Hersbach et al., 2020), which were then computed into annual attributes: net surface solar radiation (SolarRad, W/m 2 ) and volumetric water content (VWC, %).

Analysis Method
β between MAAT and MAST was calculated as the slope of the ordinary least squares regression, in which MAST was regressed against MAAT for a specified period (Lenoir et al., 2017). A slope value close to one corresponds to high coupling, whereas low coupling has a slope value close to zero. First, we calculated β for each 15-year running window from 1984 to 2013. Then, interannual variation of β at each station was investigated using the linear trend of β s from the running windows (in total, 16 windows over the period). Before the regression analyses, all variables were detrended by removing linear trends for each window.
We also analyzed Pearson's correlation (R) between MAAT and MAST to test the robustness of the results (Benesty et al., 2009).
To examine the potential mechanisms of change in β, we first divided each 20-year running window from 1984 to 2013 into two 10-year groups based on the magnitude of the investigated variables (AMSD, FDD, TDD, etc.): one group consists of 10 years with greater values, while another group includes the remaining 10 years with smaller values. As such, each of the two groups has the same set of stations, and values of β s from the two groups are fully comparable. Then, the values of β s from the two groups were compared to gain insights into the potential mechanisms of decreased β throughout the region. Likewise, to examine the potential impacts of climate extremes, we compared β and its change in the two 10-year groups divided based on hot and cold, rainfall, and snow cover extremes. Following Seneviratne et al. (2014) method for computing annual extreme weather frequency (days per year), the days with extreme snow cover (Snow 90P , the number of days in a year with snow depth higher than the 90th percentile of the daily snow depth during 1984-2013), hot (Hot 90P ) and cold (Cold 90P ) weather, and rainfall (Rain 90P ) extremes were computed for analysis.
Over the period, in situ observations at some stations and depths were inconsistent. Thus, we focused on the depths of 0.2 and 1.6 m and excluded the stations with less than 24 years of qualified data from the analyses. Figure S1 shows that an average of 27.6 and 28.1 years of data are available across 181 and 183 stations pairing all variables at 0.2 and 1.6 m, respectively. In total, the study involves 4,998 and 5,148 annual observations at 0.2 and 1.6 m, respectively, during the period (Table S1) (1984-1998) and last (1999-2013) 15 years are reported accordingly. The 5th and 95th percentiles of variables are depicted with colored bands. Permafrost zonation (continuous permafrost region (extent of permafrost ≥ 90%), discontinuous permafrost region (0% < extent of permafrost < 90%) is derived from the International Permafrost Association map (Brown et al., 1997).

Decreased Mean Annual Air and Soil Temperature Coupling
The majority of stations had positive β during 1984-2013, except for two and six negative stations at depths of 0.2 and 1.6 m, respectively (Figure 1a). The spatially averaged β 0.2 , which was 0.48 ± 0.18 (mean ± standard deviation) for 1984-1998 (the first 15 years), decreased to 0.36 ± 0.25 for 1999-2013 (the last 15 years) ( Figure 1b). Likewise, at deeper ground, β 1.6 decreased from 0.30 ± 0.16 for the first 15 years to 0.19 ± 0.20 for the last 15 years. Across the stations, β decreased at more than 70% of the stations (Figure 2a). The trends in β varied spatially, and a significant (P < 0.05) decrease in β was uniformly observed in the western part of the region affected by the seasonal frost. However, compared with the seasonal frost area, there was a less CHEN ET AL.  clear pattern of β changes for the stations located in the permafrost area, which may be attributable to the more complex and significant land-atmosphere-snow cover interactions (Grundstein et al., 2005;Henderson et al., 2018;Jia et al., 2019).
In the region as a whole, β derived from the total 16 running windows significantly decreased at 0.009 and 0.007 per window at depths of 0.2 and 1.6 m, respectively (Figure 2b). Notably, 80% of stations had significant β 0.2 during 1984-1998, while β 0.2 remained significant at 53% of stations during 1999-2013. At 1.6 m, over half of significant stations for the first 15 years (67%) became insignificant for the last 15 years (33%). To test the robustness of decreased β, we alternatively examined R between MAAT and MAST. A decrease in R was also found, and the loss of significant stations from the first to the last 15 years was substantial ( Figure S2). Hence, all these results confirm that the recent weakening relationship between soil and air temperatures at the continental scale is not an artifact but a real phenomenon.

Key Mechanisms of Decreased β
The factors associated with decreased β may vary at different spatiotemporal scales and cannot be elucidated from statistical analyses alone. Here, we mainly focus on mechanisms of the compound effects of altered FDD, TDD, and snow cover conditions throughout the period. It was found that β for the group with greater TDD was systematically lower than that for the smaller TDD group (Figure 3). Notably, the decrease in β 0.2 for the greater TDD group (−0.0048/year) was more significant and the decrease was accelerating faster than for smaller TDD years (−0.0002/year), while the situation for β 1.6 was similar (Figures 3b and 3c). Although FDD decreased across most stations, there was less difference between βs for two groups divided based on FDD compared with TDD, indicating the insulating function of seasonal snow cover.
Snow has been widely addressed as a key in mediating soil-air temperature tracking behaviors at different spatiotemporal scales (e.g., Tingjun, 2005;Zhang et al., 2018;Way et al., 2018). β 0.2 for the group with greater AMSD (Group 1) was systematically lower than for the group with smaller AMSD (Group 2), while the difference in β 1.6 between the two groups was smaller compared with β 0.2 (Figures 3b and 3c). Unlike AMSD, there was no significant difference in β between the two groups divided based on SCD. This is partly because SCD is less informative than AMSD, which is jointly determined by SCD, timing, and daily depth variations. Furthermore, no major difference in β was found between two groups divided based on Rain and VWC ( Figure S3). From 1984 to 2013, changes in AMSD and SCD were inconsistent at some stations (Figure 3a), suggesting more frequent snow cover extremes at these stations. It is also noted that the change in snow cover extreme days (Snow 90P ) is spatially heterogeneous (Figure 4a). Hence, we also tested whether β varied with Snow 90P and found that β was, in general, lower for the greater Snow 90P group (Group 1) than for the smaller (Group 2) (Figures 4b and 4c). Although Hot 90P increased at most stations, there was no significant difference in β between groups divided based on Hot 90P . Likely, the differences in β between the two groups divided based on Rain 90P and Cold 90P were mostly insignificant.

Discussion
The results highlight that coupling of MAAT-MAST was not stable at an interannual timescale and decreased from 1984 to 2013 in northern Eurasia. Meanwhile, the number of stations with significant β reduced substantially during the period, which was related to the environmental changes. Among the investigated attributes, the changes in TDD, FDD, and AMSD played vital roles in the decreased β.
During the period, the air temperature in the snow-free season (TDD) significantly increased, which masked β due to the soil's thermal inertia and filtering functions with depth, such that MAST was unable to fully track the changes in MAAT (Figure 3) (Cheruy et al., 2017). Moreover, β for warmer snow-free seasons (TDD, Group 1) showed a more notable trend than for the cooler seasons (−0.0048 and −0.0002 per year, respectively), which may imply a decoupling consequence of the MAAT-MAST relationship if climate warming in spring and summer continues in the future (Blunden & Arndt, 2019). This sensitivity of β to the air temperature change during the thermal growing season is of great relevance for temperature-sensitive processes near the land-atmosphere interface, such as plant growth, soil carbon decomposition, and the hydrological cycle ( However, compared with TDD, FDD changes were less influential on the interannual variability of β due to the insulating function of seasonal snow cover in winter attenuating the sensitivity of β to FDD variations (Figures 3b and 3c) (Decharme et al., 2016;Tingjun, 2005). Nevertheless, it is found that β 1.6 tended to be significantly lower for the groups with smaller FDD than the groups with greater FDD after 2011 (representing the running window of 1992-2011 and afterward), which potentially implies the enhanced sensitivity of β to FDD. This rising sensitivity of β to FDD variations is closely related to the dynamics of snow cover characteristics (e.g., onset, duration, and thickness) and may continue with the ongoing snow cover retreat CHEN ET AL.
10.1029/2021GL092500 6 of 10  1999-2013 and 1984-1998. The sub-density plots depict the distribution of the variables' differences, while the black dashed and solid orange lines represent the null (meaning no change) and mean values, respectively. (b and c). Comparison of the changes in β in different climatic conditions. The x-axis shows the last year of each 20year running window (e.g., 2013 stands for a running window from 1994 to 2013). For each 20-year window, we divided the 20 years into two groups based on the magnitudes of TDD, FDD, AMSD, SCD, and SolarRad: one group consists of 10 years with greater values (red, Group 1), while the other group (blue, Group 2) includes the remaining 10 years with smaller values. The difference between the two groups was assessed by the paired t-test and indicated with asterisks (P < 0.05) at the top. The numbers in b and c are linear regression slopes, and the asterisks indicate the significance level of the regressions (one for P < 0.05, two for P < 0.001).
The frequency and intensity of extreme weather events (such as extreme precipitation, hot temperature) have increased at some stations ( Figure 4a) (IPCC, 2013;Stott, 2016). The extreme snow cover events showed significant impacts on β, while other extreme attributes were less influential (Figures 4b and 4c). Of note, the rate of decrease in β for the greater Snow 90P group (Group 1) was faster than for the smaller group (Group 2), implying a potential loss of the significant relationship between MAST and MAAT if more frequent snow cover extreme events emerge in the future.
The differences in β between groups divided based on other investigated factors were relatively minor, suggesting that the decrease in β at interannual timescales was mainly attributable to the compound influences of spring and summer air temperature warming, and snow cover changes (Figures 3 and S3). Nevertheless, CHEN ET AL.   1999-2013 and 1984-1998. The sub-density plots depict the distribution of the differences, while the black dashed and solid orange lines represent the null (meaning no change) and mean values, respectively. (b and c). Comparison of the changes in β in different extreme weather conditions. The x-axis shows the last year of each 20-year running window (e.g., 2013 stands for a running window from 1994 to 2013). For each 20-year window, we divided 20 years into two groups based on the magnitudes of Hot 90P , Cold 90P , Snow 90P , and Rain 90P : one group consists of 10 years with greater values (red, Group 1), while the other group (blue, Group 2) includes the remaining 10 years with smaller values. The difference between the two groups was assessed by the paired t-test and indicated with asterisks (P < 0.05) at the top. The numbers in b and c are linear regression slopes, and the asterisks indicate the significance level of the regressions (one for P < 0.05, two for P < 0.001).
it is worth being aware that other factors (such as vegetation, rainfall, and soil moisture) can be essential in reducing the soil-air temperature relationship at finer spatial and temporal scales Fisher et al., 2016;Luo et al., 2020). However, we note that the use of less than 200 stations may hinder the mechanisms that underlie the observed weakening in β over regions with heterogeneous biogeophysical and climatology conditions, such as landscape, vegetation, and soil characteristics (Groisman et al., 2017;Monier et al., 2017). Nevertheless, snow cover retreat and warming in winter have been demonstrated to enhance β but warming in spring and summer to attenuate β. Since the mean air temperature is projected to increase both during snow-free and winter seasons incorporating snow cover retreat (Callaghan et al., 2011;Meehl & Tebaldi, 2004;Russo et al., 2014), the future dynamics of β are largely dependent on the net impacts of these factors. As such, compound effects of the individual factors on the variation of β will be uncertain under a changing climate. Thus, a continuous and effective observation network, more robust LSMs, and further in-depth studies on the dynamics of land-atmosphere thermal interactions are needed.

Summary and Conclusions
Multilayered soil temperature and paired air temperature observations provide evidence of an overall weakening interannual variability of air-soil temperature coupling during 1984-2013 over northern Eurasia. Such a decline mainly relates to an air temperature increase in spring and summer, and snow cover retreat, but we also note that winter air temperature warming appears to play a positive role in β incorporating snow cover retreat since the latter part of the study period.
The distribution of observation stations is irregular, especially over the high Arctic permafrost area, where the pattern of β change is not as clear as over the seasonal frost area, which may cause uncertainties in the results. Meanwhile, it remains uncertain whether the observed weakening in the interannual relationship of soil-air temperature reflects decadal variations of the relation to climate warming and other environmental shifts because of the relatively short run of the analyses. Further, some studies concluded that no simulated soil temperature outperforms across LSMs and time spans due to the distinct structures and parameterization schemes of the models (Rodell et al., 2004;Wang et al., 2016). Hence, this work also can be referred to as an indirect way of evaluating simulated soil temperature and the soil-air temperature relationship at an interannual scale, which will help in assessing long-term land-climate thermal interactions, which will be relevant to various biogeochemical and biophysical processes near the land-atmosphere interface in the future climate. Moreover, with the wide-ranging decline in usable weather stations at regional and global scales since the 1980s (Lawrimore et al., 2011;Sun et al., 2018), all our results demonstrate the value of longterm in situ soil temperature observations for understanding land-climate thermal interactions.