Trends of freezing period and its main cause on the Qinghai-Tibetan Plateau from 1961 to 2018

The ecosystems of Qinghai-Tibetan Plateau (QTP) are very sensitive to climate change because of their unique structure and function. However, little attention has been paid to variations in cold non-growing season. In this study, based on daily mean temperature from 63 meteorological stations throughout the QTP during the period 1961 − 2018, the spatial and temporal variations in the freezing period (FP) were investigated. The FP was defined as the period between the date of the first autumn freeze and the date of the first spring thaw in the second year. Understanding how the FP changes are imperative in predicting future climate change and decision-making for implementing ecological conservation on the plateau. The results showed that the start of freezing period (SFP) exhibited a pronounced increasing trend with a rate of 0.0704 days year−1 and the end of freezing period (EFP) showed an obviously decreasing trend with a rate of − 0.2537 days year−1 at the regional scale. The length of freezing period (LFP) presented a significant negative trend at a rate of − 0.3256 days year−1 for regional scale, which was mainly attributed to the earlier EFP. Spatially, later mean SFP, earlier mean EFP, and shorter mean LFP mainly occurred in the south of the QTP, covered the plateau temperate semi-arid (HIIC2) and plateau temperate humid/sub-humid (HIIAB1). For interannual trends, greater delayed SFP and greater advanced EFP were mainly observed in the south and east of the plateau. Furthermore, this study found that the variations in the SFP, EFP, and LFP were highly dependent on the elevation with EFP and LFP are positively correlated with elevation, while SFP is negatively correlated with elevation. At the regional scale, mean annual temperature was positively correlated with SFP and negatively correlated with EFP and LFP. Increasing temperature dominated interannual variation in FP on the plateau.


