The characteristics of gravelly soil physical properties and 1 their effects on permafrost dynamics : A case study on the 2 central Qinghai-Tibetan Plateau

Soils on the Qinghai-Tibetan Plateau (QTP) have distinct physical properties from 27 agricultural soils due to weak weathering and strong erosion. These properties might affect 28 permafrost dynamics. However, few studies have investigated both quantitatively. In this 29 study, we selected a permafrost site on the central region of the QTP and excavated soil 30 samples from 20 cm to 200 cm. We measured soil porosity, thermal conductivity, saturated 31 hydraulic conductivity and matric potential in the laboratory. Finally, we ran a simulation 32 model replacing default sand or silty clay parameters with different combinations of these 33 measured parameters. Results showed that gravel content (diameter >2 mm) was ~55% on 34 average in soil profile; soil porosity was less than 0.3; saturated hydraulic conductivity ranged 35 from 0.004-0.03 mm s -1 ; saturated matric potential ranged from -14 to -604 mm. When 36 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-11 Manuscript under review for journal The Cryosphere Discussion started: 12 February 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Permafrost covers 25% of the earth surface.Degradation of permafrost has been reported extensively in Alaska, Siberia and the Qinghai-Tibetan Plateau (QTP; Boike et al., 2013;Jorgenson et al., 2006;Wu and Zhang, 2010).It has global impacts by releasing large quantities of soil carbon previously preserved in a frozen state and enhancing concentrations of atmospheric greenhouse gases, which will promote further atmospheric warming and degradation of permafrost (Anisimov, 2007;McGuire et al., 2009).Permafrost dynamics also have local to regional impacts on ecosystems by altering soil thermal and hydrological regimes (Salmon et al., 2015;Wang et al., 2008;Wright et al., 2009;Ye et al., 2009;Yi et al., 2014a).In addition, degradation of permafrost affects infrastructure, e.g.QTP railways and roads (Wu et al., 2004), and the Trans-Alaska Pipeline System in Alaska (Nelson et al., 2001).
Therefore, it is critical to develop mitigation and adaptation strategies in permafrost regions for ongoing climate change.Accurate projection of the degree of permafrost degradation is a prerequisite for developing these strategies.
Significant effort has been made to improve modeling accuracy and efficiency of permafrost dynamics along two primary lines of inquiry.One is to create suitable freezing and The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-11Manuscript under review for journal The Cryosphere Discussion started: 12 February 2018 c Author(s) 2018.CC BY 4.0 License.
thawing algorithms for different applications, including land surface models (Chen et al., 2015;Oleson et al., 2010;Wang et al., 2017), permafrost models (Goodrich, 1978;Langer et al., 2013;Qin et al., 2017) and other related models (Fox, 1992;Woo et al., 2004).The other line of inquiry is focused on schemes of soil physical properties (Chen et al., 2012;Zhang et al., 2011), which play a critical role in permafrost dynamics.For example, thermal diffusivity (thermal conductivity/heat capacity) directly determines how quickly energy can be conducted into and out of permafrost from the top and from the bottom of the permafrost horizon.Porosity determines the maximum amount of water that can be contained in a soil layer, and hydraulic properties determine the exchange of soil water between soil layers.The amount of water then affects not only soil thermal properties, but also determines the large amount of latent heat loss/gain for freezing/thawing.On the QTP, soil is coarse due to weak weathering and strong erosion (Arocena et al., 2012).Soils with gravel content (particle diameter >2 mm) has been reported in several studies (Wang et al., 2011;Wu et al., 2016;Yang et al., 2009;Qin et al., 2015;Chen et al., 2017;Du et al., 2017).These gravelly soil properties are different from those used in current modeling studies (Wang et al., 2013).For example, Soil properties in Community Land Model are calculated from fractions of sand, silt and clay based on measurements of agriculture soils (Oleson et al., 2010).However, those of gravelly soil on the QTP and their effects on permafrost dynamics are under studied (Pan et al., 2017).
In this study we investigated the characteristics of soil physical properties at a site on the central QTP and its effects on permafrost dynamics.We first measured soil physical properties of excavated soil samples in laboratory.We then conducted sensitivity analyses with an ecosystem model by substituting the default soil physical properties by those that we measured.We aimed to emphasize the effects of gravel content on soil physical properties and on permafrost dynamics.It is not our purpose to develop general schemes of soil physical properties for using in modeling studies on the QTP.
The site is located in the continuous permafrost region of the central QTP (Figure 1a al., 2017).Based on the soil map of Li et al. (2015), soil of this region belongs to Gelisols and Inceptisols, which occupy 34% and 28% of the total area of permafrost region of the QTP, respectively.Meteorological variables (air temperature, radiation, precipitation and humidity), soil variables (soil temperature and moisture) down to 1.6 m, and borehole temperatures down to 60 m have been measured at this site since 2002.The meteorological station and borehole are located on a gentle slope with sparse vegetation (Figure 1b).Based on measurements, the mean annual precipitation and air temperature are 366 mm yr -1 and -3.6 o C, respectively.
Active layer depth is ~3.3 m and the lower boundary of permafrost is at a depth of ~20 m.Details of meteorological and borehole variables can be found in Qin et al. (2013).
Above 120 cm in the soil pit, coarse soil material was small enough to be fitted in cut rings.
Below 150 cm, there exists weathered mudstone, which could also be sampled with our cut rings.
We used the KD2 Pro (Decagon, US) to measure thermal conductivity of soil samples.The steps were: 1) soil samples were dried in oven and weighed to calculate bulk density; 2) soil samples were exposed to a constant temperature (20 o C) over 24 h, samples were then saturated with water (20 o C) weighed (0.001g precision), and the KD2 probe (SH-1) was then inserted into soil samples to measure thermal conductivity; 3) samples and the KD2 probe were then put into a refrigerator (0~-26 o C) at -5 o C over 12 h, at which time thermal conductivity was measured.Steps 2 and 3 were repeated at different levels of soil water content.4) Finally, soil samples were immersed into water over 24 h and weighed to calculate porosity; and the saturated unfrozen and frozen thermal conductivity were then measured, accordingly.
We used pressure membrane instruments (1500F1, Soilmoisture Equipment Corp, US) to measure matric potential of soil samples (Azam et al., 2014;Wang et al., 2007).In this study we used both 15 bar and 5 bar pressure chamber.We used soil permeability meter (TST-70, China) to measure saturated hydraulic conductivity of soil samples (Gwenzi et al., 2011).Finally, soil samples were sieved through meshes with diameters of 2.0, 1.0, 0.5 and 0.25 mm in a sequence and weighted to calculate fractions.

