The physical properties of coarse-fragment soils and their effects on permafrost dynamics: a case study on the central Qinghai–Tibetan Plateau

Soils on the Qinghai–Tibetan Plateau (QTP) have distinct physical properties from agricultural soils due to weak weathering and strong erosion. These properties might affect permafrost dynamics. However, few studies have investigated both quantitatively. In this study, we selected a permafrost site on the central region of the QTP and excavated soil samples down to 200 cm. We measured soil porosity, thermal conductivity, saturated hydraulic conductivity, and matric potential in the laboratory. Finally, we ran a simulation model replacing default sand or loam parameters with different combinations of these measured parameters. Our results showed that the mass of coarse fragments in the soil samples (diameter >2 mm) was ∼55 % on average, soil porosity was less than 0.3 m³ m⁻³, saturated hydraulic conductivity ranged from 0.004 to 0.03 mm s⁻¹, and saturated matric potential ranged from −14 to −604 mm. When default sand or loam parameters in the model were substituted with these measured values, the errors of soil temperature, soil liquid water content, active layer depth, and permafrost lower boundary depth were reduced (e.g., the root mean square errors of active layer depths simulated using measured parameters versus the default sand or loam parameters were about 0.28, 1.06, and 1.83 m). Among the measured parameters, porosity played a dominant role in reducing model errors and was typically much smaller than for soil textures used in land surface models. We also demonstrated that soil water dynamic processes should be considered, rather than using static properties under frozen and unfrozen soil states as in most permafrost models. We conclude that it is necessary to consider the distinct physical properties of coarse-fragment soils and water dynamics when simulating permafrost dynamics of the QTP. Thus it is important to develop methods for systematic measurement of physical properties of coarse-fragment soils and to develop a related spatial data set for porosity.

1 Introduction

Permafrost underlies 25 % of Earth's 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). Permafrost thaw 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, such as QTP railways and roads (Wu et al., 2004) or 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 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 and Ward, 2011), which play a critical role in permafrost dynamics. For example, porosity determines the maximum amount of water that can be contained in a soil layer, thermal properties determine the heat conduction within soil layers, and hydraulic properties determine the exchange of soil water between soil layers. The soil water content also determines the large amount of latent heat lost or gained by freezing or 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) have been reported in several studies (Chen et al., 2017; Du et al., 2017; Qin et al., 2015; Wang et al., 2011; Wu et al., 2016; Yang et al., 2009). These soil properties are likely different from those used in current modeling studies (Wang et al., 2013). For example, soil properties in the Community Land Model are calculated from fractions of sand, silt and clay based on measurements of agriculture soils (Oleson et al., 2010). However, the physical properties of coarse-fragment soils within the QTP and their effects on permafrost dynamics are understudied (Pan et al., 2017).

In this case study we investigated the characteristics of soil physical properties at a site on the central QTP and their effects on permafrost dynamics. We first measured soil physical properties of excavated soil samples in a laboratory. We then conducted a sensitivity analysis with an ecosystem model by substituting the default soil physical properties with those that we measured. We aimed to emphasize the effects of coarse-fragment content on soil physical properties and on permafrost dynamics, rather than develop general schemes of soil physical properties for use in modeling studies on the QTP.

Figure 1 Locations of (a) Beiluhe permafrost station on the Qinghai–Tibetan Plateau, and (b) the weather station and the surrounding environment (map data: Google, DigitalGlobe).

2 Methods

2.1 Site description

The site (34^∘49^′46.2${}^{\prime \prime }$ N, 92^∘55^′56.58${}^{\prime \prime }$ E, 4628 m a.s.l.) is located in the Beiluhe Basin, in the continuous permafrost region of the central QTP (Fig. 1a, Zou et al., 2017). Based on the map of Li et al. (2015), soils of this region belong to Gelisols and Inceptisols, which occupy 34 % and 28 % of the total area of permafrost region of the QTP. Land surface types include alpine meadow, alpine steppe, barren surface, and thermokarst lakes (Fig. 1b; Lin et al., 2011).

The site is on top of upland plain landforms, which are formed from fluvial and deluvial sediments. The surficial sediments are dominated by fine to gravelly sands and stones (Fig. 2; Yin et al., 2017). Soils at this site are Inceptisols (Wangping Li, Lanzhou University of Technology, personal communication, 18 April 2018) that are commonly underlain by mudstone. The plant community type is mainly alpine meadow, which is dominated by monocotyledonous species, primarily Poaceae and Cyperaceae. The dominant species are Kobresia pygmaea, accompanied by Elymus nutans, Carex moorcroftii, Oxytropis pusilla, Tibetia himalaica, Leontopodium nanum, and Androsace tapete (Fig. 2c–e).

Figure 2 Images of site conditions: (a) the aerial view of the weather station and the excavated soil pit (the borehole is located in the lower-left corner of white fence), (b) the detailed view of the excavated soil pit, and (c–e) examples of vegetation, gravel, and stones (iron frame is about 0.5 m × 0.5 m).

A weather station was set up in 2002 (Fig. 2a) to measure air temperature and relative humidity (2.2 m, HMP45C-L11/L36, Campbell Scientific Inc., USA), solar radiation (MS-102, EKO, Japan), and precipitation (QMR102, Vaisala Company, Finland). Soil temperatures were measured at depths of 5, 10, 20, 40, 80, and 160 cm using a PT-100 (EKO, Japan); soil moisture were measured at depths of 20, 40, 80, and 160 cm using a CS616-L50 (EKO, Japan). A CR3000 data logger (Campbell Scientific Inc., USA) was used to store these data at 30 min intervals. These readings were averaged or summed (e.g., precipitation) into monthly values to drive and validate the model. Based on measurements, multi-year mean annual air temperature, precipitation, downward solar radiation, and relative humidity were −3.61 ^∘C, 365.7 mm, 206.3 W m⁻² and 51.1 % (Fig. 3). The multi-year mean summer (June to August) air temperature and precipitation were 5.27 ^∘C and 248.3 mm. The multi-year mean winter (December to February) air temperature and precipitation were −12.44 ^∘C and 5.3 mm. The multi-year mean annual, summer, and winter soil temperatures at 40 cm were 0.17, 6.65, and −7.15 ^∘C. Those at 80 cm were 0.11, 4.32, and −4.86 ^∘C.