Introduction
Global average surface temperature increased approximately 0.85 °C during the period 1880 − 2012, and the warming rate over the past 50 years (0.12 °C/10 yr) was nearly twice that since 1880 (IPCC 2013). Moreover, regional differences in temperature variability were remarkable, and climate change was considered to be more significance at high latitudes and high altitudes (Post et al. 2018;Wang et al. 2014). In the Northern Hemisphere, springtime events were advanced more rapidly at higher latitudes due to the enhancement of warming with latitude (Post et al. 2018). The high-altitude regions exhibited greater warming than their low altitude counterparts in same latitudes (Wang et al. 2014). Owing to temperature controls vegetation growth by affecting photosynthesis and respiration (Drake et al. 2016), such climate warming may have profoundly affected the composition, structure, and function of terrestrial ecosystems (Gao et al. 2013;Shen et al. 2021;Wang et al. 2015). Remote sensing monitoring showed continued glacier shrinkage and considerable permafrost degradation over the high latitudes and high altitudes in recent years driven by increasing temperature (Bolch et al. 2010;Zheng et al. 2020), which have induced a series of ecological and environmental effects, such as the release of large amounts of deep permafrost carbon and declining water table (Wang et al. 2020e). Furthermore, the area and aboveground biomass of alpine marshlands have been significantly changed. Marshland loss largely reduces carbon sequestration and increase greenhouse gas emissions (Shen et al. 2020;Wang et al. 2021b).
The freeze-thaw cycle is a key feature of the frozen soil ecosystem, and the frozen soil ecosystem is very sensitive to climate change (Man et al. 2019). To study the effects of climate warming in alpine ecosystems, most previous studies have focused on soil freeze-thaw processes (Wang et al. , 2020d. Changes in soil freeze-thaw processes may alter hydrothermal properties of the soil, which in turn may alter the regional climate and directly or indirectly affect the growth of alpine vegetation, such as changing the activity of vegetation roots (Cleavitt et al. 2008;Liu et al. 2020a;Niu et al. 2019). With the increasing temperature, annual soil freezing days has decreased over the past half century (Henry 2008). Moreover, the start of vegetation growing season has been advanced and the full phenological period has extended Liu et al. 2020b), and ecosystem productivity and terrestrial carbon sinks have been increased significantly (Piao et al. 2007;Williams et al. 2019). Others studies have observed a delay in budburst when the chilling requirements of the species are not fulfilled (Laube et al. 2014). Although increased spring temperatures advance leaf unfolding and flowering, increased winter temperatures may counteract this effect , and previous studies may overestimate the advancement of spring phenology with winter waring. Different opinions existed in studies on effect of changes in freezing period. To quantify this effect, it is essential to examine variations in the end of freezing period (EFP) and the start of freezing period (SFP), and identify contribution of these two variations to increasing length of freezing period (LFP).
Nevertheless, soil freeze-thaw processes require consideration of heat transfer in the soil, and glacial melt requires consideration of air temperature lapse rates, whereas using climate data to analyze the freezing period (FP) is straightforward, highly indicative, and strongly correlated with soil freeze-thaw period and glacial melt period (Odland et al. 2017;Peng et al. 2016). Previous studies mainly focused on the pattern of increasing temperature, and found that warming in winter was greater than that in other seasons (Kreyling 2010;You et al. 2013). Changes in temperature mainly affect vegetation growth by altering some key biotemperature thresholds, such as 0 °C, 5 °C and 10 °C (Zhao and Wu 2016), which associated with phenology and enzymatic activity in plants. Among them, 0 °C is a key limit temperature that can transform water from liquid to solid. When temperature is below 0 °C, the surface soil begins freeze and vegetation stalls growth. In recent decades, most attentions were primarily paid on thermal growing season Zhu et al. 2020). Majority of studies on cold non-growing season mainly focused on the decrease of the number of frost days, with rare attention paid to changes in the SFP, EFP and LFP (Meng et al. 2021;Park et al. 2021;Wang et al. 2021a). Park et al. (2021) found that changes in frost days altered vegetation phenology, which reduce vegetation adaptability. (Laube et al. 2014) suggested that lack of chilling in spring could lead to a delay in plant budburst and cause a change in the chronological order of budburst. Therefore, understand the spatial and temporal variability of non-growing season with temperature below 0 °C and investigate the relationship with the thermal growing season and soil freeze-thaw period is important for predicting future climate change.
The Qinghai-Tibetan Plateau (QTP), which is located at low and middle latitudes and has an average altitude higher than 4000 m, is known as the "Third Pole" of the earth (Qiu 2008). Due to the complex topographic and extremely low temperature, the QTP is the largest distribution region of glaciers and permafrost outside the polar regions (Zou et al. 2017). To adapt the unique alpine grassland ecosystem, the QTP has developed corresponding production patterns (grazing) which are vulnerable to climate change (Luo et al. 2015). Observations indicated that the surface temperatures on the QTP increased by about 1.8 °C over the past 50 years (Wang et al. 2008). Concurrent with climate warming, the FP was also changed, which further altered ecosystems and hydrological processes. Several studies showed that the number of frost days on the QTP were significantly decreased in recent years (Ning et al. 2012;Shi et al. 2018). However, little was known whether regional pattern of the LFP was changed as that of growing season, the SFP and the EFP contributed equally to LFP changes, and a relationship existed between the FP (LFP, SFP and EFP) and topography.
The objective of this study was to investigate the spatial and temporal trends in FP across the QTP with climate warming. 0 °C was selected as a threshold, and the FP is defined as the period between the date of the first autumn freeze and the date of the first spring thaw in the second year (Watkins 1991). Based on the FP, we analyzed the spatial and temporal trends of the SFP, the EFP, the LFP, and their relationship with elevation and mean annual temperature on the QTP over the period 1961 − 2018 to reveal FP response to climate warming. This study will be helpful for predicting future climate change and decision-making for implementing ecological conservation on the plateau.