Model description
The model used in this study is dynamic organic soil version of Terrestrial Ecosystem Model (DOS-TEM).Models of TEM family simulate the carbon and nitrogen pools of vegetation and soil, and their fluxes among atmosphere, vegetation and soil (McGuire et al., 1992).They have been widely used in studies of cold region ecosystems (e.g.McGuire et al., 2000;Yuan et al., 2012;Zhuang et al., 2004;2010) The DOS-TEM consists of four modules, these being the environmental, ecological, fire disturbance, and dynamic organic soil modules (Yi et al., 2010).The environmental module operates on a daily time interval using mean daily air temperature, surface solar radiation, precipitation, and vapor pressure, which are downscaled from monthly input data (Yi et al., 2009b).The module takes into account radiation and water fluxes among the atmosphere, canopy, snow pack, and soil.

Implementation of soil thermal processes
Earlier versions of TEM did not simulate soil temperature (McGuire et al., 1992).Zhuang et al. (2001) incorporated Goodrich permafrost model into TEM.Yi et al. (2009a) incorporated a two-directional Stefan algorithm to simulate soil freezing and thawing for complex soil situation with changes of organic soil and moisture.Soil temperatures of all soil layers in the DOS-TEM are updated daily.Phase change is calculated first before heat conduction.A twodirectional Stefan algorithm is used to predict the depths of freezing or thawing fronts within the soil (Woo et al., 2004).It first simulates the depth of the front in the soil column from the top downward, using soil surface temperature as the driving temperature.It then simulates the front from the bottom upward using the soil temperature at a specified depth beneath a front as the driving temperature (bottom-up forcing).The latent heat used for phase change is recorded for each soil layer.If a layer contains n freezing or thawing fronts, this layer is then explicitly divided into n+1 soil layers All soil layers are grouped into 3 parts: 1) the soil layers above the uppermost freezing or thawing front; 2) the soil layers below the lowermost freezing or thawing front; and 3) the soil layers between the uppermost and lowermost fronts.
Soil temperatures are then updated by solving finite difference equations of each part with phase change latent heat as energy source or sink (Yi et al., 2014a).Soil surface temperature, which is used as a boundary condition, is calculated using daily air maximum, air minimum, radiation, and leaf area index (Yi et al., 2013).
The version of the DOS-TEM in this study uses the Côté and Konard (2005) scheme to calculate thermal conductivity (Yi et al., 2013;Pan et al., 2017), which is also used by other studies on the QTP (e.g.Chen et al., 2012, Luo et al., 2009).