A borehole was drilled in 2002, and thermistors made by the State Key Laboratory of Frozen Soil Engineering, Chinese Academy of Sciences were installed at 0.5 m intervals from 0.5 to 10 m, at 2 m intervals from 12 to 30 m, at 4 m intervals from 34 to 50 m,, and at 55 and 60 m. Temperature accuracy of this type of thermistor is ±0.05 ^∘C (Wu et al., 2016). The temperatures were recorded on the 5th and 20th days of each month using CR3000 data logger (Campbell Scientific Inc., USA). Based on our measurements, the active layer depth is ∼3.3 m, the depth of zero annual amplitude is ∼6.2 m, and the lower boundary depth of permafrost is at a depth of ∼20 m. The multi-year mean ground temperatures at 0.5, 6, and 60 m are about −0.52, −0.30, and 1.81 ^∘C.

Figure 3 Time series of data measured at the Beiluhe weather station, Qinghai–Tibetan Plateau, 2003 to 2011: (a) air temperature (TA, ^∘C), (b) downward solar radiation (R, W m⁻²), (c) precipitation (PREC, mm), and (d) relative humidity (RH, %).

2.2 Soil sampling and measurement

Permafrost dynamics are affected by atmosphere, vegetation, and soil textures; therefore, we excavated soil close to the weather station and borehole (Fig. 2a) down to 2 m (Fig. 2b) in August 2014. We used cut rings (10 cm diameter, 6.37 cm height and 500 cm³) to take soil samples at depth ranges of 0–10, 10–20, 20–30, 40–50, 70–80, 110–120, 150–160, and 190–200 cm. Three replicates were sampled from the top of each depth range and sealed for analysis in the laboratory. Above 120 cm in the soil pit, coarse-soil material was small enough in the cut rings. Below 150 cm, the material is weathered mudstone, which could also be sampled with our cut rings. Based on the excavated soil pit and measured soil temperature, this site is dominated by Inceptisols with the suborder Gelept (soil taxonomy, ST, Soil Survey Staff, 2014). The soil pit consists of A horizon (∼20 cm), Bw horizon (∼20–80 cm) and C material dominated by fractured bedrock.

We used the KD2 Pro (Decagon, US) to measure thermal conductivity of soil samples. The steps we took to determine soil properties for each sample were as follows: (1) the soil sample was dried in an oven and weighed (0.001 g precision) to calculate bulk density; then (2) the soil sample was exposed to a constant temperature (20 ^∘C) for 24 h, after which a certain volume of water was injected into the soil samples and a KD2 Pro (Decagon, USA) was used to measure the thermal conductivity; next (3) the sample and the KD2 probe were put into a refrigerator at −15 ^∘C for 12 h and thermal conductivity was measured again; (4) steps 2 and 3 were repeated at increasing levels of soil volumetric water content until soil samples were up to the point of saturation; finally, (5) the soil sample was saturated by immersion in water under atmospheric pressure for 24 h and then it was weighed to calculate porosity, and the saturated unfrozen and frozen thermal conductivity were measured accordingly. The bulk density (ρ_b, g cm⁻³), porosity (φ_m, m³ m⁻³) and volumetric water content (θ_liq, m³ m⁻³) were calculated with the following equations:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{b}}={\displaystyle \frac{m_{\mathrm{dry}}-m_{\mathrm{cr}}}{V_{\mathrm{cr}}}}\end{array}$$

$$\begin{array}{}\text{(2)}& {\mathrm{\varnothing }}_{\mathrm{m}}={\displaystyle \frac{m_{\mathrm{sat}}-m_{\mathrm{dry}}}{V_{\mathrm{cr}}}}/{\mathit{\rho }}_{\mathrm{w}}\end{array}$$

$$\begin{array}{}\text{(3)}& {\mathit{\theta }}_{\mathrm{liq}}={\displaystyle \frac{{\mathrm{Wm}}_{\mathrm{all}}-m_{\mathrm{dry}}}{V_{\mathrm{cr}}}}/{\mathit{\rho }}_{\mathrm{w}},\end{array}$$

where m_dry, m_sat, m_all, and m_cr are masses of oven-dried samples, saturated samples, samples with some water with cut ring, and empty cut ring (g), respectively. V_cr is the volume of the cut ring (cm³). ρ_w is the density of water (1 g cm⁻³). We also calculated porosity from bulk density (φ_c, g m⁻³):

$$\begin{array}{}\text{(4)}& {\mathrm{\varnothing }}_{\mathrm{c}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{p}}}},\end{array}$$

where ρ_p is particle density (2.65 g cm⁻³).

We used pressure membrane instruments (1500F1, Soilmoisture Equipment Corp, US) to measure the matric potential of soil samples (Azam et al., 2014; Wang et al., 2007), using both 15 bar and 5 bar pressure chambers. Pressure values were set to 0, 10, 20, 40, 60, 80, 100, 150, 200, 300, and 400 kpa. It usually took 3–4 days to finish one measurement at one pressure level. We used a soil permeability meter (TST-70, Nanjing T-Bota Scietech Instruments & Equipment Co., Ltd. China) to measure saturated hydraulic conductivity of soil samples (Gwenzi et al., 2011). Finally, soil samples were sieved through a 2.0 mm mesh, and soil particle size distribution was determined with a laser diffraction analyzer (Malvern-2000, Worcestershire, UK).