Data
Daily mean temperature dataset on the QTP was provided by the China Meteorological Administration for the period of 1961 − 2018. Quality check and homogeneity test were undergone to ensure high data quality ( Fig. 1) (China Meteorological Administration 2003). To maintain data completeness and continuity, the missing records of observation stations used for this study were less than 2% and the missing records were interpolated using linear regression from neighboring stations to achieve complete time series over the same period. The interpolated daily records in the dataset concentrated in 1968 and 1969, no more than 1% of the total records. The homogeneity test was undertaken using the RHtestV4 software (http:// etccdi. pacifi ccli mate. org/ softw are. shtml) to eliminate effects of station migrations and equipment updates (Vincent et al. 2012;Wang et al. 2010). Finally, 63 stations on the QTP were adopted, mainly located in the center and east and sparsely in the west.

Methods
To reduce effect of extreme values, the 10-day is commonly adopted as a period to identify the threshold of growing season, non-growing season or frost period (Dong et al. 2012;Zhao et al. 2021;Deng et al. 2020;Kamyar et al. 2020). In this study, the FP was defined as the period between the date of the first autumn freeze and the date of the first spring thaw in the second year. The SFP was defined as the first day of the first 10-day period with a daily mean temperature less than 0 °C in autumn and the EFP was defined as the first day of the first 10-day period with a daily mean temperature greater than 0 °C in spring of the second year. The LFP was the number of days from SFP to EFP, in which 5 days break was permitted.
The simple linear regression method was used to analyze the relationships between the regional mean SFP, EFP, LFP and elevation, and the relationships between the interannual SFP, EFP, LFP and mean annual temperature. Meanwhile, the Mann-Kendall test was adopted to assess the significance of the observed trend (Kendall 1990;Mann 1945). Mann-Kendall test is a nonparametric test, which does not require the sample to follow a standard distribution pattern and is not disturbed by sporadic outliers. Mann-Kendall test often is applied in non-normally distributed data, such as hydrology and meteorology (Sneyers 1990).
With vast area, the climatic conditions of the QTP are complex and the ecosystem are diverse so the FP across the QTP has huge variety. In order to explain this variety, the QTP was divided into 11 eco-regions with unique characteristics based on the ecological regionalization scheme of Zheng (2008) (Fig. 1, Table 1). The division of the ecoregions considers the difference in temperature and humidity over the QTP.

Changes in SFP
Over the period 1961 − 2018, the mean SFP on the QTP ranged from 269.4 Julian days to 363.7 Julian days with Fig. 1 Distribution of 63 meteorological stations in this study across the QTP and location of the QTP (stations that occurred in the text have been shown on the figure) a regional average of 310.9 Julian days ( Fig. 2(a)). Spatially, there is obvious heterogeneity in the mean SFP. The relatively early SFP primarily occurred in HIC1 and HIB1 regions, which are located in the middle region of the QTP with the average elevation greater than 4000 m. Of these stations, the earliest SFP was observed at the Wudaoliang station (269.4 Julian days). A later SFP was mainly found in HIIAB1, mostly greater than 330 Julian days. The latest SFP was observed at the Nyingchi station (363.7 Julian days). Furthermore, the analysis of correlation indicates a strongly negative correlation between SFP and elevation (R 2 = 0.2881, P < 0.01), and earlier SFP with high elevation ( Fig. 3(a)). At the regional scale, temperature and SFP showed a significant positive correlation (R 2 = 0.1104, P < 0.01), with the higher mean annual temperature contributed to the later SFP ( Fig. 6(a)).
Figures 4(a) and 5(a) illustrate the interannual variations in SFP on the QTP during 1961 − 2018. At the regional scale, SFP increased at a rate of 0.0704 days year −1 (R 2 = 0.0868, P < 0.01), with a variation range of − 0.075 to 0.574 days year −1 . For the spatial distribution of trends in SFP, most of stations (90.5%) showed a positive trend, of which 40 stations (63.5%) passed significance at the 0.05 level. Larger increasing trend in SFP (> 0.2 days year −1 ) was mainly located in HIIC1, HIIC2, and HIIAB1 regions, and the largest rate (0.574 days year −1 ) was observed at the Lhasa station. Additionally, only seven stations exhibited a negative trend and one (Lenghu station) passed the significance test at the 0.05 level. Such changes in SFP may be influenced by urbanization effects. For example, the increasing trend at the Lhasa station (0.574 days year −1 ) was significantly larger than at the Shigatse station (0.227 days year −1 ) and Zedang station (0.217 days year −1 ), which are close in latitude, longitude, and altitude.

Changes in EFP
Generally, EFP accompanied by snow and ice melt and frozen soil thaw, and is considered the beginning of the growing season (Dong et al. 2012). During the whole study period, the variation range in the mean EFP across the QTP was from 5.4 Julian days to 133.7 Julian days with a regional average of 71.9 Julian days ( Fig. 2(b)). Spatially, an earlier EFP (< 40 Julian days) was primarily observed in HIIC2 and HIIAB1 regions with lower elevation, of which the earliest EFP was found at the Jiulong station (5.4 Julian days). A relatively late EFP (> 100 Julian days) mainly occurred in HIC1, HIC2, and HIB1 regions with higher elevation. On the other hand, the mean EFP exhibited a pronounced positive correlation with elevation (R 2 = 0.3574, P < 0.01), and later EFP with high elevation (Fig. 3(b)). A closely negative correlation was examined between temperature and EFP (R 2 = 0.2907, P < 0.01), suggested earlier EFP associated with higher mean annual temperature ( Fig. 6(b)). At the regional scale, the EFP on the QTP exhibited a decreasing trend at a rate of − 0.2537 days year −1 during the last 58 years (R 2 = 0.5502, P < 0.01) (Fig. 4(b)), with a variation range of − 0.857 to 0.005 days year −1 (Fig. 5(b)). Among the 63 stations considered in this study, 62 stations (98.4%) experienced a decreasing trend in EFP except for Tuotuo River station, of which 59 stations (93.7%) were found to be significant at the 0.05 level. Moreover, the rate of this decrease varied considerably across different eco-regions. For example, EFP experienced a greater decreasing trend (< − 0.3 days year −1 ) mainly in HIIAB1 and HIIC1 regions, where elevations are relatively low. The largest decreasing rate was observed at the Deqin station (− 0.857 days year −1 ).

Changes in LFP
Over the period 1961 − 2018, the mean LFP ranged from 7.9 days to 228.3 days with a regional average of 127.3 days (Fig. 2(c)). Spatially, the longer LFP (> 160 days) was mainly found in HIC1, H1C2 and H1B1 regions, while the shorter LFP (< 80 days) was mostly observed in HIIC2 and HIIAB1 regions. Figure 3(c) illustrates the relationships between the mean LFP and elevation. Similar to the mean SFP and the mean EFP, the mean LFP was also closely depend on the elevation, and positively correlated with elevation (R 2 = 0.3256, P < 0.01). For interannual variation, LFP exhibited a pronounced negative correlation with mean annual temperature (R 2 = 0.3277, P < 0.01), and shorter LFP with higher temperature (Fig. 6(c)).
Changes in LFP are controlled by changes in SFP and EFP. We further analyzed the interannual trends in LFP and found that LFP across the QTP exhibited a decrease at a rate of − 0.3256 days year −1 (R 2 = 0.5498, P < 0.01) at the regional scale (Fig. 4(c)), with a variation range of − 1.222 to − 0.033 days year −1 (Fig. 5(c)). Spatially, there were also regional differences in LFP trends, although 62 stations (98.4%) experienced a decreasing trend excluding Minhe station, and 56 stations (88.9%) of them passed 0.05 significance level test. In addition, the greater decrease trend (< − 0.5 days year −1 ) in LFP was mainly recorded in HIIAB1 and HIIC1 regions.

Discussion
In the present study, we analyzed spatial and temporal variations in SFP, EFP, and LFP derived from 63 meteorological stations on the QTP during the period 1961 − 2018. SFP,

Variations in SFP
Overall, SFP exhibited a pronounced increasing trend with a regional average rate of 0.0704 days year −1 for the past half century. Spatially, larger increase in SFP (> 0.2 days year −1 ) was mainly located in HIIC1, HIIC2, and HIIAB1 regions with lower elevations. Most urban stations throughout the QTP were located below 3500 m, thus localized urbanization and land use change impact on surface warming cannot be ignored (You et al. 2008). Human-induced greenhouse gas emissions were significant, and the urban heat island effect exacerbates the local temperature increase (Dong et al. 2012).
As the temperature below 0 °C, the photosynthetic rate decreases, thus inhibiting the growth of vegetation (Ensminger et al. 2006). Alpine meadow, as the dominant ecosystem on the QTP, is highly sensitive to climate change due to its unique structure and function (Wang et al. 2020c). Usually, SFP relatively coincides with the end of the thermal growing season and the first frost day in autumn (Dong et al. 2012;Zhang et al. 2014). Dong et al. (2012 indicated that the end of the thermal growing season has been delayed by 1.47 days/decade from 1960 to 2009. Through a study of frost-free season, Zhang et al. (2014) suggested that the first frost day in autumn occurred later and delayed by 1.5 days/decade during 1960 − 2010. Moreover, many studies reported a delayed in autumn phenology on the QTP in recent years, which can slow down the rate at vegetation degradation in autumn, Cong et al. 2017;Liu et al. 2016). An et al. (2020) found the daily minimum air temperature can be a key factor for controlling autumn phenology in alpine grasslands, which a 1 °C increase in the mean autumn minimum temperature during the optimum length period may induce a delay of 1.6 to 9.3 days in the middle senescence date across the alpine grasslands.
Although a certain time lag exists, it is generally accepted that soil temperature and air temperature are highly correlated (Odland et al. 2017). Wang et al. (2020a) attributed the change in seasonally frozen ground to the rising of minimum air temperatures in winter. Numerous studies suggested that the onset date of soil freeze was clearly lagged on the QTP, which further supports the result of the delayed SFP (Jia et al. 2020;Li et al. 2012).

Variations in EFP
EFP showed a significant decreasing trend with a rate of − 0.2537 days year −1 at the regional scale, which indicated an earlier EFP. EFP usually consists with snow and ice melt and seasonally frozen ground thaw, and is close to the start of the thermal growing season and the last frost day in spring (Dong et al. 2012;Zhang et al. 2014). Recently, earlier spring soil thaw and advanced spring phenology on the QTP were widely reported (Dong et al. 2012;Li et al. 2012). For example, Li et al. (2012) found a trend toward earlier onset date of soil thaw by approximately 14 days, and later onset date of soil freeze by approximately 10 days during 1988 − 2007. Dong et al. (2012) found an advanced start of growing season at the rate of − 1.82 days/decade from 1960 to 2009, which is a relatively consistent with the results of this study of − 2.537 days/decade.
In temperature-limited areas, plants need a critical level of forcing to trigger spring phenology, which could be satisfied earlier as the advance of EFP (Wang et al. 2019). The early of the EFP advances leaf unfolding and flowering time in perennials. However, the advanced EFP and delayed SFP in warming winter could reduce chilling accumulation, increasing the heat requirement and attenuate the advance of spring phenology ). Therefore, current phenology models with effects of only spring warming can overestimate the advance of spring phenology.

