Increasing Snow–Soil Interface Temperature in Farmland of Northeast China from 1979 to 2018

The presence of seasonal snow cover in the cold season can significantly affect the thermal conditions of the ground. Understanding the change of the snow–soil interface temperature (TSS) and its environmental impact factors is essential for predicting subnivean species changes and carbon balance in future climatic conditions. An improved Snow Thermal Model (SNTHERM) is employed to quantify TSS in farmland of Northeast China (NEC) in a 39-year period (1979–2018) firstly. This study also explored the variation tendency of TSS and its main influencing factors on grid scale. The result shows that annual average TSS and the difference between TSS and air temperature (TDSSA) increased rapidly between 1979 and 2018 in the farmland of NEC, and we used the Mann–Kendall test to further verify the increasing trends of TSS and TDSSA on aggregated farmland of NEC. The correlation analysis showed that mean snow depth (MSD) is the most pivotal control factor in 95% of pixels and TDSSA increases as MSD increases. Snow depth can better predict the change of TSS in deep–snow regions than average winter temperature (TSA). The results of this study are of great significance for understanding the impact of snow cover on the energy exchange between the ground and the atmosphere in the cold climate.


Introduction
Snow cover has a great impact on the thermal conditions of the soil due to its low thermal conductivity and high albedo [1,2]. Snow cover maintains warmer soils allowing infiltration and gas exchange during winter but also leads to permafrost thawing and increased greenhouse gas emissions and further alters the continental global carbon budget [3,4]. The presence of snow cover causes the reduction of soil heat loss and the mitigation of soil-frost barrier formation, which significantly affects the shallow soil temperature and protects plant roots and soil microbial populations [5][6][7]. However, climate change has altered the duration of snow cover and snow depth [8,9]. Snow cover in the Northern Hemisphere has significantly changed [10,11].
Surface ground temperature (T SG ) is the result of physical and biological couplings between the atmosphere and the ground, and snow cover is an important factor affecting the change of T SG [12,13]. Some authors have defined the snow-soil interface as "subnivean space" or "subnivean environment" [14,15]. The thermal insulation effect of snow cover can prevent species from being affected by extreme winter temperatures. The change of the snow-soil interface temperature (T SS ) affects the severity of soil freezing and the stability of the soil carbon pool [7,16]. Understanding the response of T SS is essential for climate change. Snow cover affects the coupling between the air temperature (T NSA ) and T SG due to its insulating effect [17]. T SG is often higher than T NSA ; for example, a maximum difference of 36.3 K was observed in the winter of 1988 at the Franklin Bluffs area in northern Alaska [7]. Goncharova et al. [18] found that winter T SG has a great spatial variability and a significant negative correlation with altitude. Furthermore, altitude also affects the spatial distribution and thickness of snow cover [19,20].
The response of snow cover to climate change in China is complicated, and the number of snow cover days and snow depth have large geographical heterogeneities [21,22]. It is difficult to carry out field observations to obtain long-term and large-area observational data. T SS has only been observed at Chinese meteorological stations since 2005, and it is thus difficult to analyze the changes in T SS against the background of climate change with global warming. Therefore, an effective alternative is to perform numerical simulations [23,24]. At present, some authors use numerical simulations to simulate the long-term trend of T SG and T NSA in the context of climate change. For example, 32 general circulation models (GCMs) are used to simulate the past and future relationship between T SG and T NSA in North America. In the CMIP5 (Coupled Model Intercomparison Project Phase 5) model, the relationship between T SG and T NSA shows great variability and inconsistency in GCM simulations on annual and seasonal timescales by the comparison of simulated and observed values [25]. Bartlett [17] established a numerical model of snow-ground thermal to reveal the effect of seasonal snow cover on the seasonal and annual variations in T SG . It was found that snow cover caused T SG to be higher than T NSA . Additionally, Zhang et al. [26] developed a multi-layer snow-soil scheme at EALCO (Ecological Assimilation of Land and Climate Observations) to study the interaction between snow and soil in the cold season and found that snow cover significantly affected surface energy flux, soil temperature, and soil freezing depth. These studies have determined the response of snow cover on the ground by analyzing changes in the depth of soil freezing and soil temperature by model simulation or in situ observation at a point scale; however, it is not clear how T SS changes at a regional scale in the context of climate change. Snow Thermal Model (SNTHERM) is one of the most advanced multi-layer thermodynamic models available at present and can be used to predict energy conversion between snow and soil. The model can automatically identify and process snow melting and accumulation processes, and some studies have proved its excellent simulation accuracy [27,28].
Northeast China (NEC) is regarded as the second-largest stable snow area in China, with the number of snow cover days being 120-180 d. In this region, the snow depth and the number of snow days are high and the average winter snow depth and maximum snow depth are increasing [9,22]. In the present study, to understand the change of T SS in the farmland of NEC, T SS with a 10 km grid resolution in farmland of NEC from 1979-2018 is acquired through the improved SNTHERM model. The spatiotemporal variation of T SS is analyzed and its main influencing factors are analyzed. NEC is the largest commercial grain production base in China. The thermal insulation effect of snow cover creates a suitable environment for overwintering crops and various soil microbes; therefore, it is of great significance to understand the variation of T SS in farmland of NEC, and our study can provide support for agriculture and animal husbandry production in NEC.

