Analysis of soil hydraulic and thermal properties for land surface modeling over the Tibetan Plateau

Soil information (e.g., soil texture and porosity) from existing soil datasets over the Tibetan Plateau (TP) is claimed to be inadequate and even inaccurate for determining soil hydraulic properties (SHP) and soil thermal properties (STP), hampering the understanding of the land surface process over TP. As the soil varies across three dominant climate zones (i.e., arid, semi-arid and subhumid) over the TP, the associated SHP and STP are expected to vary correspondingly. To obtain an explicit insight into the soil hydrothermal properties over the TP, in situ and laboratory measurements of over 30 soil property profiles were obtained across the climate zones. Results show that porosity and SHP and STP differ across the climate zones and strongly depend on soil texture. In particular, it is proposed that gravel impact on porosity and SHP and STP are both considered in the arid zone and in deep layers of the semi-arid zone. Parameterization schemes for porosity, SHP and STP are investigated and compared with measurements taken. To determine the SHP, including soil water retention curves (SWRCs) and hydraulic conductivities, the pedotransfer functions (PTFs) developed by Cosby et al. (1984) (for the Clapp–Hornberger model) and the continuous PTFs given by Wösten et al. (1999) (for the Van Genuchten–Mualem model) are recommended. The STP parameterization scheme proposed by Farouki (1981) based on the model of De Vries (1963) performed better across the TP than other schemes. Using the parameterization schemes mentioned above, the uncertainties of five existing regional and global soil datasets and their derived SHP and STP over the TP are quantified through comparison with in situ and laboratory measurements. The measured soil physical properties dataset is available at https://data.4tu.nl/repository/uuid:c712717c-6ac0-47ff-9d58-97f88082ddc0.


Introduction
As the highest plateau in the world, the Tibetan Plateau (TP) exerts a significant influence on the Earth's climate system and plays a prominent role in the evolution of the Asian monsoon system (Yao et al., 2012;Qiu, 2008;Ma et al., 2017;Kang et al., 2010).Studying this influence can advance our understanding of climate change (Ma et al., 2017).Soil moisture (hereafter referred to as SM) -one of the lower boundary conditions of the atmosphere -is a crucial land surface state (Koster et al., 2004) and therefore of high interest to investigate the land-atmosphere interactions, reflecting the trend and the variability of feedback between the water cycle and climate over the TP (Su et al., 2013(Su et al., , 2011)).
Accurate SM information is claimed as a necessity for improving precipitation and hydrology forecasts (Drusch, 2007;Dirmeyer, 2000;Robinson et al., 2008), especially on the TP, where it undergoes evident climate change (Ma et al., 2017;Douville et al., 2001;Yang et al., 2014Yang et al., , 2011)).Consistent spatial-temporal SM data can be obtained by using land surface models (LSMs) assimilating in situ and satellite observations.In these models, the specification of soil hydraulic properties (SHP) (i.e., soil water retention curve, SWRC; hydraulic conductivities) and soil thermal properties (STP) (i.e., thermal conductivities and heat capacity) are considered more decisive for SM simulation than atmospheric forcing and land surface characteristics (Shellito et al., 2016; Published by Copernicus Publications.Livneh et al., 2015;Kishné et al., 2017;Gutmann and Small, 2005), because SHP govern the partitioning of SM between infiltration and evaporation flux and STP regulate water heat transport processes (Zeng et al., 2009a, b;Garcia Gonzalez et al., 2012).
In situ measurements of basic soil properties and SHP and STP are crucial for soil moisture and heat flux simulation by LSMs.LSMs frequently use the Clapp and Hornberger (1978) model and the Van Genuchten (1980) and Mualem (1976) model for SHP, as well as the Farouki (1981) and Johansen (1975) schemes for STP.Since direct measurements of SHP and STP are always too time, labor and cost consuming, pedotransfer functions (PTFs) (Bouma, 1989;Van Looy et al., 2017) using basic soil property information have been developed to estimate parameters in the above SHP and STP schemes.Examples are the Cosby et al. (1984) PTF (e.g., based on sand fraction) for the Clapp and Hornberger (1978) (CH) scheme estimate in the Noah and Community Land Model (Chen and Dudhia, 2001;Oleson et al., 2008) and the soil class PTF (e.g., based on soil texture types) for the Van Genuchten (1980) and Mualem (1976) (VG) scheme (Balsamo et al., 2009) in the Hydrology Tiled European Centre for Medium-Range Weather Forecasts (ECMWF) Scheme for Surface Exchanges over Land (H-TESSEL).However, these aforementioned PTFs sometimes cannot predict well the SHP and STP, especially when soils contain organic matter or gravels (particle diameter ≥ 2 mm), because gravels and organic matter have different hydraulic and thermal properties than other fine mineral soils, and this suggests the need to obtain comprehensive soil property information (e.g., not only soil texture and porosity information, but also soil organic matter content, SOC, and gravel fraction).
Furthermore, studies using information on state variables (e.g., near-surface soil moisture or brightness temperature) can retrieve effective SHP and STP directly or indirectly through PTFs and LSMs (Ines and Mohanty, 2008a;Han et al., 2014;Dimitrov et al., 2014Dimitrov et al., , 2015;;Yang et al., 2016).Nevertheless, most of such retrievals only focus on the basic surface soil properties and SHP, based on the assumption of a homogenous soil column.If the system is highly heterogeneous (e.g., along the vertical profile), retrieval may be problematic (Ines and Mohanty, 2008b) and in situ measurements of soil property profiles can shed light on the soil property retrieval of the vertical profile.
Many global and local efforts have been made to compile and develop soil databases, but uncertainties in soil datasets might also cause bias in predicting SHP and STP and hence introduce uncertainties in representing the land surface states by LSMs.It has been reported that the overestimations of ECMWF SM analyses in the central TP could be partly attributed to the unrepresentative soil information from the Food and Agriculture Organization (FAO) Digital Soil Map (2003) as used in the H-TESSEL (Su et al., 2013).Currently, there is only soil texture information and few soil organic matter content profiles available over the TP when extracted from the published global in situ soil profiles (Batjes et al., 2017).The profiles of other vital soil properties, such as dry bulk density (BD) and porosity, are not provided (e.g., no in situ BD or porosity profiles available).Moreover, there are no comprehensive in situ measurements of basic soil properties, SHP and STP for land surface modeling over the TP.
In this study, we implemented the in situ and laboratory measurements of soil physical property profiles across the three climate zones of the TP and compiled the Tibet-Obs soil properties dataset.Based on the dataset, the variations in basic soil properties and SHP and STP across the three climate zones were investigated.Applications of the Tibet-Obs dataset were demonstrated through two cases: (1) we discussed the appropriate parameterization schemes of porosity and SHP and STP for their applicability in land surface modeling over the TP and (2) we evaluated the uncertainties of the five existing regional and global soil datasets and their derived SHP and STP over the TP.In Sect. 2 of this paper, the field campaign and laboratory experiments are described as well as the parameterization schemes for porosity and SHP and STP estimates.The specification of the Tibet-Obs dataset with data availability is documented in Sect.3. Results on the application of this dataset are presented in Sect. 4. Conclusions are presented in Sect. 5.This paper is expected to contribute to land surface modeling and hydro-climatology communities for their studies of the third pole environment, as well as to soil community in terms of filling geographic gaps of the published existing global soil databases.

