A Review and Evaluation of Predictive Models for Thermal Conductivity of Sands at Full Water Content Range

The effective thermal conductivity (λeff) of sands is a critical parameter required by applications in geothermal energy resources, geo-technique and geo-environment and in science disciplines. However, the availability of the reliable λeff data is not sufficient and predictive models are usually used in practice to estimate λeff. These predictive models may vary in complexity, flexibility, accuracy and applications. There is no universal model that can be applied to all soil types and full water content range. The choice of different models may result in distinctive estimates of λeff. The objectives of this study were to conduct an extensive review of the thermal conductivity models of sands and evaluate their performance with a large dataset consisting of various sand types from dry to saturation. A total of 14 models to predict λeff of sands were evaluated with a large compiled dataset consisting of 1025 measurements on 62 sands from 20 studies. The results show that the models of Chen 2008 (CS2008) and Zhang et al. 2016 (ZN2016) give the best estimates of thermal conductivity of sands, with Nash–Sutcliffe efficiency = 0.9 and RMSE = 0.3 W m−1 ◦C−1. These two models are potentially applied to accurately estimate thermal conductivity of sands of different types.


Introduction
The effective thermal conductivity (λ eff , W m −1 • C −1 ) of sands is of great interest because it is required by petroleum or natural gas production from underground oil-or methane-bearing sands [1,2]. It is also the critical parameter to determine heat exchange between the surrounding environments and the borehole heat exchanger [3][4][5][6], nuclear waste disposal [7] or energy storage [8], buried power cables and steam/hot water pipelines [9], and other thermo-active ground structures (e.g., energy pile and tunnel linings). The value of λ eff is often an input for calculating the ground thermal regimes or land-atmosphere energy balance and associated soil physical, chemical and biological processes [10,11].
The values of λ eff can be either measured or indirectly predicted with mathematic models. The available λ eff data of sands are not sufficient to meet the demands, although great progress has been achieved in the measurement technologies of λ eff [12][13][14]. In addition, reliable λ eff data can only be measured by the transient method (e.g., heat pulse method) at full water content range, because the transient method is less likely to be affected by the water redistribution or phase change compared to the steady-state method (e.g., guarded hot plate). On the other hand, modeling of λ eff is a highly thriving field, because there is no single model for universal applications. Previous studies generally applied models of other applications to estimate λ eff of sands, while some developed predictive models for sands but based on a small size dataset [15,16] or unreliable measurements [17]. The complexity, flexibility, accuracy and applications of the predictive models largely determine the accuracy of their applications [18,19]. Therefore, it is critical to evaluate the performance of the predictive models for sands for accurate applications in engineering and science.
The objectives of this study were: (1) to conduct an extensive review of the predictive thermal conductivity models of sands; (2) to compile a large and reliable dataset consisting of measurements on various sands at full water content range with the transient heat pulse methods; and (3) to evaluate the performance of the collated models with the compiled dataset.

Woodside and Messmer 1961 Model (WM1961II)
Woodside and Messmer [1] also presented another model similar to the geometric mean model [31][32][33]] Kiyohashi and Deguchi [30] were among the first to extend Equation (1) of Woodside and Messmer [1] to unsaturated solid rocks (e.g., sandstones, zeolites, tuff, etc.) with 0.05 < n < 0.51 (with accuracy of λ eff = ±0.20 W m −1 • C −1 ) Tarnawski and Leong [20] found Equation (3) overestimates λ eff of unsaturated soils due to the model is lack of physical characteristics for unsaturated soil applications. They introduced the inter-particle thermal contact resistance factor (α) where α is expressed as [20,34,35] where S r−cr is the degree of saturation of a miniscule pore space, S r−cr = (0.22f clay + 0.01)/P [20,36], ε is a dimensionless interface contact coefficient (0 ≤ ε ≤ 1, ε = 1 indicates an ideal contact between soil particles). ε depends on soil compaction, soil specific surface area, and grain size distribution. Tarnawski and Leong [20] obtained ε value by reverse modeling of experimental λ eff data of 40 Canadian soils; the resulting ε varied between 0.988 and 0.994 for coarse and fine soils, respectively. For dry soils, αǎ1 due to the highest λ reduction arising from the thermal contact resistance.