Meteorological Station Data in NEC
NEC is located between 38.42 and 53.36 • N and between 115.24 and 135.12 • E, with a total area of 1.24 million km 2 . The land cover is dominated by farmland, forest, and grassland. The three major mountain ranges from west to east are the Greater Khingan Range, the Lesser Khingan Range, and the Changbai Mountains. The farmland in NEC is mainly distributed on the Northeast China Plain. This plain, one of only three black soil regions in the world, consists of the Sanjiang Plain, Songnen Plain, and lower Liaohe River Plain. The Northeast China Plain has a large, cultivated area and wide distribution, which is an important grain, soybean, and animal husbandry production base in China. The land area of NEC accounts for 8.2% of China, with the cultivated land area accounting for 13.3% of the national cultivated land area and the grain output accounting for 20.0% of China's total grain output [29]. The climate of NEC is continental monsoon, and the average annual rainfall is 460 mm. Winter is long and cold-arid, while summer is short and warm. The annual average temperature of NEC is below 0 degrees. Additionally, NEC is one of the three stable snow cover regions in China. The snow cover period lasts from November to April of the following year. The region contains abundant snow resources and is of great significance for agricultural irrigation.
The daily datasets of China's surface climate provided by the National Meteorological Science Data Center (air pressure, ground temperature, snow depth, snow pressure, wind speed) (http://data.cma.cn/, accessed on 1 August 2021) were used in our study. Platinum resistance thermometers were used to observe T SS at the Meteorological stations. The observation surface condition is bare soil: half of the temperature sensor is buried in the soil and closely bonded to the soil and half is exposed on the ground. Finally, the measured minute data is averaged to obtain the daily average T SS [30]. Data between 2005 and 2018 from 36 stations whose underlying surface is farmland were selected from among 98 stations in NEC for this study (Figure 1). important grain, soybean, and animal husbandry production base in China. The l of NEC accounts for 8.2% of China, with the cultivated land area accounting for the national cultivated land area and the grain output accounting for 20.0% o total grain output [29]. The climate of NEC is continental monsoon, and the ave nual rainfall is 460 mm. Winter is long and cold-arid, while summer is short an The annual average temperature of NEC is below 0 degrees. Additionally, NEC the three stable snow cover regions in China. The snow cover period lasts from N to April of the following year. The region contains abundant snow resources great significance for agricultural irrigation.
The daily datasets of China's surface climate provided by the National Met cal Science Data Center (air pressure, ground temperature, snow depth, snow p wind speed) (http://data.cma.cn/, accessed on 1 August 2021) were used in ou Platinum resistance thermometers were used to observe TSS at the Meteorological The observation surface condition is bare soil: half of the temperature sensor is the soil and closely bonded to the soil and half is exposed on the ground. Fin measured minute data is averaged to obtain the daily average TSS [30]. Data betw and 2018 from 36 stations whose underlying surface is farmland were selec among 98 stations in NEC for this study (Figure 1).

Forcing Datasets
Due to the lack of long-term and large-scale meteorological observation da good alternative to use a ground meteorological forcing dataset as the driving d model. The China Meteorological Forcing Dataset (CMFD) was made through remote sensing products, reanalysis dataset, and in situ observation data at wea tions, with high time-resolution (3-hour step, [31]). CMFD from 1979 to 2018 in

