On the freeze–thaw cycles of shallow soil and connections with environmental factors over the Tibetan Plateau

Changes in the freeze–thaw cycles of shallow soil have important consequences for surface and subsurface hydrology, land–atmosphere energy and moisture interaction, carbon exchange, and ecosystem diversity and productivity. This work examines the shallow soil freeze–thaw cycle on the Tibetan Plateau (TP) using in–situ soil temperature observations in 0–20 cm soil layer during July 1982–June 2017. The domain and layer averaged beginning frozen day is November 18 and delays by 2.2 days per decade; the ending frozen day is March 9 and advances by 3.2 days per decade; the number of frozen days is 109 and shortens by 5.2 days per decade. Altitude and latitude combined could explain the spatial patterns of annual mean freeze–thaw status well. Stations located near 0 °C contour line experienced dramatic changes in freeze–thaw cycles as seen from subtropical mountain coniferous forest in the southern TP. Soil completely freezes from surface to 20-cm depth in 15 days while completely thaws in 10 days on average. Near-surface soil displays more pronounced changes than deeper soil. Surface air temperature strongly influences the shallow soil freeze–thaw status but snow exerts limited effects. Different thresholds in freeze–thaw status definition lead to differences in the shallow soil freeze–thaw status and multiple-consecutive-day approach appears to be more robust and reliable. Gridded soil temperature products could resolve the spatial pattern of the observed shallow soil freeze–thaw status to some extent but further improvement is needed.