Implementation of soil hydrological processes
Surface runoff, infiltration, and water redistribution among soil layers are simulated in a similar way as Community Land Model 4 (Oleson et al 2010).Soil matric potential (Ψ) determines the direction of water movement.And hydraulic conductivity determines the rate of water movement.
where Ψ sat is saturated soil matric potential (mm), θ liq is volumetric liquid water content (m 3 m -3 ), and B is pore size distribution parameter.
Several important features relating to permafrost have been considered in the DOS-TEM (see Yi et al., 2014b), e.g.runoff from perched saturated zone and exchanges of water between soil and a water reservoir.Runoff from the perched saturated zone above the permafrost is implemented following Swenson et al. (2013), Where αis an adjustable parameter (0.6 m -1 ), K p is the mean saturated hydraulic conductivity within perched saturated zone (mm s -1 ), z frost and z perched are the depths to permafrost table and perched water table (m), respectively, and θis slope ( o ).
The DOS-TEM has been verified against the Neumann Equation for water, mineral and organic soil under an idealized condition (Yi et al., 2014b), and validated against field measurements for various locations in Alaska, the Arctic, and the QTP (Yi et al., 2009b, Yi et al., 2013, Yi et al., 2014a).

Model inputs and initialization
We used the measured air temperature, downward radiation, precipitation and humidity (monthly) as input to drive the DOS-TEM.Leaf area index, one half of the total green leaf area per unit ground surface area, was specified to be 0.6 m 2 m -2 in July and August, 0.1 m 2 m -2 in April and October, and interpolated linearly in other months.It is used in the DOS-TEM to calculate ground surface temperature in combination with other meteorological variables (Yi et al., 2013).
Soil temperature and moisture were initialized at -1 o C and saturation.Sand and silty clay were used in testing to represent coarse and fine soil textures, respectively, on the QTP (FAO/IIASA/ISRIC/ISSCAS/JRC, 2009).The temperature gradient at the bottom of bedrock was set to be 0.06 o C cm -1 based on borehole observations.Volumetric unfrozen liquid water in winter was set to be 0.1 based on observations.The DOS-TEM was the spun up using the driving data from 2003-2012 for 100 yr.