Forcing Datasets
Due to the lack of long-term and large-scale meteorological observation data, it is a good alternative to use a ground meteorological forcing dataset as the driving data of the model. The China Meteorological Forcing Dataset (CMFD) was made through fusion of remote sensing products, reanalysis dataset, and in situ observation data at weather stations, with high time-resolution (3-hour step, [31]). CMFD from 1979 to 2018 include air temperature, surface pressure, specific humidity, wind speed, downward shortwave radiation, downward longwave radiation, and precipitation rate (http://www.geodata.cn, accessed on 1 August 2021).
GLDAS-2.0 Data and GLDAS-2.1 Data of 0.25-degree spatial resolution and 3-h time resolution were also used, which provides a temporally consistent series from 1948 to 2014 of GLDAS-2.0 Data and 2000 to the present of GLDAS-2.1 Data. Soil temperature and humidity data of GLDAS-2.0 Data from 1979 to 2013 and of GLDAS-2.1 Data from 2013 to 2018 were used in this paper.

Snow Depth Dataset
In this study, we used a snow depth dataset (http://data.tpdc.ac.cn/zh-hans/, accessed on 1 August 2021) spanning the period from 1979 to 2018 to retrieve the change of the average snow depth (M SD ) in NEC. Daily passive microwave brightness data (SMMR, SSM/I, SSMI/S) were used for the inversion of snow depth in this dataset. The brightness temperature of different sensors was inter-calibrated before the derivation of snow depth to ensure the temporal consistency of the brightness temperature [32]. The dataset was based on the modified Chang algorithm to retrieve snow depth to obtain the daily snow depth in China with a spatial resolution of 0.25 degrees [20].

Land Cover
The land cover between 2000 and 2015 was used in this paper, which was retrieved from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/, accessed on 1 August 2021) with a 1-km spatial resolution. Land cover was mainly interpreted from Landsat-TM/ETM and Landsat 8 remote sensing images to distinguish farmland, forest, grassland, and water. In order to match the regional reanalysis data, land cover were regridded into 0.1 • cells in this study.

SNTHERM Model Description and Calibration
SNTHERM is a one-dimensional mass and energy balance model developed by the US Army Cold Regions Research and Engineering Laboratory (CRREL) to predict the temperature distribution of snow and frozen soil [33]. The model is intended for seasonal snow covers and addresses conditions found throughout the winter, from initial ground freezing in the fall to snow ablation in the spring. This model is suitable for the simulation of various winter meteorological conditions such as snowfall, sleet, rainfall, freeze-thaw cycles, and snow cover. Patankar's volume control method is adopted in SNTHERM to subdivide the snow layer and the soil layer into horizontal infinite layers (each layer is considered as a node) to obtain numerical solutions (Figure 2), each of which conforms to energy-conservation equation and mass-conservation equation [34,35].  The sum of the constituent bulk densities is the density  of the total mediu ten as: Snow and soil are both porous media, the unfrozen water coexists in equilibrium with ice at temperatures below 0 • C. In the SNTHERM code, mass liquid water fraction ( f l [kg·m −3 /kg·m −3 ]) within a given medium is a single-valued function of temperature and has a freezing curve characteristic of the snow-soil properties in SNTHERM. Thus, the semi-empirical function is introduced for each node in SNTHERM for f l [33]: where γ l is the bulk density of liquid water (kg·m −3 ); γ w is the bulk density of the combine liquid water and ice constituents (kg·m −3 ); J p is the plasticity index; γ d is the soil density (kg·m −3 ); a 1 and a 2 are defined as 0.2/ 0.01 + J p and 0.01/ 0.1 + J p , respectively; T D is the depression temperature (K), defined as 273.15-T, and T is the node temperature (K). γ l and γ w are further specified in the SNTHERM model: where k represents v, l, i, a, or d for water vapor, liquid water, ice, air, or the dry soil solids, respectively. θ k is the volume fraction of constituent k (m 3 ·m −3 ), ρ k is the density of constituent k, ρ a is the density of air at 0 • C and 1000 mb (1276 kg·m 3 ), ρ i is the density of ice (917 kg·m 3 ), and ρ l is the density of liquid water (1000 kg·m 3 ). Taken over the five possible constituents in the medium, the sum of the volume fraction is unity, or The sum of the constituent bulk densities is the density ρ of the total medium, written as: A melt function P melt is introduced in SNTHEM to obtain the node temperature, defined as the change in γ l with temperature over the time step in which the γ w is held Agriculture 2021, 11, 878 6 of 17 constant at the time of t. A melt function P melt is further derivatived as follows according to Equation (1): where F is the temporal average on ∆t of F(F = ∂ f l /∂T) in temperature evaluated for the current water content on ∆t; Noda temperature is then expressed in SNTHEM: where T t is the node temperature at time of t, g v = 1/γ w t F, and g v = T t−∆t . J p is an important feature that characterizes the physical properties of soil and directly affects the accuracy of soil and snow temperature simulation from Equations (1) and (5). The J p for snow is 0 and the ranges from 0 to 0.3 for soils. The default J p of SNTHERM is 0.2. The J p of the model is modified to suit the soil characteristics of the farmland in NEC. The calibrated model parameters are shown in Table 1. The input data required by the SNTHERM are three parts: meteorological drive data, initial soil and snow parameters, and empirical parameters. The outputs of the SNTHERM are the predicted T SS and predicted surface ground temperature (T ST ). There is no snow cover when T SS is equal to T ST in the model. Therefore, we obtained actual T SS when the tolerance between T SS and T ST is greater than 0 • C and Snow depth greater than 0 (Equation (7)).
The original SNTHERM only can simulate at a single point. In this study, multi-point simulation can be carried out by improving the model. Finally, the 10 KM resolution map of T SS is formed. The following steps describe our simulation process: (1) Grid-by-grid reading of temperature and precipitation data from September 1 to May 31 from ITPCAS forcing data (no snowfall will occur in NEC from June to August). It is considered to snowfall when the temperature is less than 0 • C and the precipitation rate is greater than 0 mhr −1 . First day of snowfall (FDS) was marked as the initial runtime of the model.
(2) At each grid, set the number of soil layers to 6 layers, initial soil temperature and humidity at −5 cm, −10 cm, −15 cm, −30 cm, −50 cm, −80 cm is obtained by reading GLDAS forcing data on the FDS. The number of snow cover layers on the FDS is set to be 2 layers, in which thicknesses are initialized at 0.01 and 0.02 and the bulk water density are 150 kg·m 3 respectively. The initial snow grain diameter is set to 0.0005 m. Furthermore, meteorological parameters with 3 h step (air temperature, relative humidity, wind speed, incident solar radiation, reflection solar radiation, precipitation) in ITPCAS between FDS and May 31 were read grid by grid, where reflected solar radiation are the default value of 9999 due to missing data.
(3) Combined with the soil texture of NEC and the verification analysis of the model, the J p of SNTHERM was set to 0.28 and other general parameters adopt the model default values. SNTHERM model is run year by year to obtain the snow and soil parameters of the farmland of NEC from the FDS to May 31 with a resolution of 10 km with 3 h step.
(4) Superimposing land cover information to extract regions that have always been farmland in NEC from 1979 to 2018, we obtained the long-term predicted temperature and other parameters in the whole farmland of NEC.
BIAS and root mean square error (RMSE) were used to assess the performance of the SNTHERM model in NEC. The calculation formula is as follows: where T mod,n represents the simulated T SS by the model on the n-th day, T mea,n represents the observed value at the meteorological station on the n-th day, and N is the total number of days that snow cover exists.

SNTHERM Model Validation
The SNTHERM model performances are shown Table 1. We firstly used simulated T SS of 36 weather stations from 2005 to 2018 with a three-hour step. In order to match the measured T SS , we process the simulated T SS into the daily average T SS and verify the accuracy of the simulated data and the measured data at the site. Figure 3a shows the scatterplots of the measured daily average Tss and model simulation results of 36 farmland weather stations from 2005 to 2018. The overall BIAS value is 0.43 • C, and RMSE value is 3.48 • C. In most stations, the BIAS of the simulated T SS is higher than that of in situ observations, with BIAS ranging from −2.32 to 1.65 • C. Our results also indicate nonsignificant relationships between snow depth and the simulation accuracy of T SS , while the accuracy of T SS is more affected by land cover type. The larger RMSE values are mainly concentrated in regions with mixed pixels, where land cover is predominantly forest and grassland with less farmland. The model has a good performance in other regions (mainly concentrated in the Sanjiang Plain, Songnen Plain, and lower Liaohe River Plain), where the underlying surface is homogeneous farmland and is less affected by other factors. There are many positive T SS in both simulations and observations in Figure 3a, which mainly concentrated during periods of snow ablation. When the snow cover is melting, the T SS will remain at about 0 • C, and more than freezing point occurs which leads to soil cooling. The calculation results of the Pearson correlation coefficient show that the trends of year-to-year fluctuations between the model and the actual measurement are consistent during 2005 to 2018 (Figure 3b). The improved SNTHERM model has a superior applicability in NEC and can simulate the trend of T SS in the Northeast.

Spatiotemporal Variations of TSS and TDSSA
The accumulative total of TSS and TDSSA (the difference between TSS and air temperature) from September 1 to May 31 is divided by the number of days when snow exists to obtain annual average TSS and TDSSA respectively.  Figure 4a). Similarly, it was found that TDSSA increased in 98% of the farmland in NEC, with 69% of the pixels passing the significance test (p < 0.05). A rapidly increasing trend of TDSSA also occurred in most pixels of the Sanjiang Plain and Songnen Plain, with a growth rate between 0.06 and 0.22 °C /a (Figure 4b). In the lower Liaohe River Plain, a lower and non-significant growth occurred in TSS (0.00-0.06 °C /a) and TDSSA (0.00-0.06 °C /a). Both TSS and TDSSA showed a non-significant decline in the southern coastal area and the west of NEC (Figure 4a,b). Additionally, the Mann-Kendall (M-K) trend analysis method further proved that an overall increasing trend of TSS and TDSSA occurred in the farmland of NEC at a confidence level of p < 0.05, with the linear trend indicating rapid increases at a rate of 0.10 and 0.09 °C /a, respectively, (Figure 4c,d).