Introduction
As a sensitive indicator of climate change, shallow soil (usually within 1 m depth below surface) freeze-thaw cycles (SSFTC for short) affect the water and energy flux interactions between land surface and atmosphere (Zhang and Armstrong 2001). Changes in SSFTC have important consequences for land surface and subsurface hydrological processes (e.g., Wang et al. 2009;Cheng and Jin 2013;Cuo et al. 2015;Zheng et al. 2018;Smith et al. 2019), soil microbial processes (Urakawa 2014), carbon exchange (Turetsky 2019), vegetation growth and ecological processes (Van Wijk et al. 2003;Beer et al. 2007;Kreyling et al. 2008;Cuo et al. 2016). Thus, understanding the changes in SSFTC and the environmental factors behind the changes is essential for understanding climate change impacts on surface water and energy balance.
SSFTC is closely related to hydrological processes. The seasonal freeze-thaw layer acts as an aquiclude or aquitard layer, which affects the water storage and infiltration capacity of the soil, resulting in the redistribution of the timing and magnitude of precipitation and melt-water traveling into rivers (Wang et al. 2009;Zheng et al. 2018). SSFTC also adds complexity to the hydrological cycles due to the intrinsic connections of both processes to groundwater (Bookhagen 2012).
Climate warming induced changes in SSFTC have a significant influence on soil carbon and nitrogen balance 1 3 through two mechanisms. Firstly, the shortened freeze or the prolonged thaw due to soil temperature rising above freezing under climate warming accelerates soil microorganism's activities to break down organic matter in the soil (Liang 2015;Bracho 2016;Cheng 2017). Secondly, the warming permafrost and seasonally frozen soil decompose more quickly and release more soil organic carbon and nitrogen into the atmosphere (Zimov et al. 2006;Schuur et al. 2009;Natali et al. 2019), potentially further accelerating global warming.
Vegetation growth is limited not only by heat and nutrient availability, but also by the availability of liquid soil water (Beer et al. 2007), which is closely related to SSFTC. Given this importance, the correct prediction of SSFTC's timing is essential for most ecosystem models (Van Wijk et al. 2003;Wania et al. 2009). SSFTC is also in phase with the plant life cycle, with the Gross Primary Productivity almost initiated in the beginning of soil thaw and terminated when soil freezes (Kreyling et al. 2008). Furthermore, climate warming increases the frequency of SSFTC which in turn enhances nitrogen uptake by plants (Urakawa et al. 2014;Sanders-DeMott et al. 2018).
Over the past decades, several approaches including in-situ observations, numerical modeling and satellite monitoring have been utilized to investigate SSFTC. In-situ observations can reliably describe SSFTC at different depths and therefore have been widely used (e.g., Henry 2008; Sinha and Cherkauer 2008;Anandhi 2013). However, in-situ observations are limited in space and time and are especially limited in high-elevation and high-latitude areas (Zhu et al. 2018). Numerical modeling can be and has been applied at regional and global scales (e.g., Slater et al. 1998;Schaefer et al. 2009;Rawlins et al. 2013;Guo et al. 2018). However, the inherent uncertainties in model structure, forcings and parameterizations bring uncertainties to model simulations and thus limit the application of modeling results (Williams et al. 2015).
Microwave remote sensing technology is a new approach for detecting soil freeze-thaw status. Microwave at low frequency (< 10 GHz) is very sensitive to water phase change in soil because liquid water has a much higher permittivity than ice (Entekhabi 2010). Many studies investigated surface soil freeze-thaw status using microwave remote sensing over the past decades (e.g., Zhang and Armstrong, 2001;Smith et al. 2004;Kim et al. 2011). Although microwave remote sensing can collect high-resolution and continuous observations, it is difficult for microwave remote sensing to detect soil freeze-thaw status below 10 cm (Chen et al. 2019). Also, to improve the accuracy of freeze-thaw detection by microwave remote sensing, the inversion algorithms need to be optimized using in situ measurements and adjusted for different land cover and vegetation types (Park et al. 2011;Derksen 2017). In addition, there do not exist long-term microwave remote sensing observations of soil freeze-thaw status, rendering it impossible to investigate the long-term changes of soil freeze-thaw cycles and their responses to climate change. Clearly, each aforementioned approach has its strengths and weaknesses in the studies of SSFTC.
About 57% of the exposed land in the Northern Hemisphere is subject to seasonal freeze-thaw processes (Zhang et al. 2003), of which the Tibetan Plateau (TP,(25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40) has the largest coverage of permafrost and seasonally frozen ground in the low-to-middle latitudes (Yang et al. 2010;Zou 2017). The TP features an average altitude of 4500 m ASL (above sea level) and an area of approximately 2.57 × 10 6 km 2 with annual mean air temperature of less than 5 °C (Fig. 1a). The permafrost and seasonally frozen ground over the TP are currently degrading due to warming (Cheng and Wu 2007;Cheng and Jin 2013). Given the TP's relatively remote location from China's industrialized areas and relatively pristine environment, changes in SSFTC on the TP could well serve as an indicator of high-altitude regional climate change impacts.
Many earlier studies investigated SSFTC on the TP at locations mainly along the Qinghai-Tibetan Railway (Cheng and Wu 2007;Wu et al. 2010Wu et al. , 2015, in the center Yuan et al. 2020) and the Three-River Headwater Region (Wang et al. 2019a), which cannot provide a plateau-wide picture of SSFTC. In recent years, a few studies examined TP's SSFTC at the regional scale based on in-situ observational data obtained from the China Meteorological Administration (CMA) (e.g., Peng 2017; Wang 2019b; Luo et al. 2020). However, these studies provided only limited information about the freeze-thaw status at various soil depths.
Remote sensing measurements have also been used to detect near-surface freeze-thaw cycles on the TP (e.g., Li et al. 2012;Kou et al. 2017;Zhao 2017). For example, Li et al. (2012) examined the Special Sensor Microwave/Imager (SSM/I) data and reported that the TP experienced a trend toward an earlier thawing date, a later freezing date and decreased frozen days during 1988-2007. Similar results were also obtained from numerical simulations (Guo and Wang 2014), although there were apparent biases especially in deep soil layers (Guo and Wang 2014;Yang 2020) due to uncertainties in model forcings and shortcomings of freeze-thaw parameterization schemes (Niu and Yang 2006;Hu 2019).
In this study, we use the long-term and up-to-date in-situ soil temperature observations at various depths in 0-20 cm and over different vegetation types to investigate the spatiotemporal patterns and variations of SSFTC on the TP. For the first time we compare the differences in freeze-thaw status calculated using several widely used freeze-thaw threshold schemes. This comparison could not only test the robustness of our method against the other widely used methods but also provide reference for the selection of freeze-thaw threshold schemes for future studies. We also investigate the possible mechanisms behind SSFTC changes in relation to environmental factors such as latitude, altitude, vegetation, air temperature and snow cover. Furthermore, we evaluate several popular global gridded soil temperature datasets in resolving the observed SSFTC on the TP. This study is the first to study SSFTC in full 0-20 cm depths using the most up-to-date observations, to investigate the relationships between SSFTC and multiple environmental factors,  (Zou et al. 2017) and distribution of stations relative to frozen soil (c). Vegetation types on the TP and distribution of stations relative to vegetation types (d) to examine the uncertainty of various thresholds related to SSFTC. SSFTC at various depths derived from first-hand observations provide essential information of regional warming impacts and can be used for optimizing remote sensing algorithms as well as evaluating and improving numerical model simulations. The findings from this study could also help with designing mitigation strategies for alleviating climate change impacts on ecosystems over the TP.

Observations
The observations used in this study include soil temperature of shallow five soil layers (0, 5, 10, 15 and 20 cm), 2-m daily mean surface air temperature and snow depth at 88 meteorological stations on the TP during the 36-year period of 1982-2017 obtained from the Meteorological Bureaus in Qinghai Province and Tibet Autonomous Region. The soil temperature data is used to describe the spatiotemporal variations of the soil freeze-thaw cycles at various depths. Air temperature and snow depth data are used to investigated their effects on SSFTC.

Gridded soil temperature products
This study is the first so far in that it investigates the spatiotemporal variations of the soil freeze-thaw cycles at various depths using the longest and most up-to-date in-situ soil temperature records on the TP. Admittedly, sparse stations restrict the representativeness of our results and global soil temperature products should have the potential to fill in the gap given that these products show reasonable performance in resolving the observed SSFTC. Here we also collect four popular gridded soil temperature products: two reanalysis products (CFS/CFSv2, hereafter CFSR; ERA-interim, hereafter ERA) and two modelling products (GLDAS-CLM, hereafter CLM; GLDAS-Noah, hereafter Noah). CFS is a global high-resolution reanalysis product covering 1979-2010 developed by the National Centers for Environmental Prediction (NCEP) Climate Forecast System (Saha 2010; http:// cfs. ncep. noaa. gov). CFSv2 is an extension of CFS covering 2011 to the present (Saha 2014). ERA-interim global analysis product is provided by ECMWF (European Center for Medium-Range Weather Forecasts) covering January 1979 to August 2019 (Dee 2011; http:// www. ecmwf. int/ en/ resea rch/ clima te-reana lysis/). Noah and CLM soil temperature products are from the Global Land Data Assimilation System (GLDAS) and driven by the Noah model and CLM model, respectively (Rodell et al. 2004; http:// disc. sci. gsfc. nasa. gov/ hydro logy/ data-holdi ngs). In this study, Noah soil temperature for 1982-1999 and 2000-2017 come from GLDAS version 2.0 and 2.1, respectively. All gridded products output soil temperature at much deeper depths but their vertical resolution in the soil is much coarser than the in-situ observations. To compare with the observations, we selected soil temperatures in 0-10 cm for CFSR and Noah, 0-7 cm for ERA, and 0-1.8, 1.8-4.5 and 4.5-9.1 cm for CLM. The detailed information of these products is shown in Table S1.

Functional vegetation types and elevation
The functional vegetation types used in this study come from the 1:1,000,000 vegetation map generated by the Institute of Geographic Science and Natural Resources Research of the Chinese Academy of Sciences (CAS) and obtained from CAS's Data Center for Resources and Environmental Sciences (http:// www. resdc. cn/). The spatial pattern and descriptions of the 11 functional vegetation types on the TP are shown in Fig. 1d and Table 1, respectively. Elevation is from the NASA Shuttle Radar Topography Mission (STRM, http:// earth explo rer. usgs. gov/) digital elevation models (Farr 2007).

Data processing procedures
Firstly, a thorough data quality control on three observational variables is conducted following these three steps: (1) removing outliers with ≥ 3 standard deviations (3σ) around the long-term means; (2) filling the missing values (including outliers) if the days with missing values are less than or equal to 5 consecutive days through averaging the values two days before and after the missing values. If the missing period is greater than 5 consecutive days, the missing values are left unchanged; (3) filling in the missing values at a station with the adjacent stations that have good linear correlations if the station has more than 300 valid records in a year. Otherwise, the records for this year are set as missing; and (4) after above 3 steps, 69 out of 88 stations with 27 years of freeze-thaw indices or more (i.e., ≥ 75% out of 36 years) (Jones and Hulme 1996) are selected for subsequent analysis.
Secondly, we design three indices in a frozen year, i.e., the period from July 1 of current calendar year to June 30 of next calendar year for representing the entire freeze-thaw cycle on the TP: beginning frozen day (BFD), ending frozen day (EFD), and number of frozen days (NFD). BFD (EFD) is defined as the first (last) day of first (last) 5 consecutive days when soil temperature is less than or equal to 0 °C. NFD is the number of days with soil temperature less than or equal to 0 °C between BFD and EFD with the days with soil temperature greater than 0 °C excluded (Fig. S1). NFD represents actual frozen soil days. These three indices are calculated separately at each layer. BFD of 0-20 cm is the average of the BFDs of 0, 5, 10 15 and 20 cm layers. EFD and NFD of 0-20 cm are calculated similarly. The adoption of five consecutive days is to reduce the impact of day-to-day variations in soil temperature that frequently occur in the shallow layer on the TP (Cheng and Wu 2007). For the ease of computation, all of the dates are expressed as the number of days in a frozen year, e.g., July 1 of 1982 is Day 1 of frozen year 1982, July 2 of 1982 is Day 2 of frozen year 1982, June 30 of 1983 is Day 365 of frozen year 1982, etc. Due to high soil temperature at some stations in some years, BFD and EFD could be missing in those years at those stations. In total, 67 stations for at least 27 freeze-thaw cycles in five soil layers are selected for further analyses. The elevation of these 67 stations ranges from 1814 to 4800 m ASL (Fig. 1b). Also, 64 out of the 67 stations are located on seasonally frozen ground while the remaining 3 are situated on permafrost ground (Fig. 1c). Table S2 contains the detailed information of these 67 stations.
To compare with the observations, we obtained the equivalent 0-10 cm soil temperature from the observations and the gridded products by averaging the available data within the 0-10 cm layer. More specifically, we selected the freeze-thaw indices calculated in 0-10 cm for CFSR and Noah, 0-7 cm for ERA, and the average of 0-1.8, 1.8-4.5 and 4.5-9.1 cm for CLM for comparison. All gridded soil temperature products were averaged to daily time scale to match the in-situ observations. We then computed BFD, EFD and NFD as defined before.
Given that conclusions based on a single approach can sometimes have large uncertainties when examining the characteristics of SSFTC (Brinkmann 1979), here, we use four different freeze-thaw thresholds applied in previous studies (S1-S4 in Table 2) and compare the freeze-thaw cycles derived from those four approaches with current results using the same observations.
To investigate the effects of seasonal snow on SSFTC, here we use the in-situ snow depth observations to examine the means and changes of four snow indices and their relations to SSFTC on the TP. The four snow indices are: snow depth (SD), beginning snow day (BSD) defined as the first snow day starting from 1 July, end snow day (ESD) defined as the last snow day ending on 30 June. and number of snow days (NSD) defined as the actual snow days between BSD and ESD.

Statistical analysis
The temporal trends of the freeze-thaw indices, air temperature and snow indices are examined using the Mann-Kendall test (Mann 1945;Kendall 1975) and Sen's slope approach (Sen 1968). The relationship of SSFTC with surface air temperature and snow indices is investigated through correlation and partial correlation analysis and multiple linear regression analysis. We also use root mean square errors (RMSE) to evaluate the performance of four gridded products in describing SSFTC on the TP. The detailed information of these statistical methods is provided in supplementary material. And all the statistical analysis mentioned above are conducted under the R3.6.0 programming environment. Altitude is the mean altitude averaged over the same vegetation types (before '/') and averaged over the stations located in the same vegetation types (after '/') NDVI is the mean NDVI averaged over the same vegetation types (before '/') and averaged over the stations located in the same vegetation types (after '/'

Climatology and changes of shallow soil freezethaw cycles
The annual mean BFD averaged over the 0-20 cm depth ranges from the earliest Day 113 (October 21) at Wudaoliang to the latest Day 175 (December 22) at Luolong with a regional average of Day 141 (November 18) (Fig. 2a, b). Early BFDs occur primarily in the sparsely vegetated central TP and the Qilian Mountains with high elevation, while late BFDs are found over the dry and low Qaidam Basin, and the relatively densely vegetated northeast and southern TP (Fig. 2b). This can also be seen from the distributions of annual mean BFD against vegetation types presented in Fig. 2c which shows that low-biomass types such as AM (alpine meadow), AG (alpine grassland), AS (alpine shrubland), and TD (temperate desert) correspond to early BFDs while the high-biomass SMCF (subtropical mountain coniferous forest) corresponds to late BFDs, with TFG (temperate forest grassland) lying in between. AWG (alpine warm grassland) is associated with late BFDs due to its low Over the 0-20 cm soil depth, the earliest (latest) BFDs always occur at 0 cm (20 cm) for all vegetation types (Fig. 2c), indicating that freezing happens first at the surface and then progresses downward. The difference in BFDs between the surface and 20 cm depth for most vegetation types is about 15 days. Furthermore, the vegetation types with earlier BFDs display smaller inter-layer differences in BFDs than those with later BFDs, e.g., a standard deviation of only 3 days in 0-20 cm for AM vs. 7 days for SMCF, implying that the entire soil layer experiences a faster freezing process in colder and less-vegetated areas compared to warmer and more-vegetated areas. Notably, for the forest types of TFG and SMCF that are situated in the relatively low and humid eastern and southern TP where the environment is favorable for vegetation growth, there is a small difference (less than 2 days) in BFD between the 0-and 5-cm depths (Fig. 2c). An explanation for this is that compared to the other vegetation types on the TP, TFG and SMCF are filled with thicker humus and denser root in the top 5 cm that permit relatively easy flow of air and water, thus leading to more or less synchronized heat and moisture cycles between the surface and the 5-cm depth.
A significantly positive (i.e., delay) BFD trend of regionally averaged 2.2 day dec -1 over the 0-20 cm depth is seen during 1982-2017 (Fig. 2a). Out of the 67 stations, 66 show delayed BFD trends (37 are statistically significant at p < 0.05), although there is also clear spatial heterogeneity in the trends (Fig. 2d). Only Xining exhibits an insignificantly negative (i.e., advance) trend of -1.2 day dec -1 , due to the decreasing trend in soil temperature at this station (not show). Averaged over the same vegetation types, BFDs show positive trends for the entire soil layer and for all vegetation types (Fig. 2e), but there are differences in the magnitudes and vertical distributions in the BFD trends amongst the vegetation types. SMCF clearly displays the largest delay (4.0 ± 1.4 day dec -1 ), which is followed by the alpine vegetation types including AM, AG and AS (2.5 ± 1.5 day dec -1 ), with TD, TFG and AWG corresponding to the smallest delay (1.6 ± 1.5 day dec -1 ). Overall, the delayed trends of BFD at surface are greater than those in the lower layers, and for most vegetation types there is a tendency for the magnitudes of the trends to decrease first with depth and then slightly increase with depth.
The annual mean EFD averaged over the 0-20 cm depth ranges from the earliest Day 212 (January 28) at Luolong to the latest Day 285 (April 11) at Wudaoliang with a regional average of Day 252 (Fig. 3a, b). Luolong and Wudaoliang also exhibit the latest and earliest BFD, respectively. EFD shows an opposite spatial pattern compared to BFD (Fig. 3b vs. Fig. 2b), i.e., the southern and northeastern TP and the Qaidam Basin experience early EFD (late BFD) while the central TP experiences late EFD (early BFD) (Fig. 3b). This also holds for the vegetation types (Fig. 3c vs. Fig. 2c), where the earliest EFDs and latest BFDs correspond to AWG and SMCF while the opposite is true for AM. For all vegetation types except for AG, the earliest and latest EFDs occur at surface and at 20 cm depth, respectively, with a mean difference of about 10 days (Fig. 3c), indicating that thawing, similar to freezing, starts at surface and then progresses downward. However, on average, from surface to 20 cm soil completely freezes in 15 days and completely thaws in 10 days.
During 1982-2017, EFD shows a significantly negative (i.e., advance) trend of -3.2 day dec -1 when averaged over the region and over the 0-20 cm depth (Fig. 3a), which is 1 day greater in magnitude than the delayed trend of BFD (Fig. 2a). Spatially, only 4 stations display slightly positive EFD trends while all other stations correspond to negative EFD trends (Fig. 3d). Similar to BFD but opposite in sign, the greatest EFD advance trends are seen for SMCF throughout the 0-20 cm layer, and the EFD advance trends are generally larger at surface than deep down for most vegetation types (Fig. 3e). Also, the EFD trends display the tendency of decreasing in magnitude with depth between the 0-10 cm depth but stays roughly constant further down (Fig. 3e).
The annual mean NFD averaged over the 0-20 cm depth ranges from 33 days at Luolong to 170 days at Wudaoliang with a regional average of 109 days (Fig. 4a, b). As mentioned before, Luolong and Wudaoliang display the latest (earlies) and earliest (latest) BFD (EFD), respectively. Among all vegetation types, AM that is located in the central TP with high elevation corresponds to the longest NFD (142 ± 12 days), while SMCF (69 ± 26 days) and AWG (69 ± 18 days) that are situated in the southern TP with relatively warm temperature are associated with the shortest NFD (Fig. 4c). Compared to BFD and EFD, the depth-to-depth differences in NFD for each vegetation type are small (Fig. 4c) and this is because both BFD and EFD are synchronized throughout the layer in that the earliest and the latest EFD and BFD always occur at surface and at 20-cm depth, respectively.
Due to the opposite trends in BFD and EFD, the TP as a whole experiences a significantly negative (i.e., decrease) NFD trend of -5.2 day dec -1 averaged over 0-20 cm during 1982-2017 (Fig. 4a). Most stations show greater trends in NFD than in BFD and EFD. The greatest decrease trend in NFD is -12.6 day dec -1 at Tongde, followed by -10 day dec -1 at both Leiwuqi and Changdu (Fig. 4d). Among the seven vegetation types, SMCF experiences the greatest NFD decrease with layer averaged -7.7 ± 5.5 day dec -1 (Fig. 4e), due to the largest BFD delays but the largest EFD advances for SMCF. AM has the second highest layer-averaged NFD decrease rate of -5.8 ± 1.7 day dec -1 (Fig. 4e). Within the layer, the NFD trends tend to decrease significantly in magnitude from surface to 5 cm but change relatively little deep down for all vegetation types except for SMCF and AWG.
To understand why SMCF experiences the most pronounced changes in the freeze-thaw status, we generate annual mean 0-cm soil temperature grid with XX spatial resolution for freezing season (October to December) and thawing season (January to March) using inverse distance weight approach and located 0 °C contour line in the two  (Fig. 5). SMCF is located in the southeast and its southern part corresponds to soil temperature well above 0 °C in the two seasons but its northern part where the majority of the stations for this vegetation zone are situated is in the vicinity of the 0 °C contour line, especially for freezing season. Thus, soil temperature for SMCF is sensitive to warming which leads to great changes in the freeze-thaw status.

Dependence of the freeze-thaw cycles on altitude and latitude
Despite large spatial variations of the three indices over the TP, clear altitude-latitude dependence in these indices can be detected as presented in the following. Annual mean and layer averaged BFD, EFD and NFD all show statistically significant correlation coefficients (CC) with altitude: negative for BFD (CC = -0.29, p < 0.05), positive for both  (Fig. 6b). This indicates that overall BFD advances, EFD delays and NFD increases with both altitude and latitude on the TP. An exception is noticed within the 3500-4500 m range and 34-38°N band for which BFD, EFD, and NFD changes with altitude and latitude differ from the surrounding areas (Fig. 6c, d). Further investigation reveals that most stations in 3500-4500 m zone is located in relatively low latitudes, whereas 34-38°N corresponds to relatively low elevation. The combination of high altitude and low latitude or high latitude and low elevation confounds the respective effects of latitude and altitude that would have been expected individually. Although altitude or latitude alone could account for the spatial characteristics of the freeze-thaw status to some extent, the combined effects of altitude and latitude could explain the spatial heterogeneity of BFD by 89%, EFD by 81%, and NFD by 80%, respectively, based on our multiple Fig. 7 Spatial distributions of the soil minus air in annual mean BFD (a), EFD (c), and NFD (e). The positive values in (a) indicate that soil experiences later BFDs than air. The negative values in (b) indicate that soil experiences earlier EFDs than air. The negative values in (c) indicate that soil has shorter NFDs than air. The histograms of the annual mean BFD/EFD/ NFD are also shown at the bottom of (a), (c) and (e). And the background colors represent the vegetation types. Scatter plot of the soil minus air differences in annual mean BFD (b), EFD (d) and NFD (e) vs. altitude (yellow) and latitude (black) regression analysis (not show). This appears to be consistent with Lin and Wu (1981) who showed that altitude and latitude together dictated the spatial distributions of heat and moist on the TP. The altitude-latitude dependence could also explain the differences in the freeze-thaw status among the vegetation types. For example, TD and TFG lie within the similar latitudes of 36-37°N but have an approximate 600-m difference in altitude, and this difference in altitude may be attributed to an earlier BFD by 13 days, a later EFD by 11 days and a longer NFD by more than one month for TD compared to TFG (Table 1, Figs. 2c, 3c, 4c).
Next we apply the same calculation of soil BFD, EFD and NFD to surface air temperature to obtain the equivalent air BFD, EFD and NFD although we recognize that air never freezes or thaws. Through investigating the differences in the air-soil freeze-thaw cycles, we aim to shed further light on the altitude-latitude dependence. Clearly, air and soil are a coupled system within which heat is transferred and modulated by local air and soil states and physical properties such as specific heat capacity and heat conductivity. Heat moves along the heat gradient from high to low in this air-soil system. The differences in the three indices between air and soil are showed in Fig. 7. Here we focus on surface air temperature and soil surface (0 cm layer).
Owing to high terrain accompanied by strong solar radiation, soil surface temperature on the TP is higher than 2 m air temperature all year round (Zhang et al. 2006;Zhu et al. 2018). Indeed, soil BFD at 0 cm lags air BFD by 10 days on average (Fig. 7a). The lag time in BFD by soil increases with altitude (R 2 = 0.53, p < 0.001) but decreases with latitude (R 2 = 0.20, p < 0.01) (Fig. 7b), indicating that soil freezing occurs much later than air "freezing" in higher altitude (lower latitude) than lower altitude (higher latitude). The largest lag in mean BFD is seen for AM (15 ± 6 days) that has altitude exceeding 4700 m ASL, followed by SMCF (14 ± 5 days) that lies at a relatively low latitude of around 30°N. For EFD, surface soil leads surface air by 21 days on average (Fig. 7c). The lead time in EFD by soil increases with altitude (R 2 = 0.67, p < 0.001) but decreases with latitude (R 2 = 0.19, p < 0.001) (Fig. 7d). The greatest lead in mean EFD occurs for AM (31 ± 8 days), followed by AG (28 ± 9 days), both of which have average altitude above 4700 m ASL. Due to the lag in BFD and the lead in EFD by surface soil relative to surface air, soil NFD is shorter than air NFD by an averaged 27 days (Fig. 7e). The soil-air differences in NFD increases with altitude (R 2 = 0.62, p < 0.001) but decreases with latitude (R 2 = 0.17, p < 0.001) (Fig. 7f). Similar to BFD and EFD, high altitude AM displays the greatest difference in NFD (45 ± 15 days). We also find that altitude (> 53%) can explain more spatial variability of the soil-air differences in the freeze-thaw status than latitude (≤ 20%). We also examine the trends of the SSFTC in relation to altitude and latitude (Table S3) but we identify only weak correlations, which suggests that altitude and latitude exert limited influence on the temporal variations of the SSFTC on the TP.  BFD/EFD is defined as the first/last day from 1 July with temperatures below -2.2 °C Peterson and Abatzoglou (2014), Schwartz et al. (2006) S3 BFD (EFD) is defined as the first (last) day of first (last) 3 or more consecutive days when soil temperature is less than or equal to 0 °C from 1 July Li et al. (2012) S4 BFD (EFD) is defined as the first (last) day of first (last) 5 or more consecutive days when soil temperature is less than or equal to 0 °C from 1 July Jin et al. (2009)

Relationship between surface air temperature and shallow soil freeze-thaw cycles
SSFTC is the result of water and heat flux exchanges between the atmosphere and land surface. Surface air temperature, as the main embodiment of heat, is the primary meteorological factor influencing SSFTC (Smith et al. 2004;Wang et al. 2015). On the TP, annual mean air temperature ranges from -4.8 °C at Wudaoliang to 9.2 °C at Xunhua with an average of 2.7 °C. Among all vegetation types, AM has the lowest seasonal and annual average air temperature, followed by AG and AS. The stations corresponding to AM, AG and AS have an average altitude of about 4000 m ASL and are located in the central TP (Table 1; Fig. 1d). Seasonal and annual mean air temperature displays statistically significant warming trends for all 7 vegetation types during 1982-2017 (Table S4). The greatest warming trends occur in winter for all vegetation types except for TD and TFG for which the greatest warming trends are in summer. Annually, AM exhibits the greatest warming trend (0.63 ± 0.05 °C dec -1 ), followed by AS (0.62 ± 0.24 °C dec -1 ). The spatial patterns of mean BFD, EFD and NFD correlate strongly with spatial patterns of mean surface air temperature within the entire 0-20 cm layer at both seasonal and annual time scales, although the correlation weakens slightly with depth due to weakened heat flux exchange between the atmosphere and soil. The CC reach more than 0.85 in winter and autumn when averaged over the 0-20 cm layer. It is worth noting that both BFD and EFD show the strongest correlation with surface air temperature in winter. It is also interesting to note that all three indices display statistically significant correlation with surface air temperature even in summer (Table 3).
Next we examine the spatial patterns of partial correlation coefficients (PCC) of BFD/EFD/NFD at 0 cm with leading, concurrent and lagged seasonal air temperature (Figs. 8,9,  A black " + " sign at the center of a circle indicates that the partial correlation coefficient is significant at p < 0.05. The same notation applies hereafter 10). BFD displays significantly positive correlation with concurrent autumn air temperature at more than 83% of the stations (Fig. 8d), illustrating the strong influence of concurrent air temperature on BFD. The exception is over SMCF and AWG where most stations exhibit insignificant correlation in autumn and this is because BFD at these stations generally occurrs in winter (December 10 on average) instead of autumn. BFD has rather weak correlations with air temperature in previous seasons (Fig. 8a-c), implying the weak lag effects of air temperature on soil freezing.
For EFD, more than 86% of the stations show negative correlation with air temperature in spring and the 55 significantly negative correlations are associated with vegetation types other than SMCF and AMG (Fig. 9d). The correlation is weak for the previous summer and autumn (Fig. 9a,  b). Stations in SMCF and AMG generally experience EFDs in early spring (March 11 on average) and strong negative correlations between EFD and air temperature are found for previous winter (Fig. 9c). Accordingly, for NFD, most stations located in vegetation types other than SMCF and AMG display significantly negative correlations with air temperature in autumn and spring (Fig. 10a, c). For SMCF and AMG, significantly negative correlation is seen mostly in winter (Fig. 10b). Clearly, concurrent seasonal air temperature exerts a profound impact on soil freeze-thaw status but the exact timing varies spatially.

Relationship between snow and shallow soil freeze-thaw status
Besides altitude, latitude and surface air temperature, seasonal snow also affects SSFTC through modulating the atmosphere-soil heat transfer due to its high albedo, absorptivity and emissivity but low thermal conductivity. The net effect of snow cover mainly depends on the timing, duration, depth, accumulation and melting processes (Groisman For all vegetation types, annual mean snow depth (SD) spans within 2-3 cm with an average of 2.3 cm, and decreases slightly at an average rate of -0.08 cm dec -1 during 1982-2017. Beginning snow day (BSD) ranges from Day 40 (August 9, Wudaoliang) to Day 233 (February 18, Jiangzi) with an average of Day 133 (November 10). Among all vegetation types, AM displays the earliest snow day (Day 77 or September 15), followed by AS (Day 100 or October 8). BSD trends are all positive (i.e., delayed) and the average is 6 day dec -1 , but only AG and AM correspond to statistically significant delay trends. End snow day (ESD) ranges from Day 207 (March 9, Guide) to Day 353 (June 18, Wudaoliang) with an average of Day 294 (April 20). ESD trends are all negative (i.e., advanced) with an average of -4 day dec -1 but none of the trends are statistically significant. Number of snow day (NSD) ranges from 2 at Xunhua to 101 at Qingshuihe and the average is 25 days. All NSD trends are negative (i.e., shortened) and spans from -5 to -1 day dec -1 (Table S5, Figs. S2-4).
The scatter plot of snow indices vs. soil freeze-thaw indices at surface reveals significantly positive correlations for BSD vs. BFD (R 2 = 0.42, p < 0.01), ESD vs. EFD (R 2 = 0.61, p < 0.01), and NSD vs. NFD (R 2 = 0.31, p < 0.01) (Fig. 11). This indicates that snow may exert strong influence on the timing of SSFTC on the TP. Assuming snow falls within 5 days prior to BFD/EFD has an impact on the soil freeze-thaw status, we divide all years during 1982-2017 into snow years and snow-free years using the 5-day window. For example, if a station whose BFD was on the 150th day in frozen year 2000 recorded any snow during 146th-150th days, then 2000 was regarded as a snow year for that station. Here we compare the differences in average BFD and EFD during snow and snow-free years at 0 cm to avoid the complexity of heat and water flux transfer in deeper layers that might offset or strengthen snow effects. The average BFD (EFD) of snow years is 1.9 days earlier (5.3 days later) than that of snow-free years (Fig. 12) which is perhaps due to the fact that snowfall is often accompanied by low air temperature.
Although the stations with earlier BSDs/ESDs also show earlier BFDs/EFDs and subsequently longer NFDs (Fig. 11), we speculate that the effects of snow by itself on SSFTC could be relatively limited on the TP due to the involvement of other factors such as air temperature. We obtain significant correlations between annual mean surface air temperature and BSD (CC = 0.79, p < 0.001), ESD (CC = -0.84, p < 0.001) and NSD (CC = -0.69, p < 0.001), and after excluding the influence of air temperature, PCC (-0.06) between BSD and BFD becomes statistically insignificant. Similarly, PCC for ESD vs. EFD and NSD vs. NFD is reduced to insignificant -0.002 and -0.135, respectively. It appears that the advanced BFD and delayed EFD revealed in snow year compared to snow-free year (Fig. 12) are likely resulted from colder air temperature accompanied with snow that consequently cools shallow soil, rather than from cooling effects of snow itself. The limited effects of snow by itself are perhaps due to the fact that snow on the TP is relatively thin and the duration of snow is short in both the freezing and thawing seasons. To demonstrate the limited snow effect on SSFTC on the TP, we investigate the timing, duration and depth characteristics of seasonal snow during soil freeze/thaw period on the TP. Before 5 days of BFDs, more than 50 stations have > 20 years with no snow records. And only 11 stations experience no more than 4 snow years with 5 consecutive snow records (Fig. 13a). Before 5 days of EFD, more than 30 stations have > 20 years with no snow records (Fig. 13b). We also investigate the characteristics of seasonal snow between BFD and EFD on the TP. For all station, the averaged ratio of snow days to NFD is 0.13 and more than 75% of these stations have the ratio < 0.2 (Fig. 14a). The averaged SD between BDF and EFD ranges from 1.2 to 8.2 cm with an average of 2.1 cm. Among all the stations, more than 80% of stations have the SD < 2.5 cm (Fig. 14b). We also detect the maximum consecutive snow days during BFD to EFD. Only 12 stations experience more than two weeks' consecutive snow cover during BFD to EFD (Fig. 14c). The scarcity of snow indicates the limited snow effects on SSFTC on the TP. Hence, the advanced BFD and delayed EFD revealed in snow year compared to snow-free year are likely resulted from colder air temperature accompanied with snow that consequently cools shallow soil, rather than from cooling effects of snow itself. This is obviously related to the fact that more than 70% of annual precipitation falls in May-September at 95% of climate stations on the TP due to summer monsoon (Cuo and Zhang 2017). Some previous investigations also stated that snow may not be a major factor in influencing SSFTC on the TP (Luo et al. 2017(Luo et al. , 2020Peng et al. 2017). Still more detailed studies about the effects of snow properties (e.g., timing, duration, depth, density and snow water equivalent) on the soil freeze-thaw status are needed.

Comparisons of the freeze-thaw status using various approaches
In order to assess the uncertainty associated with single approach, various approaches are applied to examine the characteristics of SSFTC on the plateau. Over the conterminous United States using -2.2 °C and 0 °C as thresholds for freezing and thawing generated similar results (McCabe et al. 2015). However, on the TP, using 0 °C as the threshold (S1) would produce an earlier BFD by 13 days, a later EFD by 17 days and an extended NFD by 33 days than using -2.2 °C as the threshold (S2) when averaged over the 0-20 cm depth (Table 4). The differences in the freeze-thaw indices among the various thresholds decrease with soil depth. The high sensitivity of the freeze-thaw indices to single-day temperature thresholds (S1, S2), especially at surface is perhaps due to fact that on the TP frozen ground temperature is relatively warm (Cheng and Wu 2007;Wu et al. 2010) and the freeze-thaw transition of shallow soil is frequent, especially for surface soil (Liu and Chen 2000). There is small but negligible difference in the shallow soil freeze-thaw status derived using three consecutive days (S3) or five consecutive days (S4). On the other hand, the trends of the freeze-thaw indices derived using the four thresholds are consistent in sign across the soil layer, i.e., delayed  (Table 5). Also, S3 and S4 show similar magnitudes in trend in general that are smaller compared to S2 and S1. This comparison appears to lend credence to the robustness of the S3 and S4 approaches in the soil freeze-thaw analyses which is similar to current approach.
In addition to temperature thresholds, observed shallow soil frozen depth has also been used to study the freeze-thaw status (Wang et al. 2019b;Luo et al. 2020) although with generally similar results to S1. However, compared to their studies, BFD from S1 is delayed by about 10 days and EFD from S1 is advanced by about 25 days. The trend of BFD is 1.8 day dec -1 (Luo et al. 2020) and 2.0 day dec -1 (Wang et al. 2019b), close to what S1 generates for BFD at 0 cm (2.1 day dec -1 ). Over the TP, Li et al. (2012) also examined the changes in near-surface soil freeze-thaw cycles using the SMM/I microwave remote sensing data during 1988-2007 and stated that BFD was delayed by 5.1 day dec -1 , EFD advanced by -7.2 day dec -1 and NFD shortened by -16.8 day dec -1 , all of which are greater than what we obtain for 1982-2017. It is worth noting though that their results are average grid values with statistical significance at p ≤ 0.01. In order to compare with the SMM/I results, we use the in-situ observations for the same period and we obtain the trend of 7.2, -8.8 and -10.7 day dec -1 for BFD, EFD and NFD, respectively. Guo and Wang (2014) conducted a modeling study and obtained the averaged trend of 1.7, -4.7 and -6.4 day dec -1 for BFD, EFD and NFD, respectively, during 1981-2010 on the TP. For the same period, this study finds the averaged trend of 2.5, -3.8 and -6.7 day dec -1 for BFD, EFD and NFD, respectively, similar to what the modeling results showed.
Also based on in situ observations, Wang et al. (2015) reported that the ground surface experienced advanced BFD by 7.2 day dec -1 , delayed EFD by -6.0 day dec -1 and shortened NFD by -8.9 day dec -1 during 1990-2006 over China. However, that study used limited stations on the TP and the analysis results may not be representative of the TP.

Evaluation of the freeze-thaw cycles derived from soil temperature products
As stations are limited on the plateau, it is imperative to evaluate the performance of gridded products in representing the freeze-thaw cycles so that the observations-void areas may also be examined in terms of the freeze-thaw cycles.
Here we evaluate several popular products as mentioned in Sect. 2.2. The spatial distributions of annual average BFD/ EFD/NFD from the gridded products show consistent spatial patterns, i.e., later BFD, earlier EFD and shortened NFD in the eastern and southern TP and the Qaidam Basin; and earlier BFD, later EFD and prolonged NFD in the central and western TP and the Qilian Mountains (Fig. 15). The SMM/I microwave remote sensing data (Li et al. 2012) and numerical simulations (Guo and Wang 2014) also showed similar spatial patterns. Further, values from the gridded cells covering ground stations are extracted to compare the freeze-thaw indices between the in situ observations and gridded products. CFSR, ERA and Noah display earlier BFD, later EFD and prolonged NFD, whereas CLM exhibits later BFD, earlier EFD and shortened NFD than observations (Fig. 16). Compared to other products, CLM shows better quality. Spatial scale and soil depth mismatches between the gridded products and in-situ observations account for the discrepancies to some extent. However, the underestimation of soil temperature for CFSR, ERA and Noah (Yang and Zhang 2018;Hu et al. 2019) and overestimation of soil temperature for CLM (Wang 2016) should also be the reasons for the discrepancies. The underestimation of soil temperature by CFSR, ERA and Noah might be the result of underestimated model forcing such as air temperature and downward longwave radiation at most sites on the TP (Wang et al. 2016;Yang et al. 2020). The sharp mutation around 2000 in BFD/EFD/NFD for CLM and Noah (Fig. 16) The characteristics of snow cover between the beginning and ending frozen days. The black bars represent the number of stations in each group. The orange dotted line represents the cumulative probability, a the ratio of snow days to NFD, b the averaged snow depth during beginning and ending frozen days, c the maximum continuous snow days is due to the switch of the forcing data from GLDAS-2.0 to GLDAS-2.1. The same land surface parameters derived from MODIS were used in GLDAS-2.0 and GLDAS-2.1. However, GLDAS-2.0 was forced by the global meteorological forcing dataset from Princeton University (Sheffield et al. 2006), while GLDAS-2.1 was forced by National Oceanic and Atmospheric Administration (NOAA)/Global Data Assimilation System (GDAS) atmospheric analysis fields (Derber et al. 1991), the disaggregated Global Precipitation Climatology Project (GPCP) V1.3 daily analysis precipitation fields (Huffman 2001;Adler 2003), and the Air Force   Weather Agency's agricultural meteorological modeling system (AGRMET) radiation fields. Besides forcing data, different model parameterization schemes also contribute to differences in simulations. The Noah model was also used in CFSR (Rodell et al. 2004;Saha et al. 2014), and TESSEL was employed in ERA. Among the numerous model parameters, soil thermal conductivity plays a great influence on soil thermal regime. For example, CLM employed the Johansen (1975) scheme, while Noah used the Farouki (1981) scheme. Some studies reported that these two schemes overestimated soil thermal conductivity over the TP (Chen et al. 2012;Du 2020) especially by the Farouki scheme (Wang et al. 2016;Du et al. 2020). , and NFD (c) derived from four gridded soil temperature products and in-situ soil temperature observations (red lines) at depth 0-10 cm on the TP. The correlation coefficients (CC) and root mean square errors (RMSE) between gridded soil temperature products and in-situ soil temperature in 2000-2016 are also shown in the figure We further evaluate the four gridded products for 2000-2017 after sharp mutation (Fig. 16). CLM and Noah show the highest correlations and smallest root mean square errors (RMSE) among the four products. Although CLM and Noah are relatively superior to the other two products, their RMSEs are still large. In the future, more effort is needed to develop reliable soil temperature products on the TP.

Conclusions
By using various approaches for defining the soil freeze-thaw status on the TP, we suggest that the approach of applying the threshold of multiple (e.g., 5) consecutive days be used in future as it is more robust and reliable than a single day value. During July 1982-June 2017, domain and 0-20 cm layer averaged BFD is 141 (November 18) and delays by 2.2 days per decade; EFD is 252 (March 9) and advances by 3.2 days per decade; NFD is 109 and shortens by 5.2 days per decade. Among 7 vegetation types, AWG located in the south, SMCF located in the southeast and TFG located in the northeast of TP exhibit later BFDs, earlier EFDs and shorter NFDs than the other vegetation types. In particular, SMCF experiences the most pronounced changes in the freeze-thaw status throughout the layer because most of its stations are located near the 0 °C contour line. Within the 0-20 cm layer, on average, soil completely freezes from surface to the bottom in 15 days while completely thaws in 10 days. Near-surface soil experiences more pronounced changes than deeper soil. There is strong dependence of the freeze-thaw status on altitude and latitude. The differences in surface soil-air freeze-thaw status increase with altitude and decrease with latitude.
Air temperature is the primary factor influencing SSFTC on the TP. The spatial pattern of annual average shallow soil freeze-thaw status is consistent with that of surface air temperature. Although there are significant linear correlations between snow depth and surface BFD/EFD/NFD, partial correlation reveals that the impacts of snow on soil freeze-thaw status are limited on the TP. The four gridded soil temperature products could resolve SSFTC spatial pattern on the TP to some extent but not completely. It appears that CLM is relatively superior over ERA, CFSR and Noah in reproducing temporal change of the shallow soil freeze-thaw cycles on the TP.
Notably, for the relationship between vegetation and soil freeze-thaw cycles, most studies focus on the impact of soil freeze-thaw cycle on vegetation. By categorizing SSFTC based on vegetation, we show that vegetation indeed affects SSFTC. However, we still know little about how vegetation impact SSFTC. This study is only a preliminary attempt in this respect, which needs in-depth study in the future.

Appendix
See Table 6.