Sensitivity analyses
We considered two soil textures in the DOS-TEM, i.e. silty clay and sand.The former is considered as the major soil type of this region (Lin et al., 2011), and the latter is the coarsest soil type considered in land surface modeling (Oleson et al., 2010).We first ran the DOS-TEM using the default porosity, soil thermal conductivity (Equation 1), hydraulic conductivity (Equation 5) and matric potential schemes of these two soil types (Equation 4).
The parameters Φ,Ψ sat , K sat and B were calculated based on soil texture used in Community Land Model (Equation 4 and 5; Oleson et al., 2010).We then substituted the original values of Φ,Ψ sat , K sat and B based on laboratory measurements and calibration.We did not calibrate soil thermal conductivity to retrieve parameters of Equation 2 and 3. Instead, we interpolated measured thermal conductivity over a range of the degree of saturation (0 to 1), which was used as a lookup table by the DOS-TEM.Therefore, our sensitivity analyses considered a set of 4 factors, i.e. porosity, matric potential (Ψ sat and B), hydraulic conductivity (K sat and B) and thermal conductivity.We also analyzed 3 different slopes ( 0 In summary, our sensitivity analyses with the DOS-TEM involved 288 different combinations of parameter values. We did not measure the heat capacity of gravelly soil.The maximum and minimum heat capacities of mineral soil types considered in land surface model are 2.355 and 2.136 MJ m -3 , respectively.The relative difference is less than 10%.Therefore, in this study, we did not make sensitivity tests using thermal diffusivity (the ratio between thermal conductivity and heat capacity).

Soil porosity, particle size and bulk density
The mean weight fraction of gravel (particle size diameter > 2 mm) of different soil layers ranged from 0.38 to 0.65 with a mean of 0.55 (Table 1).The weight fraction of soil with particle size diameter > 0.25 mm ranged from 0.77 to 0.86 with a mean of 0.84 among layers.
The default porosities of sand and silty clay were 37.3% and 48.1%, respectively.The mean porosity of 2 m depth ranged from 21% to 30% with a mean of 27%.The mean bulk density ranged from 1.61 to 1.86 g cm -3 with a mean of 1.74 g cm -3 .No significant relationships were found among soil porosity, bulk density and the fraction of gravel.

Thermal conductivity
The mean unfrozen dry soil thermal conductivity of different soil layers ranged from 0.24 to 0.40 W m -2 with a mean of 0.36 W m -2 (Table 2).The mean frozen dry soil thermal conductivity ranged from 0.25 to 0.41 W m -2 with a mean of 0.35 W m -2 .The difference of dry thermal conductivity between frozen and unfrozen states was small.The mean unfrozen saturated soil thermal conductivity of different soil layers ranged from 2.15 to 2.74 W m -2 with a mean of 2.48 W m -2 (Table 2).The mean frozen saturated soil thermal conductivity ranged from 3.06 to 3.72 W m -2 with a mean of 3.33 W m -2 .The difference of saturated thermal conductivity between frozen and unfrozen states was about 0.85 W m -2 .There existed a threshold of soil wetness, below which frozen soil thermal conductivity was slightly smaller than unfrozen soil (Figure 2a).
The default dry frozen and unfrozen thermal conductivities using Côté and Konard (2005) scheme of sand and silty clay were about 0.42 and 0.22 W m -2 , respectively.The saturated frozen and unfrozen thermal conductivities of sand were 3.11 and 1.90 W m -2 , respectively.
Those of silty clay were about 2.35 and 1.24 W m -2 , respectively (Figure 2b).The default dry frozen and unfrozen thermal conductivities using Farouki scheme of sand and silty clay were about 0.97 and 0.33 W m -2 , respectively.The saturated frozen and unfrozen thermal conductivities of sand were 5.21 and 3.18 W m -2 , respectively.Those of silty clay were about 2.87 and 1.52 W m -2 , respectively (Figure 2c).

Saturated hydraulic conductivity
The mean saturated hydraulic conductivity of soil layers ranged from 0.0036 to 0.0315 mm s -1 .
The maximum saturated hydraulic conductivity was about 8.7 times larger than the minimum (Table 3).The saturated hydraulic conductivity tended to be larger with increasing proportion of gravel content in the soil samples (Figure 3a), and was about 0.03-0.06mm s -1 for some samples with gravel content greater than 70%.The default saturated hydraulic conductivities of sand and silty clay were 0.024 and 0.0011 mm s -1 , respectively.

Matric potential
Saturated matric potential and B were fitted using measured matric potential values.The correlation coefficients between calculated and fitted matric potential were all greater than 0.96.The mean absolute value of saturated matric potential of soil layers ranged from 27.02 to 603.7 mm, and those of B ranged from 5.22 to 1.89 (Table 3 and Figure 3b).The default absolute value of saturated matric potential of sand and silty clay were 47.29 and 632.99 mm, respectively, and the B values 3.39 and 10.38, respectively.

Soil temperature
The mean root mean squared errors (RMSEs) between monthly measured soil temperatures and model runs with measured parameters using different combination of soil thicknesses (3.25, 4.25 and 5.25 m) and slopes (0, 5 and 10 o ) were about 1.07 o C at 20 cm (Figure 4c).
The mean RMSEs for all model runs with default sand and silty clay parameters were about 0.97 and 1.37 o C, respectively.For other soil layers, the RMSEs of model runs with measured parameters were much smaller than those with default sand and silty clay parameters (Figure 4d-l).The simulated soil temperatures using default sand and silty clay parameters were all lower than measured ones in summer at 100 and 200 cm; and in winter at 400 cm.The RMSEs can be as large as 2.91 o C (Figure 4e).
The standard deviations of soil temperatures among different slopes and soil thicknesses using measured parameters were larger than those using the default parameters (Figure 4 deviations using default silty clay parameters were smaller (<0.06 o C at all depths) than those using default sand parameters.

Soil liquid water
The mean RMSEs between monthly measured liquid soil volumetric water content (VWC) and model simulations with measured parameters ranged from 0.03 to 0.09, which were smaller than RMSEs for sand and silty clay parameters (Figure 5).The model simulations for silty clay parameters have larger RMSEs than those for sand parameters.VWCs were always overestimated in warm seasons at depths of 10, 40 and 80 cm.VWCs were underestimated at a depth of 160 cm, where the simulated soil was frozen.All model simulations overestimated VWC at 40 cm, where the maximum measured VWCs were about 0.1 (Figure 5d-f).
The standard deviations of VWC among different slopes and soil thicknesses using sand parameters were about 0.077, which were larger than those using measured parameters (~0.062).The standard deviations of VWC using silty clay parameters (<0.011) were less than those using measured parameters.

Active layer depth (ALD)
The mean RMSEs between measured ALDs (derived from linear interpolation of soil temperatures) and modelled ALDs (simulated explicitly) were about 1.06, 1.83 and 0.28 m for model runs with sand, silty clay and measured parameters (Figure 6a).The mean standard deviations were about 0.088, 0.007 and 0.28 m.All simulations using sand and silty clay parameters underestimated ALDs.

Permafrost lower boundary (PLB)
The mean RMSEs between measured PLBs (derived from linear interpolation of temperatures) and modelled PLBs (derived from linear interpolation of simulated bed rock temperatures) were about 10.25, 7.96 and 6.71 m for model runs with sand, silty clay and measured parameters (Figure 6b).The mean standard deviations were about 1.89, 0.29 and 6.62 m.All simulations using sand and silty clay parameters overestimated PLBs.

Model sensitivity analyses
Deep soil layers used in models are usually specified as being thick.For example, a 1 m thick soil layer was used in our simulations starting around 3 m soil depth.Soil temperatures at this depth are usually close to 0 o C. Therefore, the RMSEs of deep soil layers were small and did not facilitate evaluation of model sensitivities.In the following subsections, we used 20 and 100 cm soil temperatures, ALDs and PLBs for sensitivity analysis.

Hydraulic conductivity/Matric potential
Replacing default sand or silty clay hydraulic conductivity with measured parameters had very small effects on mean RMSEs of soil temperatures and ALDs (Figure 7 and 8).The same was true for matric potential.When hydraulic conductivity of default sand or silty clay was substituted, mean RMSEs of PLB were decreased or increased, however, when matric potential was substituted, mean RMSEs of PLBs were increased or decreased, respectively.
When hydraulic conductivity or matric potential parameters were substituted in default sand or silty clay parameters, mean RMSEs of VWC changed slightly (Figure 9 and 10).

Effects of combined parameters
We compared model simulations with different combinations of measured parameters (porosity, thermal conductivity, hydraulic conductivity and matric potential) with those with one substituted measured parameter.We selected those model runs with less RMSEs than any of model runs with one substituted measured parameter (Table 4 and 5).We didn't consider the 10 cm soil temperature, which were similar among all model runs.
For sand, model simulations with porosity and thermal conductivity or hydraulic conductivity substituted had 4 outcomes with lower RMSEs (Table 4 and Figures 7 and 9).
Only 2 out of 7 outcomes had lower RMSEs with all 4 parameters were substituted.Among all the 18 cases with RMSEs less than the individual "best" RMSE, porosity was included 18 times, followed by thermal conductivity and hydraulic conductivity with 10 times.
For silty clay, model simulations with porosity and thermal conductivity and/or matric potential substituted had 5 outcomes with lower RMSEs (Table 5 and Figures 8 and 10).
Among all the 29 cases with RMSEs less than the individual "best" RMSE, porosity was included 29 times, followed by thermal conductivity with 20 times and matric potential with 16 times.

Effects of slope and soil thickness
Changes of slope alone had small effects on simulated soil temperatures and ALDs (Figures 7    and 8).An increase of slope generally reduced RMSEs of VWCs (Figures 9 and 10).Model simulations with porosity substituted had smaller difference of VWC RMSE between different cases of slopes.For example, the mean RMSEs of model simulations with slope of 0 o or 5 o and porosity substituted in default sand parameters were 0.078 or 0.048, respectively.
While those with porosity not substituted were 0.141 or 0.055, respectively.Similarly, the mean RMSEs of model simulations using default silty clay parameters with porosity substituted were 0.081 or 0.06 for slope of 0 o or 5 o , respectively.The mean RMSEs were 0.21 or 0.138 with porosity not substituted, respectively.For a further increase of slope to 10 o , changes of RMSEs of VWCs at depths of 10-160 cm were small.
Soil thickness had small effects on 20 and 100 cm soil temperatures and 10-160 cm VWCs, and it had prominent effects on PLB for a few cases only with a slope of 10 o (Figures 7 and 8).

Characteristics of soil physical properties
Although the effects of gravelly soil on permafrost dynamics have been considered in a few modelling studies, the thermal and hydraulic properties of gravelly soil were calculated without validation or calibrated (Pan et al., 2017;Wu et al., 2018).To our knowledge, this is the first study measuring physical properties of gravelly soil samples from permafrost region of the QTP.
The weight fraction of gravel (diameter > 2mm) in the soil samples we analysed was greater than 55% on average.While the typical soil types considered in land surface models and other models usually have much smaller diameter.For comparison, the fractions of gravel considered in Pan et al. (2017) ranges from 5% to 33% and from 10% to 28% for the Madoi and Naqu sites, respectively.The Beiluhe site and the aforementioned sites are located in regions with Gelisols and Inceptisols, which occupy ~62% of the permafrost regions of the QTP (Li et al., 2015).It is possible that gravelly soil commonly exists on the QTP.The saturated hydraulic conductivity and matric potential of soil samples measured in this study were more similar to sand than to silty clay (see Section 3.1).It is consistent with the study of Wang et al. (2013) that coarse soil material has poor water holding capability.
The measured saturated thermal conductivities of soil samples were relatively close to those estimated by the Côté and Konard (2005) scheme.But they were much less than those estimated by the Farouki scheme (Figure 2).Several other studies also found that Farouki scheme overestimated soil thermal conductivity (Chen et al. 2012;Luo et al., 2009).
One important finding of this study is the relatively small value of porosity.The measured porosity ranged from 0.206 to 0.302, which is less than those of soil types considered in land surface models.For example, the porosities of mineral soil types considered in Community Land Model range from 0.37 to 0.48 (Oleson et al., 2010).Porosity determines the maximum water stored in a soil layer, and affects soil thermal conductivity, hydraulic conductivity and matric potential (Equation 2-5).It plays a more important role than other parameters in simulated soil thermal and hydrological dynamics (Table 4 and 5

Effects of soil water on permafrost dynamics
Soil water not only affects soil thermal properties, e.g.thermal conductivity and heat capacity, but also affects the amount of latent heat lost or gained, for freezing or thawing, respectively (Goodrich, 1978;Farouki, 1986).Soil water is determined by infiltration, evapotranspiration, water movement among soil layers, subsurface runoff and exchange with a water reservoir.
Therefore, processes or parameters that affect soil water dynamics will also affect permafrost dynamics.This study quantitatively assessed the effects of soil water on permafrost dynamics.
For example, when default silty clay parameters with high porosity and low saturated hydraulic conductivity were used, soil layers were almost saturated (Figure 5).The simulated ALDs were about 1.47 m, which was less than half of measured ALDs (Figure 6a).When the slope was 0 o , subsurface runoff didn't happen in saturated zone above the bottom of the active layer.The simulated soil water content was generally higher in the active layer.However, when the slope was 5 o , the simulated soil water content was less and the RMSE was smaller (Figure 9 and 10).These patterns were especially obvious when both porosity and saturated hydraulic conductivity were large (Equation 6; Figure 9 and 10).Other studies have also emphasized the importance of subsurface runoff above the bottom of the active layer (Frey and McClelland, 2009;Walvoord and Striegl, 2007).The effects of soil water content on soil thermal dynamics increased with soil and rock depth (Figure 7 and 8).The biggest effects were on PLB, which became manifest during long-term spinup procedures.
Land surface models generally represent soil water dynamics (e.g.Chen et al., 2015;Oleson et al., 2010;Wang et al., 2017).However the permafrost processes in these models usually use specified thermal properties, which were static during model simulations (Li et al., 2009;Nan et al., 2005;Qin et al., 2017;Zou et al., 2017).As shown in this study, when permafrost degraded, the thermal and hydrological regimes of soil also changed.It is critical to simulate soil water dynamics to properly project permafrost dynamics in the future.

Sampling and laboratory measurement
We used cut rings with 10 cm diameter to take soil samples.There are weathered mudstones in our study site, which can be sampled in cut rings.However, it is very likely that there are The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-11Manuscript under review for journal The Cryosphere Discussion started: 12 February 2018 c Author(s) 2018.CC BY 4.0 License.
soil samples with much bigger gravels.Therefore, larger containers should be used to take samples for further laboratory analysis in the future.
During our laboratory work, we found two phenomena.First, we originally used the QL-30 thermophysical instrument to measure thermal conductivity.It worked properly under unfrozen condition.However, when frozen, surface of soil samples was uneven due to frost heave.The contact between plate of QL-30 and soil sample surface was not ideal.The measured frozen thermal conductivities were smaller than unfrozen thermal conductivity even for the case of saturation, which were definitely wrong.The second phenomenon was that there seems to be a threshold of soil wetness, below which unfrozen soil thermal conductivity is greater than frozen soil thermal conductivity (Figure 2a).This pattern was somewhat exhibited in estimates of the Côté and Konard (2005) scheme, but not in the estimates of the Farouki scheme (Figure 2c).More measurements using instruments with higher accuracy should be made in the future.

Model simulation
Although the DOS-TEM using measured parameters provided satisfactory results, there are some aspects requiring further improvement in the future.For example, the measured soil moistures at 40 cm depth were less than 0.1 m 3 /m 3 .However, the simulated soil moistures were always much greater (Figure 5f).There were spikes of measured soil moistures at 80 and 160 cm depths, which were not presented in simulation (Figure 5 i and l).In the DOS-TEM, the unfrozen soil water content, or supercold water, was prescribed to be 0.1 m 3 /m 3 .When soil is freezing, if soil liquid water content is less than this value, no phase change will happen (Figure 5k).It is ideal to simulate the dynamics of unfrozen soil water content (Romanovsky and Osterkamp, 2000).
Field studies have shown that gravel content in root zone affects vegetation growth (Qin et al., 2015), which affects ground surface temperature (Yi et al., 2013).In the current study, we used specified leaf area index.The fractions of gravel content in soil are also dynamic.For example, Chen et al. (2017) found that plateau pika excavated subsurface soil with gravel on to surface.Fine soil particles were carried away by wind and water erosion, which resulted in gravel remaining at the surface.Our ongoing research is working towards representing the coupling of vegetation growth, small mammal disturbances, and soil erosion on permafrost dynamics of the QTP in the future.

Regional applications
Soil texture playes an important role in permafrost dynamics (Figure 6).However, the dominant soil texture on the QTP from FAO/IIASA/ISRIC/ISSCAS/JRC ( 2009) is silty clay.
The specification of silty clay in simulations results in estimates of ALD that are much smaller than meansurements (Yi et al., 2014a).To properly simulate the distribution and dynamics of permafrost on the QTP under climate change scenarios, it is important to develop proper schemes of soil physical properties in relation to gravel content and to develop regional datasets of soil texture for input.
Gravel content affects soil physical properties.For example, soil porosity and saturated hydraulic conductivity are determined by the fraction of gravel, diameter and degree of mixture (Zhang et al., 2011).Organic soil carbon content in mineral soil on the QTP affects soil porosity and thermal conductivity (Chen et al., 2012).We only measured physical properties of soil samples from one site on the QTP in this study.The mean gravel fraction was about 0.55, which is definitely not suitable for direct use in regional simulations.More laboratory work is needed to develop proper schemes for representing mixed soil with fine mineral, gravel and organic carbon in permafrost models.It is the first priority to develop schemes that make use of porosity data sets, due to its importance and simplicity of measurement.
The development of a spatially explicit dataset of soil texture is also required for regional applications of projecting permafrost changes on the QTP.One way is to collect relevant data through extensive field campaigns (e.g., Li et al., 2015).Currently, gravelly soil has only been mentioned in scientific literature on the QTP (Chen et al., 2015;Wang et al., 2011;Yang et al., 2009), and no systematic analysis of gravelly soils exists across the QTP.Ground penetrating radar is a feasible tool to retrieve soil thickness above gravel layer (Han et al., 2016).
Unmanned aerial vehicles has been used recently (Yi, 2017), and gravel on the ground surface can be identified easily in aerial photos (Chen et al., 2017).In combination with ancillary datasets, e.g.geomorphology, topography, vegetation, it is possible to generate spatial datasets of soil texture (Li et al., 2015;Wu et al., 2016).Another way is to retrieve soil physical properties using data assimilation technology, e.g.Yang et al. (2016) assimilated porosity using a land surface model and microwave data.

Conclusions
In this study, we excavated soil samples from a permafrost site on the central QTP and measured soil physical properties in laboratory.Gravel was common in the soil profile and porosity was much smaller than the typical soil types used in land surface models.We then performed sensitivity analysis of these parameters on soil thermal and hydrological processes within a terrestrial ecosystem model.When default sand or silty clay parameters were substituted with measured soil properties, the model errors of soil temperature, soil liquid water content, active layer depth and permafrost low boundary were generally reduced.
Sensitivity analyses showed that porosity played a more important role in reducing model errors than other soil properties examined.Though it is unclear how representative this soil is in the QTP, it is clear that soil physical properties specific to the QTP should be used to properly project permafrost dynamics into the future.
); and they increased from 0.40 o C at 100 cm to 0.61 o C at 200 cm (Figure 4f and i).The standard The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-11Manuscript under review for journal The Cryosphere Discussion started: 12 February 2018 c Author(s) 2018.CC BY 4.0 License.
Replacing default sand or silty clay porosity with measured porosities changed mean RMSEs of soil temperatures (model runs with 3 different slopes and 3 different soil thicknesses at 2 different soil depths) from 1.18 or 2.11 o C to 1.25 or 1.08 o C, respectively (Figure7 and 8).Mean RMSEs of ALD were reduced from 1.06 or 1.84 m to 0.22 or 0.83 m, respectively.Mean RMSEs of PLB were changed from 10.26 or 7.96 m to 6.61 or 10.26 m.Mean RMSEs of VWC were reduced from 0.074 or 0.20 to 0.06 or 0.073 when measured porosities were used for replacing default sand or silty clay porosity, respectively (Figure9and 10).Thermal conductivityReplacing default sand or silty clay thermal conductivity with measured thermal conductivity reduced mean RMSEs of soil temperatures from 1.18 or 2.11 o C to 1.02 or 1.33 o C, respectively (Figure7 and 8).Mean RMSEs of ALD were reduced from 1.06 or 1.84 m to 0.56 or 1.18 m, respectively.Mean RMSEs of PLB were changed from 10.26 or 7.96 m to 4.18 or 2.54 m, respectively.Mean RMSEs of VWC changed very slightly (Figure9and 10).

Figure 1 .
Figure 1.a) The location of Beiluhe permafrost station on the Qinghai-Tibetan Plateau 1 (Permafrost type is from Li and Cheng, 1996); b) the aerial view of the meteorological station 2 and the excavated soil pit; and c) the detailed view of the excavated soil pit. 3