Variations in LFP
As important indicators of the FP, changes in SFP and EFP together determine the change in LFP. Over the last 58 years, LFP presented a significant negative trend at a rate of − 0.3256 days year −1 for regional scale on the QTP, and the effect of advanced EFP (− 0.2537 days year −1 ) greatly contributed to shortened LFP relative to the delayed SFP (0.0704 days year −1 ), especially in HIIAB1 and HIIC1 regions with lower elevations. Zou et al. (2014) found that the winter warming of the southern QTP was greater than other regions at a rate of larger than 0.6 °C per decade from 1979 to 2010. It is consistent with the greater decrease in LFP over HIIAB1 region. Changes in LFP can directly affect soil freeze-thaw cycles, inducing permafrost degradation and early seasonally frozen ground thawing, which alters soil water storage causes the release of large amounts of soil carbon (Song et al. 2020;Wang et al. 2009).
The decrease in LFP contributes an extended thermal growing season. Dong et al. (2012) reported that the regional average growing season length on the QTP extended at the rate of 3.29 days/decade. Over the past half century, the growing season duration has been significantly lengthened for the majority of the Northern Hemisphere due to climate warming, which can stimulate vegetation growth (Piao et al. 2007;Sun et al. 2020a). Using remote sensing data, Sun et al. (2020b) found that advanced spring phenology more controlled the length of the growing season than autumn phenology on the QTP, which further confirms that the shorter LFP was mainly attributed to the earlier EFP. Notably, Piao et al. (2007) discovered that a 1-day extension in length of growing season led to an increase in annual gross primary productivity (GPP) of 5.8 gC m −2 yr −1 , and an increase in net primary productivity (NPP) of 2.8 gC m −2 yr −1 in the Northern Hemisphere over the past 2 decades. Yang et al. (2011) found that changes in cropping systems have resulted a significant increase in grain yield across Qinghai Province during the past 50 years with climate warming. It was the result of the shortened FP and the increased heat accumulation during the growing season. A shortening of the LFP also contributes to the elevational advance of alpine tree line, which eventually increases biomass accumulation in terrestrial ecosystems (Baker and Moseley 2007).

