Impact of Initial Soil Conditions on Soil Hydrothermal and Surface Energy Fluxes in the Permafrost Region of the Tibetan Plateau

The initial state of land surface models (LSMs) and regional climate models plays an important role in regional and global climate predictions due to their “memories” (Koster and Suarez, 2001; Hagemann and Stacke, 2015). Since the atmosphere memory is too short, the seasonal to decadal predictability of the atmosphere depends on slowly varying underlying surface components, including oceans, glaciers, sea ice, soil, snow cover, and the terrestrial biosphere (Koren et al., 1999; Rodell et al., 2018). Representative of the storage of soil heat and water, the initial soil temperature (ST) and soil moisture (SM) have long-term influences on land-atmosphere interaction on different time and spatial scales. The ST integrates the effects of energy exchange at the air-ground interface, continuously recording the energy balance at the surface (Beltrami, 2002). The initial ST memory can last several months to one year or more and it is stronger in deep soil layers than in subsurface soil layers (Hu and Feng, 2004; Yang and Zhang, 2016; Qiu et al., 2021; Xue et al., 2021). Winter and spring ST can influence summer precipitation (Tang et al., 1987; Hu and Feng, 2004; Xue et al., 2012, 2022), further noting that arid and semiarid areas in the northwestern part of China have a longer ST memory than the humid and semi-humid areas of the southeastern part of China (Yang and Zhang, 2016). The role of initial SM in the land–atmosphere interaction and the climate system also has been widely investigated by observations and simulations. Soil can “remember” the wet or dry conditions that control the partitioning of the available energy into latent and sensible heat fluxes and link the energy, water, and carbon fluxes (Entin et al., 2000; Guo et al., 2011; Koster et al., 2011; Liu and Mishra, 2017; Song et al., 2019; Zhang et al., 2021). The memory of the initial SM can last 2-3 months and demonstrate spatial and seasonal variations. (Vinnikov et al., 1996; Entin et al., 2000; Wu and Dickinson, 2004). Compared to atmospheric conditions, the initial SM anomaly has a greater influence on the decay rate of SM anomalies (Song et al., 2019). Therefore, SM is considered to be an important predictor in monthly to interannual climate prediction (Seneviratne et al., 2006; Song et al., 2019).

Permafrost represents a complex, integrated response to the water and energy balance at Earth’s surface (Nelson, 2003). As permafrost underlies about one-quarter of the exposed land in the Northern Hemisphere (Zhang et al., 2008), it is imperative to study and accurately model its behavior under current and future climate conditions (Elshamy et al., 2020). In the permafrost region, SM consists of two parts, soil liquid water (SLW) and soil ice (SI). SI, together with freeze-thaw cycles, can also act as agents of long-term climate memory (Matsumura and Yamazaki, 2012). Liquid water and ice have very different physical and hydraulic properties. For example, the heat capacities of liquid water and ice are 4.2 MJ m^–3 K^–1 and 1.9 MJ m^–3 K^–1, respectively, while their thermal conductivities are 0.57 W m^–1 K^–1 and 2.29 W m^–1 K^–1, respectively. As the active layer changes and freeze-thaw processes occur, liquid water and ice can frequently transition and strongly affect the thermal and hydraulic properties of the soil (Johansen, 1975; Côté and Konrad, 2005; Luo et al., 2009b, 2017; Dai et al., 2019b). A large amount of latent heat of SI phase change is released or absorbed, which has a specific influence on the ground heat flux and also links the surface energy budget during the freeze-thaw processes (Luo et al., 2009a; Lawrence et al., 2012; Chen et al., 2014; Dobiński, 2020; Wang et al., 2021; Cuesta-Valero et al., 2023). Moreover, due to the “heat reservoir” effect of deep soils (Alexeev et al., 2007; Stevens et al., 2007), the memories of initial ST and SM have been further amplified by ice-rich soils in permafrost regions (Zhang et al., 2000; Elshamy et al., 2020; Ji et al., 2022). In contrast to SLW, studies of the SI memory together with freeze-thaw processes and its impacts on climate are relatively scarce due to the observational and model uncertainties in permafrost region.

The Tibetan Plateau (TP) hosts the largest and thickest frozen soil at the middle and low latitudes, including permafrost and seasonally frozen ground (Cheng et al., 2019; Wang et al., 2019; Yang et al., 2019; Luo et al., 2020), which play important roles in the surface energy budget and in the local and global climate through the memory of frozen soil (Wang et al., 2003; Chen et al., 2014; Yang and Wang, 2019). In recent years, to better simulate the hydrothermal characteristics and climate effects of frozen soil on the TP, the input of atmospheric forcing data, soil hydrothermal parameterization schemes, and initial values of soil have been widely investigated and tested in LSMs. With the establishment of observation networks and the emergence of various assimilation and other new techniques, the accuracy of the input atmospheric forcing data in time and space has been continuously improved upon in the TP (Huang et al., 2014; He et al., 2020; Ma et al., 2020; Zhao et al., 2021). Much work has been done to improve the soil hydrothermal parameterization schemes of permafrost and seasonally frozen ground, including the development of a soil thermal conductivity parameterization scheme, gravel parameterization scheme, organic matter parameterization scheme, etc. (Luo et al., 2009a, b, 2017, 2018; Pan et al., 2017; Yang et al., 2018; Dai et al., 2019a; Gao et al., 2019; Li et al., 2020; Deng et al., 2021; Yang et al., 2021). Furthermore, many soil characteristics, including soil texture, gravel content, and organic matter content, were also investigated and applied to land surface modeling (Luo et al., 2009b; Chen et al., 2012b; Dai et al., 2013; Shangguan et al., 2013; Fang et al., 2016; Liu et al., 2021b). These studies dramatically improve the accuracy of the simulation results. However, there are still significant gaps in the model's accuracy in simulating the hydrothermal state of the soil, especially in the permafrost region (Li et al., 2020; Yang et al., 2021; Ma et al., 2023). The inaccuracy of the initial soil conditions is still an important factor regarding the large deviations in land surface modeling (Song et al., 2019; Ji et al., 2022). In general, SI cannot be observed in field measurements. To avoid the biases of the initial SI, the simulation of seasonally frozen soil is usually started in summer (Luo et al., 2009a, b, 2017). However, deep soil ice in permafrost exists year-round on the TP, making it particularly important to set the initial SI conditions when simulating permafrost.