Spatiotemporal Variations of T SS and T DSSA
The accumulative total of T SS and T DSSA (the difference between T SS and air temperature) from September 1 to May 31 is divided by the number of days when snow exists to obtain annual average T SS and T DSSA respectively. T SA and M SD were defined as the average value from December 1st to March 31st. The linear trend of the annual average T SS and T DSSA in NEC from 1979 to 2018 were showed in Figure 4. 99% of T SS increased in the farmland of NEC, 80% of which passed the significance test (p < 0.05). The Sanjiang Plain and Songnen Plain have a significant and rapid positive increase of T SS (0.06-0.21 • C/annual [ • C/a], Figure 4a). Similarly, it was found that T DSSA increased in 98% of the farmland in NEC, with 69% of the pixels passing the significance test (p < 0.05). A rapidly increasing trend of T DSSA also occurred in most pixels of the Sanjiang Plain and Songnen Plain, with a growth rate between 0.06 and 0.22 • C/a (Figure 4b). In the lower Liaohe River Plain, a lower and non-significant growth occurred in T SS (0.00-0.06 • C/a) and T DSSA (0.00-0.06 • C/a). Both T SS and T DSSA showed a non-significant decline in the southern coastal area and the west of NEC (Figure 4a,b). Additionally, the Mann-Kendall (M-K) trend analysis method further proved that an overall increasing trend of T SS and T DSSA occurred in the farmland of NEC at a confidence level of p < 0.05, with the linear trend indicating rapid increases at a rate of 0.10 and 0.09 • C/a, respectively, (Figure 4c,d).  Deep snow has better thermal insulation effects on the ground in the Sanjiang Plain (R 2 = 0.35, p < 0.01, Figure 4b). Therefore, although TSA non-significant decreased from 1979 to 2018 (R 2 = 0.01, p = 0.92), TSS still exhibited a rapid increasing trend in the Sanjiang Plain (R 2 = 0.46, p < 0.01, Figure 4a). Although TSS and TSA of the Songnen Plain are increasing rapidly, TDSSA still increased rapidly (R 2 = 0.35, p < 0.01, Figure 4d), the changes of TSS are less affected by TSA accordingly. In the lower Liaohe River Plain, TSS and TDSSA both increased at a slow rate (TSS: R 2 = 0.29, p < 0.01; TDSSA: R 2 = 0.13, p < 0.01; Figure 5e,f).We also found snow cover has a very weak thermal insulation effect on the ground and even has a cooling effect at certain times due to the shallow snow. Table 2 further reveals the impacts of TSA and MSD on TSS in the Sanjiang Plain, Songnen Plain, and lower Liaohe River Plain. The main reason for the change in TSS in the Sanjiang Plain is the MSD. Although TSA decreased slowly, TSS underwent a significant increase due to the large contribution of MSD to TSS between 1979 and 2018. Additionally, both TSA and MSD made a great contribution to the change in TSS in the Songnen Plain, where TSS experienced the fastest increase under the impact of increasing TSA and MSD. A rapid increasing trend of TDSSA also occurred in areas with deep snow in the Sanjiang Plain and Songnen Plain (Figures 4b and 5b,d). In the lower Liaohe River Plain, TSA increases Deep snow has better thermal insulation effects on the ground in the Sanjiang Plain (R 2 = 0.35, p < 0.01, Figure 4b). Therefore, although T SA non-significant decreased from 1979 to 2018 (R 2 = 0.01, p = 0.92), T SS still exhibited a rapid increasing trend in the Sanjiang Plain (R 2 = 0.46, p < 0.01, Figure 4a). Although T SS and T SA of the Songnen Plain are increasing rapidly, T DSSA still increased rapidly (R 2 = 0.35, p < 0.01, Figure 4d), the changes of T SS are less affected by T SA accordingly. In the lower Liaohe River Plain, T SS and T DSSA both increased at a slow rate (T SS : R 2 = 0.29, p < 0.01; T DSSA : R 2 = 0.13, p < 0.01; Figure 5e,f).We also found snow cover has a very weak thermal insulation effect on the ground and even has a cooling effect at certain times due to the shallow snow. Table 2 further reveals the impacts of T SA and M SD on T SS in the Sanjiang Plain, Songnen Plain, and lower Liaohe River Plain. The main reason for the change in T SS in the Sanjiang Plain is the M SD . Although T SA decreased slowly, T SS underwent a significant increase due to the large contribution of M SD to T SS between 1979 and 2018. Additionally, both T SA and M SD made a great contribution to the change in T SS in the Songnen Plain, where T SS experienced the fastest increase under the impact of increasing T SA and M SD . A rapid increasing trend of T DSSA also occurred in areas with deep snow in the Sanjiang Plain and Songnen Plain (Figures 4b and 5b,d). In the lower Liaohe River Plain, T SA increases faster than the Sanjiang Plain and Songnen Plain while M SD rises slowly. Under the influence of T SA and M SD , the T SS and T DSSA in the lower Liaohe River Plain increased slowly between 1979 and 2018 (Figure 5e,f, Table 2).  Table 2).