Correlation of FP with elevation and temperature
This study found that the variations in the SFP, EFP, and LFP were highly dependent on the elevation. EFP and LFP are positively correlated with elevation, while SFP is negatively correlated with elevation. This phenomenon could be explained by the temperature lapse rate that the surface temperature decreases with the elevation (Fang 1992). With rising elevation, agricultural production and ecosystems on the QTP are more strongly constrained by low temperature. However, a shortened LFP with an extended growing season may shift the plateau cropping boundary to higher elevations, thus increasing the agricultural production area (Dong et al. 2012). At higher elevations, ecosystems are more sensitive and vulnerable, and smaller changes in FP may induce greater effects (Hao et al. 2021). Furthermore, regional differences in geomorphological conditions, ecological environments and urbanization effects may also contribute to the spatial heterogeneity of SFP, EFP and LFP (Dong et al. 2012). Meanwhile, we further analyzed the relationship between FP and the interannual variation of temperature. This study found that increasing annual mean temperature was positively correlated with SFP and negatively correlated with EFP and LFP, which indicated the delay of SFP and the advance of EFP. Therefore, we can infer that the increasing temperature may have contributed to the rapid decrease in LFP throughout the QTP.

Implications and uncertainties
Vegetation response to climate change existed seasonal differences. Different trends for SFP and EFP resulted in varied responses of ecosystems. Understanding the spatial and temporal patterns and main causes of FP changes is important for predicting future ecosystem responses in the context of climate warming. However, only 10-day period temperature determines the SFP and EFP in this study, while vegetation growth requires a prolonged and stable temperature. For example, when the growing season begins, a sudden drop in temperature below 0 °C could have a more serious impact on vegetation growth and even lead to death. Therefore, calculation methods need to be improved by considering more factors to define the FP in future research.

Conclusions
At the regional scale, EFP and LFP experienced a decreasing trend, while SFP exhibited an increasing trend across QTP over the period 1961 − 2018. Spatially, the trends in SFP, EFP and LFP showed the large heterogeneous. Greater decreasing trends in EFP and LFP were observed in HIIAB1 and HIIC1 regions. More obvious delay trends in SFP were primarily recorded in HIIC1, HIIC2, and HIIAB1 regions. Furthermore, changes in SFP, EFP and LFP are influenced by elevation, positively correlated with EFP and LFP, and negatively correlated with SFP. For interannual variation, increasing annual mean temperature was positively correlated with SFP and negatively correlated with EFP and LFP.