Haigh 2012 Model (HS2012)
Based on dataset of Chen [15], Haigh [16] developed an analytical model based on unidirectional heat flow through a three-phase soil system [17,37], where F = 1.58 is a correction factor, λ w and λ air are thermal conductivity of water (λ w , W m −1 • C −1 ) and air (λ air , W m −1 • C −1 ), respectively, that are divided by λ s where coefficients a and b are related to the thickness of water film and degree of saturation where e is the void ratio, e = P/(1 − P). The left side of Equation (10) where a and b are fitting parameters. Lu et al. [21] recommended a = 0.76 and b = 1 for unfrozen soils and a = 0.72 and b = 2.06 for frozen soils.

Rubin and Ho 2018 Model (RH2018)
Rubin and Ho [22] showed that the Haigh [16] model underestimates λ eff and they provided an alternative empirical equation for the correction factor F [38] where x is the estimated λ eff with the Haigh [16] model when F = 1, multiplying the λ eff calculating with Haigh [16] by the F gives the corrected λ eff . The R 2 for Equation (12a,b) are 0.65 and 0.87, respectively.

Ewen and Thomas 1987 Model (ET1987)
Ewen and Thomas [23] proposed a modified Johansen [39] model where λ sat and λ dry are thermal conductivities when soil is saturated or dry, respectively. K e is an exponential function that is dependent on degrees of saturation. K e is calculated by [23] where ξ is a fitted parameter, ξ = 8.9 was recommended by Ewen and Thomas [23] which returns K e ≈ 1 at S r = 1. S r = θ/P. Values of λ sat and λ dry are calculated by a two-phase system (solid-air or solid-water) [39] Unfrozen soil Frozen soil (a) (b) (15) where θ l is the unfrozen water content (cm 3 cm −3 ), λ w is the thermal conductivity of water, λ w = 0.56 W m −1 • C −1 at above-zero temperature, λ w = 0.269 W m −1 • C −1 at sub-zero temperature [40], λ ice is the thermal conductivity of ice (2.2 W m −1 • C −1 ), ρ b is the bulk density (g cm −3 ), λ s is the thermal conductivity of the solid that can be calculated based on volumetric fraction and thermal conductivity of quartz (f quartz , λ quartz ) and other minerals (λ other ), where

Zhang Et Al. 2016 Model (Z2016)
Zhang et al. [25] further modified the Côté and Konrad [41] model with experimental data measured on sand-kaolin clay mixtures with thermo-TDR method, the resulting formula is where λ kaolin is thermal conductivity of kaolin (2.9 W m −1 • C −1 ), Kaolin is assumed to be zero in this study. Equation (16) is used to calculate λ s .

Zhao Et Al. 2019 Model (ZY2019)
Zhao et al. [26] presented a closed-form thermal conductivity model for mixtures of sands and peat materials where A and B are parameters related to soil properties (composition, texture and structure)-they can be directly derived from the lower and upper water saturations. λ eff = λ dry = A when S r = 0 and B = (λ sat − λ dry )/Ln2 when S r = 1. Equation (20) can be converted to the normalized form as Ke = Log 2 (1 + S r ) where Equations (15a), (16) and (17a) are retained for calculation of λ sat , λ s and λ dry , respectively.

Becker 1992 Model (BB1992)
Becker et al. [27] developed an unified methodology for predicting λ eff of gravel, sand, silt, clay and peat in both frozen and unfrozen states using literature data [42] S r = 0.01a 1 sinh 6.9348a 2 · λ e f f + a 3 − sinh(a 4 ) (22) or where for consistence with the original study. a 1 to a 4 are coefficients depending on soil type. For unfrozen sand, a 1 to a 4 are 6.8, 0.4, −2.9 and −1.5, respectively.

Chen 2008 Model (CS2008)
Chen [15] proposed a simple empirical relationship for sandy soils where b = 0.0022 and c = 0.78 are the fitting parameters obtained using DPS code with R 2 = 0.996, λ s is calculated with Equation (16).

Li Et Al. 2012 Model (LY2012)
Li et al. [28] developed an empirical relationship for powdery sand while test soil thermal conductivity along a crude oil pipeline in tropical desert and grassland in southwest Sahara Desert in West Africa,

Experimental Data of Thermal Conductivity of Sand
A large dataset was collated by digitaltizing literature data or by collecting through personal communications following several important criteria: (1) soil must be sands with zero or a negligibly small amount of organic matter. Full dataset of Chen [15] was included because it has been used by a few researchers to compare, develop, or test new models [16,18,22,38], additional soils contain sand content >90%; (2) λ eff were measured on soil samples with the transient method (e.g., single/dual-probe heat pulse method, thermo-TDR) at room temperature; (3) sufficient sample size and a wide range of water contents or degrees of saturation (e.g., >5) and bulk densities should be provided for useful evaluation; and (4) detailed description of the sample preparation and complete soil information such as grain size distribution (f sand , f silt , and f clay ), P, ρ b , and ρ s should be available. The values of λ s and f q that were measured or recommended by original research were used, otherwise f q = f sand is assumed and missed λ s was calculated by Equation (16).