The Response of TSS to TSA and MSD
We selected MSD and TSA as explanatory variables and calculated the correlation co efficient between TSS and TSA as well as MSD to analyze the factors that cause TSS to chang at the grid-scale ( Figure 6). The majority of pixels (83%) exhibited a positive trend of MS in farmland, with a faster growth trend (0.06-0.31 cm/a) in the northeast of NEC (mainl the Songnen Plain and Sanjiang Plain). We observed a decreasing trend in MSD in th southeast (including Parts of the Lower Liaohe River Plain), north (mainly forest), an west (mainly grassland) of NEC (Figure 6a). The long-term change of TSA in farmland o

The Response of T SS to T SA and M SD
We selected M SD and T SA as explanatory variables and calculated the correlation coefficient between T SS and T SA as well as M SD to analyze the factors that cause T SS to change at the grid-scale ( Figure 6). The majority of pixels (83%) exhibited a positive trend of M SD in farmland, with a faster growth trend (0.06-0.31 cm/a) in the northeast of NEC (mainly the Songnen Plain and Sanjiang Plain). We observed a decreasing trend in M SD in the southeast (including Parts of the Lower Liaohe River Plain), north (mainly forest), and west (mainly grassland) of NEC (Figure 6a). The long-term change of T SA in farmland of NEC is complicated; however, the majority of pixels (78%) show an increasing trend. As shown in Figure 6b, the areas with a rapid increase of T SA are mainly concentrated in the southern coastal area of NEC and the southeast and southwest of the Songnen Plain (0.04-0.09 • C/a). In the middle of the Songnen Plain and the northern Sanjiang Plain, T SA had a decreasing trend, especially in parts of the Sanjiang Plain, where a rapid decline of T SA was observed (−0.08-−0.02 • C/a). NEC is complicated; however, the majority of pixels (78%) show an increasing trend. As shown in Figure 6b, the areas with a rapid increase of TSA are mainly concentrated in the southern coastal area of NEC and the southeast and southwest of the Songnen Plain (0.04-0.09 °C /a). In the middle of the Songnen Plain and the northern Sanjiang Plain, TSA had a decreasing trend, especially in parts of the Sanjiang Plain, where a rapid decline of TSA was observed (−0.08-−0.02 °C /a). The results show that TSA and TSS have a widespread negative correlation. For example, in the Songnen Plain and some parts of the Sanjiang Plain, TSS unexpectedly exhibits an increasing trend when TSA decreases (Figures 4a and 6b). Therefore, we analyzed the relative importance of TSA and MSD to the change in TSS in each pixel (Figure 7, Table 3) according to the approach of Zhu et al. [36]. The results show that MSD was more important to the change of TSS; with increasing MSD, TSS also increased, especially in the Sanjiang Plain (95% of increasing TSS) and Songnen Plain (99% of increasing TSS). At the same time, TSS and TDSSA in these regions underwent rapid growth (0.06-0.22 °C /a, Figure 4a,b). The main reason for this is that the significant influence of the rapid increase of MSD (0.06-0.31 cm/a) counteracts the negative effect of the decrease in TSA on TSS, which makes TSS and TDSSA experience a rapid increasing trend (Figure 4 and 7c). There is a stronger correlation between TSS and TSA in part of the lower Liaohe River Plain (Table 3). When a decrease in MSD cooccurred with an increase in TSA, an increase of TSS was observed. The pixels for which increasing TSA had a stronger effect on TSS offset the slow decrease of MSD, causing the TSS to increase at a slower rate in these areas (0.00-0.06 °C /a, Figures 4a, 6, and 7c). Moreover, for the pixels in the northwest of NEC and the southern coast of NEC, approximately 1% of TSS decreased between 1979 and 2018 (Figure 4a), where the decrease of MSD and TSA simultaneously appears. The results show that T SA and T SS have a widespread negative correlation. For example, in the Songnen Plain and some parts of the Sanjiang Plain, T SS unexpectedly exhibits an increasing trend when T SA decreases (Figures 4a and 6b). Therefore, we analyzed the relative importance of T SA and M SD to the change in T SS in each pixel (Figure 7, Table 3) according to the approach of Zhu et al. [36]. The results show that M SD was more important to the change of T SS ; with increasing M SD , T SS also increased, especially in the Sanjiang Plain (95% of increasing T SS ) and Songnen Plain (99% of increasing T SS ). At the same time, T SS and T DSSA in these regions underwent rapid growth (0.06-0.22 • C/a, Figure 4a,b). The main reason for this is that the significant influence of the rapid increase of M SD (0.06-0.31 cm/a) counteracts the negative effect of the decrease in T SA on T SS , which makes T SS and T DSSA experience a rapid increasing trend (Figures 4 and 7c). There is a stronger correlation between T SS and T SA in part of the lower Liaohe River Plain (Table 3). When a decrease in M SD co-occurred with an increase in T SA , an increase of T SS was observed. The pixels for which increasing T SA had a stronger effect on T SS offset the slow decrease of M SD , causing the T SS to increase at a slower rate in these areas (0.00-0.06 • C/a, Figures 4a, 6, and 7c). Moreover, for the pixels in the northwest of NEC and the southern coast of NEC, approximately 1% of T SS decreased between 1979 and 2018 (Figure 4a), where the decrease of M SD and T SA simultaneously appears.