2.3 Model description

To simulate soil temperatures, soil liquid water content, temperature in rock layers, active layer depth (ALD) and permafrost low boundary depth (PLB) dynamics we used a dynamic organic soil version of the Terrestrial Ecosystem Model (DOS-TEM). Models from the 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: environmental, ecological, fire disturbance, and dynamic organic soil (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., 2009a). The module takes into account radiation and water fluxes among the atmosphere, canopy, snowpack, and soil.

2.3.1 Implementation of soil thermal processes

Earlier versions of TEM did not simulate soil temperature (McGuire et al., 1992). Zhuang et al. (2001) incorporated Goodrich (1978) permafrost model into TEM. Yi et al. (2009b) incorporated a two-directional Stefan algorithm to simulate soil freezing and thawing for complex soils with changes in the soil organic and moisture content. Temperatures of all soil layers in the DOS-TEM are updated daily. Phase change is calculated first before heat conduction. A two-directional 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 three parts: (1) those above the uppermost freezing or thawing front; (2) those below the lowermost freezing or thawing front; and (3) those between the uppermost and lowermost fronts. Soil temperatures are then updated by solving finite difference equations of each part with latent heat from phase change as an 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 Konrad (2005) scheme to calculate thermal conductivity (Yi et al., 2013; Pan et al., 2017), which is also been used by other studies on the QTP (e.g., Chen et al., 2012; Luo et al., 2009), and is as follows:

$$\begin{array}{}\text{(5)}& \mathit{\lambda }=\left\{\begin{array}{l}{k}_{\mathrm{e}}{\mathit{\lambda }}_{\mathrm{sat}}+\left(\mathrm{1}-k_{\mathrm{e}}\right){\mathit{\lambda }}_{\mathrm{dry}}s>{\mathrm{10}}^{-\mathrm{5}}\\ {\mathit{\lambda }}_{\mathrm{dry}}s\le {\mathrm{10}}^{-\mathrm{5}}\end{array}\right.,\end{array}$$

where λ, λ_sat, λ_dry are soil thermal conductivity, saturated soil thermal conductivity, and dry soil thermal conductivity (W m⁻¹ K⁻¹), and k_e is the Kersten number (Côté and Konrad, 2005). Dry thermal conductivity varies with soil properties according to the following:

$$\begin{array}{}\text{(6)}& {\mathit{\lambda }}_{\mathrm{dry}}=\mathit{\chi }{\mathrm{10}}^{-\mathit{\eta }\mathit{\varphi }},\end{array}$$

where χ (W m⁻¹ K⁻¹) and η (no unit) are parameters accounting for particle shape effects, which are specified for gravel, fine mineral and organic soil (Côté and Konrad, 2005), and φ is porosity. Saturated thermal conductivity varies with water content and phase state according to the following:

$$\begin{array}{}\text{(7)}& {\mathit{\lambda }}_{\mathrm{sat}}=\left\{\begin{array}{l}{\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{1}-\mathit{\varphi }}{\mathit{\lambda }}_{\mathrm{liq}}^{\mathit{\varphi }}T\le T_{\mathrm{f}}\\ {\mathit{\lambda }}_{\mathrm{s}}^{\mathrm{1}-\mathit{\varphi }}{\mathit{\lambda }}_{\mathrm{ice}}^{\mathit{\varphi }}T>T_{\mathrm{f}}\end{array}\right.,\end{array}$$

where λ_liq, λ_ice, λ_s are thermal conductivities of liquid water, ice, and soil solid (W m⁻¹ K⁻¹), which are all constant values. T is soil temperature (^∘C) and T_f is the soil freezing point temperature that is assumed to be 0 ^∘C in DOS-TEM, which is consistent with most land surface models (e.g., Oleson et al., 2010).

2.3.2 Implementation of soil hydrological processes

Surface runoff, infiltration, and water redistribution among soil layers are simulated in a similar way to the Community Land Model 4 (Oleson et al., 2010). Soil matric potential (Ψ) determines the direction of water movement, and hydraulic conductivity describes the ease with which water can move through the soil.

$$\begin{array}{}\text{(8)}& \mathrm{\Psi }={\mathrm{\Psi }}_{\mathrm{sat}}({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{liq}}}{\mathit{\varphi }}}{)}^{-B},\end{array}$$

where Ψ_sat is the saturated soil matric potential (mm H₂O, hereafter mm), and B is the pore size distribution parameter. The soil hydraulic conductivity (K, mm s⁻¹) is a function of the saturated soil hydraulic conductivity (K_sat) as follows:

$$\begin{array}{}\text{(9)}& K=K_{\mathrm{sat}}({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{liq}}}{\mathit{\varphi }}}{)}^{\mathrm{2}B+\mathrm{3}}.\end{array}$$

Several important features relating to permafrost have been considered in the DOS-TEM (see Yi et al., 2014b), including runoff from a perched saturated zone or exchanges of water between the soil and a water reservoir. Runoff from a perched saturated zone above permafrost is implemented following Swenson et al. (2013):

$$\begin{array}{}\text{(10)}& Q_{\mathrm{perch}}=\mathit{\alpha }{k}_{\mathrm{p}}(z_{\mathrm{frost}}-z_{\mathrm{perched}})\sin \mathrm{?}\left({\displaystyle \frac{\mathrm{\Theta }}{\mathrm{180}}}\mathit{\pi }\right),\end{array}$$

