Sensitivity of soil freeze/thaw dynamics to environmental conditions at different spatial scales in the central Tibetan Plateau

Anenhancedunderstandingofenvironmentalcontrolsonsoilfreeze/thaw(F/T) dynamicsatdifferentspatialscales iscriticalforprojectingpermafrostresponsestofutureclimateconditions.Inthisstudy,a1-Dsoilthermalmodeland multi-scaleobservationnetworkswereusedtoinvestigatethesensitivityofsoilF/TdynamicsinthecentralTibetan Plateau(TP)toenvironmentalconditionsatlocal(~10km)-,medium-(~30km),andlarge(~100km)-scales.Model simulated soil temperature pro ﬁ le generally agrees well with the observations, with root-mean-square errors (RMSE) lower than 1.3 °C and 2.0 °C for two in-situ networks, respectively. Model simulated maximum frozen depths (MFD) closely related to elevation (R 2 = 0.23, p b 0.01), soil moisture content (R 2 = 0.25, p b 0.01), and soil organic carbon (SOC) content (R 2 = 0.18, p b 0.01); however, the impact of SOC on MFD may be due to the close correlation between SOC and soil moisture. The main factors affecting MFD vary with scale. Among the environmental factors examined, topography (especially elevation) is the ﬁ rst-order factor controlling the MFD at the large-scale, indicating the dominance of thermal control. Aspect shows sizeable impacts at the medium-scale, whilesoilmoistureplays animportantroleatthelocal-scale.Soilthaw onsetshowsa closecorrelationwiththeex-aminedenvironmentalfactorsincludingsoilmoisture,whilefreezeonsetseemstobein ﬂ uencedmorebyotherfac-tors.Besidesthewell-knownthermaleffect,ourstudyhighlightstheimportanceofsoilmoistureinaffectingsoilF/T dynamics at different scales in the central TP region, and reliable soil moisture products are critical to better project the response of the TP frozen ground to future warming at ﬁ ner scale.


Introduction
Frozen ground, including permafrost and seasonally frozen ground (SFG), is an important component of the cryosphere.Permafrost is defined as the lithospheric material where temperatures have remained at or below 0 °C for a period of at least two consecutive years, while SFG refers to the soil that freezes and thaws annually (Wang and French, 1994;Yuan et al., 2020).Permafrost and SFG cover approximately 24% and 56% of the exposed land surface of the Northern Hemisphere, respectively (Wu and Zhang, 2008;Zhang et al., 1999;Zhang et al., 2003), exerting an important influence on energy exchanges, hydrological processes, natural hazards and carbon budgets, with potentially large consequences on the global climate system (Riseborough et al., 2008;Walter Anthony et al., 2018;Walvoord and Kurylyk, 2016).As the result of energy and mass exchange between the ground surface and atmosphere, frozen ground is sensitive to climate change (Wu et al., 2016).The recent warming trend has resulted in significant changes in the frozen ground, and widespread soil thawing has occurred in both the high-latitude and high-altitude regions (Chen et al., 2015;Park et al., 2016;Shakhova et al., 2017;Wu and Liu, 2004;Yao et al., 2019;Zhao et al., 2004).
Long-term observation records have shown that permafrost is warming at the global scale.In the northern high latitudes, continuous-permafrost temperature has increased by 0.39 ± 0.15 °C from 2007 to 2016, while discontinuous permafrost warmed by 0.20 ± 0.10 °C during the same period (Biskaborn et al., 2019).The mountain permafrost temperature has also increased by 0.19 ± 0.05 °C during the past decade (Biskaborn et al., 2019).There is also evidence showing that MFD in the mid-latitude SFG region decreased at a rate of 1.2 cm yr −1 for 1950-2009 (Peng et al., 2017).These changes in frozen ground are closely linked with cold-region hydrological and ecological processes.Soil thawing and permafrost degradation can alter soil water storage, aggravate surface water deficit, and further impact runoff process (Liljedahl et al., 2016;Wang et al., 2009).Soil freeze/thaw (F/T) cycles are also closely linked to vegetation growth due to their influences on soil water condition and nutrient release (Kreyling et al., 2008;Zhang et al., 2016).Therefore, an improved understanding of environmental control on soil F/T dynamics is crucial for projecting the hydrological and ecological responses to future climate changes in frozenground regions.
Soil F/T dynamics are predominantly controlled by ground thermal conditions, soil moisture content, local topography and thermal buffering of underlying soils, surface vegetation and snow covers (Cao et al., 2017;Way and Lewkowicz, 2018).Previous studies have highlighted the impacts of soil moisture, snow cover and organic soils on soil F/T dynamics, especially in the Arctic region (Fisher et al., 2016;Yi et al., 2018b;Yi et al., 2015).Soil moisture content impacts the soil F/T processes and heat transfer through regulating the soil hydraulic and thermal properties (Hinkel et al., 2001;Pellet et al., 2016;Romanovsky and Osterkamp, 2000).Active layer in the Arctic continuous permafrost areas is generally shallow with higher saturation level (Jafarov et al., 2017).Organic soil generally acts as an insulation layer due to its low thermal conductivity in summer and large hydraulic conductivity, which can have a significant impact on surface energy fluxes, soil temperature and soil moisture conditions (Lawrence and Slater, 2008).Winter snow cover also plays an important role in determining the response of soil temperature to surface warming due to its strong insulation effect and its impact on surface albedo and the energy budget (Zhang, 2005).
The Tibetan Plateau (TP) has the largest areas of permafrost in the mid-and low-latitude regions of the world (Zou et al., 2017).Compared with high latitude area, the TP region is characterized by overall shallow snow cover, relatively low soil organic matter concentration, prevalent coarse soil and low soil moisture content (Cao et al., 2018;Chen et al., 2012;Li et al., 2008;Zhou et al., 2013), and permafrost in the TP is relatively warm and thin (Wang and French, 1994).With those differences, soil F/T dynamics in the TP region may respond differently to environmental conditions and recent climate trends compared with the Arctic region.Previous studies found that elevation in the TP plays an important role in permafrost distribution, and other topographic factors such as slope and aspect also have a non-negligible influence (Luo et al., 2014;Ma et al., 2015;Zhang et al., 2008).Organic soil layer, often accompanied by high soil water content in summer and ice content in winter, is also a key factor controlling soil thermal state through affecting the soil thermal conductivity as well as latent heat (Cao et al., 2018;Zhou et al., 2013).However, most of the existing studies focused on permafrost region and during the thaw season, while environmental control on soil F/T dynamic in SFG area was relatively less investigated.In addition, due to the strong heterogeneity of local terrain and large temporal and spatial dynamic range of soil moisture in the TP (Gruber et al., 2017;Yang, 2013), soil F/T dynamics may respond differently to environmental conditions at different spatial scales.Indeed, simulation results at the coarse scale may not represent significantly variability of environmental control on soil F/T dynamics (Westermann et al., 2015).
Therefore, in this study, we aimed to determine the underlying environmental controls including elevation, slope, aspect, soil moisture, soil organic carbon content and enhanced vegetation index (EVI) on soil F/T dynamics at different spatial scales in the central TP region.Two in-situ soil temperature and moisture observation networks were used in this study.Considering the lack of deep-soil observation, we also used a numerical model (Rawlins et al., 2013;Yi et al., 2018b;Yi et al., 2015) to extrapolate the soil temperature profile.Then an integrated analysis combining both the in-situ observations and model simulations was conducted to investigate the sensitivity of soil F/T dynamics to the environmental conditions at different spatial scales.This study offers an enhanced understanding of the soil F/T sensitivity over the TP, which is beneficial to future model parameterization.

