The freeze–thaw changes of seasonally frozen ground (SFG) are an important indicator of climate change. Based on observed daily freeze depth of SFG from meteorological stations on the Tibetan Plateau (TP) from 1960 to 2014, the spatial–temporal characteristics and trends in SFG were analyzed, and the relationships between them and climatic and geographical factors were explored. Freeze–thaw changes of SFG on a regional scale were assessed by multiple regression functions. Results showed multiyear mean maximum freeze depth, freeze–thaw duration, freeze start date, and thaw end date that demonstrate obvious distribution characteristics of climatic zones. A decreasing trend in maximum freeze depth and freeze–thaw duration occurred on the TP from 1960 to 2014. The freeze start date has been later, and the thaw end date has been significantly earlier. The freeze–thaw changes of SFG significantly affected by soil hydrothermal conditions on the TP could be assessed by elevation and latitude or by air temperature and precipitation, due to their high correlations. The regional average of maximum freeze depth and freeze–thaw duration caused by climatic and geographical factors were larger than those averaged using meteorological station data because most stations are located at lower altitudes. Maximum freeze depth and freeze–thaw duration have decreased sharply since 2000 on the entire TP. Warming and wetting conditions of the soil resulted in a significant decrease in maximum freeze depth and freeze–thaw duration in the most area of the TP, while drying soil results in a slight increase of them in the southeast of the TP.
The Tibetan Plateau (TP) hosts the largest and thickest frozen ground at the middle and low latitudes, including permafrost and seasonally frozen ground (SFG) (Cheng et al. 2019; Yang et al. 2019). Frozen ground plays important roles in local and global circulation, climate, hydrology, and terrestrial ecosystems through its influence on surface energy, water, and carbon cycles in the Earth system (Derksen et al. 2012; Rong et al. 2005; Schmidt et al. 2011; Serreze et al. 2000). The annual freeze–thaw of frozen ground directly alters soil hydrology and thermal properties, enhances the energy exchange between land surface and atmosphere, and subsequently influences the upper–lower circulation that contributes to the South Asia high and East Asia climate (B. Chen et al. 2014; Cheng and Wu 2007; Li et al. 2002; Luo et al. 2009a,b; Wang et al. 2003; J. Wang et al. 2019). Soil carbon stored in frozen ground can be released into the atmosphere by these freeze–thaw cycles (Baumann et al. 2009; Ding et al. 2017; Tarnocai et al. 2009). Regional hydrological processes in the TP such as runoff generation, groundwater–surface water interaction, and soil moisture conditions, among others, are affected by frozen ground (Jin et al. 2009; T. Wang et al. 2019; Zhao et al. 2014). This influences water resources security downstream, where more than 1.4 billion people rely on the rivers originating from the TP (Immerzeel et al. 2010).
Recent research reveals that the area of SFG is larger than that of permafrost, accounting for 50%–56% of the total TP area (Chang et al. 2018; Shi et al. 2018; T. Wang et al. 2019; Zou et al. 2017). With further warming, nearly half the permafrost on the TP will degrade to SFG during the twenty-first century under the representative concentration pathway 4.5 scenario (Chang et al. 2018; T. Wang et al. 2019). To assess the freeze–thaw changes of SFG over the TP, it is crucial to delineate spatial–temporal characteristics and their changes in freeze depth and freeze–thaw duration. Recent assessments indicated that the maximum freeze depth and freeze–thaw duration of SFG in the TP have decreased steadily over the past decades (Guo and Wang 2013; Jin et al. 2009; Luo et al. 2017a, 2016; Zhao et al. 2004). This reflects the combined effects of changes in soil temperature and moisture conditions, soil hydrology and thermal state, and can lead to a redistribution of energy, water and carbon budget between land surface and the atmosphere (B. Chen et al. 2014; Luo et al. 2009a, 2017b; S. Luo et al. 2018; Luo et al. 2009b; Yi et al. 2013). Despite its vast extent and importance, the historical distribution and changes of SFG, as well as interactions between the soil and the atmosphere at regional scales have received little attention (B. Chen et al. 2014; Shiklomanov 2012). One of the major obstacles in assessing changes in SFG is the lack of long-term observations. For this reason, other climatic indicators, such as air temperature, soil temperature, precipitation, snow depth, and freeze–thaw index, as well as remote sensing data, numerical simulation, and machine learning algorithms are exploited to characterize freeze depth and freeze–thaw duration (Frauenfeld and Zhang 2011; Guo and Wang 2013, 2014; X. Li et al. 2012; Luo et al. 2017a; Peng et al. 2017; Qin et al. 2018; T. Wang et al. 2019). These results are important for understanding the freeze–thaw changes of SFG and the main factor of influence. It should be noted that although the proxy climatic data are close to the true value, it still has a large deviation because of other nonclimatic factors (Shiklomanov 2012). Though the remote sensing data can be analyzed to obtain spatial variations, the long-term historical variation fractures cannot be shown due to data shortage (Luo et al. 2017a). Although numerical simulations can provide higher spatial–temporal resolution and long duration, they are still prone to errors because of the large uncertainty in parameter and driving variables (Williams et al. 2015). Therefore, more in-depth studies are needed to project the freeze–thaw changes of SFG on the TP based on long-term observations.
Under the background of global warming, the TP shows an accelerated warming trend in recent years (Duan and Xiao 2015; Fang et al. 2019). Our previous study, found that a continuous, accelerated decreasing trend in freeze depth and freeze–thaw duration of SFG appeared in the Three Rivers Source Region (TRSR) of the TP (Luo et al. 2017a). But how do the SFG freeze–thaw changes over the entire TP? And are there any differences in their internal regions? If so, what are the possible reasons? This work addresses these issues by comparing maximum freeze depth and freeze–thaw duration from 1960 to 2014 over the TP, based on historical observations and empirical regression equations. Observed soil freeze depth data from 75 meteorological stations were analyzed for freeze–thaw changes in maximum freeze depth, freeze–thaw duration, freeze start date, and thaw end date of SFG. Freeze–thaw changes of SFG on a regional scale were assessed using the most important climatic and geographical factors. Potential causes and possible feedback of SFG change to climate and environment have also been discussed.
2. Data and methods
a. Freeze depth data, climate data, and map data
The basic data used in this study includes observed daily freeze depth of SFG. The freeze depth data were obtained from China Meteorological Administration (CMA) and the China Meteorological Data Service Center (CMDC) (http://data.cma.cn/) for a total of 87 meteorological stations on the TP (Fig. 1). The data were collected daily at these stations, and those from 1960 to 2014 were used in this study. The freeze depth of the SFG was observed once per day (0800 Beijing time, UTC + 8 h) using a frost tube when the ground surface temperature was below 0°C (CMA 2003; Luo et al. 2017a). A frost tube is a common instrument for monitoring soil freeze and thaw depths in cold regions (Phukan 1985), defined as the distance from the ground surface to the freezing front in autumn and winter, and from the ground surface to the thawing front in spring and summer, respectively. Study shows that frost tubes are excellent tools to monitor the accurate soil freeze depth conveniently and inexpensively in SFG areas (Iwata et al. 2012). In China, frost tube was buried in soils under natural land cover in the observation field (CMA 2003). A frost tube consists of an outer tube and an inner tube. The inner tube was a rubber tube containing clean water. Soil freeze depth was defined by the water freeze depth in the inner tube. The upper and lower boundary of freeze depth correspond to the upper and lower ice values (CMA 2003; Luo et al. 2017a).
Monthly air temperature and precipitation at the collocated meteorological stations which were from CMDC (http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_MON.html) were used to investigate the correlation between main climatic factors and SFG. The mean monthly gridded air temperature and precipitation data during 1960–2014 with 0.5° spatial resolution, which were also from the CMDC (http://data.cma.cn/data/cdcindex/cid/00f8a0e6c590ac15. html), were used to calculate SFG on a regional scale. The air temperature and precipitation lapse rates changed with increasing elevation. Annual air temperature and precipitation were downscaled to 1-km resolution for this study by adding an elevation adjustment term at each grid. The mean annual air temperature lapse rate was −4.16°C km−1 calculated by a linear regression method and the observations at the meteorological stations across the TP (X. Guo et al. 2016; T. Wang et al. 2019). The mean annual precipitation lapse rate was 168.2, 261.4, and 120.8 mm km−1, respectively, at 2–3, 3–4, and more than 4 km (X. Guo et al. 2016), respectively, and were used to downscale the precipitation analysis to 1-km resolution.
In this study, a new scheme for climate regionalization in China was designed to provide a climatic zones map on the TP (Fig. 1) using daily observations at the meteorological stations during the period 1971–2000. This scheme mainly relies on five principles, including zonal and nonzonal integration, genetic unity and regional relative consistent climate integration, comprehensiveness and leading factors integration, bottom-up and top-down integration, and spatial continuity and small patch omission (Zheng et al. 2010). The 1-km resolution elevation is resampled by the 90-m Shuttle Rader Topography Mission Digital Elevation Model (SRTMDEM) (http://www.gscloud.cn) using bilinear interpolation, which provides elevation to calculate the air temperature, precipitation, maximum freeze depth and freeze–thaw duration at 1-km resolution over the TP. In the meantime, a new map of frozen ground type (SFG, permafrost and unfrozen ground) on the TP (Zou et al. 2017) (https://doi.org/10.5194/tc-11-2527-2017-supplement) was selected to calculate maximum freeze depth and freeze–thaw duration on a regional scale (Fig. 1). This map is derived from the Temperature at the Top of Permafrost model, combining both modified Moderate Resolution Imaging Spectroradiometer land surface temperatures data and ground observation datasets (Zou et al. 2017). As shown in Fig. 1, SFG is mainly distributed at the southern and eastern TP and Qaidam basin. Some SFG also exists on the continuous permafrost edge and in valleys of mountains (Zou et al. 2017).
To cover the entire period with potential freeze–thaw events, annual values were calculated for each year beginning on 1 September of the preceding year and ending on 31 August of the current year. Maximum freeze depth was selected from all daily lower boundary freeze depth data for each year to represent annual maximum freeze depth. Annual freeze start date was defined as the first day after September when the lower boundary freeze depth was larger than zero. Annual thaw end date was defined as the first day when the lower boundary freeze depth and the upper boundary layer freeze depth were all equal to zero. Annual freeze–thaw duration was defined as the number of days from the freeze start date to thaw end date for each year. The stations were first classified as either permafrost or SFG, depending on the annual thaw end date. If, for the entire record, a station has no thaw end date, that station was classified as permafrost. The remaining stations were considered to be characterized as SFG. Based on the above method, all stations were located in SFG area.
Although the qualities of the freeze depth of SFG have been carefully controlled before releasing, some false and missing data are obtained inevitable. The data were flagged when freeze depth exceeded the maximum range of the frost tube in the original data (CMA 2003). In this study, annual values of maximum freeze depth of individual station were retained if the readings had less than 5 flagged and missing datapoints in the month of maximum freeze depth and the two months adjacent to it, while the annual values of freeze–thaw duration were retained if the readings were not missing data of freeze start date and thaw end date in each year. Based on the above quality control, not all of the meteorological stations in this study have continuous data over this 55-yr period. To expand the number of stations, annual values of individual stations were accepted while there were fully records of more than 30 years during these 55 years. Based on this criterion, 75 meteorological stations’ data were selected to analyze the freeze–thaw changes of SFG on the TP from 1960 to 2014 in this paper. As shown in Fig. 1, there were 36 stations (48%) with maximum soil freeze depth and freeze–thaw duration recorded over 50 years, 17 stations (23%) with maximum soil freeze depth and freeze–thaw duration recorded between 40 and 50 years, and 22 stations (29%) with maximum soil freeze depth and freeze–thaw duration recorded between 30 and 40 years. The remaining 12 stations’ data, which were between 10 and 30 years of record during these 55 years, were used to validate the SFG algorithm in section 3e.
The TP spans a very large elevation from 2000 to 8000 m, and the spatial distribution of climatic and environmental factors varies sharply resulting in a variety of climatic zones. To investigate the internal variation and influence factors in different climate zones, SFG were further discussed in each climatic zone. Considering the zonal air temperature and moisture as an index of climate, it was divided into 12 climatic subregions, consisting of three temperature zones ranging from a plateau subfrigid zone (HI) and a plateau temperate zone (HII) to a plateau subtropical zone (HIII), and four moisture regions, ranging from a humid zone (A), subhumid zone (B), and subarid zone (C) to arid zone (D) (Zheng et al. 2010) (Fig. 1). Table 1 provides details regarding temperature zones, moisture regions, climatic subregion classifications and number of stations for SFG. Three climatic subregions (HIB, HIIC1, and HIIC2) had more than 10 stations, with a maximum of 21 in HIIC1. The number of stations in HIA, HIIA, HIIB, and HIID1 was 4, 3, 8, and 9, respectively. The representativeness of the SFG stations chosen for this paper was reasonable, although there was only one station in HIC1, HIC2, HIID2, and HIIIA, respectively, and no station in HID, because these zones are mainly permafrost region (HIC1, HIC2, and HID) and unfrozen ground region (HIIIA) (Ran et al. 2012; Zou et al. 2017). A total of 12 stations with shorter period in six climatic subregions were used to validate the SFG algorithm.
In this study, linear trends were used to detect trends in maximum freeze depth, freeze–thaw duration, freeze start date, thaw start date and thaw end date of SFG. The Student’s t test was used to assess the significance of trends, within which the >95% confidence level was used to judge the significance. Linear trends were also used to detect trends in annual air temperature and annual precipitation over the entire TP. Multiple linear regression was conducted to explore the possible empirical relationship between the SFG and the influence factors. A p value ≤ 0.05 was considered significant.
The calculated results from the regression function were compared with the observed data based on two statistical measures: mean bias (MB) and correlation coefficient R. The methods used to calculate MB and R are shown as follows:
where Ci is calculated data, Oi is observed data, n is the total number of data, is the average calculation, and is the average observation.
a. Spatial distribution and variations of maximum freeze depth
As shown in Fig. 2a, on average, the multiyear mean maximum freeze depth was 99.9 cm, with the maximum being 287 cm at Ando station in HIC2 zone and the minimum being just 5.5 cm at Chayu station in HIIIA zone over the past 55 years on the entire TP. For more than 2/3 of the stations, the maximum freeze depth was between 50 and 200 cm. Six sites (8%) had a multiyear mean maximum thickness of more than 200 cm and are located mainly in two areas, the subfrigid zone (HI) of the central TP at an elevation of more than 4000 m, and the temperate subarid zone (HIIC1) of the northeastern TP, at a latitude of more than 37°N. The multiyear mean maximum freeze depths were more than 100 cm at 44% of all sites, mainly located north of 31°N. There are also eight sites with multiyear mean maximum freeze depths of no more than 20 cm, located in the subtropical humid zone (HIIIA), the temperate humid zone (HIIA) and the temperate subarid zone (HIIC2) of the south TP at an elevation of less than 3500 m and a latitude of less than 30°N.
As can be seen from the climate map (Fig. 3), the maximum freeze depth shows obvious distribution characteristics of the climatic zones. The freeze depth in the subfrigid zone (HI), with an average of 141.6 cm, was deeper than those in the temperate zones (HII) (85.2 cm) and subtropical zone (HIII) (5.5 cm). Meanwhile, in the same temperature zones, the freeze depth deepened significantly with decreasing moisture. For example, in HI zone, for HIA, HIB, HIC1, and HIC2, the multiyear mean maximum freeze depths were 73.2, 145.2, 215.9, and 287.0 cm, respectively, and in HII zone, for HIIA, HIIB, HIIC1, and HIID1, the multiyear mean maximum freeze depths were 9.3, 54.5, 117.8, and 125.7 cm, respectively.
Over the past 55 years, maximum freeze depth decreased strongly (at α = 0.01) at 2/3 of all sites, with an average rate of −4.9 cm decade−1 and a net change of −26.8 cm on the entire TP (Fig. 2b). This decreasing trend is slightly greater than that of the 71-yr period from 1930 to 2000 in the Eurasian high latitudes, which is −4.5 cm decade−1 (Frauenfeld and Zhang 2011), and it is also greater than that of the same period from 1960 to 2014 in TRSR of the inner TP, which is −4.0 cm decade−1 (Luo et al. 2017a). The result is a bit less than that of the same period on the TP, which is −5.6 cm decade−1 (Fang et al. 2019), because of the different number of meteorological stations. The maximum decreasing trend was −19.3 cm decade−1, also located at Ando station in HIC2 zone. There are 20 sites (27%) where the decreasing trends were more than −8 cm decade−1. About 28% and 24% of all sites, the decreasing trends were between −8 and −4 and between −4 and 0 cm decade−1, respectively. Although 16 sites (21%) showed an increasing trend, the increase was not strong and only three sites (Lenghu, Huangyuan, and Linzhi, 4%) passed the significance test (at α = 0.01). For the Lenghu and Huangyuan sites, the main reason is that these two stations had unusual cooling trends in deep soil temperature in winter (Fang et al. 2019). For the Lizhi site, the main reason is that the precipitation and soil moisture decreased in southeast of the TP (Rodell et al. 2018), which caused the soil to be drier and the freeze depth increased.
As indicated in the climatic zone average (Fig. 4), decreasing trends occur in 10 climatic zones, while only one climatic zone (HIIIA) increases, and that only marginally. The trend of maximum freeze depth is closely related to temperature and moisture. The colder and drier the place, the greater the decreasing trend. The average changing rates in the HI, HII, and HIII zones were −7.4, −4.0, and 0.1 cm decade−1, respectively.
b. Spatial distribution and variations of freeze–thaw duration
Comparing Fig. 2c with Fig. 2a, the spatial distribution of the multiyear mean freeze–thaw duration was similar to that of the multiyear mean maximum freeze depth. As maximum freeze depths increase, the freeze–thaw durations get longer. As shown in Fig. 2c, over the past 55 years, the multiyear mean freeze–thaw duration of SFG on the TP exceeded half a year, at 187.6 days, varying between 314.6 days (about 10.5 months) at Gande station in HIB zone and 63.5 days (about 2 months) at Chayu station in HIIIA zone. For nearly half the sites, the multiyear mean freeze–thaw duration was between 150 and 200 days. For nearly 35% of all sites, the multiyear mean freeze–thaw duration was more than 200 days. There are also six stations (8%), all located in the subfrigid zone (HI) of the central TP, where the multiyear mean freeze–thaw duration was more than 250 days.
From the average of the climatic zone (Fig. 3), the multiyear mean freeze–thaw duration in HI zone, with an average of 230.1 days, was obviously longer than those in HII zone (173.1 days) and in HIII zone (64.5 days). The same temperature zones still show that the freeze–thaw duration in the arid areas is longer than that of the humid areas.
As shown in Fig. 2d, most sites (88%) have decreasing trends in freeze–thaw duration from 1960 to 2014. On average, the decreasing trend is −6.0 day decade−1, or a net change of 33 days for the 55 years. This result is similar to the result in the same period in TRSR of the inner TP, which is −7.5 days decade−1 (Luo et al. 2017a). The maximum decreasing trend was −28.4 days decade−1, located at Qumalai station in the HIB zone. There are four sites (5%) where the decreasing trend was more than −16 days decade−1. At 9%, 21%, 33%, and 19% of all sites, the decreasing trends were between −16 and −12, between −12 and −8, between −8 and −4, and between −4 and 0 days decade−1, respectively. There are still nine sites (12%) showing increasing trends, among which five sites (7%) pass the significance test (at α = 0.01). The Dangxiong, Linzhi, and Daocheng stations, where the increasing trends were more than 8 days decade−1, are all located in southeast of the TP.
Regarding the climatic zone average (Fig. 3a), freeze–thaw duration decreased in all 11 climatic zones. Temperature and moisture are closely linked with the trends of freeze–thaw duration. The average changing rate in the HI zone (−7.4 days decade−1) was higher than those in the HII (−5.5 days decade−1) and HIII zone (−1.3 days decade−1). In the same temperature zone, the rate of changing trend in the arid and subarid region was higher than those in subhumid and humid areas. For example, in HI zone, from HIA, HIB, and HIC1 to HIC2, the changing trends are −1.6, −8.4, −12.1, and −11.1 days decade−1, respectively; in the HII zone, from HIIA, HIIB, HIIC1, HIID1, and HIIC2 to HIID2, the multiyear mean maximum freeze depths are −3.0, −4.0, −4.8, −6.1, −7.6, and −12.4 days decade−1, respectively.
c. Spatial distribution and variations of freeze first date and thaw last date
As shown in Fig. 5a, the multiyear mean freeze start date of SFG occurs from early September to early December (from day 256 to 343) across the TP. On average, the freeze start date was the day 290.6 (about mid-October). Two stations (Zeku and Gande) in HIB zone of the northeastern TP first began to freeze in early September (before day 260). After that, freezing occurred in areas north of the subfrigid and temperate zone (HI and HII). By the end of September (before day 275), 16% of the sites began to freeze. Then, the ground in the temperate zone of the southern TP froze in October and nearly 87% of all sites entered the freezing period by the end of October (before the 305th day). From November, freezing started in the area south of HIII. In early December, all sites were freezing. As illustrated in Fig. 6c, the multiyear mean thaw end date occurs from mid-February to late July (day 41 to 205), with an average of day 112.9, and spans nearly 5.5 months (164 days) across the TP. The Gande station was also the last station of thaw end date. For 45% of all sites, the multiyear mean thaw end date was from late April to mid-May (from the 110th day to the 135th day).
As illustrated in Fig. 3, the freeze start date and thaw end date of SFG also show distribution characteristics of temperature zones. The freeze start date in HI, with an average of day 278 was earlier than those in HII (day 295) and HIII (day 343). Thaw end date in HI, with an average of day 142.8 was later than those in HII (day 102.4) and HIII (day 41.2). Similarly, freeze start date and thaw end date of SFG also show distribution characteristics of moisture zones. The freeze start date in the humid region is later than that in the arid area, while the thaw end date is earlier than that in the arid area.
Trends in multiyear mean freeze start date and thaw end date of SFG at the 75 stations and 11 climatic zones from 1960 to 2014 are shown in Figs. 5b, 5d, and 4. Over the past 55 years, the multiyear mean freeze start date has been slightly delayed, at an average rate of 1.83 days decade−1, and a net change of 10 days for the entire TP (Fig. 5b). At 41% of all sites, significant delays occur with the maximum of 15.32 days decade−1 at Qumalai station in the HIB zone. At 51% of the sites, the trends are not significant. There are also six sites (8%) whose freeze start dates significantly advanced. Among these six sites, Maduo and Linzhi have a larger advancing rate, which are −11.37 and −7.97 days decade−1, respectively. As shown in Fig. 5d, the thaw end date exhibits a statistically significant advancement over the past 55 years, and the advancing trend, with an average rate of −4.10 days decade−1 or a net change of 22.55 days, is significantly greater than that of freeze start date. This means that decreasing trends in the freeze–thaw duration are mainly due to advancing trends in the thaw end date. Similar results were found in TRSR of the inner TP for the same period (Luo et al. 2017a). For the spatial distribution, significant advancing trends (at α = 0.01) of thaw end date occur in 53% of all stations. The maximum advancing rate of thaw end date was −23.49 days decade−1 at Maduo station in HIB zone. There are no significant changing rates of thaw end date in 41% of the sites. In the meantime, there are four sites (Jianzha, Banma, Linzhi, and Daocheng, about 5%) experiencing significant delays in thaw end date.
Regarding climatic zone distribution (Fig. 4), temperature is also the main factor in controlling the trends in thaw end date, but not significant in freeze start date. The trends in thaw end date in the HI zone, with an average of −6.2 days decade−1 are obviously greater than those in the HII zone (average of −3.3 days decade−1) and the HIII zone (average of −1.1 days decade−1). Moisture is not a significantly correlated factor for the trends in freeze start date and thaw end date.
d. Temporal distribution and variations of SFG
From sections 3a–3c, the SFG shows marked spatial distribution with air temperature and precipitation across the TP in the past 55 years. Here the temporal changing of SFG with air temperature and precipitation is discussed for the entire TP.
As shown in Fig. 6, air temperature and precipitation at the 75 stations exhibited statistically significant increases of 0.42°C decade−1 and 8.3 mm decade−1, or approximately 2.3°C and 45.5 mm from 1960 to 2014, respectively. Maximum freeze depth and freeze–thaw duration were both significantly negatively correlated with air temperature and precipitation by the correlation coefficients of −0.85, −0.34, −0.78, and −0.46, respectively. The results indicate that as annual mean air temperature increased, maximum freeze depth and freeze–thaw duration decreased over the 1960–2014 period on the TP. The higher annual precipitation amount corresponds to a shallower maximum freeze depth and shorter freeze–thaw duration over the past 55 years. As shown in Fig. 7, the first freeze date is significantly positively correlated with air temperature and precipitation by correlation coefficients of 0.54 and 0.29, respectively, while last thaw date was significantly negatively correlated with air temperature and precipitation by correlation coefficients of −0.78 and −0.46, respectively. The results also indicate that as annual air temperature and precipitation increased, the first freeze date was later, while the thaw end date advanced, from 1960 to 2014 on the TP.
e. Development of the SFG algorithm
To show the possible empirical relationship of latitude and elevation with SFG, the coefficients of determination between them are mapped in Figs. 8a and 8b. As shown, there were generally good correlations between two variables with latitude and elevation. The multiple linear relationships between multiple-year mean SFG and latitude and elevation were statistically significant (p < 0.0001). With an increase of latitude, maximum freeze depth gradually deepens and freeze–thaw duration becomes longer. With an increase of elevation, maximum freeze depth also deepens and freeze–thaw duration becomes longer. Their fitting equations can be expressed as
where y1 is maximum freeze depth, y2 is freeze–thaw duration, x1 is latitude, and x2 is elevation (m). The coefficient of determination (R) was analyzed to show the fitting performance.
Figures 8c and 8d illustrate that the multiple-year maximum freeze depth and freeze–thaw duration vary with air temperature and precipitation. As illustrated, air temperature and precipitation were the two most important factors affecting SFG on the TP. Maximum freeze depth and freeze–thaw duration were highly correlated with air temperature and precipitation. The multiple linear relationships between them are statistically significant (p < 0.05). The best fit linear regressions to the SFG algorithms are
where y1 is maximum freeze depth, y2 is freeze–thaw duration, x1 is air temperature (°C), and x2 is precipitation (mm).
Figure 9 illustrates the trends in maximum freeze depth and freeze–thaw duration, which vary with geographical factors and climatic factors. As illustrated in Figs. 9a and 9c, the trend in maximum freeze depth is highly negatively correlated with elevation (p < 0.0001) and air temperature (p = 0.0016) but is not significantly correlated with latitude and precipitation. Based on observation and satellite data, a clear signal of elevation-dependent warming (EDW) of air temperature has been found in recent decades over the TP (Gao et al. 2015, 2019; X. Liu et al. 2009). The elevation-dependent decreasing of maximum freeze depth is a response to EDW of air temperature. As shown in Figs. 9b and 9d, the trend in freeze–thaw duration is negatively correlated with elevation and air temperature (p < 0.06) but not significantly. Simulation results have shown EDW of air temperature is only present at relatively low elevation ranges but is absent when the elevation continues to rise (Gao et al. 2018; D. Guo et al. 2016). Recent satellite data even revealed that the rates of annual mean 2 m air temperature rapid decreased above 4500 m on the TP from 2001 to 2015 (Guo et al. 2019). So, long-term trends in freeze–thaw duration are complex, as other local parameters play roles in addition to changes in air temperature and precipitation or latitude and elevation.
Annual maximum freeze depth and freeze–thaw duration were calculated by observed annual air temperature and precipitation [Eqs. (4) and (6)] at 12 stations, and the fitted results are validated against observed data in Fig. 10. The calculation of maximum freeze depth using Eq. (4) is in good agreement with the observations (Fig. 10a). The MB and R between the calculation and the observation were 22.1 cm and 0.93, respectively. Using annual mean air temperature and precipitation [Eq. (6)] the calculated freeze–thaw duration was close to the measured (Fig. 10b), with the MB and R of 23.8 days and 0.87, respectively.
f. Freeze–thaw changes of SFG on a regional scale
Based on 1-km resolution elevation resampled from SRTMDEM and 1-km air temperature and precipitation derived from elevation adjustments from CMA gridded analysis, the multiyear mean maximum freeze depth of SFG from 1960 to 2014 on the TP was assessed by Eqs. (3) and (5). As shown in Figs. 11a and 11b, the spatial distribution of multiyear mean maximum freeze depth calculated by Eq. (3) is very close to the distribution of that by Eq. (5), with positive R of 0.78 and MB of 29.3 cm over the entire TP. Part of the high spatial correlation is due to the high correlation of temperature and precipitation with elevation and latitude. The maximum freeze depth ranges from 0 to 383.5 cm, with an average of 143.5–153.6 cm and varies with climatic zones, elevation, and latitude. The mean regionally accessed result is larger than that from the stations (99.9 cm) because most stations are located at lower altitudes. For the two assessments, half of all SFG region (49%–54%), the maximum freeze depth was great than 150 cm. The largest values for maximum freeze depth occur at the edge of permafrost, where they can exceed 250 cm for 0.5%–3% of all SFG area. The maximum freeze depth is from 200 to 250 cm (11%–18% of all SFG area) in the north, northeastern, center of the TP and TRSR. The maximum freeze depth in the Qaidam basin is lower than in the surrounding area. In the southwest, east and Qaidam basin, maximum freeze depth ranged from 100 to 200 cm in 61%–70% of all SFG area. In the south and southeast regions, maximum freeze depth is less than 100 cm in 16%–18% of all SFG areas.
Figure 12 shows the spatial variability of the maximum freeze depth anomaly for the each decade from 1960 to 2010, with respect to the 1960–2014 mean across the TP. Results show that spatial changes of the maximum freeze depth anomaly ranged from larger than 30 cm to less than 25 cm. For each decade from the 1960s to 2010s, in most areas of the TP, the maximum freeze depth anomaly changed from positive to negative, with an average of 8.5, 4.3, 2.9, −2.6, −12.5, and −14.7 cm, respectively. This means that maximum freeze depth has a decreasing trend during this period, especially in the 2000s and 2010s, although there are only five years of data in the 2010s. In very few areas of the southeast, the maximum freeze depth anomaly changed from negative to positive, which is consistent with the results from stations. The main reason for this is that precipitation and soil moisture increased in the inner TP and decreased in the southeast of the TP (Rodell et al. 2018). A larger decrease of maximum freeze depth occurs in the southwest of the TP and Qaidam basin from the 1960s to 2010s. After 2000, especially since 2010, in TRSR and the southwest and eastern areas of the TP, maximum freeze depth decreases sharply. In addition, a smaller variability of the maximum freeze depth anomaly occurs in the southeastern TP.
Figures 13a and 13b show the spatial pattern of multiyear mean freeze–thaw duration. The multiyear mean freeze–thaw duration calculated by Eq. (4) is very close to that calculated by Eq. (6), with positive R of 0.76 and MB of 24.4 days across the TP. As shown in Figs. 13a and 13b, the freeze–thaw duration ranges from less than 50 to 360 days, with an average of 219.8 to 223.5 days and also varies with climatic zones, elevation, and latitude. The regionally accessed result is also larger than the station result (187.6 days) due to the lower elevation of the stations. In 82%–83% of all SFG regions, freeze–thaw duration is greater than 180 days (six months). The largest value of freeze–thaw duration also occurred on the edge of permafrost, where it exceeded 300 days (about 10 months) over 0.7%–3% of the SFG area. In the north, northeastern, TRSR and central part of the TP, freeze–thaw duration was from 270 to 300 days (9–10 months). In the southwest, east, and Qaidam basin, it ranged from 180 to 270 days (6–9 months). Only in the south and southeast region of the lower altitudes, was the freeze–thaw duration less than 180 days (6 months) over 16% to 17% of the SFG area.
Figure 14 demonstrate the spatial variability of the freeze–thaw duration anomaly for the decades from the 1960s to the 2010s, with respect to the 1960–2014 mean across the TP. Over this time frame, freeze–thaw duration had a clear decreasing trend with an average anomaly of 6.4, 3.6, 2.3, −2.2, −9.5, and −11.4 days in each decade, respectively. The freeze–thaw duration anomaly has been very large since 2000. A large decrease of freeze–thaw duration also occurred in the southwest of the TP, TRSR, and Qaidam basin, while a smaller variability of freeze–thaw duration anomaly occurred in the southeast of the TP.
a. Potential causes of SFG change
The freeze–thaw changes of SFG have been widely regarded as a response to global warming (Frauenfeld and Zhang 2011; Luo et al. 2017a; Yang et al. 2019). It was altered significantly by soil heat and moisture conditions, which were mainly decided by energy balance at the land surface interface, based on the flux of surface net radiation, sensible heat flux, latent heat flux and ground heat flux (Fisher et al. 2016; Guo and Wang 2014; Luo et al. 2009a,b). Because of the model uncertainties and the lack of land–atmosphere interaction data, especially shortage of high spatial–temporal land surface heat flux data, some easily obtained climate variables and geographical factors such as air temperature (index of the sensible and latent heat flux), freezing or thawing index (index of the soil temperature and ground heat flux), precipitation (index of the soil moisture and latent heat flux), snow depth (index of albedo and net radiation), latitude and elevation (index of the net radiation, air temperature and precipitation) have been proposed to represent the main influential factors on SFG (Frauenfeld and Zhang 2011; Frauenfeld et al. 2004; Luo et al. 2017a; Wang et al. 2015). To reveal the direct causes of SFG changes, it is necessary to explore the long changes of surface energy budget on the TP.
Meteorology station data showed surface net radiation decreased slightly over the TP from 1960 to 2010, with the average trend of −0.38 W m−2 decade−1 (Gao et al. 2013). Remote sensing data revealed surface net radiation increased slightly after 2000, especially in the central part of the TP (X. Chen et al. 2014; Han et al. 2017). The latent heat flux presented a weak increasing trend, with the average of 0.70 W m−2 decade−1 during 1980–2003 (Duan and Wu 2008). The most remarkable change was the sensible heat flux that exhibited a significant decreasing trend due to the subdued surface wind speed since 1980s, especially in spring (Chen et al. 2019; Duan and Wu 2008; Han et al. 2017; Yang et al. 2011). Annual mean sensible heat flux was weakened by 2% to 7% decade−1, with the average decrease trend was about −3.4 W m−2 decade−1 from the 1980s to the 2000s (Duan and Wu 2008; Yang et al. 2011). It decreased sharply since 2000, with the trend from −5.3 to −6.8 W m−2 decade−1 on the whole TP (Han et al. 2017). From the above research results, it is clear that although the surface net radiation presented a weak decreasing trend and latent heat flux presented a slightly increasing trend, the ground heat flux over the TP might be increased due to the significant weakened sensible heat flux according to the equation of surface energy balance. Due to the limitation of data, there were few studies on the long-term change of surface heat flux over the TP. In spite of this, studies on long-term soil temperature changes and soil thermal states indirectly revealed a significant increase in surface heat flux (Fang et al. 2019; Jin et al. 2009; Luo et al. 2017a, 2016; Zhou et al. 2019). In our study time, soil temperature in surface (0 cm), shallow (5–20 cm), and deep (40–320 cm) all increased significant with rates of 0.47°, 0.36°, and 0.36°C decade−1, respectively (Fang et al. 2019). The soil thermal diffusivity appeared an overall decreasing trend from 1980 to 2001 (Zhou et al. 2017). The long-term trend of soil enthalpy was mainly increased in Eurasia from 1979 to 2010 (Zhao 2017). In some permafrost region of the TP, observation of soil heat flux presented an increasing trend, with a variation rate from 1.0 to 2.9 W m−2 decade−1 in the recent 10 years (R. Li et al. 2012; D. Luo et al. 2018). Analysis showed that, as soil flux increases 1.0 W m−2, the active layer thickness of permafrost increases 24 cm (R. Li et al. 2012). The increased soil heat flux might be the main direct reason for the decreased freeze depth and freeze–thaw duration of SFG on the TP. More ground heat flux analyses and relevant studies are needed to be further strengthened.
On the other hand, soil water conditions will also affect the soil heat conditions. As shown in section 3, spatial–temporal distributions of freeze depth and freeze–thaw duration were both significantly negatively correlated with precipitation. A significant increasing trend of precipitation and soil moisture were evident on most of the TP in the past 50 years (Fang et al. 2019; Li et al. 2011; Rodell et al. 2018; Shi et al. 2016). The soil moisture increased with the rate of 9.1 mm decade−1 from 1980 to 2012 (Shi et al. 2016). The freeze–thaw process is a buffer to the seasonal changes in soil (B. Chen et al. 2014; Luo et al. 2016). In the freeze phase, the more liquid water freezes into ice, the greater release of latent heat and the slower the cooling of the soil column is, because wetter soil takes up a much larger heat loss to facilitate the phase change of water compared to drier soil. In the thaw phase, an increase in precipitation or soil moisture resulted in an increase in thawing rate, because those temperatures were higher than that of SFG in spring and summer (Luo et al. 2016). So warming and wetting conditions within the soil result in a decrease in freeze depth and a short freeze–thaw duration. Similar results were also found in permafrost regions (Iijima et al. 2010; Li et al. 2019; Neumann et al. 2019).
Different types of land cover, such as snow, vegetation, and soil type, that directly affect the surface energy budget in land surface processes, are all associated with changes of permafrost and SFG in the Northern Hemisphere (Frauenfeld and Zhang 2011; Luo et al. 2009a,b; Shiklomanov 2012; Yang et al. 2019). It is noteworthy that although frozen ground is closely linked with snow because of its high albedo, which plays an important role in surface energy budget in high latitudes, snow may not be a main influence factor of SFG changes on the TP (Luo et al. 2017a; Peng et al. 2017). Satellite data showed that snow is dominantly distributed in high and large mountains on the TP (Che et al. 2019; Hao et al. 2019). In the central and eastern TP, which are also the common places of SFG, snow depth is generally less than 0.5 cm and the number of snow-covered days is no more than 30 days per year (Jiang et al. 2020; Luo et al. 2017a; Peng et al. 2017; Xu et al. 2017). During 1961–2014, the annual average of snow depth and the number of snow-covered days were just 0. 26 cm and 23.78 days, respectively (Jiang et al. 2020). Even though, snow depth and the number of snow-cover days showed a slightly decreasing trend, especially after the year of 2000 (Che et al. 2019; Jiang et al. 2020; Xu et al. 2017). The decreased snow cover which allows surface absorption of more incoming solar radiation might be increased the surface net radiation, which is consistent with the results that the surface net radiation increased slightly after 2000 (X. Chen et al. 2014; Han et al. 2017). The increase of surface net radiation and the decrease of sensible heat flux further increased the soil heat flux since 2000. This conclusion can well explain why the freeze depth and the freeze–thaw duration decreased rapidly since 2000 in this paper.
Vegetation and soil are heterogeneous and have different thermal and hydraulic properties, causing varied behavior in soil heat and moisture transport and freeze–thaw processes (Chang et al. 2012; Luo et al. 2017b, 2009b). Over the TP, regional changes in the structure and functioning of ecosystem have been found in response to recent globe warming and human activities (Piao et al. 2019). Some scholars found degraded grasslands accelerate the rates of change in soil temperatures and soil moisture because of declining vegetative coverage in some frozen ground areas (G. Liu et al. 2009). Since the 1970s, a sharp increase in human activities, such as overgrazing of increasing cattle herds, and imprudent utilization of wetlands, pasturelands and water resources, have been considered as another environmental factors in a general degradation of frozen ground on the TP (Jin et al. 2009). Vegetation of observation field in meteorological stations is typical alpine meadow or alpine steppe on the TP. It is relatively stable by artificial maintenance with canopy height no more than 0.20 m in summer and no more than 0.05 m in winter (CMA 2003). So vegetation may not be a main influencing factor of SFG changes in this paper.
b. Possible feedback of SFG change to climate and environment
In the process of land–atmosphere interaction, the change of frozen ground also produces important regional climatic and environmental effects (Chen et al. 2014a,,2017; Poutou et al. 2004; Viterbo et al. 1999; Yang et al. 2019). Simulation showed that the soil freeze–thaw process “buffers” the seasonal changes of soil and near-surface temperatures (B. Chen et al. 2014, 2017; Viterbo et al. 1999). Compare to no freeze–thaw process, the multiyear regional averages of shallow layer soil temperature (16.5–28.7 cm) and air temperature (2 m) increased by about 2.05° and 0.59°C in winter, while it declined by about 0.32° and 0.05°C in summer on the whole TP (B. Chen et al. 2014). Sensible heat flux and latent heat flux were both larger than simulations without soil freeze process in winter while they were both smaller than simulations without soil thaw process in summer, which leads to an increase of ground source heat in winter and a decrease of ground source heat in summer (B. Chen et al. 2014). With the decreasing of freeze depth and freeze–thaw duration on the TP, the buffering effect of soil freeze–thaw process will decrease and surface heat source will continue to reduce, resulting in important climatic feedback.
SFG change plays an important role on the ecosystem of the TP. The conditions of soil hydrothermal are curial factors for vegetation germination, growing and wilting. Soil thaw depth is an important parameter because it affects the timing of cultivation and seeding in early spring (Gao et al. 2019). During the past several decades, spring vegetation phenology has been significantly advanced and the plant-community structure has been changed (Piao et al. 2019). It is obvious that the decreasing of soil freeze–thaw duration will lead to the increasing of vegetation growing season on the TP. The spring vegetation phenology advancing and regional vegetation greening should be linked with the declining of freeze depth and freeze–thaw duration.
The change of soil carbon in frozen ground is one of potentially most significant carbon–climate feedbacks (Cheng et al. 2019; Tarnocai et al. 2009). On the whole TP, warming had significantly increased vegetation productivity, which consequently lead to an enhanced carbon sink (Piao et al. 2009, 2019; Zhang et al. 2014). While the variations of carbon sources or sinks had obvious highly uncertain mainly because of the high spatial heterogeneity of soil properties and lack of information regarding deep-layer soil processes (Piao et al. 2019). Study show that the temperature sensitivity of the ecosystem respiration rate in the SFG areas was greater than the permafrost areas, indicating that the carbon emissions in the SFG areas were more sensitive to warming (Mu et al. 2017). A comprehensive understanding of freeze–thaw changes will contribute to the in-depth study of soil carbon changes. Relevant feedback and interaction mechanism between SFG, vegetation and climate are needed to be further strengthened.
This study has focused on investigating the freeze–thaw changes of SFG on the TP from 1960 to 2014. The following conclusions can be drawn:
Meteorology data from 75 stations showed that the multiyear average for mean maximum freeze depth was 99.9 cm, the freeze–thaw duration was 187.6 days, the freeze start date was day 290.6, and thaw end date was day 112.9 with distribution characteristics along climatic zones. A decreasing trend in maximum freeze depth (−4.9 cm decade−1) and freeze–thaw duration (−6.0 day decade−1) occurred over the TP from 1960 to 2014. The freeze start date was delayed (1.8 day decade−1) and the thaw end date (−4.1 day decade−1) was advanced significantly. The average changing rate in the subfrigid zone was larger than those in the temperate and subtropical zones, while the rates were higher in the arid and subarid region than those in subhumid and humid areas.
The freeze–thaw changes of SFG significantly affected by soil hydrothermal conditions on the TP could be assessed by elevation and latitude or by air temperature and precipitation, due to their high correlations. The estimated multiyear mean maximum freeze depth, averaged from 143.5 to 153.6 cm, and the multiyear mean freeze–thaw duration, averaged 219.8–223.5 days, showing a larger result than those developed from average meteorological stations data only, because most meteorological stations are located at lower altitudes. In half of all SFG regions, the maximum freeze depth was greater than 150 cm. For more than 82% of all SFG regions, the freeze–thaw duration was greater than six months.
The estimated maximum freeze depth and freeze–thaw duration showed a clearly decreasing trend from 1960 to 2014 on a regional scale over the entire TP, especially after 2000. The largest decrease of freeze depth and freeze–thaw duration occurred in the southwest, TRSR, and Qaidam basin, while a smaller increase of them occurred in very few areas of the southeast.
The increased soil heat flux might be the main direct reason for the freeze–thaw changes of SFG on the TP. Warming and wetting conditions with the soil result in a significant decrease in maximum freeze depth and freeze–thaw duration in the most area of the TP, while drying soil result in a slightly increase of them in the southeast of TP. More ground heat flux analyses and interaction mechanism between SFG, vegetation, and climate are needed to be further strengthened.
This work was supported jointly by the Second Tibetan Plateau Scientific Expedition and Research (STEP) Project (2019QZKK0105), the National Natural Science Foundation of China (Grants 91537104, 41975096, and 41775016), the Global Water Futures Programme of the Canada First Research Excellence Fund, and the Canada Research Chairs Programme. The authors would like to acknowledge CMA and CMDC for supplying the freeze depth data for the TP. We wish to thank the editor and three anonymous reviewers for their helpful and constructive comments.
Denotes content that is immediately available upon publication as open access.