In permafrost, deeper soil profiles have larger hydraulic and thermal memory thereby requiring more effort to properly initialize an LSM (Elshamy et al., 2020). Spin-up becomes a plausible way to provide initial states for LSM simulations when they are not supplied with reliable initial conditions (Gao et al., 2015; Luo et al., 2017; Ji et al., 2022). Studies have shown that many LSMs reached a state of rough equilibrium within several months or 1 to 2 years (Takle et al., 1999; Cosgrove et al., 2003; Luo et al., 2017). Gao et al. (2015) found a slightly longer spin-up time (4 years) due to the cryospheric processes in the permafrost region in the TP. Because of the long soil memory of permafrost, the strategy of creating the spin-up forcing including the total spin-up length and cycling scheme inevitably affects subsequent simulations (Schlosser et al., 2000). A recent study reveals that the impacts of total spin-up length on the resulting initial field of ST are relatively small in the permafrost region on the TP (Ji et al., 2022). Even when the LSMs reached equilibrium, they still experienced some gaps between the real soil states and the spun-up ST and SM due to the inaccuracies of forcing variables, soil parameters, and modeled physical processes. One study shows that the spun-up ST differed from the observed value by 1°C–3°C, while the spun-up SLW differed from the observed value by 0.02–0.16 m³ m⁻³ at one permafrost site in the TP (Gao et al., 2015). Inaccurate initial ST and SM obtained by spin-up methods is another important factor related to inaccurate simulated results of LSMs on the TP.

To improve the simulation ability of LSMs in permafrost regions on the TP, accurate initial ST and SM are necessary. If the spin-up method is relatively limited, is there a better way to obtain a better initialization of ST and SM in the permafrost region on the TP? In this paper, we first apply one observation site in the TP to evaluate the performance of the Community Land Model version 5.0 (CLM5.0) in simulating ST, SLW, SI, and surface energy fluxes. To address the bias in CLM5.0, we conducted three modified initial soil schemes experiments alongside two long-term spin-up experiments to improve the accuracy of CLM5.0. Also, the difference in simulation results between these two types of methods was analyzed. The remainder of this paper is organized as follows. Section 2 describes the data, model, default and modified initial soil schemes; section 3 presents the experiment design and statistical method in this study; section 4 analyzes the impacts of different initial soil schemes and long-term spin-up methods on simulated soil hydrothermal and surface energy fluxes in permafrost areas. The conclusions and discussions are presented in section 5.

The simulations were conducted for one typical permafrost site, the Xidatan station (XDT, 94.13°E, 35.72°N, Alt.: 4538 m a.s.l) in the Three River Source Region (TRSR) on the TP (Fig. 1). It is located in a faulted valley at the northern foot of the eastern Kunlun Mountains near the northern boundary of the permafrost zone in the TP and is representative of the characteristics of discontinuous permafrost (Liu et al., 2021a; Zhao et al., 2021).

Figure 1.  The Tibetan Plateau (TP), Three River Source Region (TRSR) and the location of the Xidatan (XDT) station.

The XDT station is one of the automatic meteorological stations in the permafrost monitoring network. Meteorological data, including precipitation, air temperature, radiation, relative humidity, and wind speed, were recorded every half hour with a CR1000/CR3000 data acquisition instrument from 2004 (Zhao et al., 2021). The ST and SM data were collected from an integrated observation system in the active layer. The ST for the active layer monitoring system was measured at different depths with a 105T/109 thermocouple probe (Zhao et al., 2021). The SLW was measured by a Hydra SM sensor (Zhao et al., 2021). The above meteorological, ST, and SM data can provide basic information for both the analysis of physical processes and LSM research.

This station is mainly covered by an alpine meadow, and the vegetation coverage is about 85% (Yang et al., 2021). The annual mean air temperature varied from −4.4°C to −2.4°C during the period between years 2005 and 2016, with a multiyear average of −3.6°C (Liu et al., 2021a; Zhao et al., 2021). During the same period, the annual accumulated precipitation varied between 341 and 492 mm, with an average of 384.5 mm, of which summer precipitation accounted for 50%–70% (Liu et al., 2021a; Zhao et al., 2021). The average active layer thickness of the XDT site from year 2005 to 2017 was 1.54 m (Liu et al., 2021a).