Field experiments
Soils show spatial variation over the TP due to the varying soil formation factors (e.g., climate and parent material).The TP can be categorized into three dominant climatic zones: an arid zone (0.03 < aridity index (AI) < 0.2), a semi-arid zone (0.2 < AI < 0.5) and a subhumid zone (0.5 < AI < 1.0), according to the Food and Agriculture Organization (FAO) aridity index map (Fig. 1a) (Zeng et al., 2016).The Tibetan Plateau observatory stations of plateau-scale soil moisture and soil temperature (Tibet-Obs) (Su et al., 2011) are distributed throughout these climatic zones: (1) the Ngari network in the arid zone, located in the western part of the TP with the elevation ranging between 4200 and 6300 m above mean sea level (a.s.l.), where the annual mean temperature is 1.01 • C and the annual mean precipitation amount is 66.4 mm, the land cover is a typical desert environment dominated by bare soil surrounded with spare grass and soils are prevailed by sandy soils mixed with gravel (Fig. 1b); (2) the Naqu network in the semi-arid zone, located in a flat terrain with rolling hills at an average elevation of 4500 m a.s.l., where the annual mean temperature is −0.6 • C and the an- nual mean precipitation amount is 482 mm, the land cover is characterized as grasslands consisting of prairie grasses and mosses, and soils are dominated by loamy sand with organic matter and gravels (Fig. 1c); and (3) the Maqu network in the subhumid zone, located at the northeastern edge of the TP with elevations ranging between 3430 and 3750 m, where the annual mean temperature is 1.8 • C and the precipitation is 600 mm annually with more than 70 % falling during the monsoon season (e.g., from June until September).The land cover is dominated by short grassland, and soils are dominated by fine minerals with large silt proportions (Fig. 1d).Of these, the Naqu network is collocated with the multi-scale Soil Moisture and Temperature Monitoring Network on the central Tibetan Plateau (CTP-SMTMN) (Yang et al., 2013).
A field experiment was carried out across the TP in August 2016, taking soil core samples and measuring field saturated hydraulic conductivity (K s ) at various soil depths (Table 1, Fig. 1).Soils were vertically sampled using sample rings and augers (Eijkelkamp Soil & Water Company) in the vicinity of existing soil moisture and soil temperature stations of the Tibet-Obs (Su et al. 2011).Table 1 lists the specific sampling approach: (1) the soil was sampled (ca.200 g) with a plastic bag used to measure gravel content, soil texture and soil organic matter content; (2) the soil was sampled with standard sample rings (5 cm in height, 100 cm 3 in volume) for the determination of dry bulk density, porosity and thermal conductivity (λ); (3) for deriving the soil water retention curve, a dedicated small sample ring (1 cm in height, 20 cm 3 in volume) was used; (4) the in situ K s was measured using the Aardvark permeameter (2840 operating instructions -Eijkelkamp), a fully automated constant-head borehole permeameter.The Reynolds and Elrick solution aided with soil texture-structure category information (Elrick et al., 1989) was chosen for calculating K s .
Within the Maqu network, soil samples were collected at eight stations, located in areas to the east, west and southeast of the ELBARA-III radiometer location as well as in the southwest corner of the Maqu network (Fig. 1a).The K s was measured at three locations near the ELBARA station and at one location (CST05-near) in the southwest corner.Within the Naqu network, soil samples were taken at eight sites along the southwest branch of the CTP-SMTMN network, and the K s was measured at seven sites at BJ, Naqu_west, NQ01-04 and MS3608 (Fig. 1b).Within the Ngari network, soils were sampled at 14 stations (Fig. 1c).Eight sites at Ali02, SQ03, SQ07, SQ10, SQ17, SQ18, SQ20 and SQ21 were chosen for K s measurement.In total, 155 soils samples were taken and loaded into plastic bags, 101 samples were collected in standard rings and another 96 samples in small rings.Due to the remoteness and harsh environment on the TP, the locations chosen for soil sampling and fieldwork needed to take practical considerations into account -for example, (1) the location should be accessible by track, local road or national road and (2) the surrounding area should be flat enough to be representative of the local area.