Performance Metrics
The similarity between the predicted and actual values was assessed by plotting them on a 1:1 diagram with ±10% line and two goodness-of-fit parameters were used, including (1) root mean square error (RMSE) and (2) Nash-Sutcliffe efficiency (NSE): where Y j is the j th measured value,Ŷ j is the j th predicted value, Y is the mean of measured values, and n is the number of measured values. NSE = 1 indicates a perfect match betweenŶ j and Y j ; NSE = 0 indicates thatŶ j is as accurate as Y; while NSE < 0 indicates that Y is a better predictor. RMSE indicates the absolute difference betweenŶ j and Y j -the closer the RMSE value is to zero, the better the predictor.

Analytical or Mixing Models
The best performing model among the six thermal conductivity models of this category is TL2016 (Table 2 and Figure 1(3)) with RMSE = 0.43 W m −1 • C −1 and NSE = 0.79. The HS2012 (Figure 1(4)) model retains a physical origin but it is relatively complex [17,37] and does not have a satisfactory performance on the dataset investigated. The models of TL2016 (Figure 1(3)) and WM1961I (Figure 1(1)) generally underestimate the measured λ eff . The WM1961II (Figure 1(2)), LY2018 (Figure 1(5)) and RH2018 (Figure 1(6)) generally overestimate the measured λ eff .   The WM1961I and WM1961II were originally developed for water saturated soils with λ s /λ w ≈ 10 [1]. Although the geometric mean model was extended for frozen soil study by Cosenza et al. [31] and Endrizzi et al. [32], and incorporated into the GEOtop program [33], the extension to three-phase system may not always be valid because this model lacks a physical basis for the heat conduction in granular materials [58,59]. The LY2018 model was developed based on the basic series-parallel model using a small dataset (28 measurements) of Aeolian sands from the Tibetan Plateau [21]. The RH2018 model (Figure 1(6)) is in fact the modified form of HS2012 (Figure 1(4)) model by multiplying an F factor based on Equation (12). However, it should be noted that Equation (12a) gave an unreasonably high value at S r < 0.1 and only Equation (12b) was used for all degrees of saturation in this study. Rubin and Ho [22] found that the correction factor worked well for the dataset of Chen [15] and others (a total of 155 measurements), but a large variation was reported when applying this method to their own measurements (124 measurements). Therefore, it is not surprising to see that it showed a greater variation on a much larger dataset used in this study (a total of 1025 measurements).

Normalized Models
The four normalized thermal conductivity models generally performed better than the analytical or mixing models of Section 4.1. The ZN2016 (Table 2, Figure 2(3)) performed the best, with RMSE = 0.3 W m −1 • C −1 and NSE = 0.9, followed by ZN2015 (Table 2, Figure 2(2)), with RMSE = 0.33 W m −1 • C −1 and NSE = 0.88. It should be noted that both models, ZN2015 and ZN2016, are a modified form of the Côté and Konrad [41] model. The difference is that the ZN2015 model was developed for quartz sand while the ZN2016 was developed for a sand-clay mixture [24,25]. The normalized thermal conductivity models are among the most prevalent soil thermal conductivity models, because they generally give satisfactory estimates [19,60]. The ZN2015 and ZN2016 models extend the capability of the normalized thermal conductivity model to accurately predict λ eff of sands at a high range of λ eff (or water content range), but they still faced difficulties in λ eff prediction at low to middle ranges of λ eff . The ET1987 model (Figure 2(1)) generally overestimated the measured λ eff , which agrees with the findings of other research [18,19]. The ZY2019 model (Figure 2(4)) generally slightly underestimated measured λ eff , the same underestimates was reported by previous study [19]. The possible reason for the underestimates of the ZY2019 model is that it was developed for a sand-peat mixture, in which peat has a much smaller magnitude of thermal conductivity [26].