Figure 2 .
Figure 2. The relationship between soil wetness (solid and dotted lines represent frozen and unfrozen cases) and soil thermal conductivity (W/m o C) from a) measured values (Measured; dots and empty diamonds represent measured frozen and unfrozen soil thermal conductivities, respectively), b) using the Côté and Konard (2005) scheme (CK); and c) using the Farouki (1986) scheme(Farouki).Thick and thin lines represent relationships for sand and silty clay, respectively.

Figure 3
Figure 3. a) the relationship between saturated hydraulic conductivity (mm s -1 ) and gravel content fraction (Solid dots represent measured value; empty circle and empty triangle represent the corresponding values of sand and silty clay used in Community Land Model, respectively) ; b) the relationship between soil wetness (lines) and absolute value of matric potential (mm) at three representative depths.Solid and dashed lines represent default values of sand and silty clay, respectively(Oleson et al., 2014).

Figure 4 .
Figure 4. Comparisons of soil temperatures simulated using default parameters of sand, silty clay, and measured parameters (lines) with measured soil temperatures (dots) at 20, 100, 200 and 400 cm depths.Error bars showed the standard deviation calculated based on 9 simulations with 3 different slopes and 3 different soil thicknesses.

Figure 5 .
Figure 5. Comparisons of soil volumetric liquid water content simulated using default parameters sand, default silty clay, and measured parameters (lines) with measured soil moistures (dots) at 10, 40, 80 and 160 cm depths.Error bars showed the standard deviation calculated based on 9 simulations with 3 different slopes and 3 different soil thicknesses.

