Evaluation of Soil Thermal Conductivity Schemes for Use in Land Surface Modeling

Soil thermal conductivity (STC) is an important physical parameter in modeling land surface processes. Previous studies on evaluations of STC schemes are mostly based on direct measurements of local conditions, and their recommendations cannot be used as a reference in selecting STC schemes for land modeling use. In this work, seven typical STC schemes are incorporated into the Common Land Model to evaluate their applications in land surface modeling. Statistical analyses show that the Johansen (1975) scheme and its three derivatives, Côté and Konrad (2005, https://doi.org/10.1139/t04‐106), Balland and Arp (2005, https://doi.org/10.1139/s05‐007), and Lu et al. (2007, https://doi.org/10.2136/sssaj2006.0041), are significantly superior to the other schemes with respect to both STC estimations and their applications in modeling soil temperature and the partitioning of surface available energy in the Common Land Model. The Balland and Arp (2005, https://doi.org/10.1139/s05‐007) scheme ranks at the top among the selected schemes. Uncertainty analyses based on single‐point and global simulations both show that the differences in STC estimations can induce significant differences in simulated soil temperature in arid/semiarid and seasonally frozen regions, especially at deep layers. The hydrology‐related variables are slightly affected by STC variations, but slight changes in these variables can induce notable changes in soil temperature by altering soil thermal properties. These results emphasize the important role of STC in modeling soil thermodynamics and suggest the necessity of further developing STC schemes based on land modeling applications. According to the evaluation analyses, we recommend the Balland and Arp (2005, https://doi.org/10.1139/s05‐007) scheme as a more suitable selection for use in land surface models.


Introduction
Soil thermal properties (generally including heat capacity and thermal conductivity) are highly important in land surface process modeling. They influence a wide range of physical, biological, and chemical processes by regulating energy partitioning at the ground surface and energy distribution in subsurface soil layers (e.g., Lawrence & Slater, 2008;Luo et al., 2009;Ochsner et al., 2001;Pan et al., 2017;Peters-Lidard et al., 1998;Zhao et al., 2018). The soil thermal conductivity (STC), which quantifies the rate of heat transfer across different soil layers, directly determines the vertical profiles of soil heat flux and soil temperature and further affects soil freeze-thaw cycles, soil water movement, and other related processes such as root water uptake and transpiration (e.g., Cuntz & Haverd, 2018;Wang & Yang, 2018). Given the strong control on soil thermodynamics, STC can be considered as one of the most important physical parameters in land modeling studies (Luo et al., 2009;Peters-Lidard et al., 1998;Zhao et al., 2018).
Accurate estimation of STC is always difficult because it can be affected by a considerable number of factors, such as the mineral soil composition (particularly the quartz content), particle size distribution, porosity, dry density, soil moisture, and soil temperature (e.g., Farouki, 1981;Peters-Lidard et al., 1998). Based on these controlling factors, many theoretical and empirical models have been developed to predict STC (Dong et al., 2015;Farouki, 1981). For example, Wiener (1912) proposed parallel and series flow models, which respectively provided theoretical upper and lower limits for the prediction of STC. De Vries (1963) developed a Maxwell equation analogous model, which provided a theoretical description of STC as the weighted average of thermal conductivities of all soil components. Johansen (1975) provided a classical empirical method to estimate STC by interpolating between the values of dry and saturated soils, and the weighting coefficient of the saturated STC was obtained by fitting a logarithmic function with the degree of soil saturation using experimental measurements. Johansen's method also led to some other similar empirical models, in which the weighting coefficient of the saturated STC was determined not only by the degree of soil saturation but also by some other properties such as soil types and porosity (e.g., Côté & Konrad, 2005;Lu et al., 2007). A recent comprehensive review of STC models for unsaturated soils can be found in Dong et al. (2015). They analyzed key influential factors of STC and noted some common pitfalls among all the existing models.
For years, many studies have been conducted for the comparison of different STC models to evaluate their abilities to reproduce measured STCs. Early studies agreed that the Johansen (1975) model yielded the most accurate prediction of STC among traditional models by comparison with experimental measurements over a range of soil types and saturation levels (e.g., Farouki, 1981;Peters-Lidard et al., 1998). Barry-Macaulay et al. (2015) further compared the Johansen (1975) model to its three derivative models developed by Balland and Arp (2005), Côté and Konrad (2005), and Lu et al. (2007) and showed that all four models yielded comparably good fits to the experimental data, with the overall best estimates provided by the Côté and Konrad (2005) model. Zhang and Wang (2017) presented a thorough summary of the advantages and disadvantages of 13 typical STC models, and they demonstrated that the methods of Chen (2008), Haigh (2012), and Zhang et al. (2015) were superior to other schemes in predicting sand thermal conductivity.
Although the aforementioned work and many other related studies tried to determine which STC scheme was the best, most of them were based on direct STC measurements under specific experimental conditions or on local soil samples. Their conclusions or recommendations cannot be used as the reference to determine which STC scheme is more suitable for incorporation into land surface models (LSMs) due to the nonlinearity and uncertainty in depicting soil thermodynamics in land processes. Currently, few studies have evaluated STC schemes based on their land modeling applications, and the uncertainties of land surface modeling induced by STC estimations are often ignored. Therefore, there is a need to couple STC schemes with LSMs to reevaluate these schemes for regional or global land modeling use.
Based on early soil property data sets, most STC schemes were derived from only the mineral fractions of the soil (particle diameters <2 mm), while few have considered the effects of soil organic matter (SOM) and gravel (particle diameters >2 mm). SOM and gravel have been proven to have different hydraulic and thermal properties (e.g., relatively low thermal conductivity and high heat capacity for SOM and high thermal conductivity for gravel) than the mineral soil and thus can significantly affect soil moisture and temperature simulations and even alter the dynamics of boundary layer development in global climate models (e.g., Lawrence & Slater, 2008;Pan et al., 2017). With the increase in the amount of collected soil information, the SOM and gravel fractions have been gradually complied into several soil data sets, such as the Global Soil Dataset for Earth (GSDE) system model  and SoilGrids (Hengl et al., 2017). These data sets make it possible to add the effects of SOM and gravel into the existing STC schemes, and thus, a more complete description of STC estimation can be implemented in LSMs.
The presented paper focuses mostly on applications of STC schemes in land surface process modeling. Seven typical STC schemes are selected for analysis. In addition to direct comparison with laboratory measured STCs, these schemes are incorporated into the Common Land Model (CoLM), which is the land component of many Earth system models, for example, the Climate-Weather Research and Forecasting model (Liang et al., 2012) and Beijing Normal University Earth System Model ; hence, they can be evaluated based on the corresponding CoLM performances by comparison with in situ observed land state variables such as soil temperature and surface energy fluxes. Moreover, the uncertainty in land surface modeling introduced by STC estimations can also be quantified by analyzing the discrepancies among CoLM performances incorporating different STC schemes. This work is expected to (1) identify STC schemes that are more appropriate for land modeling use by evaluating the schemes from the aspects of both their abilities to reproduce measured STCs and their applications to simulate soil thermodynamics in LSMs and (2) clarify the importance of the STC parameterization in land surface modeling. The results could contribute to filling cognitive gaps between physical soil scientists and land process modelers.

STC Schemes
In this work, seven typical STC schemes existing in the literature are selected for analysis. Five of them are the Johansen (1975) scheme and its derivatives (Balland & Arp, 2005;Côté & Konrad, 2005;Farouki, 1981 andLu et al., 2007), and the other two are an empirical scheme from Tarnawski and Leong (2012) and a theoretical scheme from De Vries (1963). The formulations of these schemes are modified as a form of volumetric weighted combination of mineral soil, SOM, and gravel (see details in the following subsections) to accommodate the effects of all the available soil components. The key formulations of each scheme implemented in the CoLM are described below, and more details can be found in Dai et al. (2014). Johansen (1975) calculated the STC k (W/m/K) as a combination of the dry k dry and saturated k sat thermal conductivities, weighted by the coefficient K e (also called the Kersten number) as equation (1) depicts:

Johansen (1975) Scheme
where K e is expressed as a function of the degree of soil saturation S r (S r = θ/θ s , in which θ and θ s are the real and saturated soil water content [cm 3 /cm 3 ], respectively); the phase of water and the soil particle size are as follows: K e ¼ 0:7logS r þ 1; for unfrozen coarse−grained soils with S r >0:05 logS r þ 1; for unfrozen fine−grained soils with S r >0:1 S r ; for frozen medium and fine sands and fine−grained soils In equation (1), the dry thermal conductivity k dry was estimated using the volumetric weighted arithmetic mean of thermal conductivities of mineral soil, SOM, and gravel: where f minerals , f om , and f gravel are the volumetric fractions of mineral soil, SOM, and gravel among the soil solids and k minerals_dry , k om_dry , and k gravel are the corresponding thermal conductivities under dry conditions, which were empirically calculated as k minerals dry ¼ 0:135ρ d þ 0:0647 2:7−0:947ρ d ; (4) k om dry ¼ 0:05; where ρ d is the dry bulk density of the mineral soil (g/cm 3 ).

Journal of Advances in Modeling Earth Systems
The saturated thermal conductivity k sat was approximated as the volumetric weighted geometric mean of thermal conductivities of soil solids and water: where v m , v om , and v g are the volumetric fractions of mineral soil, SOM, and gravel among all soil components (including water) and k minerals_wet , k om_wet , and k gravel are the corresponding thermal conductivities under wet conditions. Here, k minerals_wet was calculated based on the quartz content: where v q is the volumetric fraction of quartz in the mineral soil and k q (7.7 W/m/K) and k o (2.0 W/m/K for v q > 0.2, and 3.0 W/m/K for otherwise) are the thermal conductivities of the quartz and nonquartz fractions of the mineral soil, respectively. v q was given by the lookup table of the 12 U.S. Department of Agriculture soil textural classes provided by Peters-Lidard et al. (1998), in which the quartz content was approximated as the median percentages of the sand fraction in each U.S. Department of Agriculture soil class. k om_wet was given as 0.25 W/m/K, and k gravel had the same value as that under dry conditions. k w is the thermal conductivity of water and has a value of 2.29 W/m/K for frozen soil and 0.57 W/m/K for unfrozen soil.
The Johansen (1975) scheme is the first scheme to introduce the Kersten number K e to calculate the bulk thermal conductivity by weighting the dry k dry and saturated k sat thermal conductivities. The K e -S r relationship includes the effects of the soil particle size distribution, the phase of water, and the degree of saturation simultaneously. Based on the Johansen (1975) scheme, a series of derivative schemes have been developed, since the effects of soil properties on STC can be studied through different K e -S r relationships and k dry and k sat parameterizations. The following four derivatives of the Johansen (1975) scheme are all based on equation (1), with only the K e -S r relationship or the formulations of k dry and k sat differing from the Johansen (1975) scheme. Farouki (1981) used the same equation proposed by Johansen (1975) to predict the thermal conductivity of dry soils (equations (3)- (6)). Under saturated conditions, however, instead of using the quartz content to calculate the thermal conductivity of the mineral soil, Farouki (1981) directly used the sand and clay content, and the bulk thermal conductivity of soil solids was estimated as the volumetric weighted arithmetic mean of thermal conductivities of all soil components. Thus, the saturated thermal conductivity k sat was empirically calculated as

Farouki (1981) Scheme
k minerals wet ¼ 8:80%sand þ 2:92%clay %sand þ %clay ; where %sand and %clay represent the gravimetric fractions of sand and clay in the mineral soil, respectively. Farouki (1981) also simplified the K e -S r relationship by ignoring the effects of the soil particle size distribution as follows: K e ¼ logS r þ 1; for unfrozen soils S r ; for frozen soils : & The Farouki (1981) scheme is easily used due to its simplicity, and the required inputs, such as the contents of sand, clay, SOM, and gravel, can be directly obtained from soil property databases. For these reasons, the Farouki (1981) scheme is always adopted as the default scheme to calculate STC in LSMs such as the CoLM and the Community Land Model. Côté and Konrad (2005) analyzed a large data set of measured thermal conductivities of dry soils, and they demonstrated that the relationship between the dry thermal conductivity and porosity depended strongly on the soil particle shape and grain size, which was not reflected by the Johansen (1975) scheme.

Côté and Konrad (2005) Scheme
Moreover, the logarithmical K e -S r relationship proposed by Johansen (1975) could lead to negative predictions of STC in dry conditions, which were out of range for realistic STC values. Due to these limitations, Côté and Konrad (2005) proposed a new method to predict the dry thermal conductivity and built a new K e -S r relationship to account for the effects of the soil saturation level on STC estimations. The new method to calculate dry thermal conductivity k dry was formulated as follows: where χ i (W/m/K) and η i are the empirical soil type parameters that account for the particle shape effects of the mineral soil, SOM and gravel. The new K e -S r relationship was proposed incorporating a parameter κ as follows: where κ is an empirical parameter that can be expressed as a function of the soil type and phase of water. The suggested values of χ i , η i , and κ were given in Côté and Konrad (2005).
The saturated thermal conductivity in the Côté and Konrad (2005) scheme was calculated in the same way as in the Johansen (1975) scheme (equations (7) and (8)).

Balland and Arp (2005) Scheme
Balland and Arp (2005) followed the formulations of Johansen (1975) to calculate the dry and saturated thermal conductivities (equations (3)- (8)), but they pointed out a fundamental flaw in the Johansen (1975) scheme that it failed to be seamless in estimating STC from dry to saturated and from fineto coarse-textured soils. To solve this issue, they introduced a new function to build a K e -S r relationship: where v s , v om , and v g are the volumetric fractions of sand, SOM, and gravel among all soil components, respectively, and α (0.24 ± 0.04) and β (18.1 ± 1.1) are adjustable parameters that are determined based on unfrozen experimental data by Kersten (1949) and Ochsner et al. (2001). This function can lead to continuous variations in STC predictions over the entire range of soil saturation levels and particle size distribution. Lu et al. (2007) conducted laboratory measurements of STC for 12 different soil samples at multiple saturation levels, and they showed that the Johansen (1975) scheme could not always account for situations at low saturations, especially for fine-grained soils. They also validated the Côté and Konrad (2005) scheme against these measurements and found that the K e -S r relationship proposed by Côté and Konrad (2005) was extremely sensitive to the parameter κ, which therefore introduced substantial uncertainty into STC predictions.

Lu et al. (2007) Scheme
To improve the K e -S r relationship of Côté and Konrad (2005), Lu et al. (2007) proposed the following new function: where α and β are soil texture (fine or coarse) dependent variables (see their fitted values in Lu et al., 2007).
The dry and saturated thermal conductivities of Lu et al. (2007) were calculated following the formulations of Johansen (1975), but a linear relationship to estimate the thermal conductivity of dry mineral soil was proposed by fitting the measurements to soil porosity:

Tarnawski and Leong (2012) Scheme
In the Johansen (1975) scheme and its derivatives, soil components were assumed to be parallel or in a way between parallel and serial to heat flow, and thus, the bulk thermal conductivity could be calculated as the arithmetic or geometric mean of the thermal conductivities of all soil components (Zhang & Wang, 2017). As

10.1029/2019MS001723
Journal of Advances in Modeling Earth Systems a more complex scheme, Tarnawski and Leong (2012) proposed a mixed serial-parallel arrangement for soil components. They assumed that heat was conducted through three pathways, namely, a solid uniform passage (Θ sb ), a series-parallel passage composed of solids connected to a parallel path of a portion of soil water (n w ) and a portion of soil air (n a ), and a path of water (Θ w ) and air (Θ a ) in a parallel arrangement. The symbols in the parentheses represent the volumetric fractions of each component in a soil. By applying a classical resistor model to each heat pathway, the bulk thermal conductivity of Tarnawski and Leong (2012) could be calculated by the following equation: where k s ¼ k vm mineralswet k vom omwet k vg gravel , k w , and k a are the thermal conductivities of the soil solids, water, and air. Θ sb and n wm could be empirically calculated by fitting to the contents of gravel and sand: where m gravel and m sand are the gravimetric fractions of gravel and sand, respectively. The portion of soil water (n w ) in the serial-parallel passage can be obtained by the following relationship: where X = 0.6 − 0.3(m gravel +m sand ) 3 is a minuscule pore water retention factor.

De Vries (1963) Scheme
The De Vries (1963) scheme is a classical theoretical scheme that was developed from the Maxwell equation.
In this scheme, the soil structure is assumed to be composed of ellipsoidal grains freely floating in a continuous pore fluid (air or water). Thus, the bulk thermal conductivity is estimated as a weighted average of thermal conductivities of soil components with consideration of their shape factors: where F s , F w , and F a are the shape factors of the soil solids, water, and air, respectively. These factors are empirically given as follows: where g i is a fitting parameter for ellipsoidal particles, calculated as

Journal of Advances in Modeling Earth Systems
taken from 40 Canadian soils and 21 soils from Italy, China, and Japan, which cover a wide diversity of soil conditions from loose to compact, organic to mineral, fine to coarse textured, and dry to wet. The database provides basic soil properties, such as the textural fractions and bulk density, for the STC schemes as the required inputs. The STC values of the Canadian soils are measured at a full range of saturation levels (S r = 0, 0.1, 0.25, 0.5, 0.7, and 1), and the other soils provide measurements under absolutely dry (S r = 0) and saturated (S r = 1) conditions. The measured STC of the Canadian soils as a function of saturation level and porosity are shown in Figure 1. It can be seen that the soils consistently display an increasing trend of STC with increasing saturation level, especially for those with lower porosity. This pattern is expected because soil water can significantly reduce the thermal resistance between the soil particle surface and pore fluid and promote interparticle heat conduction. The soils with low porosity generally have a higher STC than those with high porosity. Note that the STC values among different soils vary in a larger range with porosity at higher saturation levels, since soil water can enhance the effects of porosity on soil heat conduction. These features are essential for depicting the STC variations among different soil types and saturation levels, and thus, the STC measurements of Tarnawski et al. (2015Tarnawski et al. ( , 2018 are quite suitable for use as a reference to evaluate STC schemes.

Scheme Evaluation With Respect to CoLM Performances
The STC schemes are then incorporated into the CoLM to evaluate their applications in land surface modeling. The CoLM is derived from a community effort. The initial version of the CoLM was adopted as Community Land Model 2.0 for use with Version 2 of the Community Climate System Model (Bonan et al., 2002). Afterward, it underwent further development in China in many areas, such as the two-big-leaf model for calculating leaf temperatures and photosynthesis-stomatal resistance and the two-stream approximation model for simulating canopy radiation (Dai et al., 2004). To date, two versions of the CoLM have been released: CoLM2005 (Dai et al., 2003) and CoLM2014 . The updates of the latter version concentrate mostly on the global land surface data (e.g., basic soil properties and leaf area index), pedotransfer functions for estimating soil hydraulic and thermal parameters, and the numerical solution of the Richards equation for simulating soil water movement. Li et al. (2017) compared the two versions of the CoLM and demonstrated that CoLM2014 outperformed CoLM2005 in many aspects of energy and water budget simulations. Therefore, the version CoLM2014 is selected for our STC scheme evaluations. In the CoLM, the soil heat conduction is solved numerically via the following equation: where c is the volumetric soil heat capacity (J/m 3 /K), k is the STC (W/m/K), T is the soil temperature (K), t is time (s), z is the soil depth (m), and S h is the latent heat of phase change of soil water. Here, k is estimated by our selected seven STC schemes, and c is calculated as the volumetric weighted average of the heat capacities of the mineral soil, SOM, gravel, and soil water.
The comparison of the CoLM performances incorporating different STC schemes is conducted against in situ observations from 12 sites, which are listed in Table 1. Two of the sites, Nagqu and Arou, are from the Tibetan Plateau observatory networks. Nagqu is located in the central Tibetan Plateau, and Arou is located in the upstream portion of the Heihe River basin. The two sites feature a semiarid continental climate with seasonal freeze-thaw cycles, and their land cover types are both alpine steppes. The other sites are from the FLUXNET network (Baldocchi et al., 2001). The 12 sites are selected based on the following criteria: (1) The data records have little gap filling (less than 20%) for more than two consecutive years to ensure reasonable data quality; (2) the sites cover a wide range of land surface dryness from arid to moist; (3) the sites sample a range of climate zones (e.g., boreal, temperate, Mediterranean, and tropical), land cover types (e.g., trees, savanna, and grass), and soil types (e.g., sandy, loamy, and clay soils with or without gravel); and (4) the sites are distributed broadly around the globe. These criteria are chosen to ensure that the selected sites are The field observations to evaluate the CoLM performances with different STC schemes include the soil temperature profiles and surface turbulent heat fluxes. The observed soil temperature profiles can be obtained at the Nagqu and Arou sites. Nagqu provides observations at depths of 5, 10, 20, and 40 cm, and Arou provides measurements at 17 vertical nodes with the soil depth up to 3.2 m (Liu et al., 2011(Liu et al., , 2018. As a soil column is only divided into 10 layers in the CoLM, we only use the vertical measurements of Arou at depths of 2, 6, 10, 15, 20, 30, 60, 80, and 160 cm, which are closest to the depths at which the soil temperature is simulated in the CoLM. The simulated soil temperature is linearly interpolated to the measured depths of Nagqu and Arou for comparison. The observed sensible and latent heat fluxes are provided by all the sites. However, due to the known problem that eddy covariance systems do not have long-term energy balance closure (e.g., Kidston et al., 2010;Wilson et al., 2002), the observed fluxes cannot be directly compared to the simulated results. Following Sellers et al. (1989) and Blyth et al. (2010), the present work assumes that the measured Bowen ratio of land surfaces is correct and, thus, the Bowen ratio is used to evaluate the CoLM performances in terms of the model's ability to partition incoming radiant energy into sensible and latent heat. The observed and simulated latent heat values with magnitudes less than 5 W/m 2 are filtered out to avoid unreasonable Bowen ratio values during freezing seasons. Note that some of the sites also provide the observed ground heat flux but these measurements cannot be compared to the model output either, because the ground heat flux in the CoLM is calculated as the residual between the incoming net radiant flux and the outgoing turbulent heat fluxes to close the energy budget, and thus, it is not the real ground heat flux by definition.
The meteorological forcing data of the sites, including near-surface air temperature, relative humidity, pressure, wind speed, incident solar and infrared radiation, and precipitation, are used to drive the CoLM. Only the last year of the simulated results at each site are compared with the field observations from that year. For most of the sites, the duration of the available data is adequate to fully spin up the model, while for the FI-Kaa and ZA-Kru sites, the only 2 years of data are probably not sufficient to allow the model variables to reach equilibrium, especially for those in deep soil layers. Here, we repeatedly apply the 2-year data of FI-Kaa and ZA-Kru one more time to spin up the model so that the model can be spun up for three simulation years at the two sites. The basic soil properties for the CoLM runs at each of the sites except Nagqu are extracted from the global soil databases GSDE , which provides data on the mineral soil and SOM fractions and the bulk density, and SoilGrids (Hengl et al., 2017), which provides data on the gravel fraction. The two databases have kilometer and subkilometer spatial resolutions, and thus, the soil information at the sites can be accurately extracted. The local basic soil properties of the Nagqu site have been measured and provided by Pan et al. (2017).

Uncertainties in the CoLM Performances Introduced by the STC Estimations
To determine to what extent the STC can affect land surface modeling, the differences among the simulated results produced by the CoLM incorporating different STC schemes are analyzed to quantify the uncertainties. First, we mainly concentrate on the situations at two sites, Nagqu and AU-How, which represent dry and wet conditions, respectively. The differences among the simulated STC and soil temperature generated with different STC schemes are particularly investigated. The variations in the soil heat flux profiles, heat capacity, soil water content, and freeze-thaw cycles are analyzed as well. Note that we choose the Nagqu site to represent dry conditions instead of the driest site US-FPe for consistency with the evaluation analyses in which the soil temperature profiles at Nagqu are provided.
Then, we implement a global comparison between the CoLM performances with the identified superior STC scheme and the mean of the other schemes to see how the STC affects land surface modeling on the global scale. The two simulations are performed for 3 years (2000)(2001)(2002), with the first 2 years used for spin-up and the last year used for analysis in terms of the annual mean results. Similar to the runs at single point locations, the global runs also use the basic soil information from the GSDE and SoilGrids databases, which have a spatial resolution of 1 km. To eliminate the uncertainties introduced by upscaling processes, the global simulations follow the super high spatial resolution of the soil property databases, and hence, only two global runs instead of seven runs with each of the STC schemes are performed to avoid excessive calculations and data storage. The atmospheric forcings are from Version 7 of the Climate Research Unit-National Centers for Environmental Prediction data with QIAN-NCEP ocean fill (see https://www.earthsystemgrid.org/dataset/ ucar.cgd.ccsm4.CRUNCEP.v4.html for download information), which are interpolated to the spatial resolution of the soil property databases. In the analyses, the global patterns of the differences in the soil temperature and moisture at multiple layers, surface energy fluxes, and snow cover simulated with different STC schemes are presented. The results highlight the importance of the STC estimation on land surface modeling.

Statistical Metrics
The major metrics for the evaluation of the STC schemes are the bias, root mean squared error (RMSE), and relative error (σ) defined as where k pre and k obs are the predicted and measured STCs, N is the total number of sample data, and median (RMSE) is the median value of RMSE among all the schemes. The relative error depicts the relative performance of each STC scheme, indicating how much better or worse the performance is compared to the median level of the scheme predictions (Gleckler et al., 2008). In the evaluation of the CoLM performances, only the RMSE is used because the bias cannot reveal the absolute accuracy of the CoLM simulations due to the offset between positive and negative biases. The differences of the simulated results with different STC schemes are quantified using the standard deviations (SDs) of variables with respect to the change of the STC scheme incorporated into the CoLM. We limit our analyses to daily means of the model outputs to reduce the uncertainty of observations in each single measurement (Hagen et al., 2006).

Journal of Advances in Modeling Earth Systems
close to the measurements, while the scatter of the predictions increases as the soils become saturated. Statistical analyses show that the Johansen (1975) scheme and its three derivatives, Côté and Konrad (2005), Balland and Arp (2005), and Lu et al. (2007), produce comparably satisfactory predictions, with the overall biases of −0.03 to 0.01 W/m/K and RMSEs of 0.26 to 0.29 W/m/K. Among the four schemes, the Balland and Arp (2005) scheme provides the lowest absolute error. The two mechanistic schemes, Leong (2012) andDe Vries (1963), stay at the median level, and their predictions evidently underestimate STC for moderately and highly saturated soils. The Farouki (1981) scheme yields the least accurate predictions, by which the STC is significantly overestimated for most soils. The overall bias of the Farouki (1981) scheme is 0.25 W/m/K, and the RMSE is up to 0.47 W/m/K.
To further ascertain the relative performance of the seven schemes, the relative error based on the RMSE of each scheme is calculated at all soil saturation levels ( Figure 3). The results show that the Côté and Konrad (2005) scheme produces the most accurate predictions for absolutely dry soils (S r = 0) and its RMSE is 30% lower than the median RMSE. The Tarnawski and Leong (2012) scheme and the Balland and Arp (2005) scheme provide the best predictions at the saturation levels of 0.1 and 0.25, respectively. At saturation levels from 0.5 to 1.0, the Johansen (1975) scheme, Côté and Konrad (2005) scheme, Balland and Arp (2005) scheme, and Lu et al. (2007) scheme all perform beyond the median level of the scheme predictions, with the top estimations of STC made by the Balland and Arp (2005) scheme. The Farouki (1981) scheme provides the worst predictions at all saturation levels except the level 0 (absolutely dry soils), with the RMSE at least 30% higher than the median RMSE at each level. This is consistent with the pattern shown in Figure 2, in which the STC is significantly overestimated by the Farouki (1981) scheme for most soils. A possible reason for the relatively poor predictions of the Farouki (1981) scheme is that the effects of quartz content are ignored in the calculation of saturated thermal conductivity. Overall, it can be seen from the relative errors that the predictions made by Johansen (1975), Côté and Konrad (2005), Balland and Arp (2005), and Lu et al. (2007) are generally better than those from the other schemes and the Balland and Arp (2005) scheme ranks at or near the top at all saturation levels. The relatively better performance of the Balland and Arp (2005) Figure 3. Relative errors of the seven soil thermal conductivity schemes in terms of six saturation levels. The error measure, treating each saturation level separately, is calculated by normalizing the root mean square error of each scheme by the median root mean square error of all the schemes, with blue shading indicating performance being better, and red shading worse, than the median level of scheme predictions. k dry , k 0.1 , k 0.25 , k 0.5 , k 0.7 , and k sat represent the soil thermal conductivity measured at saturation levels 0, 0.1, 0.25, 0.5, 0.7, and 1, respectively. Dat1 represents the relative errors measured based on the 40 Canadian soils, and Dat2 represents the corresponding values from the 21 other soils.

Journal of Advances in Modeling Earth Systems
scheme can be likely attributed to its K e -S r relationship, which can depict STC variations continuously from dry to saturated and from fineto coarse-textured soils. This was also noted in the evaluation analyses by Zhang and Wang (2017).
Next, we compare the CoLM performances incorporating different STC schemes to evaluate the STC schemes with respect to their land modeling applications. Figure 4 presents the simulated and observed soil temperatures at the Nagqu site. In each of the soil layers, the CoLM can reproduce the general variations of soil temperature with all the STC schemes. The correlation coefficients between the simulations and observations are all greater than 0.95. The discrepancies among the simulations with different STC schemes are slight at shallow layers, but they tend to be notable in deep soils. The RMSE of the simulated soil temperature with each of the STC schemes is given in

Journal of Advances in Modeling Earth Systems
can produce significantly more accurate soil temperature than that with the other schemes. The vertical averaged RMSE of the simulations with the four better-performing schemes is between 1.52 and 1.57 K, which is at least 0.2 K lower than that with the other schemes. The situations at the Arou site are similar to those at the Nagqu site (Table 3). For all the soil layers other than the surface layer at Arou, the most accurate soil temperature is provided by the simulations with one of the aforementioned four better-performing schemes. This pattern is especially notable in the middle and deep layers. The vertical averaged RMSE of the simulated soil temperature with any of the four better-performing schemes is less than 2.0 K. Note that among the four better-performing schemes, the CoLM incorporating the Balland and Arp (2005) scheme provides the best soil temperature profiles in terms of the accuracy of vertical means at both the Nagqu and Arou sites. Table 4 presents the RMSE of the Bowen ratio between the observations and simulations with each of the STC schemes at all the selected sites. Similar to the results of soil temperature, the most accurate Bowen ratio is generally produced by the simulations with the Balland and Arp (2005) scheme at most of the sites. However, the RMSE values of the Bowen ratio simulated with different STC schemes are similar to each other, especially at sites with dense forests such as IT-Ro1 and US-Ha1. This result is not surprising because the partitioning of surface energy into sensible and latent heat is mostly determined by the descriptions of ground surface properties, plant physiological and ecological processes, the status of the atmosphere, and the associated parameterizations with respect to soil-plant-atmosphere interactions (e.g., Blyth et al., 2010;Dai et al., 2003), whose uncertainties can dominate the biases of simulated surface turbulent fluxes.
In summary, the Johansen (1975) scheme and its three derivatives, Côté and Konrad (2005), Balland and Arp (2005), and Lu et al. (2007), exhibit superiority over the other schemes with respect to both the STC Table 3 The Same as

Quantification of the Uncertainties in Land Surface Modeling Introduced by STC Estimations
In this section, we quantify the uncertainties of the CoLM performances introduced by the STC estimations to investigate to what extent the STC affects the modeling of soil thermodynamics and its related processes in LSMs. The SDs of the simulated variables with respect to the change of the STC scheme incorporated into the CoLM are used to represent the differences among the CoLM performances with different STC schemes. First, the situations at two sites, Nagqu and AU-How, are analyzed. The former site is characterized by dry and seasonally frozen conditions, and the latter is characterized by wet and nonfrozen conditions. Figure 5a presents the SD of the simulated vertical STC profiles at the Nagqu site. A high SD of STC mainly appears at surface layers when heavy precipitation occurs and at middle and deep layers when the saturation levels of the soil become relatively high (Figure 5b). The good correspondence between the SD of STC and the degree of soil saturation is expected because a high soil saturation level can enhance the differences among the STC schemes due to the different K e -S r relationship they adopt. During the unfrozen period (July to November), the differences among the estimated STC associated with the changes in soil heat conduction in deep layers can induce slight changes of the simulated soil water content (Figure 5c). Despite the small magnitude of the changes (~0.4%), the soil water content can induce significant changes of the simulated bulk heat capacity due to the strong ability of heat storage of water ( Figure 5e). Thus, the simulated soil temperature can be significantly affected, with the magnitude of the changes reaching 2 K in deep soil layers (Figure 5g). During the frozen period (December to June), the simulated soil temperature with different STC schemes remains mostly the same (i.e., the frozen point) due to the coexistence of liquid and solid soil water. Significant changes of the simulated soil temperature appear around the dates of fall ice freezeup and spring ice breakup across most soil layers (Figure 5g), which corresponds well to the SD pattern of the simulated soil ice content (Figure 5d). Very likely the changes in soil heat conduction due to different STC schemes can lead to different predicted dates of ice freezeup and breakup, during which the differences in the simulated soil heat flux (Figure 5f) can induce a large SD in the simulated soil temperature. Unlike the significant responses of the simulated soil temperature and other related variables to the differences in the estimated STC at the Nagqu site, the simulated soil temperature at AU-How seems to be insensitive to the STC estimations during the entire simulation period (Figure 6c). The SD of the simulated STC is high in the surface and deep soil layers due to the high saturation levels (Figures 6a and 6b). However, the extremely high soil water content at AU-How (due to a large amount of precipitation) can lead to a considerably high soil heat capacity, thereby increasing the soil heat inertia. Therefore, the simulated soil temperature exhibits little sensitivity to STC estimations.
Overall, under relatively dry conditions, the simulated soil temperature is significantly affected by the STC estimations due to the differences in the simulated soil water content, which alter the bulk heat capacity during the unfrozen season, and due to differences in the dates of ice freezeup and breakup during the frozen season. Under wet conditions, the high soil heat inertia due to the high soil water content can make the simulations insensitive to STC estimations.
Next, we conduct a global comparison between the CoLM performances with the identified superior scheme, the Balland and Arp (2005) scheme, and the mean of the other schemes to see how the STC affects land surface modeling on the global scale. Figure 7 shows the differences in the global patterns of the simulated soil temperature over multiple layers. Consistent with the simulated results at single-point locations, the most significant differences appear over arid and semiarid regions at middle and high latitudes. The magnitudes of the differences increase with increasing soil depth and are more than 1 K in northern frozen areas.
Examinations of the differences of the simulated ground surface energy components reveal that the ground heat flux is significantly affected by STC estimations (Figure 8a), with approximately 3 W/m 2 more heat transported into soils simulated with the Balland and Arp (2005) scheme. This difference explains the overall warming simulated with the Balland and Arp (2005) scheme as shown in Figure 7. The more ground heat flux is mostly balanced by the sensible heat flux (Figure 8b), while the latent heat flux and net radiative flux do not change significantly (Figures 8c and 8d). Here, we also examine the differences of some related hydrological variables. Figure 8e shows the less snow cover simulated with the Balland and Arp (2005) scheme over high latitudes, which likely results from the greater ground heat flux, which can melt the snow cover. As the melted snow cover provides more available water to infiltrate into the soil layers, the simulated soil water content at high latitudes tends to be higher in all soil layers (Figures 8f, 8g, and 8h). Although the magnitudes of the differences in snow cover and soil water content are very small, they can significantly affect soil thermal properties and thus can probably induce notable changes in soil temperature, as shown by the situation at the Nagqu site ( Figure 5).
The above analyses in this section illustrate that the STC can significantly affect the modeling of soil thermodynamics over arid/semiarid regions and seasonally frozen regions. The changes in soil temperature and freeze-thaw cycles are the responses to STC variations, especially in deep soil layers. Hydrological processes, such as soil water movement and surface evaporation, are slightly affected by STC variations, but these slight changes can probably induce notable changes in soil temperature by altering the soil thermal properties. This work is limited by the off-line simulations due to considerable amounts of calculations and data storage in online simulations, and thus, the effects of STC on climate systems are not considered. In fact, the climate over boreal dry and frozen regions has been shown to be sensitive to external forcings due to its stable conditions and small atmospheric effective heat capacity (e.g., Davy & Esau, 2014a, 2014b. Thus, even slight differences in simulated ground surface temperature and heat fluxes transported into the atmosphere can lead to remarkable changes in predicted atmospheric temperature, boundary layer height, cloud cover fractions and many other related variables (e.g., Lawrence & Slater, 2008;Wei et al., 2014). This suggests the need to couple an LSM incorporating different STC schemes into an earth system model to investigate the climatic effects of STC in future work.

Conclusion
In this work, seven typical STC schemes (with modifications of including SOM and gravel effects) are evaluated for their applications in land surface modeling. These schemes are first compared with a database of laboratory measurements to evaluate their abilities to predict STC. Then, they are incorporated into the CoLM to evaluate their applications in modeling soil thermodynamics by comparing the corresponding CoLM performances with in situ observations from 12 sites. Statistical analyses show that the Johansen (1975) scheme and its three derivatives, Côté and Konrad (2005), Balland and Arp (2005), and Lu et al. (2007), are significantly superior to the other schemes with respect to both the STC estimations and their applications in modeling soil temperature and the partitioning of surface available energy. The Balland and Arp (2005) scheme ranks at the top among all the schemes.
To determine to what extent the STC affects land surface modeling, the uncertainties of the CoLM performances introduced by STC estimations are quantified by analyzing the discrepancies among the simulated results of the CoLM incorporating different STC schemes. Simulations at single-point locations show that under relatively dry conditions, the simulated soil temperature is significantly affected by different STC estimations due to differences in the simulated soil water content, which alter the bulk heat capacity, and the changes in freeze-thaw cycles. Under wet conditions, the simulated soil temperature is insensitive to STC estimations due to the high soil heat inertia. Global comparisons between the CoLM performances with the Balland and Arp (2005) scheme and the mean of the other schemes suggest that the differences in STC can induce significant changes in soil temperature in arid/semiarid regions and seasonally frozen regions, especially in deep layers. The hydrology-related variables are slightly affected by STC variations, but their slight changes can probably induce notable changes in soil temperature by altering soil thermal properties. These results emphasize the important role of STC in modeling soil thermodynamics and suggest the necessity of further developing STC schemes based on land modeling applications.
We realize that our results for the evaluations of the STC schemes may be site dependent because only 12 sites are considered in our analyses. However, our selected sites sample a sufficiently broad range of land surface dryness, climate zones, and land cover types, and thus, the results demonstrating the superiority of the Johansen (1975) scheme, and its three derivatives are expected to hold for global applications. The highestranking scheme, in this case the Balland and Arp (2005) scheme, could vary with the site selection because the difference in superiority between the Balland and Arp (2005) scheme and the other three betterperforming schemes is not significant. Given the particular advantage of the Balland and Arp (2005) scheme that its K e -S r relationship can produce continuous variations of the predicted STC over the full range of soil saturation levels and particle size distribution, the Balland and Arp (2005) scheme can still be strongly recommended for use in land surface modeling. In addition, our evaluation analyses are only based on the CoLM performances, and users of other LSMs should take care in the selection of STC schemes. However, as the essential physics of soil thermodynamics are mostly the same among the existing LSMs, the Johansen (1975) scheme and its three derivatives, especially the Balland and Arp (2005) scheme, can be considered suitable for use in other LSMs.