The 2-m air temperature, wind speed, atmospheric pressure, humidity, precipitation, and downward longwave and shortwave radiation at the XDT station was used as the atmospheric forcing data with a time step of 30 minutes (Xiao and Qiao, 2020). The forcing data spanned 1 July 2015 to 31 August 2016, for a total of 14 months. The ST and SLW in 8 layers (5 cm, 20 cm, 40 cm, 80 cm, 120 cm, 140 cm, 180 cm, 240 cm, respectively) at the XDT site were selected as validation data for the model (Wu, 2020). The dataset was provided by National Cryosphere Desert Data Center (http://www.ncdc.ac.cn). The soil texture data were obtained from the measured soil texture data in this area (Luo et al., 2009b).

The Community Land Model (CLM) is the land component of the Community Earth System Model (CESM) and is used in several global and regional modeling systems (Lawrence et al., 2019). CLM5.0 is the latest version of the CLM. It builds upon progress made in CLM4.0 and CLM4.5. Scientific topics that have driven CLM5.0 development include the following: (1) improved understanding of carbon and nitrogen cycle interactions and their influence on the long-term trajectory of the terrestrial carbon sink; (2) assessment of the response and vulnerability of ecosystems to climate change and disturbances (human and natural) and the possibility for ecosystem management to mitigate climate change; (3) quantification of the role of terrestrial processes in diurnal to interannual weather and climate variability including influence on droughts, floods, and extremes; (4) establishment of the availability of water resources under climate variability and climate change; (5) quantification of key land feedbacks to climate change including the permafrost climate-carbon feedback and snow- and vegetation-albedo feedbacks; (6) examination of the impact of model structure and parameter uncertainty and exploration of parameter optimization techniques (Lawrence et al., 2019).

Compared to CLM4.5 and previous versions, CLM5.0 incorporates many new and updated processes and parameterizations: (1) dynamic weights of land units and plant functional type distribution; (2) an updated structure and parameterizations for soil hydrology; (3) new parameterizations for snow, glaciers, and ice sheets; (4) a plant hydraulics and hydraulic redistribution; (5) revised carbon dynamics and nitrogen cycling; (6) a demographically structured dynamic vegetation model (Lawrence et al., 2019). CLM5.0 and previous versions have been widely used to simulate land surface processes of the seasonally frozen ground over the TP, showing good applicability (Chen et al., 2012a, 2017; Guo and Wang, 2013; Fang et al., 2016; Luo et al., 2017; Deng et al., 2020). A new study shows that CLM5.0 can also be used to simulate the surface energy budget and soil hydrothermal regime in the permafrost regions of the TP (Ma et al., 2023).

The soil column in CLM5.0 can be discretized into an arbitrary number of layers. It provides 4 types of soil layer settings for users, including 15 soil layers (3.5 m), 23 soil layers (3.5 m), 53 soil layers (10 m), and 25 soil layers (8.5 m). The default vertical discretization uses 25 soil layers. The layer structure of the soil is described by the node depth ($ {z}_{i} $(m)) the thickness of each layer, $ ({\Delta z}_{i} $(m)), and the depths at the layer interfaces $ ({z}_{i} $(m)). When the vertical discretization uses 15 soil layers, 10 soil layers are hydrologically and biogeochemically active, and the deepest 5 layers are only included in the thermodynamical calculations (Lawrence et al., 2008, 2011)

In CLM5.0, the initial schemes of ST and SM are the same as in CLM4.5. Soil points are default initialized with ST $ {T}_{{\mathrm{i}}} $, for all 15 soil layers, of 274 K (0.85°C), and the default volumetric soil water content (including liquid water and ice) $ {\theta }_{{\mathrm{i}}} $, for the first 10 soil layers, is 0.15 m³ m⁻³, and for the last 5 soil layers of bedrock, is 0.0 m³ m⁻³. The mass of SLW ($ {w}_{\rm liq,i} $) and SI ($ {w}_{\rm{ice,i}} $) are initialized as:

where $ {\rho }_{{\mathrm{liq}}} $ and $ {\rho }_{{\mathrm{ice}}} $ are the densities of liquid water and ice (1.0 kg m⁻³ and 0.9 kg m⁻³), and $ {T}_{{\mathrm{f}}} $ is the freezing temperature of water (273.15 K).

It can be seen from the above description that the initial ST $ {T}_{{\mathrm{i}}} $, for each soil layer in CLM5.0, is greater than the freezing temperature, which is inconsistent with the actual situation of permafrost. Moreover, the initial states of SM in CLM5.0 is also inconsistent with the actual situation of permafrost. Since the initial ST is set higher than the freezing temperature, the initial mass of SI is zero and all soil water is SLW. From Eqs. (1) and (2), if the initial ST is set lower than the freezing temperature, the initial mass of SLW is zero; at that moment, all soil water is SI. It can be seen from the above analysis that the default initialization conditions of ST and SM in CLM5.0 is too simple and it does not adequately describe the actual conditions of permafrost.

In this paper, to assess the impact of the initial soil states on simulated ST, SM, and surface energy fluxes in permafrost areas on the TP, the initial soil conditions, including ST, SLW, and SI were modified.

In permafrost regions, ST fluctuates above or below the freezing temperature (273.15 K) in the active layers, while remaining consistently below freezing temperature in the permafrost layers throughout the year. In modified initial soil schemes, to better characterize ST in permafrost, the initial ST can be higher or lower than the freezing temperature. Here, we set the initial ST to either 274 K (0.85°C) or 272.3 K (–0.85°C). Moreover, we modified the initial mass of SLW and SI. In two of the experiments designed in this paper, the mass of SLW and SI are initialized as:

where $ {\theta }_{{\mathrm{i}}} $ is the initial volumetric soil water content (including liquid water and ice water), $ {\theta }_{{\mathrm{liqmax,i}}} $ is the maximum liquid water content when the ST is below the freezing temperature. The latter is calculated by:

where $ {L}_{{\mathrm{f}}} $ is the latent heat of fusion (3.337 × 10⁵ J kg^–1K^–1), $ g $ is the gravitational acceleration (9.806 m s^–2). $ {\theta }_{{\mathrm{sat,i}}} $ is the porosity of the soil, $ {\varPsi }_{{\mathrm{sat,i}}} $ is the saturated hydraulic potential, and $ {B}_{{\mathrm{i}}} $ is the exponent. They are defined as weighted averages based on mineral soil and organic soil (Lawrence and Slater, 2008):

In Eqs. (6–8), $ {\theta }_{{\mathrm{sat}},\mathrm{m}\mathrm{i}\mathrm{n},{\mathrm{i}}} $, $ {\Psi }_{{\mathrm{sat,min,i}}} $ and $ {B}_{{\mathrm{min,i}}} $ are the porosity, saturated hydraulic potential, and the exponent at soil layer $ {\mathrm{i}} $ for mineral soil, respectively. They are calculated by soil texture (Cosby et al., 1984):

In Eqs. (6–8), $ {\theta }_{{\mathrm{sat}},\mathrm{o}\mathrm{m}} $, $ {\varPsi }_{{\mathrm{sat,om}}}, $ and $ {B}_{{\mathrm{om}}} $ are the porosity, saturated hydraulic potential, and the exponent for organic soil, respectively. They are set to be 0.9 m³ m^–3, –10.3 mm and 2.7, respectively (Farouki, 1981; Letts et al., 2000). $ \mathrm{\%}sand $ and $ \mathrm{\%}clay $ are the percentage of sand and clay in the soil and they are obtained from the measured soil texture data in this area (Luo et al., 2009b). The soil organic matter fraction is defined as $ {f}_{{\mathrm{om,i}}}= {{\rho }_{{\mathrm{om,i}}}}/{{\rho }_{{\mathrm{om,max}}}} $, where $ {\rho }_{{\mathrm{om,i}}} $ is the soil organic density for the soil layer and $ {\rho }_{{\mathrm{om,max}}} $ is the bulk density of peat (130 kg m^–3).

The default vertical discretization uses 25 soil layers with a depth of 8.5 m in CLM5.0. Recent studies have indicated the need to have a deeper soil column in LSMs to properly capture changes in freeze-thaw cycles and active-layer dynamics in permafrost areas (Alexeev et al., 2007; Lawrence et al., 2012). However, a deeper soil column implies larger soil hydraulic memory and, more importantly, thermal memory that require proper initialization to be able to capture the evolution of past, current, and future changes (Elshamy et al., 2020). The wrong setting may increase the deviation of the simulation, as little is known about the soil hydrothermal state below 2.4 m without observed data. By testing different soil layers in CLM5.0, we found that for the simulation of permafrost on the XDT site, the simulated ST and SM were closer to the observation by selecting 15 soil layers with a depth of 3.5 m in CLM5.0. Meanwhile, this setting can be adapted to simulations of permafrost at this site since the average active layer thickness of the XDT site was about 1.5 m over the most recent 10 years (Liu et al., 2021a). Therefore, in this paper, the vertical discretization used 15 soil layers with a depth of 3.5 m.

We designed six simulation experiments in this study (Table 1). TEST1 was run by the default initial schemes of ST, SLW, and SI, while TEST2 to TEST4 was run by modified initial schemes of ST, SLW, and SI. As shown in Table 1, the initial ST was set as the default value of CLM5.0 (274 K or 0.85°C) in TEST1 and it was set as 272.3 K (–0.85°C) in TEST2. In these two experiments, the volumetric soil water content $ {\theta }_{{\mathrm{i}}} $ was set as the default value of CLM5.0 (0.15 m³ m⁻³) and the mass of SLW and SI were initialized by Eqs. (1) and (2). Since the initial ST was higher than the freezing temperature in TEST1, the initial mass of SI was zero and all soil water was SLW; in TEST2, since the initial ST was lower than the freezing temperature, the initial mass of SLW is zero and all soil water was SI. In TEST3, the initial ST was set as 272.3 K (–0.85°C), and the mass of SLW and SI were initialized by the modified soil conditions (see Eqs. (3) and (4). In TEST4, in the active layers (1st–7th layer for XDT), the initial ST was set as 274 K (0.85°C), while in the permafrost layers (8th–10th layer for XDT), the initial ST was set as 272.3 K (–0.85°C). The mass of SLW and SI were also calculated by Eqs. (3) and (4) in TEST4. It can be concluded from the calculation of Eq. (5) that in TEST3 and TEST4, the maximum liquid water content at 272.3 K can reach 0.17 m³ m⁻³ for the XDT station. Here, we set the initial volumetric soil water content, $ {\theta }_{{\mathrm{i}}}, $ to 0.25 m³ m⁻³ in TEST3 and TEST4. From TEST1 to TEST4, the model was run with 14 months of forcing data (from 1 July 2015 to 31 August 2016). In practice, we found CLM4.5 reached a state of rough equilibrium in the first two months when it was forced by observational data in the seasonally frozen soil on the TP (Luo et al., 2017). In this paper, the simulations of the last 12 months (1 September 2015 to 31 August 2016) for these four experiments were compared with the observations (equal to a two-month spin-up in the model).

In another two experiments (TEST1-SP and TEST2-SP), the model was designed for long-term spin-up. A study showed that LSMs need such a long term to achieve strict convergence of ST and SM (Ji et al., 2022). Here, we designed 100-year spin-up experiments. The forcing data of spin-up experiments which usually comes from gridded datasets or observed data has a certain effect on the simulation results. Despite many available reanalyses of gridded datasets that can provide long-term forcing data, most of them are substantially biased on the TP (Zhang et al., 2020; Zhou et al., 2021). Meanwhile, it is very difficult to obtain continuous long-term observational data as the long-term forcing data in the permafrost region on the TP. As a compromise, the long-term spin-up simulation can be driven by spin-up cycling with time and correspond to an average meteorological condition. There are three common spin-up cycling schemes for generating the forcing 1) observations from a single year, 2) observations from multiple years or their averages, and 3) all available observations (Ji et al., 2022; Li et al., 2023). Here, we use a single-year cycling scheme. In TEST1-SP and TEST2-SP, the model was repeatedly run for over 100 years by forcing data from one year (1 July 2015 to 30 June 2016), then the simulation continued for 14 months using 14 months of forcing data (from 1 July 2015 to 31 August 2016). Also, the simulations of the last 12 months (1 September 2015, to 31 August 2016) in TEST1-SP and TEST2-SP were used for comparison with the observations as well as the simulated results of TEST1 to TEST4. In short, the initial ST and SM conditions in TEST1-SP and TEST2-SP were the same as TEST1 and TEST2, but for 100 years of spin-up using one year of forcing data (1 July 2015 to 30 June 2016).