where αis an adjustable parameter (0.6 m⁻¹), K_p is the mean saturated hydraulic conductivity within the perched saturated zone (mm s⁻¹), z_frost and z_perched are the depths to the permafrost table and the perched water table (m), and Θ is slope (^∘).

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, 2013, 2014a).

2.4 Model inputs and initialization

We used the monthly averaged air temperature, downward radiation, precipitation, and humidity as input to drive the DOS-TEM. Leaf area index (LAI), leaf area per unit ground surface area, was specified to be 0.6 m² m⁻² in July and August, 0.1 m² m⁻² in April and October, 0 m² m⁻² between November and March, 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). Its value is unchanged in each month.

Soil temperature and moisture were initialized at −1 ^∘C and saturation. The temperature gradient at the bottom of bedrock was set to be 0.06 ^∘C cm⁻¹ based on borehole observations. Volumetric unfrozen liquid water in winter was set to be 0.1 based on observations. Multi-year (2003–2012) mean monthly driving data were used to spin up the model for 100 years. In this way, suitable initial values of soil moisture, temperature, and rock temperature of each layer are generated before driving DOS-TEM with monthly data over the period of 2003–2012.

2.5 Sensitivity analyses

The soil textures on the QTP mainly consist of loam, sand, and coarse fragment soils (Wu and Nan, 2016). We used a uniform sand or loam soil profile to represent coarse and fine soil textures, respectively. Sands are the coarsest texture considered in most modeling studies (e.g., Oleson et al., 2010). Therefore, we used our measured parameters to substitute the parameters of sand and loam to investigate the effects of coarse-fragment soil parameters on permafrost dynamics. We first ran DOS-TEM using the default porosity, soil thermal conductivity (Eq. 5), hydraulic conductivity (Eq. 9), and matric potential schemes of these two default soil textures (Eq. 8). The default parameters φ, Ψ_sat, K_sat and B were calculated based on soil texture used in the Community Land Model 4 (Eqs. 8 and 9; Oleson et al., 2010). We then substituted the default values of φ, Ψ_sat, K_sat and B based on our laboratory measurements and calibration. Parameters Ψ_sat and B were fitted with measured matric potential data using Matlab Isqucurvefit tools. We did not calibrate soil thermal conductivity to retrieve parameters of Eqs. (6) and (7). Instead, we interpolated measured thermal conductivities over a range of degrees of saturation (0 to 1), which was used as a look-up table by the DOS-TEM. Therefore, our sensitivity analyses considered a set of four factors, i.e., porosity, matric potential (Ψ_sat and B), hydraulic conductivity (K_sat and B) and thermal conductivity. We also analyzed three different slope gradients (0, 5, and 10^∘) and three different soil thicknesses (3.25, 4.25, and 5.25 m) above 56 m of bed rock. There were 11 soil layers with the top nine layers being 0.05, 0.1, 0.1, 0.2, 0.2, 0.2, 0.3, 0.3, and 0.3 m thick. The thicknesses of the bottom two soil layers were 0.5 and 1 m, 0.5 and 2 m, and 1.5, and 2 m for the 3.25, 4.25, and 5.25 m soil thickness cases. There were six rock layers with thicknesses of 2, 2, 4, 8, 16, and 20 m. Since the site is on the top of an upland plain landform, we did not test the effects of aspect variation. We instead considered the effects of slope on surface runoff. In summary, our sensitivity analyses with the DOS-TEM involved 288 different combinations of parameter values.

We did not measure the heat capacity. The maximum and minimum heat capacities of mineral soil types considered in land surface model are 2.355 and 2.136 MJ m⁻³, giving a relative difference of less than 10 %. Therefore, in this study, we did not carry out sensitivity tests using thermal diffusivity (the ratio between thermal conductivity and heat capacity).

Table 1 The mean (standard deviation in brackets) of measured soil bulk density (ρ_b, g cm⁻³), porosity calculated from bulk density (φ_c, m³ mm⁻³), and measured porosity (φ_m, m³ m⁻³) of different layers based on soil samples in this study.

Table 2 The particle size diameter fractions (for >2 mm this is the mass ratio between soil particles greater than 2 mm and total soil sample, while for the other fractions this is the ratio between mass of the soil in the size range and the mass of all particles <2 mm) and soil texture (based on USDA classification) of different layers based on soil samples in this study.

3 Results

3.1 Soil physical properties

3.1.1 Soil porosity, particle size and bulk density

Results from laboratory analysis of the soil samples are shown in Tables 1 and 2. The mean mass ratio of the coarse-soil fraction (particle size diameter >2 mm) of different soil layers ranged from 0.38 to 0.65 with a mean of 0.55. According to the USDA classification system (clay is <2 µm, silt is 2–50 µm, in this study 2–63 µm, and sand is 50 µm −2.0 mm, but in this study was 63 µm −2.0 mm), the major soil texture of this site was loamy sand, with the exception of sandy loam at 20–30 cm depth. The default porosities of sand and loam were 37.3 % and 43.5 %. The φ_m of samples down to 2 m depth ranged from 21 % to 30 % with a mean of 27 %, and the mean ρ_b ranged from 1.61 to 1.86 g cm⁻³ with a mean of 1.74 g cm⁻³. The φ_c (Eq. 4) ranged from 29.8 % to 39.2 %. No significant relationships were found among φ_m, ρ_b, and the coarse-soil fraction (p>0.05).

Table 3 The mean (standard deviation in brackets) of the measured frozen and unfrozen dry and saturated soil thermal conductivity (W m⁻¹ K⁻¹) of different soil layers.

3.1.2 Thermal conductivity