Laboratory experiments
Three categories of soil samples were handled.From the 155 samples (59 from Ngari, 45 from Naqu and 51 from Maqu) in plastic bags, the soils were first separated into gravels and fine minerals (size < 2 mm) by using a sieve of 2 mm diameter mesh and weighed separately to obtain the gravimetric gravel fraction (GGF).Sand (0.05 mm < size < 2 mm), silt (0.002 mm < size < 0.05 mm), and clay (size < 0.002 mm) percentages as well as the mean particle diameter of the fine minerals (FD) were determined with the Malvern Mastersizer 2000 particle size analyzer (http://www.malvern.com), and the SOC was determined by the total organic content analytical instrument Multi N/C 3100 (http://www.analytik-jena.de/).For gravels, a set of sieves with diameters of 2, 2. 5, 4, 5, 7, 10, 16, 20, 25, 31.5, 40 and 50 mm were used to obtain their particle size distribution and the mean particle diameter of gravels (GD).
The 101 undisturbed soil samples (35 from Ngari, 21 from Naqu and 45 from Maqu) in standard sample rings were saturated and then dried in the oven (105 • C) for 24 h.The difference between wet and dry weight with known volume was used to calculate porosity and BD.The KD2 Pro thermal property analyzer connected to an SH-1 sensor (Decagon Devices) was used to measure heat capacity C s and thermal conductivity λ, while the soil was drying, providing drying C s -SM and λ-SM curves.
The 96 samples in small rings were intended for use in the SWRC experiment by using the pressure-cell method, but to complete this entire task it was considered too time and labor consuming.Therefore, instead of utilizing all soil samples, only 30 out of 96 samples were used for the E-east, E-west, E-southwest, CST05-near, NST30 and NST33 sites in the Maqu network.As the structure of the samples at the Naqu and Ngari networks was so unconsolidated that the material did not remain enclosed within the rings, only 25 undisturbed samples contained in standard rings were used from Naqu_north, SQ17, SQ18 and SQ21 sites.
The quality of the measured soil property dataset was evaluated based on quality indicators (e.g., observation date, level of trust, data quality rating and accuracy) from the World Soil Information Service (WoSIS) institute (Ribeiro et al., 2015).These four indicators provide measures that allow investigators to recognize factors that may compromise the quality of certain data and hence their suitability for use (Ribeiro et al., 2015).The results show that the dataset is of trust level "C", which means the highest level of the subjective measure inferred from soil expert knowledge.The entered data (level "A") have been standardized (level "B"); i.e., data numbers were correspondingly aligned with measured soil properties involved in the GlobalSoilMap specifications (McMillan, 2009) and with the measurement method and unit (see above paragraphs in Sect.2.2).The level B dataset was further harmonized (level C) to be sorted in the reference table (Ribeiro et al., 2015).For instance, tables for the profile data (see the raw data in the data repository, Zhao et al., 2018) described a soil profile and its attributes (e.g., land cover, position), as well as its constituent layers with their respective soil properties.These collated raw data included error-checking for possible inconsistencies.Furthermore, the values of the measured soil properties and SHP and STP were compared to those available in the literature to cross-check whether they were within a reasonable range.
The collected basic soil properties and the SHP and STP datasets named the Tibet-Obs dataset will be further used to evaluate the existing soil datasets of the FAO-UNESCO Soil Map of the World (2007) (hereafter referred to as FAO-UNESCO), the Harmonized World Soil Database (hereafter referred to as HWSD) (FAO/IIASA/ISRIC/ISSCAS/JR, 2012), a Chinese dataset of soil properties (Shangguan et al., 2012(Shangguan et al., , 2013) ) and soil hydraulic parameters using PTFs (Dai et al., 2013) released by Beijing Normal University (hereafter referred to as BNU), SoilGrids1km (Hengl et al., 2014) and the updated version of SoilGrids250m (Hengl, 2017) released by the International Soil Reference and Information Center (ISRIC) -WoSIS institute, and the hydraulic parameters based on SoilGrids1km and Schaap et al. (2001) PTFs (hereafter referred to as HPSS) (Montzka et al., 2017).The description of the existing datasets is listed in Table S1 of the Supplement.All datasets were linearly interpolated to match the measured dataset at specific depths, to ensure the comparability (inter-comparability).

Parameterization schemes
Many basic-soil-property-dependent schemes have been proposed for porosity estimation.The Cosby et al. (1984) univariate PTF that uses sand percentage (hereafter the Cosby-S scheme, see Eq. (A1) in the Appendix) has been widely used, and it should be noted that there is a multivariate PTF (Cosby et al., 1984) that uses clay as well as sand (Van Looy et al., 2017).Porosity can be inversely related to soil dry bulk density (Hillel, 2003) and calculated from in situ BD (hereafter the BD scheme, see Eq. A2).In most cases, these schemes perform well.However, with SOC in soils, soil porosity tends to increase.Another factor affecting porosity is the gravel content.As gravel content increases, the porosity tends to decrease.Chen et al. (2012) parameterized the impact of SOC and gravel content into a porosity estimation scheme (hereafter the SocVg scheme, see Eqs. (A3)-(A6)).Zhang et al. (2011) proposed a mixing-coefficient model to calcu- late the porosity of binary mixture (BM) made of a coarse (gravel) and a fine component over a range of gravel content (hereafter the BM scheme, see Eqs.A7-A10).In this study, as Fig. 2 shows, the Cosby-S, BD, SocVg and BM schemes were evaluated for their applicability over the three climate zones.
For the SHP estimate, we selected the Clapp and Hornberger (1978) (hereafter CH) and the Van Genuchten (1980) and Mualem (1976) (hereafter VG) schemes.Based on measured SWRCs, we used the scaling method (see Eq. A15) (Montzka et al., 2017) to determine the three hydraulic parameters involving saturated soil moisture (θ s ), soil water potential at air entry (ϕ s ), and an empirical parameter related to the pore-size distribution of the soil matrix (b) in the CH equation (see Eqs. A11-A12), as well as four hydraulic parameters involving θ s , residual soil moisture (θ r ), a parameter corresponding approximately to the inverse of the airentry value (α), and a shape parameter (n) in the VG model (see Eqs. A13-A14).The field capacity (FC) and the permanent wilting point (PWP) (regarded as the SM at about −33 and −1500 kPa of matric pressure, respectively) were also derived as they are the main parameters for soil water budget.Furthermore, the selected PTFs (see Table A1 in the Appendix) were used to estimate hydraulic parameters of SWRCs-CH and SWRCs-VG.Given that a good θ s estimate will improve SWRCs prediction, the optimal porosity scheme will be pre-selected for predicting SWRCs-CH and SWRCs-VG.Estimated SWRCs from PTFs were further compared with the measurement-determined SWRCs to indicate the uncertainty of using different PTFs.Saturated hydraulic conductivity, K s , combined with SWRCs-CH or SWRCs-VG is used to calculate the unsaturated hydraulic conductivity (K) and diffusivity (D).The PTFs used for SWRCs-CH and SWRCs-VG estimations also have corresponding equations (see footnotes in Table A1, Appendix) to predict K s , and most PTFs were developed based on fine minerals.To estimate the K s of a mixture containing gravels, Peck and Watson (1979) used a heat-flow analogy correlating the K s of the mixture with the K s of both fine minerals and the volumetric gravel fraction (VGF) (hereafter the PTFs-VGF scheme, see Eq. A16).This PTFs-VGF scheme can be applied to soils with low gravel content (Zhang et al., 2011).It is noted that the PTFs-VGF scheme needs an input (K sat,f , see Appendix A3) from the PTFs K s estimation.Furthermore, Koltermann and Gorelick (1995) used the Kozeny-Carman equation to estimate the hydraulic conductivity for binary mixtures (BM-KC), and the suitable grain diameter estimation was declared important (Kamann et al., 2007).To improve the performance of the Kozeny-Carman equation, Zhang et al. (2011) introduced the BM scheme for estimating porosity and a poweraveraging method for calculating representative grain diameter (hereafter the BM-KC scheme, see Eqs.A17-A19).In this study, the standard PTFs (see Table A1 in the Appendix), PTFs-VGF and BM-KC schemes were employed as shown in Fig. 2.
Several (semi-) empirical models have been developed to estimate the soil thermal conductivity λ.De Vries (1963) developed a Maxwell equation analogous physics-based model to describe λ (see Eq. A22).This model can predict λ accurately, although this is complicated by the fact that at least five soil mineral components and their separate shape features need to be taken into account (Tarnawski and Wagner, 1992).Furthermore, the effect of vapor movement caused by the temperature gradient is also parameterized in the De Vries (1963) model.It should be noted that the consideration of soil vapor flow is critical to accurately investigate the simultaneous transfer of moisture and heat, particularly in semi-arid and arid environments (Zeng et al., 2011a, b;Zeng and Su, 2013).Farouki (1981) proposed an alternative method and regarded liquid water as the continuous medium and soil minerals as uniform particles in the De Vries (1963) model.In this model, the λ of soil minerals was estimated by using a geometric mean equation from the quartz content in soil minerals and the λ of quartz and other soil minerals (see Eq. A23).The λ of vapor together with the shape factor for air pore were calculated in terms of water content and porosity (see Eqs. A24-A25) (hereafter the D63F scheme).Tian et al. (2016) developed a simple and generalized De Vriesbased model, which assumed that the λ and shape features of soil minerals were determined by soil texture (sand, clay and silt) and that the effect of vapor movement was negligible (hereafter the T16 scheme, see Eqs.A26-A29).The empirical model proposed by Johansen (1975) used the Kersten (1949) number and λ in dry and saturated conditions to estimate λ (hereafter the J75 scheme, see Eqs.A30-A35).In this study, as shown in Fig. 2, the D63F, T16 and J75 schemes were adopted.For each λ scheme, a comparison was made using parameters (i.e., the λ of soil minerals) with (see Eq. A34) and without (see Eq. A23) gravel and SOC considerations.The De Vries (1963) model was used for calculating C s (see Eqs. A20-A21).The details of porosity and the SHP and STP schemes are listed in the Appendix A1-A5.
3 The Tibet-Obs dataset

Data availability
The soil physical dataset is available at the 4TU.ResearchData data center at https://data.4tu.nl/repository/uuid:c712717c-6ac0-47ff-9d58-97f88082ddc0 (Zhao et al., 2018).The data are stored in .XLSX files.A readme file describes the structure of the Excel files, the measurement devices and contact information.The download linkages of existing soil property datasets used in this paper are included in the .txtfile.The location of sampling is stored in the .kmzfile.The raw data for each sampling site are also provided.

Soil texture
Figure 3 shows the mean of sand, clay and silt percentages, gravimetric gravel fraction, soil organic matter content, and the mean diameter of fine components and gravels at different depths across the three climate zones over the TP.In the Ngari network under the arid zone (Fig. 3a), the mean sand content was around 80 %, with higher values at surface layers of 5 and 10 cm than at deep layers.Silt and gravel contents ranged between 10 and 20 % and the percentages increased with depth.Clay content and SOC were 3 and 0.8 %, respectively, and remained constant along the profile.The FD and GD ranged from 0.19 to 0.24 and 4 to 8 mm, respectively, and showed a tendency to increase from the top to a depth of 20 cm, but to decrease in the deeper layers.It can be concluded that soil texture in the arid zone consists of a high proportion of coarse sand accompanied by gravel, and that the gravel content increases until 20 cm and then decreases slightly in the deeper layer.
In the Naqu network under the semi-arid zone (Fig. 3b), the mean sand fraction ranged from 70 to 80 %, with a slight decrease with depth.The silt and clay contents ranged 15-25 and 4-8 %, respectively, and increased with depth.The GGF exceeded 50 % for soils at depths of 40 and 50 cm, while it was much lower at the shallow layers.Mean FD and GD ranged 0.18-0.22 and 4-8 mm, respectively.GD at deep layers was larger than that at shallow layers.SOC approached 10 % in the surface layers but quickly declined at deep layers.It can be summed up that soil texture in the semi-arid zone is dominated by a high percentage of sand mixed with Bottom panel: variations in GD and FD at different depths.GGF is the gravimetric gravel fraction.SOC is the soil organic matter content.FD is the mean particle diameter of fine minerals.GD is the mean particle diameter of gravels.a small proportion of gravels, but with high SOC in shallow layers, and mainly mixed with big gravels at deep layers.
In the Maqu network under the subhumid zone (Fig. 3c), mean silt and clay contents were around 60 and 10 %, respectively, with a smoothly decreasing trend along the profile.The mean sand fraction ranged from 28 to 40 % and increased with depth.No gravel was found.Mean FD ranged from 0.024 to 0.036 mm, and fine soil mineral particles in deep layers (40 and 80 cm) were larger than in shallow layers (see Fig. 3c, lower panel).Similar to the SOC profile distribution in the Naqu network, the SOC was almost 20 % in surface soil layers and declined to 2.8 % at 80 cm.Soil texture in the subhumid zone is characterized as being dominated by a high percentage of silt content with relatively large SOC in the shallow layers and with mainly fine sand mixed in the deep layers.

Dry bulk density and porosity
In the Ngari network under the arid zone (Fig. 4a), the BD varied slightly (between 1.55 and 1.65 g cm −3 ) with depth, showing a peak at 10 cm.The porosity of the surface layer was slightly higher than in deep layers, with a mean profile porosity of 0.33.The porosity at 20 cm was the lowest in the profile, which might be caused by this layer containing the greatest proportion of gravel as well as the greatest GD and FD (see Fig. 3a).In the Naqu network under the semi-arid zone (Fig. 4b), the BD increased continuously with depth, with a minimum of 1 g cm −3 in the top layer and a maximum of 2.1 g cm −3 in the bottom layer.The porosity peaked at around 0.6 at the top layer, while monotonously decreasing to 0.25 at the bottom layer.Combined with the soil texture analysis (see Fig. 3b), variations of BD and porosity in the profile were inferred relevant to the high SOC in the surface layer and the large gravel content in the bottom layer.In the Maqu network under the subhumid zone (Fig. 4c), the BD ranged from 0.8 to 1.5 g cm −3 and increased with depth, while porosity decreased with depth and ranged from 0.72 to 0.45.The profile pattern of BD and porosity might be induced by SOC layering in the surface and soil texture fraction variations in the deep layers as Fig. 3c reveals.In summary, profiles of BD and porosity differ with soil texture variation over the three climate zones, and both the SOC and gravels affect the porosity.Overall porosity at shallow layers (5, 10 and 20 cm) increases from the arid to the semi-arid and then to the subhumid zones, while at deep layers (>= 40 cm) it shows an increase from the semi-arid to the arid and then the subhumid zones.

Soil water retention curve (SWRC) and saturated hydraulic conductivity (K s )
Figure 5 shows the pressure-cell measured SWRCs (markers in the figure) differed across the three climate zones.In the Ngari network under the arid zone, soil water retention quickly reduced as the suction slightly increased (Fig. 5a).
The same situation occurs for the deep layers in the Naqu network under the semi-arid zone (Fig. 5b).In the Maqu network under the subhumid zone, soil water retention was high and gradually decreased as the suction increased (Fig. 5c).2.
In the Ngari network under the arid zone (Fig. 6a), the magnitude of mean K s was on the order of 10 −5 (m s −1 ).The K s at 20 cm was lower than at other depths, which might be due to the lowest values of porosity in this layer (see Fig. 4a).In the Naqu network under the semi-arid zone (Fig. 6b), the mean K s exhibited a variation of 1 order of magnitude with depth, namely 10 −6 (m s −1 ) at depths of 10, 20 and 50 cm Earth Syst.Sci.Data, 10, 1031Data, 10, -1061Data, 10, , 2018 www.earth-syst-sci-data.net/10/1031/2018/ Figure 6.Profiles of mean saturated hydraulic conductivity (K s ) at three different climate zones.and 10 −5 (m s −1 ) at a depth of 40 cm.In the Maqu network under the subhumid zone (Fig. 6c), K s also differed by 1 order of magnitude: 10 −6 (m s −1 ) at depths of 5, 10, 20 and 80 cm and 10 −7 (m s −1 ) at a depth of 40 cm.It should be noted that the K s profiles of both the semi-arid and subhumid zones presented a lower K s in shallow layers than in the deeper layer.This is mainly due to the negative correlation between saturated hydraulic conductivity and soil organic carbon in soils where hydrophobic functional groups might dominate with organic carbon composition and reduce soil wettability (Nemes et al., 2005;Wang et al., 2009;Ellerbrock et al., 2005).As can be seen, K s varies with soil texture over the three climatic zones, and both SOC and gravels have an effect.At a certain depth, where the soil basic properties undergo a transition (see Fig. 3), the K s reaches a minimum.The mean and the standard deviation of the soil properties of the profiles in the three climate zones are listed in the Supplement (Tables S2-S4).

Gravel impact on porosity and K s
Figure 7a and b show that porosity did not change with GGF increasing within 0.3 in shallow layers, while with a GGF > 0.4, porosity tended to decline with increasing GGF, especially in deep layers.For example, porosities for layers with a GGF of 0.6 and 0.72 at 20 and 40 cm depths were lower than those with a GGF < 0.3 at 5 and 10 cm depths (Fig. 7a).With more gravels embedded in the matrix, the flow paths in the soil would become blocked and the porosity reduced (Zhang et al., 2011).However, the porosity did not always decrease as the GGF increased.Porosity for the layer with a GGF of 0.84 in the semi-arid zone was higher than porosities with a GGF ranging between 0.4 and 0.6 at 50 cm depth (Fig. 7b).Porosity for the layer with a GGF of 0.7 at 20 cm depth in the arid zone was also higher than that with a GGF of 0.6 at 40 cm depth (Fig. 7a).Porosity tended to increase as the GGF continued to increase because, when a GGF is relatively high (> 1 -porosity of gravels), connected pores can form among the gravels and thus increase porosity (Zhang et al., 2011).Figure 7c and d show a slight decrease in K s at 10 cm with a GGF < 0.62 and a slight increase in K s at 20 and 40 cm with a GGF > 0.8, which is consistent with the changes in porosity.The observations clearly show that gravels have a distinct impact on the porosity and K s in the arid and semiarid zones.It should be noted that, although the in situ K s measurements were conducted at locations adjacent to the places where we took soil samples, heterogeneity may have had an effect on the values of soil properties and parameters throughout our sampling procedures, as with any soil field experimentation.Nevertheless, the current findings based on field experiments are in line with reported findings based on laboratory experiments (Zhang et al., 2011;Koltermann, 1995;Sakaki and Smits, 2015).For these samples, it is not easy to vertically insert needles of the KD2 Pro device (see Sect. 2.2); instead, the needles were buried with soils in the surface of the sample as the alternative for the measurement.Additionally, the KD2 needles might experience a slight skew when they touch the hard gravel.All these factors might cause the fluctuation of C s when SM increases, while the overall rising trend is still observed in Fig. 8a and b. Figure 8c shows C s almost steadily increased with SM for the Maqu network under the subhumid zone.Samples from the subhumid zone are all fine-grained soils and make the needle of the KD2 device insert easily and thus form a steady environment for the measurement.Figure 8a-c show that C s -SM varied slightly at different depths in the three climate zones.The C s ranged from 1 at oven dry state to 2.5 MJ m −3 K −1 as the soil reached saturation over the arid zone, 0.5 to 3 MJ m −3 K −1 over the semi-arid zone, and 0.5 to 2.4 MJ m −3 K −1 over the subhumid zone.Figure 8d-f show how the relationship of λ-SM varied with depth.For the arid zone (Fig. 8d), the λ-SM curves were very similar at each depth due to the nearly homogenous sandy soils across the whole profile (see Fig. 3a).The mean λ ranged from 1.8 at saturation and 0.2 (W m −1 K −1 ) as the soils reached oven-dry state.In the semi-arid zone (Fig. 8e), the λ-SM curves were stratified, and soils with gravels in deep layers (see Fig. 3b) clearly had a higher λ (> 2 W m −1 K −1 ) than in other layers and other climate zones.In the subhumid zone (Fig. 8f), the λ-SM curves also presented variations with depth, though within a much narrower range than in the semi-arid zone.Such variation is mainly caused by the sand distribution along the profile, which increased slightly with depth (see Fig. 3c).The mean λ in the subhumid zone ranged from 1.6 to 0.2 (W m −1 K −1 ) as soils dried out.Furthermore, the surface layers in the semiarid and subhumid zones had lower λ values (Fig. 8e and f) because of the SOC influence.

Porosity estimation
Using basic soil properties data of Tibet-Obs with four schemes, porosities were estimated.Comparisons against measured porosities (see Table A2 in the Appendix) indicate that the BD scheme performs the best for estimating porosity in profiles across the three climate zones, as it is a bulk estimation scheme that takes both gravel and fine minerals into consideration.The Cosby-S scheme overestimates porosities over the arid zone and provides constant porosity values over the semi-arid and subhumid zones.The SocVg scheme also overestimates porosity, because the assumed porosity of gravels with a theoretical minimum value (0.363) is higher than the observed maximum (0.31) (Wu and Wang, 2006).The BM scheme estimates porosity well for soils with more gravels especially in the deep layers over the arid and semiarid zones.

SWRC and K s estimations
Using basic soil properties data of Tibet-Obs with the selected PTFs (see Table A1), parameters of SWRCs-CH and SWRCs-VG were estimated (see Table A3 in the Appendix).Figure 9 shows comparisons of the estimated SWRCs from PTFs combined with the BD porosity scheme to the measurement-determined SWRCs at 5 cm (see Sect. 3.2).The Saxton et al. (1986) PTFs overestimated the SWRCs-CH in the arid zone (Ngari), while the PTFs given by Campbell and Shiozawa (1992) and Saxton and Rawls (2006) underestimated them (Fig. 9a), and the Cosby et al. (1984) PTFs (nos. 1 and 2) presented good SWRCs-CH predictions with smaller absolute biases compared to measurements (see Table A4 in the Appendix).In the semi-arid zone (Naqu) (Fig. 9b), all PTFs underestimated the SWRCs-CH at 5 cm, while the Cosby et al. (1984) PTFs (nos. 1 and 2), Saxton et al. (1986), and Saxton and Rawls (2006) PTFs captured them well with lower biases compared to measurements (see Table A4).In the subhumid zone (Maqu) (Fig. 9c), the Cosby et al. (1984) PTFs (no. 1) and Saxton et al. (1986) PTFs predicted SWRCs-CH well.It is noteworthy that, in combination with the BD scheme, the Cosby PTFs (no. 1) performed much better regarding the estimation of SWRCs-CH, compared with the estimates by the Cosby PTFs (no. 1) combined with the Cosby-S porosity scheme (see Sect. 3.2).On the other hand, without the BD scheme, the Saxton and Rawls (2006) PTFs were found to be performing better over the semi-arid and semi-humid zones (see Table S5 in the Supplement).
For the SWRCs-VG estimate, the Rosetta1-H3 and Rosetta3-H3 PTFs were developed based on the mixed database (Schaap et al., 2001).Figure 9 (right panel) shows they underestimated SWRCs-VG across the three climate zones, as did the Rawls and Brakensiek (1985) PTFs.The Weynants et al. (2009) underestimated the SWRCs-VG in the semi-arid zone and the Class Wösten et al. (1999) PTFs overestimated them (Fig. 9b).The Vereecken et al. (1989) PTFs, which were developed based on a database where hydraulic properties were measured for every sample with the same measurement techniques (Vereecken et al., 2010), performed well when the m was set at 1.However, these PTFs were not performing well for m = 1 − 1/n in the VG model and estimated SWRCs-VG out of range in the subhumid zone.The Continuous Wösten et al. (1999) Saxton et al., 1986Campbell and Shiosawa, 1992Saxton et al., 2006Rawls and Brakenssiek 1985Class Wösten et al., 1999Vereecken et al., 1989Continuous Wösten et al., 1999Rosetta1-H3 2001Rosetta3-H3 2017 Weynants et  Taking basic soil properties data of Tibet-Obs as the input, the K s was estimated by using the PTFs scheme (see footnotes in Table A1), the empirical PTFs-VGF scheme (see  Eq. A16 in the Appendix) and the semi-physical BM-KC scheme (see Eqs. A17-A19).Comparing them against the in situ measured K s , Fig. 10a and d show the PTFs scheme had a lower bias for predicting Log 10 K s than the PTFs-VGF and BM-KC schemes did for the arid zone (Ngari).In particular, the PTFs given by Cosby et al. (1984) (nos. 1 and 2), predicted good K s that can be used in the CH model for estimating hydraulic conductivity (K), and as the Rosetta1-H3 PTFs, Rosetta3-H3 PTFs and Rawls and Brakensiek (1985) did in the VG model.The K s derived from BM-KC scheme had lower RMSE with measurements at 40 cm depth, indicating the gravel impact on K s .
Figure 10b and e show that the BM-KC scheme predicted better K s at depths of 10, 20 and 40 cm in the semi-arid zone (Naqu) than most PTFs and PTFs-VGF did.For K s estimate used in the CH model, the Cosby et al. (1984) PTFs (no. 1) performed best at shallow depths, while the PTFs-VGF of these PTFs were better at deep layers of 40 and 50 cm.For the usage in the VG model, K s derived from PTFs and PTFs-VGF schemes were almost the same, indicating that the estimated K s used in the VG model is less affected by gravels.The Rosetta1-H3 PTFs predicted K s better than other PTFs.Figure 10c and f show most of the PTFs underestimated K s , while the selected PTFs (i.e., Cosby (no. 1) and Rosetta1-H3) in the arid zone also predicted K s close to the measurements in the subhumid zone (Maqu).To sum up, the PTFs resulting from Cosby et al. (1984) (no. 1) and Rosetta1-H3 PTFs are appropriate for estimating K s , respectively used in the CH and VG models, across the three climate zones.PTFs-VGF of the Saxton and Rawls (2006) scheme should be applied in deep layers in the semi-arid zone, where gravel is abundant in the soil.

C s and λ estimations
Using basic soil properties data of Tibet-Obs, the C s was estimated through the De Vries (1963) model.Comparing to C s measured, this scheme performs well over the three climate zones.Furthermore, with the consideration of SOC impact, it improves the C s estimates for soils at top layers in the semiarid and subhumid zones (see Table A5 in the Appendix).
Based on the Tibet-Obs basic soil properties data, the D63F, T16 and J75 schemes combined with the BD porosity scheme were used to estimate the λ.For the arid (Ngari) and semi-arid (Naqu) regions, the estimation of λ considered two scenarios: with (case 1) and without (case 2) gravel impact.For the subhumid region (Maqu), λ estimations with (case 1)  and without (case 2) SOC impact were considered.Table 3 shows that the λ derived from D63F model had a lower bias in all cases compared to the measurement than other schemes across the three climate zones.The T16 scheme overestimated λ, which may be due to its ideal assumption that the λ of soil minerals is totally determined by sand, clay and silt particles.The J75 scheme generally underestimated the λ.
Table 3 also shows that the D63F scheme improved the λ estimate at surface layers in the arid zone and at a depth of 50 cm when incorporating gravel impact parameterization (lower biases in case 2).The improvement also occurred with the T16 scheme, while biases tended to be greater for the J75 scheme.In the subhumid zone, biases also became larger for all schemes when SOC impact parameterization was considered.Although parameterization of the SOC impact was demonstrated to improve the λ estimate in the top layer (SOC > 12 %) over the Eastern TP (Chen et al., 2012;Zheng et al., 2015), it should be noted that, in these studies for porosity estimate, the Cosby-S scheme was used instead of the BD scheme as adopted in this paper (see Sect. 4.1).Comparisons in Fig. 11 indicate that the D63F scheme combined with the BD porosity scheme can predict λ well across the three climate zones.It should be noted that, combined with the Cosby-S scheme, the D63F scheme also performs well (see Fig. S2 in the Supplement).

Basic soil properties
Figure 11 shows that all datasets underestimated both the sand fraction and BD in the arid and semi-arid regions, while they overestimated them in the subhumid region.For the silt fraction, the pattern was reversed.Almost all datasets overestimated the silt fraction in the arid and semi-arid regions (only FAO-UNESCO underestimated silt very slightly in the semi-arid region) and underestimated the silt fraction in the subhumid region.All datasets overestimated the clay fraction throughout the three climate zones.
The estimates of SOC from all the datasets were within a 1 % range of the measurements across the arid and semi-arid zones and within 10 % across the subhumid zone, apart from the FAO-UNESCO data, which underestimated the SOC heavily in this region.Most of the GGF estimates for the arid zone were within 10 %, with the FAO-UNESCO data underestimating by 20 %.For the semi-arid and subhumid regions, all datasets consistently underestimated and overestimated the GGF, respectively.
The BD scheme was used to derive porosity from the existing datasets.Figure 12a shows that the estimations of porosity were higher than the in situ measurement for the arid zone, with the SoilGrids1km and HWSD providing the closest approximations.In the semi-arid zone (Fig. 12b), all datasets underestimated porosity at the top layer, but overestimated it at other depths.It should be noted that Soil-Grids1km and SoilGrids250m presented porosity almost as a constant figure in each profile, which is not representative for conditions in the field.The porosity estimations from FAO-UNESCO, HWSD and BNU did show profile variation, although much less than the in situ measurements did.In the subhumid region (Fig. 12c), all datasets underestimated porosity in the surface layers 5, 10 and 20 cm, and either underestimated or overestimated porosity in the deep layers.

SWRC and K s
As previous analysis of PTFs (see Sect. 4.1) suggested, the Cosby et al. (1984) and continuous Wösten et al. (1999) PTFs were used with basic soil properties (i.e., only texture, BD and SOC) from the independent datasets (e.g., SoilGrids) to estimate, respectively, SWRCs-CH and SWRCs-VG.Given the relatively homogenous soil profile derived from existing Earth Syst.Sci.Data, 10, 1031Data, 10, -1061Data, 10, , 2018 www.earth-syst-sci-data.net/10/1031/2018/ Figure 11.Average bias in soil basic properties between the existing products and the laboratory measurements taken for the three climate zones.To enable the comparison of BD with the same order of magnitude as other properties, the original BD was multiplied by 100 (with unit × 100 g cm −3 ).Likewise, a multiplication (% × 10) is applied to SOC data on the semi-arid zone.FAO-UNESCO is the FAO-UNESCO Soil Map of the World (2007).HWSD is the Harmonized World Soil Database (FAO/IIASA/ISRIC/ISSCAS/JR, 2012).BNU is a Chinese dataset of soil properties (Shangguan et al., 2012(Shangguan et al., , 2013) ) and soil hydraulic parameters using PTFs (Dai et al., 2013) released by Beijing Normal University.SoilGrids 1km (Hengl et al., 2014) and the updated version of SoilGrids250 m (Hengl, 2017) datasets are released by the International Soil Reference and Information Center (ISRIC) -WoSIS institute.HPSS is the hydraulic parameters of the Van Genuchten (1980) and Mualem (1976) model based on SoilGrids1km and Schaap et al. (2001) PTFs (Montzka et al., 2017).products (see Fig. 13), the averaged SWRCs derived from existing datasets over different depths were used for comparison with the laboratory measurements.
Figure 13a shows all datasets overestimated SWRCs in the arid zone, in the order of FAO-UNESCO > BNU > HWSD > SoilGrids250m > HPSS for VG model > SoilGrids1km > Tibet-Obs.In the semi-arid zone (Fig. 13b), all datasets underestimated SWRCs at the surface layers of 5 and 10 cm, while they overestimated them at deep layers.and BNU presented the closest estimations for deep layers.Regarding SWRCs-VG, SoilGrids250m and HWSD, respectively, matched the measurements at the surface and deep layers well.In the subhumid zone (Fig. 13c), all datasets showed similar SWRCs-CH, slightly underestimating at low suction (< 100 kPa) but then becoming consistent with the measurements.The results for SWRCs-VG were quite diverse.The HWSD and HPSS showed consistent underestimation.
The FAO-UNESCO and BNU closely matched the measurements in deep layers.The SoilGrids1km and SoilGrids250m were within the range of the measurements across the whole profile, although their mean values were larger at high suction range (> 300 kPa).Furthermore, it should be noted that the averaged profile SWRCs derived from Tibet-Obs tended to reflect SWRCs at deep layers over the three climate zones.Additionally, the SoilGrids1km-, HWSD-and SoilGrids1km-derived FC (0.37, 0.41, 0.51 cm 3 cm −3 ) and PWP (0.16, 0.20, 0.27 cm 3 cm −3 ) were found close to the mean measured values over the three respective climate zones (see Table A6 in the Appendix).Using basic soil properties (i.e., only texture, BD and SOC) from the independent datasets (e.g., SoilGrids) with the Cosby et al. (1984) (no. 1) and Rosetta1-H3 PTFs, Table 4 shows the mean predicted K s (10 −6 m s −1 ) for all existing datasets across the three climate zones.They were of a smaller order than the field measurements in the arid and semi-arid zones but of a larger order than some of the field measurements in the subhumid zone.The Tibet-Obs dataset as input in the applicable PTFs predicted K s well.The existing datasets for estimating SWRCs -SoilGrids1km, HWSD and SoilGrids1km -also performed well estimating K s in the three climate zones, respectively.

SHP and STP in LSMs
Most of the LSMs use Richards' equation for soil water flow modeling (see Sect.A6 in the Appendix) with the hydraulic conductivity, while in some LSMs (e.g., Noah and H-TESSEL) soil diffusivity (D) is used.When the soil dries down, with the largest pores in the soil draining, the K and D are reduced by many orders of magnitude from saturation to dryness (Bittelli et al., 2015).Lower K (higher D) will result in slower water transport and thereby a higher SM derived from LSMs compared to soil moisture measurements, and vice versa.
LSMs use the thermal diffusion equation for soil heat transport modeling (see Sect.A6 in the Appendix).Soil heat capacity (C s ) and thermal conductivity (λ) are the important thermal parameters in the equation.Lower λ with higher C s will lead to reduced soil heat fluxes and thereby the higher soil temperature derived from LSMs compared to soil temperature measurements, and vice versa.The curves of K, D, C s and λ derived by using basic soil properties from the independent datasets (SoilGrids etc.) with recommended parameterization schemes were compared to the measurements (see Figs. A1-A3 in the Appendix) for quantifying the LSMs' uncertainty inherited from soil dataset.A special case is formed by the FAO-UNESCO dataset, which slightly overestimated VG-K at surface layers and heavily underestimated it at deep layers, while the dataset heavily overestimated VG-D at surface layers and slightly underestimated it at deep layers.These would lead to the overestimation of derived SM, as the ECMWF SM analyses do in this region (Su et al., 2013).The uncertainty from the soil dataset also propagates to soil temperature estimation.The FAO-UNESCO dataset underestimated Cs-SM at surface layers, but overestimated λ-SM, while at other depths it estimated Cs-SM well, but underestimated λ-SM.These would lead to the underestimation of the simulated soil temperature, which is also consistent with the findings of ECMWF soil temperature analyses (Su et al., 2013).

Data availability
Detailed information on data availability can be found in Sect.3.1.

Conclusions
For this study an in situ measurement dataset of soil physical properties was set up across the arid (Ngari), semiarid (Naqu) and subhumid (Maqu) climate zones across the Tibetan Plateau.The dataset can fill geographical gaps of global profiles on the third pole region.Analyzing this in situ dataset shows that soil texture in the Ngari network under the arid zone consists of a high proportion of coarse sand accompanied by gravel, and that the gravel content increases until 20 cm and then decreases slightly in the deeper layer.Dry bulk density (BD) and porosity vary with depths slightly.Soil texture in the Naqu network under the semi-arid zone is dominated by a high percentage of sand mixed with a small proportion of gravels, but with high SOC in shallow layers, and mainly mixed with big gravels at deep layers.The BD has a minimum in the top layer and a maximum in the bottom layer, and the porosity presents the opposite.Soil texture in the Maqu network under the subhumid zone is dominated by a high percentage of silt content with relatively large SOC in the shallow layers and with mainly fine sand mixed in the deep layers.The BD increases with depth and the porosity decreases.Depending on basic soil properties varying over three climate zones, soil hydraulic properties (SHP; i.e., soil water retention curve, hydraulic conductivity) and thermal properties (STP; i.e., heat capacity and thermal conductivity) differ for each climate zone and vary within each profile (e.g., presenting layering in the semi-arid and subhumid zones), and gravels were found affecting porosity and SHP and STP in the arid zone and in deep layers of the semi-arid zone.
Various schemes for estimating the porosity and SHP and STP over the TP were examined.The Cosby et al. (1984) PTFs proved more applicable for SHP estimation by the Clapp and Hornberger (1978) (CH) model, and the continuous Wösten et al. (1999) PTFs for SHP estimation by the Van Genuchten (1980) and Mualem (1976) (VG) model.The original formulation of the De Vries (1963) model can be deployed successfully for estimating the heat capacity of a profile.Furthermore, the De Vries (1963) model combined with the Farouki (1981) scheme (D63F) and with the implementation of the BD porosity scheme proved superior for estimating thermal conductivity.
Referenced by the measurements, uncertainties of the existing soil basic property datasets and their derived SHP and STP were quantified across the TP.This information is of significance in assessing the LSMs' uncertainty inherited from soil datasets and, moreover, in screening the proper soil datasets for LSMs over the TP.Furthermore, the existing soil property datasets can also be used as the ancillary data for SM retrieval.For example, the composited datasets of FAO and HWSD were used in the Soil Moisture and Ocean Salinity (SMOS) and Soil Moisture Active Passive (SMAP) SM product generation.Therefore, the information also became valuable for understanding uncertainties in satellite SM products inherited from soil maps.Based on the dataset comparison, this paper indicates that SoilGrids1km can reduce such uncertainty and is therefore recommended for use in the arid and subhumid zones, while the combination of FAO-UNESCO at shallow layers and HWSD at deep layers is recommended for the semi-arid zone over the TP.
In summary, this paper provides a comprehensive in situ measured dataset of soil physical properties over the TP and presents the applicable schemes to use for porosity and SHP and STP estimation in the LSM across the TP.The dataset contributes significantly for generating spatiotemporally consistent soil moisture and temperature estimates by LSMs.Furthermore, the evaluation of the existing soil property datasets is crucial for quantifying the uncertainty arising from soil data used in the LSMs and in soil moisture retrieval from microwave remote sensing.where K sm is the K s of soil mixtures.K sat,f is the K s of fine minerals and was calculated using PTFs in Table A1.VGF shares the same definition as for Eq.(A8).
A4.2 BM Kozeny-Carman equation (BM-KC scheme) The Kozeny-Carman equation (Eq.A17), originally developed to quantitatively describe hydraulic conductivity vs. the mean grain size in capillary flow, was used to estimate the K s of binary mixtures.The porosity was obtained by using the BM scheme in Sect.A1.The representative grain diameter was estimated using the power-averaging method (Eq.A18) proposed by Zhang et al. (2011).This method introduced a coefficient (Eq.A19) with the critical fraction of gravels taking into account.
where φ m has the same definition as in Eq. (A7).d m is the representative grain diameter of soil mixture.ρ is the fluid density.g is gravitational acceleration, and µ is the dynamic viscosity.
where VGF, VFF, GD and FD have the same definition as in the BM scheme in Sect.A1 displays.p is a coefficient that varies sigmoidally from −1 to 0, with VGF increasing from 0 to 1. p is estimated empirically by where VGF c is the critical fraction of gravels and is approximated by VGF c = 1 − φ g (φ g from Eq. A9).α is a shape factor set at a value of 20 as in Zhang et al. (2011).
A5 Heat capacity and thermal conductivity

A5.1 Heat capacity
Soil heat capacity C s depends on the heat capacities of all constituents and is calculated using Eq.(A20) given by De Vries (1963), where θ and θ s share the same meaning as in Eq. (A11).
C represents the heat capacity (MJ m −3 K −1 ), and the subscripts "w", "soil" and "air" refer to water, solid soil and air, respectively.C w , C soil and C air are taken as 4.2, 2.0 and 0.001 MJ m −3 K −1 , respectively.If taking SOC impact into consideration, C s is calculated as Eq.(A21) shows as follows: where V soc shares the same definition as in Eq. (A4).
C soc is the heat capacity of organic matter and taken as 2.5 MJ m −3 K −1 .
A5.2 Thermal conductivity by the De Vries (1963) model revised by Farouki (1981) (D63F) The De Vries (1963) model was developed from the Maxwell equation for electrical conductivity of a mixture of granular materials dispersed in a continuous fluid (Eucken, 1932).Farouki (1981) set liquid water as the continuous medium and regarded soil minerals as uniform particles.Considering soils as the binary mixture of fine minerals and coarse gravels, λ is estimated as follows: where w is the weighting factor, x is the volume fraction, λ is the thermal conductivity, and the subscripts "w", "a", "v", "m", "g" and "soc" refer to water, air, vapor, fine minerals, gravels and SOC composed of soil, respectively.
λ m was calculated using Eq.(A23):  Table A4.Biases of estimated SWRCs from PTFs combined with the BD scheme and the measurements at 5 cm at three climate zones.
The saturated thermal conductivity is calculated using Eq.(A33): where λ m has the same definition as in Eq. (A23).If considering the SOC impact, λ m was calculated using Eq.(A34): where ρ b is the dry bulk density (g cm −3 ).

A6 Soil water flow and heat transport
The vertical movement of water in the unsaturated zone of the soil matrix obeys the following Richards equation (Richards, 1931) for the volumetric water content θ : where D(θ ) (m 2 s −1 ) and K(θ ) (m s −1 ) are the hydraulic diffusivity and hydraulic conductivity, respectively, and S θ is a volumetric sink term associated with root uptake (m 3 m −3 s −1 ), which depends on the surface energy balance and the root profile.
The soil heat transfer is assumed to obey the following Fourier law of diffusion: where C s is soil thermal heat capacity (J m −3 K −1 ) and λ is thermal conductivity (W m −1 K −1 ).T is soil temperature ( • C).

Figure 1 .
Figure 1.The location of Tibet-Obs and the spatial distribution of soil sampling across three climate zones.(a) Tibet-Obs networks are distributed under the three different climatic zones where the zones were classified based on the FAO aridity index map.The dark blue color represents the area around Tibetan Plateau, with elevation lower than 3000 m above sea level (a.s.l.) (Zeng et al., 2016).(b), (c) and (d) are sampling distributions in the Maqu network under the subhumid zone, Naqu network under the semi-arid zone and Ngari network under the arid zone, respectively, in the KML image from Google Earth.It should be noted that the image acquisition times were in August, February and December, respectively.The triangle in ginger pink represents each sampling site.

Figure 3 .
Figure3.Profiles of mean soil basic properties at three climate zones.Top panel: variations in sand, clay, silt, GGF and SOC at various depths.Bottom panel: variations in GD and FD at different depths.GGF is the gravimetric gravel fraction.SOC is the soil organic matter content.FD is the mean particle diameter of fine minerals.GD is the mean particle diameter of gravels.
1038 H. Zhao et al.: Analysis of soil hydraulic and thermal properties

Figure 4 .Figure 5 .
Figure 4. Profiles of mean dry bulk density (BD) and porosity at three climate zones.
Figure5shows the pressure-cell measured SWRCs (markers in the figure) differed across the three climate zones.In the Ngari network under the arid zone, soil water retention quickly reduced as the suction slightly increased (Fig.5a).The same situation occurs for the deep layers in the Naqu network under the semi-arid zone (Fig.5b).In the Maqu network under the subhumid zone, soil water retention was high and gradually decreased as the suction increased (Fig.5c).Figure5also shows the CH and VG models captured the re-

Figure 7 .Figure 8 .
Figure 7. Scatter points of measured porosities (a, b) and K s (c, d) with GGF varying at different depths in the arid and semi-arid zones.
Figure8aand b show that the heat capacity C s went up and down with SM increasing for the Ngari and Naqu networks under the arid and semi-arid zones.Some samples from these two networks are fine-grained soils well mixed with gravels.For these samples, it is not easy to vertically insert needles of the KD2 Pro device (see Sect. 2.2); instead, the needles were buried with soils in the surface of the sample as the alternative for the measurement.Additionally, the KD2 needles might experience a slight skew when they touch the hard gravel.All these factors might cause the fluctuation of C s when SM increases, while the overall rising trend is still observed in Fig.8a and b.Figure8cshows C s almost steadily increased with SM for the Maqu network under the subhumid zone.Samples from the subhumid zone are all fine-grained soils and make the needle of the KD2 device insert easily and thus form a steady environment for the measurement.Figure8a-cshow that C s -SM varied slightly at different depths in the three climate zones.The C s ranged from 1 at oven dry state to 2.5 MJ m −3 K −1 as the soil reached saturation over the arid zone, 0.5 to 3 MJ m −3 K −1 over the semi-arid zone, and 0.5 to 2.4 MJ m −3 K −1 over the subhumid zone.Figure8d-f show how the relationship of λ-SM varied with depth.For the arid zone (Fig.8d), the λ-SM curves were very similar at each depth due to the nearly homogenous sandy soils across the whole profile (see Fig.3a).The mean λ ranged from 1.8 at saturation and 0.2 (W m −1 K −1 ) as the soils reached oven-dry state.In the semi-arid zone (Fig.8e), the λ-SM curves were stratified, and soils with gravels in deep layers (see Fig.3b) clearly had a higher λ (> 2 W m −1 K −1 ) than in other layers and other climate zones.In the subhumid zone (Fig.8f), the λ-SM curves also presented variations with depth, though within a much narrower range than in the semi-arid zone.Such variation is mainly caused by the sand distribution along the profile, which increased slightly with depth (see Fig.3c).The mean λ in the subhumid zone ranged from 1.6 to 0.2 (W m −1 K −1 ) as soils dried out.Furthermore, the surface layers in the semiarid and subhumid zones had lower λ values (Fig.8e and f) because of the SOC influence.

Figure 9 .
Figure9.Comparisons between estimated SWRCs from PTFs combined with the BD scheme and the measurement-determined SWRCs at 5 cm for three climate zones.It should be noted that the SWRC estimated fromVereecken et al. (1989) PTFs was out of range over the subhumid zone and removed (right figure in Fig.9c).

Figure 10 .
Figure 10.Comparisons of K s , derived from PTFs, PTFs-VGF and BM-KC schemes in the CH and VG models, with field measurements in the profile over three climate zones.

Figure 12 .
Figure 12.Comparisons between the porosity estimated from various existing datasets based on BD scheme and the in situ measurements.

Figure 13 .
Figure 13.Comparisons of SWRCs, derived from the applicable PTFs based on various datasets, with laboratory measurements.(a), (c) and (e) represent the SWRCs by the CH model based on six datasets.(b), (d) and (f) represent the SWRCs by the VG model based on seven datasets, in which HPSS only provides hydraulic parameters for the VG model.

Figure A1 .
Figure A1.Comparisons of derived soil conductivity (K) and soil diffusivity (D) by the CH model based on six soil datasets with those derived from the laboratory measurements.Given the relatively homogenous soil profile derived from existing datasets (see Fig. 12 in the text), the averaged K and D derived from existing datasets over different depths were illustrated.

Figure A2 .
Figure A2.Comparisons of derived soil conductivity (K) and soil diffusivity (D) by the VG model based on seven soil datasets with those derived from the laboratory measurements.Given the relatively homogenous soil profile derived from existing datasets (see Fig.12in the text), the averaged K and D derived from existing datasets over different depths were illustrated.It should be noted that HPSS only provides hydraulic parameters for the VG model.

Figure A3 .
Figure A3.Comparisons of derived C s -SM (a, c, e) and λ-SM (b, d, f) by the D63F model based on various datasets, with the measurements.Given the relatively homogenous soil profile derived from existing datasets (see Fig.12in the text), the averaged C s -SM and λ-SM derived from existing datasets over different depths were illustrated.

Table 1 .
Sampling approach for soil basic properties, SHP and STP over the Tibet-Obs.

Table 2 .
Pressure-cell determined parameters of the CH and VG models for the three climate zones.The scaling method used for the determination is Eq.(A15) in the Appendix.
PTFs were developed based on the database of Hydraulic Properties of European Soils (HYPRES) and as such were more affiliated www.earth-syst-sci-data.net/10/1031/2018/ Earth Syst.Sci.Data, 10, 1031-1061, 2018 1042 H. Zhao et al.: Analysis of soil hydraulic and thermal properties

Table 3 .
Biases of λ estimates based on D63F, T16 and J75 schemes combined with the BD scheme in the profiles across the three climate zones and the measurements.Case 1 is the bias (listed in the upper part of the table) derived from schemes not considering gravel impact parameterization for the arid and semi-arid zone or SOC impact parameterization for the subhumid zone.Case 2 is the bias (listed in the lower part of the table) with these parameterizations taken into consideration.The unit of the listed values is W m −1 K −1 .

Table 4 .
Comparisons of mean derived K s from the applicable PTFs for the CH and VG models based on various soil datasets, with the measurements.The unit of listed value is m s −1 .

Table A2 .
Bias and RMSE between simulated porosities with measurements at three climate zones.The unit of listed value is cm 3 cm −3 .

Table A3 .
Estimated parameters of the CH model by using PTFs over the three climate zones of the TP.It should be noted that θ s values aligned with each PTFs were estimated by their own original porosity estimate.θ s is calculated from in situ BD measurement.