In this paper, the active layer of permafrost was divided into two periods, the freezing period and the thawing period. To avoid the influence of random fluctuations on ST, the freezing and thawing onsets were determined based on the average of the daily STs taken for five consecutive days; they were below and above 0°C, respectively, and varied depending on the soil layer (Guo and Wang, 2014; Yang et al., 2021).

The simulation results from six tests were compared with the observed data using the root-mean-square error (RMSE), mean bias error (MBE), and the Nash–Sutcliffe efficiency (NSE). The methods used to calculate RMSE, MBE, and NSE are given as follows:

where $ N $ is the time series length, $ {M}_{i}(i=\mathrm{1,2},3,\cdots ,N) $ is the simulation, $ {O}_{i}(i=\mathrm{1,2},3,\cdots ,N) $ is the observation, and $ \overline{O} $ represents the average of observations. Normally, MBE (smaller values are desireable) has a minimum value of zero under the hypothetical situation that the model is capable of perfect (long-term) simulations of the system. An NSE = 1 indicates a perfect correspondence between simulations and observations; an NSE = 0 indicates that the model simulations have the same explanatory power as the mean of the observations; and an NSE < 0 indicates that the model is a worse predictor than the mean of the observations (Schaefli and Gupta, 2007; Knoben et al., 2019).

The equilibrium conditions in TEST1-SP and TEST2-SP were defined as follows. If ST (year N+1) minus ST (year N) is within 1°C and SLW (year N+1) minus SLW (year N) is within of 0.1 m³ m⁻³, the soil is said to reach a rough equilibrium condition. If ST (year N+1) minus ST (year N) is within 0.1°C and SLW (year N+1) minus SLW (year N) is within of 0.01 m³ m⁻³, the soil is said to reach a strict equilibrium condition. Figure 2 shows the convergence in ST and SLW during the first 10 years of 100-year spin-up time for 10 soil layers for TEST1-SP and TEST2-SP. We found that the equilibrium conditions were established relatively quickly in these two long-term spin-up experiments. As shown in Figs. 2a and b in TEST1-SP, the equilibrium conditions of ST and SLW were established in the second year for all soil layers regardless of whether the constraint of a rough equilibrium condition or a strict equilibrium condition is used. As shown in Figs. 2c and d, given a rough equilibrium condition in TEST2-SP, the equilibrium conditions of ST and SLW for all soil layers were also quickly established in the second year; however, there were slight differences when given a strict equilibrium condition; ST needed 3 years and SLW needed 7 years. From Figs. 2a and c, though ST reached the strict threshold of 0.1°C in the second or third year, it did not change consistently. As shown in Fig. 2a, after 6 years of spin-up, the ST in layers 1 to 6 was not stable, and ST (year N+1) minus ST (year N) was larger than the threshold of 0.1°C in some years. The possible reasons for this instability are that the summer of 2016 was the end of the spin-up cycle, and it was used to reinitialize the model in the summer of 2015. This approach makes the whole study highly dependent on the climate within this year. But ST has a reasonable fluctuation within 1°C until the end of the 100-year spin-up. To summarize, in both experiments, CLM5.0 only took one year to reach a rough equilibrium in ST and SLW at the XDT site.