The results of the thermal conductivity determinations are shown in Table 3. The unfrozen λ_dry of different soil layers ranged from 0.24 to 0.40 W m⁻¹ K⁻¹ with a mean of 0.36 W m⁻¹ K⁻¹, and the frozen λ_dry ranged from 0.25 to 0.41 W m⁻¹ K⁻¹ with a mean of 0.35 W m⁻¹ K⁻¹. The difference in λ_dry between frozen and unfrozen states was small. The unfrozen λ_sat of different soil layers ranged from 2.15 to 2.74 W m⁻¹ K⁻¹ with a mean of 2.48 W m⁻¹ K⁻¹. The frozen λ_sat ranged from 3.06 to 3.72 W m⁻¹ K⁻¹ with a mean of 3.33 W m⁻¹ K⁻¹. The difference in λ_sat between frozen and unfrozen states was about 0.85 W m⁻¹ K⁻¹. There existed a threshold of soil saturation (i.e., ∼0.28 m³ m⁻³), below which frozen soil thermal conductivity was slightly smaller than unfrozen soil (Fig. 4a).

Figure 4 The relationship between soil saturation (solid and dotted lines represent frozen and unfrozen cases) and soil thermal conductivity (λ, W m⁻¹ K⁻¹) from (a) measured values (measured; dots and empty diamonds represent measured frozen and unfrozen soil thermal conductivities), (b) using the Côté and Konrad (2005) scheme (CK) and (c) using the Farouki (1986) scheme (Farouki).

Results from determining thermal conductivities using the Côté and Konrad (2005) scheme are shown in Fig. 4b. The default frozen and unfrozen λ_dry for sand and loam were about 0.42 and 0.24 W m⁻¹ K⁻¹. The frozen and unfrozen λ_sat of sand were 3.11 and 1.90 W m⁻¹ K⁻¹. Those of loam were about 2.36 and 1.33 W m⁻¹ K⁻¹. Results from determining thermal conductivities using the Farouki (1986) scheme are shown in Fig. 4c. The default frozen and unfrozen λ_dry for sand and loam were about 0.97 and 0.63 W m⁻¹ K⁻¹. The frozen and unfrozen λ_sat of sand were 5.21 and 3.18 W m⁻¹ K⁻¹. Those of loam were about 4.49 and 2.52 W m⁻¹ K⁻¹.

Figure 5 The relations between (a) saturated hydraulic conductivity (K_sat, mm s⁻¹) and coarse-fragment fraction (solid dots represent measured values; empty circles and empty triangles represent the corresponding values of sand and loam used in the Community Land Model), and (b) soil saturation (m³ m⁻³, lines) and absolute value of matric potential (Ψ, mm H₂O) at three representative depths (solid and dashed lines represent the default values (Oleson et al., 2010) of sand and loam).

Table 4 The mean (standard deviation) of measured saturated hydraulic conductivity (K_sat; mm s⁻¹) and fitted absolute value of saturated matric potential (Ψ_sat; mm), fitted pore size distribution parameter (B), and the correlation coefficients (R²) between calculated matric potential using fitted and measured equations.

3.1.3 Saturated hydraulic conductivity

The mean K_sat of soil layers, shown in Table 4, ranged from 0.0036 to 0.0315 mm s⁻¹. The maximum K_sat was about 8.7 times larger than the minimum. The K_sat tended to be larger with an increasing proportion of coarse fragment in the soil samples (Fig. 5a), and was about 0.03–0.06 mm s⁻¹ for some samples with coarse fragment greater than 70 %. The default K_sat of sand and loam were 0.024 and 0.0042 mm s⁻¹.

3.1.4 Matric potential

The correlation coefficients between calculated and fitted Ψ, shown in Table 4, were all greater than 0.96. The mean absolute value of Ψ_sat of soil layers ranged from 14.47 to 603.7 mm, and those of B ranged from 1.89 to 5.22 (Table 4 and Fig. 5b). The default absolute values of Ψ_sat of sand and loam were 47.29 and 207.34 mm, and the B values were 3.39 and 5.77.

Figure 6 Comparisons of soil temperatures (T, ^∘C) simulated using default parameters for sand, loam, and our measured parameters (lines) with measured soil temperatures (dots) at 20, 100, 200, and 400 cm depths. Error bars show the standard deviations calculated based on nine simulations with three different slopes and three different soil thicknesses (measured porosities were used in the simulation).

3.2 Comparisons between simulations using default vs. measured parameters

3.2.1 Soil temperature

The mean root mean square 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^∘) were about 1.07 ^∘C at 20 cm (Fig. 6c). The mean RMSEs for all model runs with default sand and loam parameters were about 0.97 and 1.18 ^∘C. For other soil layers, the RMSEs of model runs with measured parameters were much smaller than those with default sand and loam parameters (Fig. 6d–l). The simulated soil temperatures using default sand and loam 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.53 ^∘C (Fig. 6e).

Figure 7 Comparisons of soil volumetric liquid water content (θ_liq, m³ m⁻³) simulated using default parameters sand, default loam, and measured parameters (lines) with measured soil moisture (dots) at 10, 40, 80, and 160 cm depths. Error bars showed the standard deviation calculated based on nine simulations with three different slopes and three different soil thicknesses (measured porosities were used in the simulation).

The standard deviations of soil temperatures among different slopes and soil thicknesses using measured parameters were larger than those using the default parameters (Fig. 6), and they increased from 0.40 ^∘C at 100 cm to 0.61 ^∘C at 200 cm (Fig. 6f and i). The standard deviations using default loam parameters were smaller (<0.15 ^∘C at all depths) than those using default sand parameters.

3.2.2 Soil liquid water