The Environmental Impact of T SS Changes
We simulated T SS from 1979 to 2018 using SNTHERM and analyzed its relationship with several environmental factors. Some studies considered that winter air temperatures are the primary effects on soil frost. In addition, these studies are conducted at several locations over the short term [37][38][39]. Our analysis indicates that M SD and T SA have an important control on the change T SS in the farmland of NEC and that the impact of M SD on T SS is more significant at larger scales. The results of this study show that T SA is not a good measure of climate change in deep-snow regions, where the M SD affected T SS more. Due to the effect of snow cover on the ground-atmosphere transmission, it changes the surface energy balance [18]. In deep-snow regions, the high reflectivity of snow cover to solar radiation causes the ground to absorb less heat, and the low thermal conductivity of snow hinders surface heat dissipation, causing T SS to be more affected by snow cover. Similarly, in shallow snow regions, snow cover reduces the absorption of solar shortwave radiation, allowing snow cover to cool the underlying soil ( Figures 5 and 7).
Snow cover is a natural insulator and reduces the heat exchange between the soil and the atmosphere [18]. Gold [40] concluded that snow cover is the principal reason that the mean annual ground temperature can be many degrees warmer than the mean annual temperature in cold regions. Additionally, the difference between ground temperature and air temperature during the frozen period with snow cover is greater than that for the period of frozen ground without snow [36]. The thickness of the snow cover greatly affect the energy exchange on the ground [41]. Some studies have shown an increasing trend of the annual cumulative snow depth and maximum snow depth in NEC [22,42]. Therefore, it is very essential to understand the changes in T SS .
Our study found that the M SD in the farmland of NEC increased from1979 to 2018 in the majority of pixels (Figure 6a). This will alter the radiation balance of the earthatmosphere system and, which can explain the increasing trend of T SS in this region and may exacerbate the degradation process of permafrost in NEC ( Figure 5, Table 2). Furthermore, T SS will affect the evolution of snow parameters in turn, which will lead to the continuous change of snow thermal conductivity, thus affecting the exchange of material and energy between atmosphere and snow cover and soil [43].
Zhang [7] reported that the heat preservation of snow cover gradually increases with increasing snow depth before the snow cover reaches the optimal thickness. The presence of snow cover has a significant impact on soil frost penetration [44]. We also found that snow cover has a good insulation effect when the M SD of Songnen Plain and Sanjiang Plain is 5.02 cm and 8.54 cm, respectively (Figure 5b,d). In shallow snow areas such as lower Liaohe River Plain, the presence of snow has a certain cooling effect. Higher temperatures in these regions lead to faster snow melting, and the existence of shallow snow reduces the absorption of solar shortwave radiation by the surface. At the same time, the melting of snow cover requires a certain amount of heat to be absorbed, leaving the T SS below the T SA (Figure 5f).
In the Northern Hemisphere, snow cover provides a refuge for many terrestrial organisms during winter [15]. Increased T SS creates a more suitable environment for subnivean species and plants to overwinter. Increasing T SS can also promote plant activity, especially during late winter or early spring. Starr and Oberbauer [45] found that a higher T SS makes the air temperature much higher than the critical point of photosynthesis, thereby promoting plant photosynthesis and the metabolic activity of soil microbes in early spring. However, the increase of T SS may have a negative impact on some microbes. For example, a higher T SS will increase the metabolic rate of overwintering caterpillars living above the snow cover by altering the frequency of freezing, thus reducing the creatures' adaptability to the environment [46]. Moreover, the change of T SS has a great impact on low temperature-adapted microbial flora beneath the snow cover. This microbial flora plays a role in a certain soil temperature range and is highly sensitive to warming or cooling outside this range. Increased T SS leads to changes in the soil respiration rate of the microbial flora, thus affecting the rate of soil carbon sequestration [47].
Climate change can lead to complex changes in snow cover, which can alter the exchange of net radiation energy between snow and soil and the loss of energy on the soil surface [26]. For example, the increase of M SD in the farmland of NEC leads to a rapid increase in T SS , which may change the freeze-thaw cycle and conductive heat flow to the atmosphere [48][49][50]. Changes in these factors also affect patterns of CO 2 fluxes, as they are controlled by many biological, chemical, and physical factors such as snow accumulation, freeze-thaw cycles, soil temperature, and soil moisture [51].