Figure 2.  Convergence in soil temperature (a and c) (ST, units: °C) and soil liquid water (b and d) (SLW, units: m³ m⁻³) at 10 soil layers in TEST1-SP (a and b) and TEST2-SP (c and d).

Table 2 shows the values of ST, SLW, and SI after 100 years of spin-up in the TEST1-SP and TEST2-SP as the initialization of the simulation. It can be seen that after the CLM5.0 reaches equilibrium, the difference between ST, SLW and SI at each layer in these two experiments was really small. Regardless of whether the simulation started with an above-freezing temperature of 0.85°C or a below-freezing temperature of –0.85°C, after 100 years of spin-up, the ST of 8 layers was higher than the freezing temperature at the initialization of the simulation, and the ST in the deepest two layers was equal to the freezing temperature which was different from the actual situation of permafrost layers. After reaching equilibrium, the difference between the ST at each layer was really small, while the SLW and SI at each layer were almost the same in these two spin-up experiments. Since the initialization ST was higher than the freezing temperature from the surface to 140 cm, the initial SI was zero for these layers. As can be seen from the above, the initial fields of ST, SLW, and SI obtained by the long-term spin-up method was relatively limited at the XDT site, especially in the permafrost layers.

Figure 3 and Table 3 show a comparison of the observed and six experimental results of ST in different soil layers at the XDT site. It can be seen from the observation (red line in Fig. 3) that in summer and autumn, the STs at 0 to 120 cm were higher than 0°C (the shaded areas in Figs. 3a–d), and the soil at these depths were the active layers of permafrost. The STs at 120 cm to 240 cm (Figs. 3e–h) were below 0°C for the whole year, which represents the permafrost layers. The default CLM5.0 simulation (TEST1, black solid line in Fig. 3) greatly amplified the annual amplitude of ST, especially in the active layers, that is there was a significant underestimation during the freezing period, and conversely a significant overestimation during the thawing period in the active layers. This is also a common problem in the CLM4.5 and CLM5.0 simulations regarding the seasonally frozen ground and permafrost on the TP (Fang et al., 2016; Luo et al., 2017; Deng et al., 2020; Ma et al., 2023). In the deep layers, the STs at 140 to 240 cm were significantly higher than the observation in the whole year. The RMSE of ST in these eight layers was between 2.68°C and 3.99°C. The MBE increased with depth, from 0.85°C at 5 cm to 2.42°C at 240 cm while the NSE decreased with depth, from 0.74 in the topsoil to –12.43 in the deep soil. It can be inferred through the MBE and NSE values, that the simulation effect of ST gradually decreased with increasing depth, in TEST1.

Figure 3.  The observed and simulated daily soil temperature (ST, units: °C) at (a) 5 cm, (b) 20 cm, (c) 40 cm, (d) 80 cm, (e) 120 cm, (f) 140 cm, (g) 180 cm, and (h) 240 cm at the XDT site. Shaded areas in the figures denote the thawing period on observed soil temperature.

Compared with TEST1, the annual amplitude of ST in TEST2 (orange solid line in Fig. 3) was smaller and it was closer to the observed values, especially in four permafrost layers. As shown in Table 3, the RMSE (between 0.47°C and 3.70°C) and MBE (between 0.19°C and 1.10°C) of each layer were all smaller than that in TEST1 while the NSE (between 0.75 and –1.49) of each layer were all higher than that in TEST1. In the four permafrost layers, as soil depth increased, the RMSE and MBE decreased while the NSE increased, which was quite different from TEST1. The simulation of ST showed significant improvement, with the average RMSE and MBE reduced to only 1.44°C and 0.49°C, respectively, in these four permafrost layers. In contrast, these values were as high as 3.08°C and 2.33°C in TEST1. The resultant ST in the permafrost layers of TEST2 were significantly better than those in TEST1, evidenced by a reduction of RMSE and MBE by 53% and 79%, respectively, indicating that the simulated ST of permafrost is more accurate when the initial ST is below freezing temperature. This is because when the initial ST is below freezing temperature (–0.85°C), according to the default initial soil conditions (Eqs. (1) and (2)), all soil water becomes ice, which is more consistent with the ice-rich characteristics of permafrost. The results presented in Table 3 demonstrate that improving the simulation of soil hydrothermal processes in permafrost layers also has a significant impact on the accuracy of simulating soil hydrothermal characteristics in active layers. In TEST2, compared to TEST1, there was a reduction in both RMSE and MBE for ST, while NSE increased in the active layers.

The simulated results of ST in TEST3 (purple solid line in Fig. 3) and TEST4 (blue solid line in Fig. 3) were very similar. Compared with TEST1 and TEST2, more substantial improvements occurred in TEST3 and TEST4 both in active layers and permafrost layers. As shown in Table 3, TEST3 and TEST4 exhibited significantly lower RMSE and MBE values, along with higher NSE values in seven soil layers. In the four active layers, TEST3 had the smallest RMSE and MBE with an average of 2.78°C and 0.53°C which is consistent with a reduction of 21% and 61% compared to TEST1 and it had the largest NSE with an average of 0.51. In four permafrost layers, TEST4 had the smallest RMSE with an average of 1.21°C which was consistent with a reduction of 62% compared to TEST1 and it also had the largest NSE which was close to 0. In TEST4, the NSE of ST in six layers was higher than that in TEST3. It can be seen from the above analysis that modification of initial SM schemes (Eqs. (3) and (4)) will greatly impact the ST simulation. Due to the more reasonable distribution of the initial SLW and SI, the soil heat transfer was also closer to the actual state, which greatly improved the simulation of ST.