The mean RMSEs between monthly measured θ_liq and model simulations with measured parameters ranged from 0.03 to 0.09, which were smaller than RMSEs for sand and loam parameters (Fig. 7). The model simulations for loam parameters have larger RMSEs than those for sand parameters. θ_liq was always overestimated in warm seasons at depths of 10, 40, and 80 cm. θ_liq was underestimated at a depth of 160 cm, where the simulated soil was frozen. All model simulations overestimated θ_liq at 40 cm, where the maximum measured θ_liq were about 0.1 (Fig. 7d–f).

The standard deviations of θ_liq 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 θ_liq using loam parameters (<0.032) were less than those using measured parameters.

3.2.3 Active layer depth (ALD)

The mean RMSEs between measured ALDs (derived from linear interpolation of soil temperatures) and modeled ALDs (simulated explicitly) were about 1.06, 1.72, and 0.28 m for model runs with sand, loam, and measured parameters (Fig. 8a). The mean standard deviations were about 0.088, 0.026, and 0.28 m. All simulations using sand and loam parameters underestimated ALDs. When φ_m was replaced with φ_c, the mean RMSEs and standard deviations were about 0.55 and 0.12 m.

3.2.4 Permafrost lower boundary (PLB)

The mean RMSEs between measured PLBs (derived from linear interpolation of temperatures) and modeled PLBs (derived from linear interpolation of simulated bed rock temperatures) were about 10.25, 10.23, and 6.71 m for model runs with sand, loam, and measured parameters (Fig. 8b). The mean standard deviations were about 1.89, 1.51, and 6.62 m. All simulations using sand and loam parameters overestimated PLBs. When φ_m was replaced with φ_c, the mean RMSEs and standard deviations were about 4.78 and 2.82 m.

3.3 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 ^∘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.

Figure 8 Contour plots showing (a) soil temperature (^∘C) from borehole measurements down to 5 m superimposed with simulated active layer depths over the period of 2003–2011; and (b) ground temperature down to 50 m superimposed with the simulated permafrost low boundary. Black, blue, and magenta represent simulations with loam, sand, and measured parameters. Error bars show the standard deviation calculated based on nine simulations with three different slopes and three different soil thicknesses (measured porosities were used in the simulation; white zones in the contour plots indicate missing borehole data).

Figure 9 Root mean square errors between measurements and model simulations (with different combinations of measured porosity (I), thermal conductivity (II), hydraulic conductivity (III), and matric potential (IV) substituted for default sand parameters) for 20 and 100 cm soil temperatures (^∘C), active layer depth (ALD, m), and permafrost low boundary (PLB, m). O and All represent model runs without substitution of default parameters and with all four parameters substituted, respectively. Mean and standard deviation of model simulations with three different soil thicknesses at each slope (0, 5, and 10^∘) are shown.

Figure 10 Root mean square errors between measurements and model simulations (with different combinations of measured porosity (I), thermal conductivity (II), hydraulic conductivity (III), and matric potential (IV) substituted for default loam parameters) for 20 and 100 cm soil temperatures (^∘C), active layer depth (ALD, m), and permafrost low boundary (PLB, m). O and All represent model runs without substitution of default parameters and with all four parameters substituted. Mean and standard deviation of model simulations with three different soil thicknesses at each slope (0, 5, and 10^∘) are shown.

3.3.1 Effects of single-parameter sensitivity analyses

Porosity

Replacing default sand or loam porosity with φ_m changed mean RMSEs of soil temperatures (model runs with three different slopes and three different soil thicknesses at two different soil depths) from 1.18 or 1.84 ^∘C to 1.25 or 1.09 ^∘C (Figs. 9 and 10). Mean RMSEs of ALD were reduced from 1.06 or 1.72 m to 0.22 or 0.85 m. Mean RMSEs of PLB were changed from 10.26 or 10.24 m to 6.61 or 10.97 m. Mean RMSEs of θ_liq were reduced from 0.074 or 0.14 to 0.06 or 0.062 when φ_m was used for replacing default sand or loam porosity (Figs. 11 and 12).

Figure 11 Root mean square errors between measurements and model simulations (with different combinations of measured porosity (I), thermal conductivity (II), hydraulic conductivity (III), and matric potential (IV) substituted for default sand parameters) for 10, 40, 80, and 160 cm soil volumetric liquid water content. O and All represent model runs without substitution of default parameters and with all four parameters substituted. Mean and standard deviation of model simulations with three different soil thicknesses at each slope (0, 5, and 10^∘) are shown.

Figure 12 Root mean square errors between measurements and model simulations (with different combinations of measured porosity (I), thermal conductivity (II), hydraulic conductivity (III), and matric potential (IV) substituted for default loam parameters) for 10, 40, 80, and 160 cm soil volumetric liquid water content. O and All represent model runs without substitution of default parameters and with all four parameters substituted. Mean and standard deviation of model simulations with three different soil thicknesses at each slope (0, 5, and 10^∘) are shown.

Thermal conductivity

Replacing default sand or loam thermal conductivity with measured parameters reduced mean RMSEs of soil temperatures from 1.18 or 1.84 ^∘C to 1.02 or 1.15 ^∘C (Figs. 9 and 10). Mean RMSEs of ALD were reduced from 1.06 or 1.72 m to 0.56 or 1.04 m. Mean RMSEs of PLB were changed from 10.26 or 10.24 m to 4.18 or 1.27 m. Mean RMSEs of θ_liq changed very slightly (Figs. 11 and 12).

Hydraulic conductivity and matric potential

Replacing default sand or loam hydraulic conductivity with measured parameters had very small effects on mean RMSEs of soil temperatures and ALDs (Figs. 9 and 10). The same was true for matric potential. When hydraulic conductivity of default sand or loam was substituted, mean RMSEs of PLB decreased or increased, respectively. However, when matric potential was substituted, mean RMSEs of PLBs increased or decreased. When hydraulic conductivity or matric potential parameters were substituted in default sand or loam parameters, mean RMSEs of θ_liq changed slightly (Figs. 11 and 12).

