Characteristics of near‐surface soil freeze–thaw status using high resolution CLM5.0 simulations on the Tibetan Plateau

Soil freeze–thaw alternation is a natural characteristic of the Tibetan Plateau (TP), and plays an important role in surface energy balance and eco‐hydrological processes. The soil freeze–thaw process on the TP has changed significantly owing to global warming, affecting the alpine ecosystem structure and function. This study used high‐resolution atmospheric forcing datasets to drive the Community Land Model version 5.0 (CLM5.0) to simulate the near‐surface soil freeze–thaw status between 1979 and 2020. The simulated results were compared with in situ observations, and then the spatiotemporal distribution of the freeze start‐date (FSD), freeze end‐date (FED), freeze duration (FD), and thaw duration (TD) at a depth of 0.1 m were analyzed. The Nash–Sutcliffe efficiency coefficients (NSEs) of FSD, FED, FD, and TD between simulations and in situ observations were 0.77, 0.90, 0.98 and 0.92, and the correlation coefficients of FSD, FED, FD, TD were 0.97, 0.99, 0.99 and 0.98, respectively. The spatial distribution of FSD and TD was characterized by gradually increasing from northwest to southeast while FED and FD exhibited the opposite characteristics. FSD, FED, FD, and TD changed at an area‐mean rate of 1.1, −1.4, −2.5, and 2.5 days decade−1, respectively. This study provides an important reference for analyzing and predicting the changes in near surface soil freeze–thaw status on the TP under the warming climate.


| INTRODUCTION
The Tibetan Plateau (TP) is the largest plateau in China and the highest in the world, with an average altitude over 4000 m; it is known as the "Third Pole" (Qiu, 2008). It has formed a unique atmospheric circulation pattern due to its unique geographical position, topographic conditions, and dynamic and thermal cycles, which have an important influence on the climate of China and the world (C. G. Wang et al., 2017;M. X. Yang et al., 2010), and is a sensitive and key area of global climate change (X. Wang et al., 2018). The TP has the largest and most extensive frozen soil area in the middle and low latitudes of the world, covering an area of approximately 1.06 Â 10 6 km 2 (Liu et al., 2022;Zou et al., 2017). Studies have proven the magnifying effect of TP warming, that is, the warming range of surface temperature in the TP is higher than that in China, at the same latitude, and the global average (You et al., , 2021. Because the TP is sensitive to climate change, the soil freezethaw process has changed significantly with global warming. The soil freeze-thaw process plays an important role in regulating land-atmosphere energy exchange, hydrological processes, atmospheric circulation, the carbon and nitrogen cycle, and ecosystem functions on the TP (Guo et al., 2011;Shen et al., 2015;C. H. Wang et al., 2020). Because the exchange of energy and water between land and atmosphere mainly occurs at the near-surface (M. X. Yang et al., 2006), the freeze-thaw cycle of near-surface soil is extremely sensitive to climate change (Li et al., 2012). Therefore, accurate quantification of spatiotemporal distribution characteristics of near-surface soil freeze-thaw status over the TP would lead to the better understanding of the response of the freeze-thaw process to climate change and provide scientific evidence for exploring the mechanism of regional climate change in the context of climate warming.
In this study, we adopted high-resolution atmospheric forcing datasets to drive the Community Land Model version 5.0 (CLM5.0) to simulate the soil temperature data at a depth of 0.1 m of the near-surface on the TP from 1979 to 2020. The simulation results were evaluated, and the spatiotemporal distribution characteristics of the freeze start-date (FSD), freeze end-date (FED), freeze duration (FD), and thaw duration (TD) of the nearsurface soil on the TP were analyzed. This study provides new insights into the complex freeze-thaw cycle of the TP driven by the almost certain future intensified climate warming.
2 | DATA AND METHODS

| Model
The CLM, released by National Center for Atmospheric Research, belongs to the land surface module of the Community Earth System Model (Dickinson et al., 2006;Swenson et al., 2012), and has been developed to the latest version of CLM5.0. CLM5 includes several components related to land biogeophysics, hydrology, biogeochemistry, human dimensions, and ecosystem dynamics (Lawrence et al., 2019). High spatial land surface heterogeneity in CLM is represented as a nested subgrid hierarchy. Each grid cell is composed of multiple land units, each land unit can have a different of columns, and each column has multiple plant functional types (PFTs). In this study, simulation experiments were forced by high-resolution atmospheric forcing datasets, including a daily precipitation product with 3 km Â 3 km (Jiang et al., 2023) and other variables (temperature, pressure, dew point temperature, wind speed, incident shortwave radiation, and downward longwave radiation) from ERA5-Land (Muñoz-Sabater et al., 2021). The precipitation data were produced by dynamic downscaling derived from combing a short-term high-resolution Weather Research and Forecasting Model (WRF) simulation (Zhou et al., 2021) with ERA5 reanalysis (ERA5-WRF) and then merging gauge observations from more than 9000 rain gauges using the Climatology Aided Interpolation and Random Forest (RF) methods (Jiang et al., 2023). Compared with the data of all rain gauges on the TP, the precipitation product is typically unbiased, with a root mean square error of 4.5 mm day À1 , and a correlation coefficient of 0.84, showing that this dataset was better than the widely-used global datasets (Jiang et al., 2023). The spatial and temporal resolutions of ERA5-Land are 9 km and hourly, respectively. Before the historical simulation, a long-term cyclically spin-up run (1200 years) was prepared to allow the simulated soil water and temperature to reach their equilibrium state. After that, a normal transient simulation was initiated and run from 1979 to 2020. The simulation experiments have a spatial resolution of 0.1 Â 0.1 and temporal resolution of 1800 s covering the study area (73 -105 E, 25 -41 N). All the atmospheric datasets have been interpolated into the same spatiotemporal resolution using model configuration before the simulations. It is noted that precipitation and solar radiation data are interpolated into the model time step using the nearest and coszen methods, respectively, while adopting the linear method for the other forcing fields.

| In situ observations
To validate the suitability of the simulations, a longterm dataset of integrated land-atmosphere interaction observations (Ma et al., 2020) and a time-lapse observation dataset of soil temperature and humidity on the TP (Dente et al., 2012;Su et al., 2011Su et al., , 2013van der Velde et al., 2012), including four sites provided by the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn), were employed: BJ (91.9 E, 31.37 N), NAMORS (90.98 E,30.77 N), QOMS (86.95 E,28.36 N),and Naqu (92.31 E,31.11 N). The periods of available observations for the four sites