Since the initialization of ST and SM by long-term spin-up were very close (Table 2), the simulated ST in TEST1-SP (black dotted line in Fig. 3) and TEST2-SP (orange dotted line in Fig. 3) were basically consistent. Compared with TEST2, TEST3, and TEST4, the simulated results in these two long-term spin-up experiments were relatively limited. The MBE of ST was between 0.79°C and 2.02°C, and it was slightly smaller than that in TEST1 but much larger than that in TEST2, TEST3, and TEST4 in all eight layers. The NSE of each soil layer in these two long-term spin-up experiments was also smaller than that in the simulations of the three modified initial schemes.

Overall, in the active layers, the default CLM5.0 simulation greatly underestimated ST during the freezing period while overestimating it during the thawing period. In the permafrost layers, the default CLM5.0 overestimated ST for the entire year. When modified, the initial ST and SM schemes can effectively reduce the bias of ST. This improvement has proven to be significantly more profound than long-term spin-up methods. The simulated ST of permafrost was more accurate when the initial ST is sub-freezing. Modified initial SM schemes provide a reasonable distribution of the initial SLW and SI, which greatly improves the simulation of ST, especially in the four permafrost layers. On average, the RMSE was reduced from 3.29°C in TEST1 to 2.30°C in TEST2 to 2.17°C in TEST2 to 2.00°C in TEST4; MBE was reduced from 1.92°C in TEST1 to 0.70°C in TEST2 to 0.24°C in TEST3 to 0.45°C in TEST4. Meanwhile, the average NSEs in TEST2, TEST3, and TEST4 were closer to zero, especially in TEST4, within which it was 0.23, while it was only about –3 in the default and two long-term spin-up experiments. Compared to the default CLM5.0 simulation, the TEST2, TEST3, and TEST4 experiments improved the ST simulation, effectively reducing the average RMSE by 30%, 34%, and 39%. Compared to the results of TEST1, the average MBE was reduced by 64%, 88%, and 77% in the TEST2, TEST3 and TEST4 experiments, respectively, whereas it only experienced an 11% reduction in both the TEST1-SP and TEST2-SP experiments.

Figure 4 shows the spatial and temporal variation in the daily mean ST based on observations and simulations from TEST1 to TEST4. From the observations (Fig. 4a), the active layer at the XDT site began to thaw in early April and thawed completely in December, with a maximum thawing depth of approximately 130 cm. Before the active layer soil thawed completely, the surface soil entered a new freeze-thaw cycle in early November. At the depth of 130 cm to 320 cm in permafrost, the annual variation of ST was small, between –4°C and 0°C. As can be seen from Fig. 4b, TEST1 could not adequately simulate the characteristics of ST in permafrost. The results of TEST1 showed that the characteristics of seasonally frozen ground, and the ST below 300 cm were greater than 0°C for the whole year. The ST profiles in TEST1-SP and TEST2-SP were very similar to the results of TEST1 and also showed the characteristics of seasonally frozen ground (figures not shown). After modifying the initial ST as –0.85°C (TEST2), the model could simulate the characteristics of permafrost, but the simulated active layer thickness (about 280 cm) was obviously larger than the observation in 2016. After further modification of the initial SLW and SI conditions (TEST3), the model could simulate the characteristics of permafrost and also significantly improve the thickness of the active layer (about 170 cm), and the annual variation of the ST further decreases, which was more consistent with the observed value. In the TEST4 experiment, ST in permafrost was further increased compared with the TEST3 experiment and was closer to the observed value.

Figure 4.  The spatial-temporal variation in the daily soil temperature (ST, units:°C) for the (a) observation, (b) TEST1, (c) TEST2, (d) TEST3, and (e) TEST4.

Figure 5 shows the simulated and observed daily SLW in different soil layers at the XDT site. It can be seen from the observation (red solid line) that the SLW in the shallowest 20 cm corresponded well with precipitation, while it had a great relationship with the soil freeze-thaw processes in the soil below 40 cm. In the default CLM5.0 simulation (black solid line), the initial SLW between 20 cm and 120 cm was much lower than in the observation, which allowed for a substantial underestimation of SLW during the thawing period in active layers. The SLW between 5 cm and 140 cm decreased earlier than the observed values during the freezing period while increasing earlier than the observed values during the thawing period. The main reason was that the ST in the active layers in TEST1 was lower than the observed value during the freezing period, so the soil froze earlier; while the ST was higher during the thawing period, so it thawed earlier. The RMSE ranged from 0.032 m³ m⁻³ to 0.162 m³ m⁻³ while the MBE ranged from –0.098 m³ m⁻³ to 0.050 m³ m⁻³. The results of SLW in active layers were generally better than those in permafrost layers, with a smaller RMSE and greater NSE. The simulated results of SLW at 5 cm and 140 cm were better than other layers, with NSEs greater than 0.32. However, the simulated result at 240 cm was poor, with an NSE of up to –46.72. The simulations of the other five soil layers had the same explanatory power as the mean of the observations.

Figure 5.  The observed and simulated daily soil liquid water (SLW, units: m³ m⁻³) at (a) 5 cm, (b) 20 cm, (c) 40 cm, (d) 80 cm, (e) 120 cm, (f) 140 cm, (g) 180 cm, and (h) 240 cm at the XDT site. The blue bar is observed daily precipitation (units: mm).

In the active layer, the results in TEST2 (orange solid line) were generally consistent with those in TEST1, while in the permafrost layers, the results in TEST2 were better than those in TEST1, with a smaller RMSE and MBE, and a larger NSE, especially at 240 cm. Compared with TEST1, the average RMSE and MBE of SLW in four permafrost layers were reduced by 20% and 21%, respectively. The results in TEST2 also indicated that the simulated SLW of permafrost is more accurate when the initial ST is below freezing.

From Table 4, TEST3 and TEST4 provide a better fit for the observed SLW in all eight layers compared to the results of TEST1 and TEST2. The SLW in TEST3 and TEST4 was larger than that in TEST1 and TEST2 during the thawing period in the active layers, effectively reducing the dry biases of SLW. On average, in these two modified experiments, the RMSE of these four active layers was 0.0121 m³ m⁻³ and 0.0122 m³ m⁻³ respectively, while the average MBEs were –0.028 m³ m⁻³ and –0.029 m³ m⁻³, representing a reduction of nearly 10% compared to the default experiment. Meanwhile, the average NSE in the active layers of these two experiments reached 0.26, which was higher than the results in TEST1 (0.09) and TEST2 (0.06). Moreover, the first date of soil thawing in active layers occurred later than those in TEST1 and TEST2, aligning more closely with the observation. In the four permafrost layers, TEST3 and TEST4 still demonstrated a simulation performance that was better than TEST1, with average RMSEs of 0.052 m³ m⁻³ and 0.054 m³ m⁻³, respectively, consistent with a reduction of approximately 15%. Additionally, the average MBE for TEST3 and TEST4 were –0.020 m³ m⁻³ and –0.011 m³ m⁻³, indicating a decrease of around 30%.