Model Limitations and Uncertainties
There are some important uncertainties and limitations in our modeling analysis due to the model physics, inputs, and parameterization. It is necessary to set the initial number of soil layers as well as the soil moisture content and soil temperature of each layer when deriving the model. Because we cannot obtain the in situ soil parameters on each node layer of farmland in NEC in the gridding process, our research only used existing products, such as GLDAS data. Although we regard the first ten days of the model as the stable period and discard these data, the accuracy of GLDAS data will inevitably affect our results. Additionally, the initialization of the SNTHERM model requires the input of snow parameters, such as the number of snow layers, the snow grain diameter, the snow density, and the snow temperature of each layer. These data are currently without gridded products and in situ data. Therefore, the empirical value can only be set according to the characteristics of snow cover in NEC and model principle. Furthermore, T SS is highly sensitive to J p in the SNTHERM model. Because the soil properties of each pixel cannot be obtained, we adopted empirical values based on the soil texture data. These uncertain factors may have a greater impact on the simulation precision of the model. The underlying surface type also affects the model accuracy. The inversion accuracy of SNTHERM decreases in mixed pixels, forests, and grassland.
Additionally, although we analyzed the impact of several environmental factors on changes in T SS , other environmental factors may also have a strong influence. For example, the decrease in snow depth or increase of snow density reduces T SS [15]. However, due to the lack of snow density data, we did not discuss the changes in snow density and their impact on T SS . Future developments may enable the use of data assimilation methods to merge in situ data with satellite-based snow products and snow inversion models. Finally, in our research, snow depth products are retrieved using passive microwave data. However, due to the properties of passive microwaves, a low inversion accuracy is obtained when the snow depth is shallow [19].