| Methods
Based on the definition of key parameters for the freezethaw cycle in previous studies (Guo et al., 2011;Guo et al., 2018;Guo & Wang, 2014), FSD was defined as the date when the daily soil temperature changed from above 0 C to below 0 C, and FED was defined as the date when the daily soil temperature changed from below 0 C to above 0 C. To avoid the influence of random changes in soil temperature and small disturbance on the calculation of FSD and FED, when the soil temperature was below 0 C for three consecutive days, the first day was defined as the FSD. When the soil temperature was above 0 C for three consecutive days, the first day was defined as the FED. FD was the number of days between the FSD and FED. TD was the number of days between the FED and FSD. In this study, the FD of the current year was defined as the time from the FSD of the current year to the FED of the next year, and the TD of the current year was defined as the time from the FED of the current year to the FSD of the current year.
The Nash-Sutcliffe efficiency coefficient (NSE) (Nash & Sutcliffe, 1970), which reflects the effect and accuracy of model simulations, was adopted to measure the model prediction effect (Jia et al., 2021). It is calculated as follows (Mbaye et al., 2015): where x obs is the observed value, x sim is the simulated value, x obs is the average of the observed value, and n is the number of samples.

| Validation of the model
To validate the suitability of the simulations, the diurnal variations of soil temperature at a depth of 0.1 m were compared with observations from these four stations (Figure 1a-d). According to available data at the four Based on the comparison of FSD, FED, FD and TD at Naqu, BJ, NAMORS and QOMS stations, the simulated values were close to the observed values (Table 1), and the four freeze-thawing indexes were practically distributed around y = x (Figure 1e). As shown in Table 1, the NSEs of FSD, FED, FD and TD were 0.77, 0.90, 0.98 and 0.92, respectively. The correlation coefficients of FSD, FED, FD and TD were 0.97, 0.99, 0.99, and 0.98, respectively, and significance all exceeded the α = 0.01 level. The mean biases between simulations and in situ observations were À6.4, À5.2, 1.2, and 4.6, respectively. In addition, we used the Taylor diagram to evaluate the accuracy of model simulations (Figure 1f). The results demonstrated that the simulated FSD, FED, FD and TD values were good, and the correlation coefficients were all over 0.9 and close to the REF point (Figure 1f). Thus, CLM5.0 based simulations can capture the temporal dynamics of surface soil temperature well and be used to explore the characteristics of near-surface soil freeze-thaw status over the TP.