The simulated SLW in TEST1-SP (black dotted line) and TEST2-SP (orange dotted line) were also consistent since the initialization of ST and SM using a long-term spin-up are very close (Table 4). Compared with TEST2, TEST3, and TEST4, the simulated results in these two long-term spin-up experiments were relatively limited. During the thawing period, the simulated SLW was found to be even lower than the results from TEST1, resulting in a larger RMSE and a smaller NSE as shown in Table 4.

In general, the default CLM5.0 simulation greatly dried the SLW during the thawing period at the XDT site. The SLW simulation was significantly improved and the average RMSE was effectively reduced by 13%, 21%, and 19% in the TEST2, TEST3, and TEST4 experiments compared to the default CLM5.0 simulation. However, only a 5% reduction was found in the TEST1-SP and TEST2-SP experiments. Meanwhile, the average NSEs in TEST2, TEST3, and TEST4 were also closer to zero, especially in TEST3, which was –0.07. In contrast, the NSE was only about –5.85 in the default experiment and –1 in the two long-term spin-up experiments. Overall, TEST2 significantly improved the SLW in the permafrost, while both TEST3 and TEST4 showed significant improvements in SLW throughout the whole soil column.

From the spatial-temporal variation of observed SLW (Fig. 6a), it can be seen that during summer and autumn, there is an increase in SLW within the active layers due to permafrost thawing, reaching a maximum value of 0.35 m³ m⁻³. In the vertical direction, two high-value centers were observed at approximately 20 cm and 80 cm, respectively. Below 120 cm, SLW drastically decreased due to the reduced permeability of ice-rich soils near the permafrost table that impeded SLW migration.

Figure 6.  The spatial-temporal variation in the daily soil liquid water (SLW, units: m³ m⁻³) for the (a) observation, (b) TEST1, (c) TEST2, (d) TEST3, and (e) TEST4.

In TEST1, the initial value of SLW throughout the profile was evidently underestimated, resulting in lower simulated SLW levels during the study period compared to observations across the entire profile (Fig. 6b). Upon modifying the initial ST to –0.85 °C (TEST2), a large dry bias also existed throughout the entire profile. After further modification of the initial SLW and SI schemes (TEST3 and TEST4), the simulated SLW in the active layer increased during the thawing period, thereby demonstrating that the model was capable of accurately simulating the spatial distribution characteristics of two prominent high-value centers of SLW in 2016. The long-term memory of SLW significantly impacted the ST simulation. During the thawing period, TEST3 and TEST4 had more SLW in the active layer compared to TEST1 and TEST2, resulting in a lower ST from 5 cm to 180 cm (Fig. 3).

The simulated spatial-temporal variation of daily SI is presented in Fig. 7. In TEST1 (Fig. 7a), the model exhibited limitations in accurately simulating year-round permafrost with continuous ice, indicating a constraint on its simulation capability. The results of TEST2 (Fig. 7b), TEST3 (Fig. 7c), and TEST4 (Fig. 7d) differed significantly from TEST1. In TEST2, SI was consistently present in the soil below 1.1 m throughout the entire year, while in TEST3 and TEST4, SI was found at depths below 1.3 m and 1.4 m, respectively, for the entire year duration of the study. The high center (SI ≥ 0.10 m³ m⁻³) in TEST1 was located at a depth of 1.5–2.0 m in the soil from February 2016 to July 2016, while the high center (SI ≥ 0.13 m³ m⁻³) in TEST2 was situated at a depth of deeper than 1.5 m from February 2016 to June 2016. The SI values in TEST3 and TEST4 were higher compared to those in TEST1 and TEST2. The profiles of TEST3 and TEST4 exhibited two prominent centers of SI. The first center with SI ≥ 0.10 m³ m⁻³ was located in the active layer, while the second center with SI ≥ 0.24 m³ m⁻³ was found in the permafrost layer. In TEST2, it is evident that the adjustment of the initial ST to –0.85°C resulted in the detection of multi-year ice within the permafrost soil layers at various depths. In TEST3 and TEST4, subsequent modifications to the initial SLW and SI schemes led to a further increase in SI throughout the entire soil column. SI also has a great influence on the ST simulation due to its long-term memory. Throughout the simulation period, three modified experiments demonstrated higher SI values across the entire soil column, particularly in the permafrost layers, and resulted in a reduced annual amplitude of ST fluctuations in TEST2, TEST3, and TEST4 (Figs. 3, 4). The SI profiles in TEST1-SP and TEST2-SP were similar to the results of TEST1, and they also showed the characteristic features associated with seasonally frozen ground (figures not shown).

Figure 7.  The spatial-temporal variation in the daily soil ice (SI, units: m³ m⁻³) for (a) TEST1, (b) TEST2, (c) TEST3, and (d) TEST4.

Simulated differences in topsoil ST, SLW, and SI resulted in significant variations in simulated net radiation and surface energy fluxes, as shown in Fig. 8. The average and statistical results of the simulations and observations are presented in Tables 5 and 6. Due to the lack of observations regarding sensible and latent heat fluxes, we only compared the net radiation and ground heat flux in this study. The results of TEST1-SP and TEST2-SP have been omitted in Fig. 8 as they exhibited negligible differences when compared to the results of TEST1 and TEST2.

Figure 8.  The observed and simulated daily surface energy flux (units: W m⁻²), including net radiation (Rn), sensible heat flux (H), latent heat flux (LE), and ground heat flux (G). Shaded areas in the figure denote the thawing period based on observed ST at a depth of 5 cm.

Compared to the observation, six experiments effectively simulated the annual characteristics of net radiation. The observed annual mean net radiation was 81.43 W m^–2, while it ranged from 74.02 to 75.68 W m^–2 in the six simulated results. The RMSE ranged from 21.01 to 22.74 W m^–2, while the MBE ranged from –7.42 to –5.76 W m^–2 across these experiments. Moreover, all simulations showed NSE values close to 1, ranging from 0.83 to 0.85. During autumn and winter, a clear underestimation of net radiation was found in six experiments, resulting in significantly higher STs in the active layers (Fig. 3). Among these experiments, TEST3 and TEST4 demonstrated superior simulation results for net radiation with lower RMSE and MBE values and a higher NSE.