Figure 6
Figure 6.a) Contours of measured soil temperature ( o C) from borehole measurements down to 5 m and simulated active layer depth over the period of 2003-2011; and b) same as a) but down to 50 m and for simulated permafrost low boundary.Black, blue and magenta represent simulations with silty clay, sand and measured parameters, respectively.Error bars show the standard deviation calculated based on 9 simulations with 3 different slopes and 3 different soil thicknesses.

Figure 7 .
Figure 7. Root mean squared errors between measurements and model simulations (with 1 different combinations of measured porosity (I), thermal conductivity (II), hydraulic 2 conductivity (III) and matric potential (IV) of default sand parameters) for a) 20 and b) 100 3 cm soil temperatures ( o C), c) active layer depth (ALD, m) and d) permafrost low boundary 4 (PLB, m).O and All represent model runs without substitution of default parameters and with 5 all 4 parameters substituted, respectively.Mean and standard deviation of model simulations 6 with 3 different soil thicknesses at each slope are shown.7 8

Figure 8 .
Figure 8. Same as Figure 7 but for default silty clay parameters.

Figure 9 .
Figure 9. Same as Figure 7 but for a) 10 cm, b) 40 cm, c) 80 cm and d) 160 cm soil liquid water content.

Table 1 .
The mean (standard deviation) of measured soil bulk density, porosity, and particle 1 size diameter fractions of different layers based on soil samples in this study.The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-11Manuscript under review for journal The Cryosphere Discussion started: 12 February 2018 c Author(s) 2018.CC BY 4.0 License.

Table 2 .
The mean (standard deviation) of the measured frozen and unfrozen dry and 1 saturated soil thermal conductivity (W m -2 ) of different soil layers.The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-11Manuscript under review for journal The Cryosphere Discussion started: 12 February 2018 c Author(s) 2018.CC BY 4.0 License.

Table 4 .
Model performance of substituting default sand parameters with measured porosity 1

Table 5
Same as Table 4, but for silty clay.