Conclusions
In this study, the single-point SNTHERM model was improved according to snow cover and soil characteristics of farmland in NEC. The improved SNTHERM model was used to simulate T SS in each grid point from 01 September to 31 May between 1979 and 2018, while GLDAS forcing data were used to provide the initial attribute information of each node layer and ITPCAS forcing data as the meteorological driving data. In this paper, the interface temperature under snow cover, also known as T SS , was extracted, the changes in T SS and T DSSA were quantitatively analyzed, and the specific factors that caused the increase of T SS in farmland in NEC were deeply explored. The improved SNTHERM model employed observational data from 36 meteorological stations located in farmland in NEC from 2005 to 2018 to compare with the simulation results. The results showed that the SNTHERM model can adequately simulate T SS . We identified that 99% of annual average T SS and 98% of T DSSA showed an increasing trend. The results indicated that the snow depth is critical to the thermal insulation effect of snow cover. T DSSA is higher in areas with greater snow depth, such as the Sanjiang Plain and Songnen Plain. There are small gaps between T SS and T SA in areas with shallow snow, such as the lower Liaohe Plain. The calculation of the contribution of M SD and T SA to T SS further showed the influence of these two factors on T SS in different regions. The results revealed that M SD is the pivotal factor affecting the trend of T SS . T SS are more strongly correlated with M SD than T SA . For example, in the Sanjiang Plain, the rapid increase of T SS between 1979 and 2018 is mainly affected by M SD . Some areas of NEC, such as in the Songnen Plain, show significant rapid increases of T SS under the dual impact of increasing snow depth and increasing air temperature. In the areas, such as in the lower Liaohe River Plain, T SS is mainly affected by M SD and T SA from 1979 to 2018, and smaller changes in M SD causes T SS to increase slowly.
Our research provides a simulation of the trend of T SS in farmland in NEC under the background of climate change. The results show that the variations of snow depth and temperature have a strong correlation with the rise of T SS .