The observed ground heat flux at a depth of 5 cm had negative values during the freezing period, with an average of –5.45 W m^–2, while it showed positive values during the thawing period, averaging 10.96 W m^–2. The offsetting of both positive and negative values throughout the year resulted in an annual mean ground heat flux of only 2.67 W m^–2 (Table 5). As shown in Fig. 8, the seasonal variation of ground heat flux was successfully captured in six experiments. However, the simulated ground heat fluxes from TSET1, TSET1-SP, and TSET2-SP were found to be underestimated compared to the observed values, with an MBE ranging from –0.96 to –0.89 W m^–2. Conversely, the simulated ground heat fluxes from three modified experiments were slightly overestimated, with an MBE ranging from 0.44 to 0.69 W m^–2. Among all experiments, TEST4 demonstrated superior simulation results for ground heat flux with the smallest MBE and largest NSE.

The simulated annual mean value of sensible heat flux ranged from 44.59 to 48.88 W m^–2, while the range for latent heat flux was between 23.31 to 28.06 W m^–2 (Table 6). TEST3 and TEST4 simulated relatively lower values for sensible heat flux and higher values for latent heat flux, whereas the remaining four experiments demonstrated relatively higher values for sensible heat flux and lower values for latent heat flux.

The modified initial soil schemes as a whole significantly improved the estimations of net radiation and ground heat flux. Comparing the result of TSET1, the average MBE of net radiation was reduced by 7%, 22%, and 21%, respectively, while it was only reduced by less than 1% in the two spin-up experiments. Furthermore, the average MBE of ground heat flux was reduced from –1.78 W m^–2 to –1.64 W m^–2.

Permafrost soil profiles possess larger hydraulic and thermal memories that necessitate greater effort to initialize in LSMs modeling. Accurate initialization of ST and SM, along with robust parameterization schemes for soil hydrothermal processes and accuracy of forcing data, are equally crucial in simulating LSMs in the permafrost region on the TP. This study emphasizes the significant impact of initial ST and SM conditions on simulated soil hydrothermal dynamics and surface energy fluxes in the permafrost region on the Tibetan Plateau. The following conclusions can be drawn:

1) The default initial schemes for ST and SM in CLM5.0 were simplistic and inaccurately represented the soil initial conditions of permafrost on the TP. This resulted in the underestimation of ST during the freezing period, while significantly overestimating ST and underestimating SLW during the thawing period at the XDT site. It also resulted in the underestimation of net radiation and soil heat fluxes. The long-term spin-up method remained challenging in its ability to generate a relatively accurate initial field of ST and SM in permafrost regions on the TP. Consequently, applying this method to obtain initial field values has only led to limited improvement in simulating hydrothermal properties of the soil, with a reduction of 11% in MBE for ST and a reduction of 5% in RMSE for SLW compared to the results obtained from default simulation.

2) The coexistence of SLW with SI, in conjunction with sub-freezing ST, is more indicative of permafrost characteristics. The modified initial soil schemes proposed in this study comprehensively considered the aforementioned characteristics of permafrost, effectively enhancing the accuracy of ST and SLW simulations by reducing RMSE and MBE, as well as improving NSE. This improvement in the statistical results greatly exceeded those achieved through the long-term spin-up method. Compared to the results of default simulation, the average RMSE of ST was reduced by 30%, 34%, and 39% in the three modified initial soil schemes experiments, respectively, while the respective average MBE of ST was reduced by 64%, 88%, and 77%. The average RMSE of SLW was also reduced by 13%, 21%, and 19% in these experiments.

3) The ST and SM profiles in the default and two long-term spin-up experiments showed the characteristics of seasonally frozen ground, while three modified initial soil schemes experiments were able to simulate the characteristics of ST and SM profiles in permafrost. In TEST3 and TEST4, the active layer displayed a shallower depth, with two high-value centers of SLW that closely resembled observation. TEST3 and TEST4 exhibited higher levels of SLW and SI compared to TEST1 and TEST2, resulting in a reduced annual amplitude of ST in these modified experiments. The simulated STs are significantly influenced by the variables SLW and SI, while conversely, ST has a profound impact on the distribution of SLW and SI.

4) The modified initial soil schemes experiments significantly improved the estimations of net radiation and ground heat flux, surpassing the efficacy of the long-term spin-up method. There was a reduction in the average MBE of net radiation by 7%, 22%, and 21%, respectively, while the average MBE of ground heat flux decreased from –1.78 W m^–2 to –1.64 W m^–2.

The modified initial ST and SM schemes proposed in this study contribute to improving the accuracy of land surface energy and water simulation in the permafrost region on the TP, thereby potentially influencing the advancement of LSMs. However, simulation experiments in this paper are based on data from a single permafrost observation site, and it is necessary to validate relevant conclusions with data from multiple sites in the permafrost region of the TP. The strong coupling effect between soil heat and water should be considered, as the inaccuracy of SLW and SI significantly contributes to the deviation of ST simulation. Additionally, the long-term memories of ST and SM are further amplified by ice-rich soils in the permafrost region. Therefore, it is imperative to enhance our comprehension of SI in the permafrost region on the TP for future research endeavors.

Acknowledgements. This work was supported jointly by the National Natural Science Foundation of China (Grant No. U20A2081), West Light Foundation of the Chinese Academy of Sciences (Grant No. xbzg-zdsys-202102), and the Second Tibetan Plateau Scientific Expedition and Research (STEP) Project (Grant No. 2019QZKK0105). We also thank the National Cryosphere Desert Data Center and Cryosphere Research Station on the Qinghai-Tibet Plateau, Chinese Academy of Sciences for providing data.

Data Availability Statement. The dataset is provided by the National Cryosphere Desert Data Center. (http://www.ncdc.ac.cn) and Cryosphere Research Station on the Qinghai-Tibet Plateau, Chinese Academy of Sciences. The input atmospheric forcing data at the Xidatan site can be assessed at https://doi.org/10.12072/ncdc.CCI.db0017.2020. (http://www.ncdc.ac.cn/portal/metadata/b274594b-e6c1-444a-8399-5767b4cef62d). Soil temperature and soil liquid water at the Xidatan site can be assessed at https://doi.org/10.12072/ncdc.CCI.db0014.2020. (http://www.ncdc.ac.cn/portal/metadata/fbd8e9c7-0ae7-4bd7-b295-7de62208be5f).