Study area
The study area is located in the central TP (Fig. 1,, along the Qinghai-Tibet Railway.It is characterized by relatively high elevation variation (4000-5500 m), and a heterogeneous subsurface condition.Influenced by the East Asian and Indian monsoon, the climate in this region is cold and semi-arid (Yao et al., 2012).During the period of 1980-2015, the annual mean air temperature in the study area ranges from −4.6 to −0.5 °C, while annual precipitation varies from 306 mm in the north to 452 mm in the south.Most of the precipitation falls during the monsoon season from May to September.The study region is mainly underlain by permafrost in the north and SFG in the south based on permafrost maps in this region (Fig. S1), with only some patchy and shallow snow cover in winter, usually b10 cm, and lasting only for a few days (Dai et al., 2018).This region is dominated by alpine meadows and grasslands with sparse vegetation cover, and soil organic content is relatively low and mainly accumulates in the top layers, when there is no peat layer (Jin et al., 2008;Wu et al., 2016).The World Reference Base for Soil Resources (WRB) soil classification system defines the dominant soil types in this region as gleysols and chernozems in the alpine meadows and arenosols in the alpine steppe (Niu et al., 2019).

Datasets
Two multi-scale in-situ soil temperature and moisture observational networks were used in this study (Table 1), i.e., the CEOP (Coordinated Enhanced Observing Period) Asia-Australia Monsoon Project sites within the central TP domain (CAMP-Tibet, 1997-2007) (Ma et al., 2003), and the Central Tibetan Plateau Soil Moisture and Temperature Monitoring Network (CTP-SMTMN, 2011-2014) (Yang, 2013).The CAMP-Tibet consists of 9 sites with soil temperature and moisture measured simultaneously.These sites are located along the Qinghai-Tibet Railway, spanning from permafrost area in the north to SFG area in the south (Fig. 1b).Soil temperature and moisture were measured at 9 depths from 4 to 250 cm during the period 1997 to 2007 with some data gaps, and only the periods with both soil temperature and moisture observations available were used in this study.
The CTP-SMTMN was established in a region spanning ~1.0°× 1.0°( ~100 × 100 km) around Naqu, with an area of ~10,000 km 2 (Fig. 1c).Both soil temperature and moisture were measured at three spatial scales ranging from 0.1°(~10 km) to 1.0°(~100 km) at 4 depths (5, 10, 20 and 40 cm).A subset of 41 sites from the total 56 CTP-SMTMN sites with longer observations (N2 years), lower gravel content and better data quality was used in this study.During the study period from 2011 to 2014, there were small annual temperature changes ranging from −0.2 to 0.0 °C.One of the objectives in designing the multi-scale network is to match different scales of land hydrological modeling (Yang, 2013).The sites in the large-scale (1.0°× 1.0°) network are deployed along a cross transect, with large altitudinal variation (4470 to 4835 m, Fig. 1c).The sites in the medium-scale (0.3°× 0.3°, ~30 × 30 km) network are distributed uniformly in a small rolling-hill area and nested in the large-scale network, with elevation ranging from 4539 to 4818 m (blue rectangle in Fig. 1c).The sites in the smallscale (0.1°× 0.1°) network are further nested in the medium grid and located in a flat area, and elevation ranges from 4647 to 4737 m (red rectangle in Fig. 1c).Based on the uniform spatial distribution, we used 7, 16, and 29 sites at the small-, medium-, and large-scales respectively with some sites overlapped at different scales.Surface layer (0-5 cm) soil organic carbon (SOC, in mass, kg kg −1 ) content was measured at each site using the Total Organic Carbon Analyzers (SHIMADZU-TOC-VCPH) after soil samples (b2 mm) were dried (at 95 °C for 12 h) and grinded into powder (Chen et al., 2012).
The soil properties of each site, such as soil porosity and soil texture (sand and clay fraction) were extracted from the Soil Database of China for Land Surface Modeling (http://globalchange.bnu.edu.cn/research/soil2).The physical and chemical soil attributes from this database were developed based on 8979 soil profiles and the Soil Map of China (1:1,000,000), using polygon linkage method to characterize the spatial distribution of soil properties.The spatial resolution is 30 arc-seconds (~1 km), and six layers were used to characterize the soil vertical variations down to 83 cm below surface.A more detailed description of this database can be found in Wei et al. (2013).The ASTER global Digital Elevation Model (DEM) with 30 m resolution was also used in this study to calculate the slope and aspect at each site through spatial analysis.Instead of using the original value of aspect from 0-360°, we used the cosine values of the aspect for the statistical analysis.In addition, the MODIS growing-season EVI (averaged from June to August) during the period 2011-2014 was used to represent the vegetation density.The EVI data are available at every 16 days and 250-m spatial resolution (MOD13Q1) (Huete et al., 2002).

Model description
A detailed soil process model (Rawlins et al., 2013;Yi et al., 2018b;Yi et al., 2015), was used to simulate the soil temperature profile in this study.The soil model can simulate the ground thermal dynamics by solving a 1-D heat transfer equation with phase change: where T(z, t) is the temperature (°C), C and λ are the soil volumetric heat capacity (J m −3 K −1 ) and thermal conductivity (W m −1 K −1 ) respectively; L is the volumetric latent heat of fusion of water (J m −3 ); ζ is the volumetric water content, and θ is the unfrozen liquid water fraction (range: 0-1) which is estimated empirically as: where the constant T * is used to represent the freezing point depression, and b is a dimensionless parameter determined by fitting the unfrozen water curve (Romanovsky and Osterkamp, 2000).A significant amount of liquid water can exist even when the soil temperature is considerably lower than T * , which is characterized by different values of b.The upper boundary can be set as at the soil surface or the snow surface.When snow is present, the snowpack is divided into multiple layers, with snow thermal properties calculated based on snow density (Yi et al., 2015).
The parameterization of soil thermal properties employed in the previous model underestimated summer soil temperature and overestimated winter soil temperature at the study sites, while the Johansen's method performs better compared to previous studies (Chen et al., 2012;Farouki, 1981).Therefore, the thermal conduction calculation was modified in this study using the Johansen's method (Johansen, 1977), which are widely used in soil thermal transfer modeling in global land models (Lawrence and Slater, 2008;Tao et al., 2017).Soil thermal conductivity was calculated as: where λ refers to the thermal conductivity, and the subscripts sat and dry stand for the saturated soil and dry soil respectively, while K e is the Kersten number (K e ≥ 0) related to soil saturation degree S r : For dry mineral soils, thermal conductivity mainly depends on the soil bulk density (γ d ): where the bulk density γ d = 2700(1 − n) is in kg/m 3 , and n is the soil porosity.For unfrozen saturated soil, thermal conductivity can be calculated as: while for frozen saturated soil, thermal conductivity is: where the subscripts i and w stand for the ice and water respectively.The λ solid is the thermal conductivity for soil solids, which can be calculated based on the soil texture: where the sand and clay fraction extracted from the Soil Database of China for Land Surface Modeling (Wei et al., 2013) were averaged across the six layers (0-83 cm) to calculate the soil properties for the entire soil column.

Model setup and parameterization
In previous studies (Rawlins et al., 2013;Yi et al., 2015), the model was coupled with a water balance module that simulates soil moisture content regulating the soil thermal properties.However, in the current study, we used in-situ soil moisture measurements to parameterize the model soil thermal properties instead, due to a lack of in-situ precipitation datasets as model inputs and large uncertainties in the soil moisture modeling in both Arctic and the TP region (Bi et al., 2016;Swenson et al., 2012).Because the ice content cannot be detected by soil moisture sensors in winter, the soil saturation during frozen period was calculated as the soil water content prior to the soil freeze onset at each year, while the unfrozen water content was determined using Eq. ( 2).Soil column depth was set to 50 m below surface and divided to 30 layers, with soil layer thickness increasing with depth (Table S1).The Dirichlet boundary conditions were used as an upper boundary (z s ), i.e., soil temperatures at 4 or 5 cm depth at CAMP-Tibet and CTP-SMTMN sites, respectively, and no heat flux was assumed at the lower boundary (z b = 50 m, geothermal heat flux = 0).In order to achieve equilibrium steady state for the top 10-m soil temperatures, the model was first run using daily averaged MODIS land surface temperature (LST) at the year before simulation period with a spin-up for 60 years, then followed by a transient run driven by in-situ observation.
The soil porosity data and soil moisture measurements were available at different depths from the model defined soil layers.Therefore, we linearly interpolated the soil porosity and daily in-situ soil moisture along the soil profile.Soil moisture content below observation layers and above the bedrock layer was defined as the same as the deepest observation layer.The soil texture below 83 cm and above bedrock was also defined as the same as 83-cm layer.Because near-surface soil temperature was used as the upper boundary, and considering the shallow snow depth, short snow duration, and low organic content (usually a few percent in mass concentration) at these sites, we did not include the parameterizations of snow and soil organic matter in the current model.
The observed soil temperatures from the two CAMP-Tibet sites (D105 in permafrost region and BJ in SFG region) were first used to calibrate the model, including the parameters related to volumetric heat capacities and soil thermal conductivities, which were listed in Table 2.The soil thermal properties were primarily determined by soil texture and soil saturation conditions.Comparing with the in-situ data near the study site from other studies, we found that soil porosity from the Soil Database of China for Land Surface Modeling was slightly overestimated.Therefore, a coefficient (α) was used to adjust the soil porosity data based on the in-situ measurements from He et al. (2017).Other parameters such as freezing point depression were determined and adjusted based on the previous studies (De Lannoy et al., 2014;Yi et al., 2018b).Then the same soil parameterization was applied to the other 7 CAMP-Tibet sites and also the 41 CTP-SMTMN sites, and the model simulated soil temperatures were validated using in-situ data.Since soil moisture observations at the CTP-SMTMN network are only available at the upper 40 cm soils, soil moisture profile at CTP-SMTMN site was extrapolated based on the soil moisture profile at the CAMP sites adjacent to the CTP-SMTMN sites.

Soil F/T dynamics analysis
We performed an integrated analysis using the model simulations and the multi-scale in-situ observations to investigate the sensitivity of soil F/T conditions to environmental variables across the central TP domain.Several F/T indices based on observed and simulated soil temperature were selected to describe the F/T conditions, including freeze onset, thaw onset, ALT and MFD.Annual soil freeze (thaw) onset was defined as the first day when the soil temperature of a 5-day moving window dropped below (or raised above) freezing point.MFD for the SFG was calculated by linearly interpolating the time series of the soil temperature profile and identifying the deepest soil layer crossing the freezing point threshold throughout the year, while ALT was defined for the permafrost area by identifying the deepest layer rising above freezing point.Most of the analysis was based on the simulations and observations at CTP-SMTMN sites located in a concentrated area underlain by SFG.Linear regression and Pearson correlation analysis were used to investigate the relationship between the soil F/T indices and different environmental variables including elevation level, SOC content, vegetation density represented by growing-season EVI, and soil moisture content averaged from July to September that can well reflect the soil-water condition at study sites.We further analyzed the regional variability of MFDs and environmental conditions at the three scales defined by the CTP-SMTMN network and in a relatively flat region (mean elevation 4525 ± 37 m, green rectangle in Fig. 1c).Considering the strong surface heterogeneity at the medium-and large-scale, we also used the partial correlation analysis separately controlled by topographic parameters (i.e., slope, aspect and elevation) to understand the impact of above environmental variables on soil F/T dynamics at different scales.The statistical relationships are considered to be significant with p level of 0.05.

Model validation
After calibrating the soil parameters at D105 and BJ sites (Section 2.4), we evaluated the model performance at the other CAMP-Tibet and CTP-SMTMN sites.The model simulated soil temperatures from 20 to 250 cm are generally in close agreements with the in-situ observations at most of the CAMP sites, although there are slight overestimations at D110 site and underestimations at MS3637 site, especially at deeper layers (Fig. 2).The root-mean-square errors (RMSE) between the observed and simulated soil temperatures at different layers are overall lower than 1.3 °C, and the mean biases range from −0.5 to 0.5 °C (Fig. 2e and f).The averaged ALT and MFD simulated during the study period also well agree with the observations (Table 3), although the model overestimates the MFD at MS3637 site, likely due to the poor quality of the soil moisture data for model parameterization at this particular site.
The model also shows overall good performance at the CTP-SMTMN sites (Fig. 3a and b).The RMSEs between simulated and observed soil temperatures from 10 to 40 cm generally increase with soil depth, but are overall lower than 2.0 °C.Soil temperature at most sites are slightly underestimated, with the mean biases ranging from −1.2 to 0.8 °C.The simulated F/T onsets among the CTP-SMTMN sites are generally consistent with the observations.The correlation coefficients range from 0.63 to 0.79 (p b 0.01), although the model shows a slight negative bias (i.e.advancing) in the estimated freeze onset (Table 4 and Fig. S2).The simulated mean freeze onsets at soil depths of 20 cm and 40 cm across all the CTP-SMTMN sites are DOY 313 ± 6, and DOY 320 ± 7 respectively, compared with mean freeze onsets of DOY 317 ± 6 and DOY 326 ± 9 derived from the in-situ temperature data.
In summary, even though the model slightly underestimates the soil temperature, the validation results from both the CAMP-Tibet and CTP-SMTMN sites indicate the applicability of this model in the central TP region.

Analysis across all CTP-SMTMN sites
Simulated MFDs show a large variability among the CTP-SMTMN sites, ranging from 1.47 ± 0.10 m at BC07 site to 6.32 ± 0.87 m at MS3488 site (Fig. 3c).This might be explained by the strong variability of environmental factors in this region as shown in Fig. 4(a, b, c).The CTP-SMTMN sites are located in a relatively small area (~1.0°× 1.0°), while the elevation among these sites ranges from 4470 to 4835 m a.s.l. (Fig. 1c).The linear regression results show that MFDs have significant relationship with elevation, soil moisture and SOC content (Fig. 4a, b, c).Higher elevation sites are generally associated with deeper MFDs (R 2 = 0.23, p b 0.01, Fig. 4a), and soil generally freezes deeper during winter at sites with higher soil moisture content (R 2 = 0.25, p b 0.01, Fig. 4c).The statistical analyses indicate the positive effect of elevation and soil moisture on MFDs.A positive relationship between the MFDs and surfacelayer (0-5 cm) SOC content (R 2 = 0.18, p b 0.01, Fig. 4b) is also observed.However, there is no significant (p N 0.1) relationship between growing-season EVI and simulated MFDs.
The results also indicate significant relationships between the examined environmental factors and soil thaw onset, while soil freeze onset is less explained by these factors (Fig. 4d and e).The significant correlations between elevation and upper layer F/T onsets indicate that soil generally thaws later but freezes earlier at high-elevation sites (Fig. 4d and e, 4-20 cm), but the influence of elevation on soil freeze onset generally diminishes with increasing soil depth (non-significant correlations at layers below 40 cm, Fig. 4d).The soil moisture content shows a positive correlation with soil thaw onset, especially for the deeper layers (below 40 cm, R = 0.50-0.67,p b 0.01), indicating its overall delaying effect on thaw onset.There is no significant correlation between soil moisture and freeze onset, except for the negative correlations from 80 to 130 cm depth (R ≈ −0.31, p b 0.05), implicating an advancing effect of soil moisture.Surface-layer SOC content and vegetation density (indicated by growing-season EVI) also positively correlated with soil thaw onset, indicating their postpone effect (Fig. 4d  and e).These influences are gradually weakened with the increasing soil depth, especially for the effects of vegetation density.

Analysis for different scales
The CTP-SMTMN soil temperature and moisture were measured at three spatial scales, which can be used to investigate the sensitivity of soil F/T dynamics to environment conditions at different scales.Results indicate that MFDs show a large variability even at the small-scale (~0.1°× 0.1°) where the changes in elevation are negligible (mean elevation 4694 ± 27 m, Fig. 5a).Soil moisture content (from July to September) at this scale varies from 0.22 to 0.30 (m 3 m −3 ), which can partially explain the changes in the MFDs (Fig. 5b, R 2 = 0.36).Within a relatively large area with small elevation variations (mean elevation 4525 ± 37 m), soil moisture also exhibits a strong control on MFDs (Fig. 5c, R 2 = 0.53).
At the medium-scale (~0.3°× 0.3°, Fig. 6a), the influence of the environmental factors on MFDs is more complicated than that at small-scale, due to stronger topographic heterogeneity.Fig. 6(b, c, d, e) shows the variation of MFDs and different environmental factors at the medium scale.The elevation variation is moderate, ranging from 4539 to 4818 m (mean: 4679 ± 81 m, Fig. 6b) and the slopes range from ~4 to 23°(Fig.6d).Statistical analysis shows the strong impact of aspect at the medium-scale (R = 0.51, p b 0.05), and this impact is also distinct after excluding the other topographic influences (Fig. 6f).The ground surface temperature at south-slope sites is usually higher than others, and therefore shows a shallower MFD, e.g., at the P7 site (MFD = 2.80 ± 0.60 m, Fig. 6b and e).Meanwhile, slope may also have some complex impacts on MFDs at the medium-scale.The influence of elevation on MFDs becomes significant when excluding the slope influence (R partial = 0.52, p b 0.05, Fig. 6f), and the correlation between elevation and MFDs is also modulated by aspect.Soil moisture has some impacts on MFD, but this correlation is affected by slope (R partial = 0.49, Fig. 6f).
At the large-scale with larger elevation range from 4470 to 4835 m (mean: 4630 ± 109 m, Fig. 7a), soil freeze dynamic also shows complicated responses to the environmental conditions (Fig. 7b, c, d, e).The variation of MFDs among these sites is generally consistent with the variation of elevation, indicating the major influence of elevation on MFDs at the large-scale (R = 0.59, p b 0.01, Fig. 7b, e and f).Soil moisture also exerts significant effects on MFDs (R = 0.54, p b 0.01, Fig. 7d, e and f).Although the vegetation cover, slope and aspect have some impacts on MFD at individual sites (Fig. 7b, c, d and e), the influences of elevation and soil moisture are still dominant (Fig. 7f).At BC02 site located in the wetland, the deep MFD may be related to the high SOC at this site (~20% in mass concentration, Fig. 7c).However, similar variations in the soil moisture and SOC among these sites imply that SOC affects the soil freeze process mainly through altering the soil moisture content at the upper soil layers (Fig. 7c and d).

Environmental sensitivity of soil F/T dynamics at different scales
We examined different soil F/T indices and their relations to several environmental factors using a multi-scale in-situ soil temperature network and model simulations.Our analysis indicates the important control of topography on soil F/T conditions across the central Tibetan Plateau, and the impacts of elevation become more significant at larger-scales (Figs. 6 and 7).Elevation and topography strongly affect regional climate through the orographic effects on precipitation, temperature lapse rates and solar radiation loading (Gruber et al., 2017).A recent study showed that the average lapse rate for this region is around 4.2 °C per km (Wang et al., 2018).Topography can affect soil F/T dynamics by altering the land surface temperature and impacting the surface energy balance.For example, our analysis also shows that aspect has obvious impacts on soil freezing especially at the medium-scale sites (Fig. 6f), and slope further influences the correlation between MFDs and factors.Higher-elevation sites usually experience later thaw onset, deeper MFD, and earlier freeze onset, especially at the upper soil layers (Fig. 4a, d and e).However, our results did not show a significant relationship between elevation and deeper-layer (below 40 cm) soil freeze onset, maybe because the soil freeze onset at deeper layer was influenced more by other factors such as heat transfer process.
Our results also highlighted the important effects of soil moisture on soil F/T dynamics, especially at the small-scale.Soil moisture affects the soil F/T processes by changing soil physical properties such as thermal conductivity, heat capacity and hydraulic conductivity, and controlling the latent heat release during the F/T periods (Boike et al., 2008;Hinkel et al., 2001;Pellet et al., 2016).The large dynamic range of soil  (a-d) are the comparisons between simulated and observed soil temperatures at different depths at one of the CAMP-Tibet sites (Amdo).Panel (e) and (f) refer to the mean bias and root mean square errors (RMSE) between simulated and in-situ soil temperatures from 20 to 250 cm at CAMP-Tibet sites.D105 and BJ sites were used to calibrate the model, while the other 7 sites were used for model validation.

Table 3
Comparisons between the simulated and in-situ ALT or MFD at the CAMP-Tibet sites.The in-situ ALT or MFD was determined using soil temperature observations.The ALT or MFD (Mean ± Standard Deviation) are averaged during the observation period of each site, which was listed in Table 1.moisture in the study area might partly explain the strong variation in MFDs at the local scale (Fig. 5b and c).Even at the large-scale where the topography showed a dominant control on MFDs, the influence of soil moisture cannot be ignored (Figs. 6 and 7).Soil moisture also exerts a strong influence on soil thaw onset.Higher soil moisture generally resulted in later soil thaw onset in the study area, and this effect increased with soil depth (Fig. 4e).Higher soil moisture content is generally associated with longer zero-curtain (when soil temperature stays around 0 °C) during the F/T periods, due to large latent heat exchange with soil water phase change (Outcalt et al., 1990).Therefore, the longer zero-curtains observed at deeper layers may partly explain the increasing influence of soil moisture on soil F/T onsets.In addition, the spring zero-curtain period is more distinct than autumn zero-curtain period in SFG areas (Jiang et al., 2018), which may explain why soil moisture exerts a stronger control on soil thaw onset than on soil freeze onset at the CTP-SMTMN sites.SOC content in the TP region is relatively low and mainly contained in the topsoil layer (Chen et al., 2012;Wang et al., 2002).Our analysis still shows a large influence of surface (0-5 cm) SOC on the soil F/T dynamics.Compared with mineral soil of very low SOC content, soil with high SOC generally has a high-water retention capacity, and is generally associated with high soil moisture content (Cao et al., 2018).Organicrich soil generally has a lower thermal conductivity under thawed condition, but shows larger thermal conductivity when saturated at frozen conditions (Cao et al., 2018;Shiklomanov and Nelson, 1999), and therefore might be associated with shallower ALTs or larger MFDs.Our results show that soil with higher organic carbon content generally froze deeper (Fig. 4b).The results also show a significant delaying impact of SOC content on soil thaw onset, although this impact decreased with increasing soil depth (Fig. 4e).The surface SOC content at most CTP-SMTMN sites is lower than 10% (except for BC02, located in the wetland); therefore, the influence of SOC on soil F/T dynamics may be mainly attributed to its impact on soil moisture content rather than its insulation effect.Soil moisture sensor at CTP-SMTMN sites was calibrated to account for the impact of SOC content by the comparison between gravimetric method measurements and sensor measurements (Yang, 2013), and this may explain similar effects of soil moisture and SOC content on soil F/T onsets in our study.

Site name
There is no significant (p N 0.1) correlation between growing-season EVI and frozen soil indices (i.e., MFD and freeze onset), while EVI is closely correlated with thaw onset especially at upper soil layers (Fig. 4e, ~10-40 cm), which likely indicates the negligible effects of surface vegetation density on ground thermal conditions during the frozen period.Surface vegetation cover and soil moisture are closely linked to soil F/T dynamics through their control on surface energy balance, soil hydraulic and thermal properties (Li et al., 2015).Thus, areas with higher soil moisture were generally linked with higher vegetation density, and showed an overall later soil thaw onset.
In summary, our results indicate that environmental conditions had significant influences on soil F/T dynamics in the central TP region, with topographic factors being the first-order factor controlling the MFD at the large-and medium-scale, while soil moisture plays an important role regulating MFD at the local scale.Soil thaw onset shows a closer association with examined environmental factors than soil freeze onset, and the major influence factors on soil freeze onset still need further investigation.Moreover, since F/T dynamics may be different between

Table 4
Comparisons between the simulated and observed F/T onsets (Mean ± Standard Deviation) across the CTP-SMTMN sites (n = 41) at different soil depths.The mean F/T onsets shown here were averaged over all the sites from 2011 to 2014, and the correlation coefficient (R) was calculated between the in-situ observation and model simulation of mean (2011)(2012)(2013)(2014)  SFG and permafrost regions (Jiang et al., 2018), further studies for SFG and permafrost regions are called to fully understand the environmental controls on soil F/T dynamics.

Model limitations and uncertainties
There are some important uncertainties and limitations in our modeling analysis due to the model physics, inputs and parameterization.Although we investigated the impacts of several environmental factors, there are other environmental factors that may have strong impacts on soil F/T dynamics, but were not included in our analysis due to a lack of observation, such as soil texture (Yi et al., 2018a).Soil texture and SOC play important roles in determining how soil responds to surface warming (Lawrence and Slater, 2008;Yi et al., 2018b).Our analysis demonstrates that higher SOC content is normally associated with deeper MFD.The effect of surface SOC on subsurface soil temperature has been considered since the upper boundary condition was set as near-surface temperature (~5 cm), but the influence of deeper-layer SOC is unclear.Although the SOC is mainly stored in the surface layer and gradually decreases with soil depth in the central TP (Chen et al., 2012;Yang, 2013), ignoring the effects of organic content in the modeling experiment may result in biases in the model simulations.On the other hand, gravel soil has a larger thermal conductivity and can result in a rapid response of soil temperature to surface temperature changes   (Pan et al., 2017), and it should be considered in the model across the TP region, where the coarse soil is prevalent due to the slow soil-forming processes and strong erosion (Arocena et al., 2012).A few studies found that considering the effect of SOC content and the distinct physical properties of coarse-fragment soils can more accurately characterize soil thermal conductivity for the alpine grasslands in the TP (Pan et al., 2017;Yi et al., 2018a;Zhou et al., 2013).However, accurately representing the effects of gravel and SOC content in the model parameterization is difficult due to the lack of data.Recent studies have provided estimates of SOC content and its vertical distribution (up to 25 m) in the TP area (though limited to the permafrost area) through synthesizing a substantial number of soil pedons (Ding et al., 2016), which has a potential to improve the model parameterization in future studies.Moreover, although we focus on a few meter depth frozen ground, assuming no heat flux at lower boundary layer in the model may introduce additional uncertainty, as previous study highlights the significant impacts of geothermal heat flux on permafrost temperature (Wu et al., 2010).Seasonal snow cover is also important in determining how soil responds to surface warming, because of its strong insulation effects on ground temperature and the influences on surface albedo (Slater et al., 2016;Yi et al., 2015;Zhang, 2005).Previous studies suggested that there is little snow accumulation and only some patchy snow cover distributed in our study area (Dai et al., 2017;Dai et al., 2018;Jiang et al., 2015Jiang et al., -2016)).The observations at Naqu station also show that snow cover is usually lower than 10 cm and only present during a short period (b10 days).Actually, there is no typical snow accumulation and ablation season in the interior of the Tibetan Plateau, which is characterized by sporadic patchy and shallow snow cover during winter (Wang et al., 2020).Taking near-surface soil temperature (~5 cm) as the model's upper boundary can reflect the effect of snow on sub-soil temperature to some extent.Partly because of this, we did not include the parameterizations of snow in the current simulation.Snow cover may have substantial effects on soil F/T processes in other parts of the TP region with more persistent snow cover (Zhou et al., 2013); snow cover effects on soil thermal dynamics may become even more important with projected increases in winter precipitation in the TP area (Kang et al., 2010).However, accurately representing the snow cover distribution in the model is difficult, due to the paucity of in-situ observations and extreme elevation gradients which limited our understanding of local phenomena of snowfall (Kapnick et al., 2014).Potential improvements may come from merging in-situ and gridded meteorological datasets with satellite-based snow products, such as optical-infrared and microwave sensor derived snow cover extent and snow depth observations (Painter et al., 2016).
Finally, uncertainties in the in-situ data will also influence our model simulations.The model was driven by ground-surface temperature (~5 cm) rather than land surface skin temperature due to the data limitation.However, the land surface temperature is a critical variable to address the impact of surface energy exchange process on soil F/T especially at a regional scale.On the other hand, although we excluded the sites with poor soil moisture data, the extrapolated soil moisture content especially at the deeper soils (N40 cm) due to data limitation cannot fully reflect the actual soil moisture profile at CTP-SMTMN sites.This may result in relatively large temperature bias in deeper soils, e.g., the RMSE at 40-cm layer larger than that at the upper layer.Therefore, better soil moisture parameterization scheme is needed in future studies.

Conclusion
In this study, we used a multi-scale observational network for soil temperature and moisture and a physical model to investigate the sensitivity of soil F/T dynamics to environmental conditions in the central TP region.The validation results show that the model simulated soil temperatures generally agreed well with the observations, with RMSE lower than 1.3 °C and 2.0 °C for the CAMP-Tibet and CTP-SMTMN sites, respectively.
The analysis on the relationships between the soil F/T dynamics and environmental conditions shows that MFDs were significantly influenced by elevation (R 2 = 0.23, p b 0.01), soil moisture content (R 2 = 0.25, p b 0.01), and SOC content (R 2 = 0.18, p b 0.01); however, the impact of SOC on MFD may be due to a close correlation between SOC and soil moisture in the study area.The main factors affecting MFDs vary with scale.Topography is the first-order factor controlling the MFD at the large-scale, especially the elevation which exerts a strong impact on the soil freezing at that scale.This is not surprising as elevation has a close relationship with temperature.Among the sites at the medium-scale, aspect shows sizable impacts on MFD, and slope further influences the correlation between MFDs and soil moisture.At the local scale, soil moisture plays an important role in regulating MFD.Soil thaw onset is more closely associated with examined environmental factors than soil freeze onset at the SFG sites.However, more analyses are still needed to evaluate the potential environmental control on soil freeze onset.
Besides the well-known thermal effect, our study highlights the importance of soil moisture in affecting soil F/T dynamics at different scales in the central TP region, which should be considered when projecting the responses of TP frozen ground to future warming.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The elevation (a) and site location (b, c) in the study area.Panel (b) shows the permafrost distribution based on the IPA permafrost map, and the location of the two in-situ networks (CAMP-Tibet and CTP-SMTMN), while panel (c) shows the location of CTP-SMTMN sites (n = 41), corresponding to the black rectangle shown in panel (b).The "high" and "low" in the legend refer to the ice content of permafrost in panel (b).The rectangles in panel (c) refer to in-situ networks at small-and middle-scales as well as a flat region used for the scaling analysis.(For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Model calibration and validation results at the CAMP-Tibet sites.Panels (a-d) are the comparisons between simulated and observed soil temperatures at different depths at one of the CAMP-Tibet sites (Amdo).Panel (e) and (f) refer to the mean bias and root mean square errors (RMSE) between simulated and in-situ soil temperatures from 20 to 250 cm at CAMP-Tibet sites.D105 and BJ sites were used to calibrate the model, while the other 7 sites were used for model validation.

Fig. 3 .
Fig. 3. Evaluation of model performance (a and b) and the simulated mean MFDs (Mean ± Standard Deviation) from 2011 to 2014 (c) at the CTP-SMTMN sites.Panels (a) and (b) refer to the mean bias and RMSE values between the simulated and observed soil temperatures from 10 to 40 cm.

Fig. 4 .
Fig. 4. Relationships between model simulated soil F/T dynamics and environmental factors for all the CTP-SMTMN sites (n = 41).Panels (a-c) show the linear regression between MFDs and elevation, surface SOC content (0-5 cm), upper-layer (b40 cm) soil moisture content (cm 3 cm −3 , averaged during the period from July to September), respectively.Panel (d) and (e) are the correlation coefficients between environmental factors and soil F/T onsets at different soil depths.Soil moisture content (0-40 cm) in panel (d) was averaged from July to September (JAS, prior to freezing period), while soil moisture content in panel (e) was averaged from March to May (MAM, prior to thawing period).The gray bars in panel (d and e) denote the significance level of the correlation (from p b 0.05 to p b 0.005).

Fig. 5 .
Fig. 5. Relationship between MFDs and upper-layer (b40 cm) soil moisture content (cm 3 cm −3 , from July to September) at the small-scale (b, 0.1°× 0.1°, mean elevation 4694 ± 27 m) and within a relatively flat region (c, mean elevation 4525 ± 37 m) in the CTP-SMTMN network, represented by the red and green rectangles in panel (a) respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Variations in environmental conditions and model simulated MFDs among the sites defined by the CTP-SMTMN network aiming at resolving the medium-scale (0.3°× 0.3°).There are 16 sites in total, indicated by the blue rectangle in panel (a).Sites shown in panel (b-e) were arranged with increasing elevation, and the error bars in panel d and e indicate the standard deviation.Correlations coefficients between MFD and environmental variables in panel (f) were derived from Pearson-correlation analysis and partial correlation analysis (controlled by topographic parameter separately), respectively.The significance level is indicated by * (p b 0.05).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1
Characteristics of the CAMP-Tibet sites and a subset of the CTP-SMTMN sites (n = 41) used in this study.Only the sites with both soil temperature and moisture observations available were used in this study.Unknown soil types are indicated by "-".

Table 2
The calibrated soil parameters in the soil thermal model.A "-" indicates unitless.
soil F/T onsets among all the sites.
pths.The F/T onsets shown here were averaged from 2011 to 2014.