Table 5 Model performance when default sand parameters are substituted with combinations of measured porosity (I), thermal conductivity (II), hydraulic conductivity (III), and matric potential (IV).

Best column shows the model simulations (individual parameter substitution) with the smallest root mean square error (RMSE) for 100 cm soil temperature (ST, ^∘C); active layer depth (ALD, m); permafrost low boundary (PLB, m); 10, 40, 80, and 160 cm soil liquid water content (SM, -). Numbers indicate the combination of parameters that have smaller RMSE than the best model run using individual parameter substitution. The column All indicates the combination of all four parameters. The smallest number indicates the smallest RMSE.

Table 6 Model performance when default loam parameters are substituted with combinations of measured porosity (I), thermal conductivity (II), hydraulic conductivity (III), and matric potential (IV).

Best column shows the model simulations (individual parameter substitution) with the smallest root mean square error (RMSE) for 100 cm soil temperature (ST, ^∘C); active layer depth (ALD, m); permafrost low boundary (PLB, m); 10, 40, 80, and 160 cm soil liquid water content (SM, -). Numbers indicate the combination of parameters that have smaller RMSE than the best model run using individual parameter substitution. The column All indicates the combination of all four parameters. The smallest number indicates the smallest RMSE.

3.3.2 Effects of combined parameters

We compared model simulations with different combinations of measured parameters (porosity, thermal conductivity, hydraulic conductivity, and matric potential) to those with one substituted measured parameter. We ranked those model runs with less RMSEs than the best of the model runs with one parameter substituted with a measurement-derived value (Tables 5 and 6). We did not consider the 10 cm soil temperature, which was similar across all model runs.

For sand, model simulations with porosity and thermal conductivity and/or hydraulic conductivity substituted had four outcomes with lower RMSEs (Table 5 and Figs. 9 and 11). Only two out of seven outcomes had lower RMSEs with all four parameters substituted. Among all the 18 cases with RMSEs less than the individual “best” RMSE, porosity was included 18 times, and thermal conductivity and hydraulic conductivity were included 10 times.

For loam, model simulations with porosity and thermal conductivity substituted had five outcomes with lower RMSEs (Table 6 and Figs. 10 and 12). Among all the 27 cases with RMSEs less than the individual best RMSE, porosity was included 27 times, thermal conductivity was included 16 times, and matric potential 14 times.

3.3.3 Effects of slope and soil thickness

Changes in slope alone had small effects on simulated soil temperatures and ALDs (Figs. 9 and 10). An increase in slope generally reduced RMSEs of θ_liq (Figs. 11 and 12). Model simulations with substituted porosity had smaller differences in θ_liq RMSE between different cases of slopes. For example, the mean RMSEs of model simulations with slopes of 0 or 5^∘ and sand parameters substituted with φ_m were 0.078 or 0.048, while those with porosity not substituted were 0.141 or 0.055. Similarly, the mean RMSEs of model simulations using default loam parameters with porosity substituted were 0.08 or 0.05 for slopes of 0 or 5^∘, respectively. The mean RMSEs were 0.18 or 0.1 with porosity not substituted. For a further increase in slope to 10^∘, changes of RMSEs of θ_liq at depths of 10–160 cm were small.

Soil thickness had small effects on 20 and 100 cm soil temperatures and 10–160 cm θ_liq, and it had prominent effects on PLB for a few cases only with a slope of 10^∘ (Figs. 9 and 10).

4 Discussion

4.1 Characteristics of soil physical properties

Although the effects of coarse-fragment soils on permafrost dynamics have been considered in a few modeling studies, the thermal and hydraulic properties of coarse-fragment soils were calculated without validation or calibration (Pan et al., 2017; Wu et al., 2018). To our knowledge, this is the first study measuring physical properties of coarse-fragment soil samples from the permafrost region of the QTP.

The weight fraction of coarse fragment (diameter >2 mm, including gravel) in the soil samples we analyzed 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) range from 5 % to 33 % and from 10 % to 28 % for the Madoi and Naqu sites. 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 coarse-fragment soils commonly exist in the QTP. The data set of Wu and Nan (2016) indicated that gravel content widely exists on the middle and western parts of the QTP. The saturated hydraulic conductivity and matric potential of soil samples measured in this study were more similar to sand than loam (see Sect. 3.1). It is consistent with the study of Wang et al. (2013) that coarse-soil material has a poor water-holding capability.

The measured thermal conductivities of saturated soil samples were relatively close to those estimated by the Côté and Konrad (2005) scheme. However, they were much less than those estimated by the Farouki scheme (Fig. 4). 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 φ_m 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 the 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 (Eqs. 6–9). It plays a more important role than other parameters in simulated soil thermal and hydrological dynamics (Tables 5 and 6; Figs. 9–12). It is noteworthy that it is easy and efficient to measure porosity.

4.2 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 loam parameters with high porosity and low saturated hydraulic conductivity were used, soil layers were almost saturated (Fig. 7). The simulated ALDs were about 1.58 m, which was less than half of measured ALDs (Fig. 8a). When the slope was 0^∘, subsurface runoff did not occur in the saturated zone above the bottom of the active layer. The simulated θ_liq was generally higher in the active layer. However, when the slope was 5^∘, the simulated θ_liq was lower and the RMSE was smaller (Figs. 11 and 12). These patterns were especially obvious when both porosity and saturated hydraulic conductivity were large (Eq. 10; Figs. 11 and 12). 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 (Figs. 9 and 10). The biggest effects were on PLB, which became manifest during long-term spin-up 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 thermal processes in permafrost 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, variation in soil water content in coarse-fragment soils strongly affects the thermal and hydrological properties; thus it is critical to simulate soil water dynamics to properly project permafrost dynamics in the future.