| Spatiotemporal distribution characteristics of soil freeze-thaw
As shown in Figure 2, the spatial distribution characteristics of FSD and TD were similar and exhibited a gradient that increases from northwest to southeast (Figure 2a,d). The amounts of FSD and TD were higher in certain parts of northeast than the surrounding areas; it is likely that the altitude of this area is 1000-2000 m lower than that of surrounding areas, resulting in higher temperature. As a result, the region began to freeze later and thaw earlier. Furthermore, the spatial distribution characteristics of FED and FD were similar, showing a gradually decreasing gradient from northwest to southeast (Figure 2b,c). The amounts of FED and FD were smaller in the corresponding regions with higher FSD and TD because the higher temperature causes thawing earlier, causing the FD to be shorter. Figure 3 shows the trends of FSD, FED, FD, and TD at a soil depth of 0.1 m on the TP from 1979 to 2020. FSD and TD show increasing trends, while FED and FD show decreasing trends. The simulated FSD for the majority of the grids displayed a delayed trend, particularly in the northwest and eastern TP (Figure 3a). The simulated FED for the majority of the grids displayed a statistically advancing trend, except for central part of the TP (Figure 3b). The simulated FD for the majority of the grids displayed a statistically significant shortening trend and the simulated TD displayed an extending trend, except for certain grids in the central and southeastern parts of the TP (Figure 3c,d). The insignificant changes over parts of the central TP were also found by other studies (Guo & Wang, 2014;Luo et al., 2020), which may be related to the insignificant trends of precipitation.
Based on the annual rates of variations of FSD, FED, FD and TD at a soil depth of 0.1 m on the TP (Figure 3e-h), FSD displayed increasing trends, FED showed decreasing trends, FD showed significant shortening trends, and TD showed significant extending trends, respectively. The trends in FSD were delayed with an area-mean value of 1.1 day decade À1 . The trends in FED advanced with an area-mean value of 1.4 day decade À1 . The trends in FD were shortened with an areamean value of 2.5 day decade À1 . The trends in TD were extended with an area-mean value of 2.5 day decade À1 . As a result, the area-mean FSD was delayed by 4.5 days, and the FED advanced by 5.7 days, both of which resulted in a shortening of the FD by 10.2 days, and TD was extended by 10.2 days during entire study period.

| Diurnal variations of soil temperature
Based on the diurnal variations of soil temperature at a depth of 0.1 m in four stages (1980-1989, 1990-1999, 2000-2009, 2010-2019) on the TP (Figure 4), it could be seen that, soil temperature showed obvious diurnal variation characteristics, showing a "unimodal," which was consistent with the periodic change of air temperature. Previous studies have shown that the TP has been the fastest warming region in China, more than twice the global warming rate in the same period, which would certainly affect the near-surface soil temperature and further affect the soil freeze-thaw status (Palazzi et al., 2017;L. B. Yan et al., 2016;Y. P. Yan et al., 2020;You et al., 2017). Consequently, the FSD were delayed and the FED were advanced, which could be seen from the F I G U R E 4 The diurnal variations of soil temperature at a depth of 0.1 m in four stages are as follows: 1980-1989 (blue), 1990-1999 (purple), 2000-2009 (green), and 2010-2019 (red), shown by solid lines. The dashed lines corresponding to these four colors are FSD (right), FED (left), and the black dashed line indicates the 0 C soil temperature. movement of the dashed lines. And the trends of advance of FED were greater than the trends of delay of FSD. This was consistent with previous findings. The reason could be that soil temperature in spring is more sensitive to climate change than that in autumn and winter, and the specific reasons need to be further studied.

| CONCLUSIONS
In this study, high-resolution atmospheric forcing datasets were used to drive the CLM5.0 to simulate the nearsurface soil freeze-thaw status between 1979 and 2020. Model simulations were first evaluated using in situ observations and then used to investigate the changes in the near-surface soil freeze-thaw status on the TP during the study period.
Comparison between the simulations and in situ observations shows that the NSEs of FSD, FED, FD, and TD are 0.77, 0.90, 0.98, and 0.92, and the correlation coefficients are as high as 0.97, 0.99, 0.99, and 0.98, respectively. Thus, the simulations can accurately capture the temporal dynamics of soil temperature. Further analysis based on the simulations reveals that the FSD and TD have similar distribution characteristics, gradually increasing from northwest to southeast. The FED and FD exhibit the opposite characteristics. The FSD and TD show increasing trends for most regions of the TP, while FED and FD show decreasing trends. During the entire study period , FSD was delayed by 4.5 days, FED was advanced by 5.7 days, FD was shortened by 10.2 days, and TD was extended by 10.2 days, respectively.
Our results indicate that climate change has a significant impact on the near-surface soil freeze-thaw, providing an important reference for analyzing and predicting the changes in near-surface soil freeze-thaw status on the TP under the warming climate. Our next step is to investigate the changes of soil freeze-thaw in deeper soil layers and their impacts on ecological processes. AUTHOR CONTRIBUTIONS Qing Peng: Conceptualization; formal analysis; investigation; methodology; visualization; writingoriginal draft. Binghao Jia: Methodology; project administration; supervision; writingreview and editing. Xin Lai: Methodology; writingreview and editing. Longhuan Wang: Methodology; writingreview and editing. Qifeng Huang: Methodology; writingreview and editing.