Linear or Non-Linear Regression Models
The CS2008 model (Figure 3(2)) is the best performing model in this category and performs equally well with model of ZN2016. The CS2008 model was developed empirically for estimating λeff of quartz sands, but it actually has a form similar to that of the geometric mean model [1,39,61]. The inclusion of λs may explain why it outperformed the other three empirical models (i.e., BB1992, LY2012 and AA2016) in this category. The BB1992 (Figure 3(1)) and LY2012 (Figure 3(3)) models have mixed performance, giving overestimates of λeff at ranges below ~1 W m −1 °C −1 and underestimates of λeff at ranges greater than ~1 W m −1 °C −1 . The LY2012 model was developed for powdery sand in tropical desert and grassland in the southwest Sahara Desert in West Africa [28]-the sand characteristics may different from sands from other parts of the world. This difference may result in the unsatisfactory performance of the LY2012 model. The AA2016 model (Figure 3(4)) generally overestimated the measured λeff as a whole, which may be attributed to the fact that AA2016 was developed based on steady-state measurements (i.e., guarded hot plates). Steady-state measurements are prone to water redistribution in unsaturated soils, which changes the thermal conductivity being measured [12,62]. This may explain why Alrtimi et al. [17] found that the models of de Vries [63], Johansen [39], Côté and Konrad [41], Lu et al. [48], Chen [15] and Haigh [16] could not be used to properly predict thermal conductivity of Tripoli sand measured with the steady-state measurements. Therefore, care should be taken when applying these models.

Linear or Non-Linear Regression Models
The CS2008 model (Figure 3(2)) is the best performing model in this category and performs equally well with model of ZN2016. The CS2008 model was developed empirically for estimating λ eff of quartz sands, but it actually has a form similar to that of the geometric mean model [1,39,61]. The inclusion of λ s may explain why it outperformed the other three empirical models (i.e., BB1992, LY2012 and AA2016) in this category. The BB1992 (Figure 3(1)) and LY2012 (Figure 3(3)) models have mixed performance, giving overestimates of λ eff at ranges below~1 W m −1 • C −1 and underestimates of λ eff at ranges greater than~1 W m −1 • C −1 . The LY2012 model was developed for powdery sand in tropical desert and grassland in the southwest Sahara Desert in West Africa [28]-the sand characteristics may different from sands from other parts of the world. This difference may result in the unsatisfactory performance of the LY2012 model. The AA2016 model (Figure 3(4)) generally overestimated the measured λ eff as a whole, which may be attributed to the fact that AA2016 was developed based on steady-state measurements (i.e., guarded hot plates). Steady-state measurements are prone to water redistribution in unsaturated soils, which changes the thermal conductivity being measured [12,62]. This may explain why Alrtimi et al. [17] found that the models of de Vries [63], Johansen [39], Côté and Konrad [41], Lu et al. [48], Chen [15] and Haigh [16] could not be used to properly predict thermal conductivity of Tripoli sand measured with the steady-state measurements. Therefore, care should be taken when applying these models.

Conclusions
In this study, 14 models used to predict the thermal conductivity of sands were reviewed and evaluated with a large compiled dataset consisting of 1025 thermal conductivity measurements measured on 62 sands from 20 studies. The results show that most of the models either underestimated or overestimated the experimental measurements. Only the models of Chen [15] (CS2008) and Zhang et al. [25] (ZN2016) gave the best estimates, with Nash-Sutcliffe efficiency = 0.9 and RMSE = 0.3 W m −1 °C −1 . These two models can potentially be used for accurately estimating the thermal conductivity of sands of different types with satisfactory performance. Although sands from various locations around the world were used in this study, care should be taken when applying the selected models for much wider applications beyond the sands tested due to the highly spatial variability in the nature of soil properties. This work would provide useful information pertaining to applications in geothermal energy resources, geo-techniques and geo-environments and in science disciplines that require the effective thermal conductivity of sands.

Conclusions
In this study, 14 models used to predict the thermal conductivity of sands were reviewed and evaluated with a large compiled dataset consisting of 1025 thermal conductivity measurements measured on 62 sands from 20 studies. The results show that most of the models either underestimated or overestimated the experimental measurements. Only the models of Chen [15] (CS2008) and Zhang et al. [25] (ZN2016) gave the best estimates, with Nash-Sutcliffe efficiency = 0.9 and RMSE = 0.3 W m −1 • C −1 . These two models can potentially be used for accurately estimating the thermal conductivity of sands of different types with satisfactory performance. Although sands from various locations around the world were used in this study, care should be taken when applying the selected models for much wider applications beyond the sands tested due to the highly spatial variability in the nature of soil properties. This work would provide useful information pertaining to applications in geothermal energy resources, geo-techniques and geo-environments and in science disciplines that require the effective thermal conductivity of sands.