4.3 Limitations and outlook

4.3.1 Sampling and laboratory measurement

We used cut rings with 10 cm diameter to sample soil and weathered mudstones. However, it is very likely that there could have been much bigger coarse-fragment soils. 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 used the QL-30 thermophysical instrument (Anter Corporation, US) to measure thermal conductivity. It worked properly in an unfrozen condition. However, when frozen, the surface of the soil sample was usually uneven due to frost heave, which reduces the contact between the QL-30 plate and the soil sample surface. The measured frozen thermal conductivities were smaller than unfrozen thermal conductivity even for the case of saturation, which was definitely wrong. Thus we used the KD2 Pro to determine thermal conductivities. The second phenomenon was that there seems to be a threshold of soil saturation, below which unfrozen soil thermal conductivity is greater than frozen soil thermal conductivity (Fig. 4a). This pattern was somewhat exhibited in estimates of the Côté and Konrad (2005) scheme (Fig. 4b), but not in the estimates of the Farouki scheme (Fig. 4c). More measurements using instruments with higher accuracy should be made in the future.

The measured porosities are generally smaller than those calculated from bulk density. We made additional model simulations using porosities calculated from bulk density in combination with other measured parameters. Our results showed that the RMSEs of ALD and PLB were 0.55 and 4.78 m (figures not shown), whereas those calculated using φ_m were 0.28 and 6.71 m. There is a variety of methods for measuring soil porosity (Stephens et al., 1998). The method used in this study is widely used for its simplicity (e.g., Chen et al., 2012), and only requires measuring weights of samples under saturation and dry conditions (Eq. 2). Though soil samples were immersed in water under atmospheric pressure for 24 h to research saturation, it is possible that some air still remained in the soil after immersion but most of our soil samples contained coarse fragments and we assumed the volume of any remaining air to be negligible.

4.3.2 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 moisture at 40 cm depth was less than 0.1 m³ m⁻³. However, the simulated soil moisture were always much greater (Fig. 7f). There were also spikes in measured soil moisture at 80 and 160 cm depths, which were not presented in the simulation (Fig. 7i and l). In the DOS-TEM, the unfrozen soil water content, or super-cold water, was prescribed to be 0.1 m³ m⁻³. When soil is freezing, if the soil liquid water content is lower than this value, no phase change will happen (Fig. 7k). Therefore, model results would improve with the capability to simulate the dynamics of unfrozen soil water content (Romanovsky and Osterkamp, 2000).

The TEM family models use monthly atmospheric data to drive both site and regional applications. In this study, 30 min and daily driving data are available. Although it is possible to lose fidelity after daily interpolations, we still decided to use monthly driving data for the following reasons: (1) Zhuang et al. (2001) performed a test with daily and monthly driving data sets, and the results showed that the RMSEs of ALD were about 3 cm; and (2) we intend to apply the model over large regions where reliable daily data sets might not be available.

4.3.3 Regional applications

The coarse-fragment content of soil affects its 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). Thus soil texture plays an important role in permafrost dynamics (Fig. 8). The dominant soil textures on the QTP from Wu and Nan (2016) are loam, sand, and gravel. The specification of loam in simulations results in estimates of ALD that are much smaller than the measurements (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 coarse-fragment content (including gravel) and to develop regional data sets of soil texture for input.

Organic soil carbon content in mineral soil on the QTP affects soil porosity and thermal conductivity (Chen et al., 2012). However, in the site considered in this study, the amount of organic soil carbon in the soil was small (Fig. 2), and we did not explicitly consider the effects of organic soil carbon on soil properties. Alpine swamp meadow, alpine meadow, alpine steppe, and alpine desert are the major vegetation types on the QTP (Wang et al., 2016; see also Fig. 1b). Alpine swamp meadow and alpine meadow usually contain fine soil particles and high organic carbon density, while the other two types usually contain coarse-soil particle and low organic carbon density (Qin et al., 2015). More laboratory work is needed to develop proper schemes for representing mixed soil with fine mineral, coarse fragment (including 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 data set of soil texture is also required for regional projections of permafrost changes on the QTP. Currently, a preliminary data set for gravel exists (Wu and Nan, 2016), though gravel soil has only been mentioned in a few papers on the QTP (Chen et al., 2015; Wang et al., 2011; Yang et al., 2009). One way to improve the regional data set is to collect relevant data through extensive field campaigns (e.g., Li et al., 2015). Ground-penetrating radar is a feasible tool that retrieves soil thickness above the coarse-fragment soil layer (Han et al., 2016), and coarse-fragment soils can be identified in photos taken with unmanned aerial vehicles (Chen et al., 2017; Yi, 2017). In combination with ancillary data sets (e.g., geomorphology, topography, vegetation), it is possible to improve the accuracy of spatial data sets of soil texture on the QTP (Li et al., 2015; Wu et al., 2016). Another way is to retrieve soil physical properties using data assimilation technology, as done by Yang et al. (2016), who assimilated porosity using a land surface model and microwave data.

5 Conclusions

In this study, we excavated soil samples from a permafrost site on the central QTP and measured soil physical properties in the laboratory. Coarse fragments were common in the soil profile (up to 65 % of soil mass) and porosity was much smaller than the typical soil types used in land surface models. We then performed a sensitivity analysis of these parameters on soil thermal and hydrological processes within a terrestrial ecosystem model. When default sand or loam parameters were substituted with measured soil properties, the model errors of active layer depth were reduced by 74 % or 84 %, whereas those of the low permafrost boundary were reduced by 35 % or 34 %. Our sensitivity analyses showed that porosity played a more important role in reducing model errors than the 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.