Empirical and theoretical models for prediction of soil thermal conductivity: a review and critical assessment

On the basis of the results of thermal conductivity measurements of 19 soils, Kersten (1949) proposed empirical equations to determine the thermal conductivity of the soil medium on the basis of moisture content w [%] and the dry density ρ_d [g/cm³]. For fine-grained soils (silt, clay, etc.) the thermal conductivity of the soil medium can be estimated using the following expression: (1)λ=0.1442(0.9logw−0.2)100.6243ρd\lambda = 0.1442\left( {0.9\log w - 0.2} \right){10^{0.6243{\rho _d}}}

The analogous formula for coarse-grained soils is as follows: (2)λ=0.1442(0.7logw+0.4)100.6243ρd\lambda = 0.1442\left( {0.7\log w - 0.4} \right){10^{0.6243{\rho _d}}}\

Johansen (1975) proposed the determination of thermal conductivity of partially saturated soil medium on the basis of thermal conductivity of dry soil λ^dry and in the state of full water saturation λ^sat. For that purpose, a normalized (dimensionless) thermal conductivity, also known as the Kersten number K_e, was introduced. According to Johansen’s (1975) approach, thermal conductivity of the soil is governed by the following equation: (3)λ=(λsat−λdry)⋅Ke+λdry\lambda = \left( {{\lambda ^{{sat}}} - {\lambda ^{dry}}} \right) \cdot {K_e} + {\lambda ^{dry}}

Based on the data contained in Kersten’s work (1949), Johansen proposed the correlation formulas between the Kersten number K_e and the degree of saturation S_r. The conductivity of the dry soil was proposed to be evaluated using an empirical equation based on the dry density of the soil. The value of λ^sat was proposed to be computed using the geometrical mean approach, which involves the conductivity of water and of soil solids. All the equations constituting Johansen method are summarized in Table 1.

The Johansen method is recognized in the literature as providing one of the best predictions of the thermal conductivity of soil medium (Farouki, 1981; Dong et al. 2015; Różański, 2018; Zhang & Wang, 2017). In the following years, various researchers modified the Johansen method. The changes mainly concern the introduction of a slightly improved form of the function describing the K_e − S_r relationship and the proposal of new approach for determining the thermal conductivity of dry soils.

Côté and Konrad (2005a; 2005b) modified the Kersten number by introducing an additional parameter κ, reflecting the type of soil. In addition, a modified formula for calculation of λ^dry value was proposed – it was assumed to be dependent not only on the porosity but also on two additional parameters reflecting the type of soil and grain shape (χ, η). Another modification of Kersten number and conductivity of dry soil was proposed by Lu et al. (2007). According to the work of Lu et al. (2007), the Kersten number is governed by the exponential function with parameter α whose value depends on whether we are dealing with coarse- or fine-grained soil. The empirical formula for calculation of λ^dry is the linear function of soil porosity n. The most recent proposal for a functional dependence between Kersten number and S_r is the one proposed in the work of He et al. (2017). All formulas for the Johansen method and its subsequent modifications are collected in Table 1. The sign ‘−’ means that the authors did not modify the parameter or they did not refer to it at all in their work.

Donazzi et al. (1979) studied the thermal and hydrological characteristics of soils surrounded by underground installations. Based on the results of the laboratory tests, they proposed the following correlation formula to determine the thermal conductivity of the soil, which can be applied to all types of soil: (4)λ=λswnλs(1−n)exp(−3.08n(1−Sr)2).\lambda = \lambda _w^n\lambda _s^{(1 - n)}\exp \left( { - 3.08n{{\left( {1 - {S_r}} \right)}^2}} \right).

Chen (2008) carried out tests on the thermal conductivity of quartz sands in different saturation states. On this basis, the empirical equation of thermal conductivity as a function of porosity and degree of saturation was proposed: (5)λ=λwnλs(1−n)[(1−b)Sr+b]cn\lambda = \lambda _w^n\lambda _s^{(1 - n)}{\left[ {(1 - b){S_r} + b} \right]^{cn}}

Based on the results of thermal conductivity measurements, the model parameters were determined as: b = 0.0022 and c = 0.78. Equation (5) has a good accuracy in estimating the conductivity of sandy soils with relatively high quartz content.

In the work of Lu et al. (2014), the following empirical equation for thermal conductivity at various saturation states was proposed: (6)λ=λdry+exp[β−(Srn)−α] Sr>0\lambda = {\lambda ^{dry}} + \exp \left[ {\beta - {{\left( {{S_r}n} \right)}^{ - \alpha }}} \right]\,\,\,\,\,\,\,\,\,\,\,\,\,\,{S_r} > 0

The conductivity of dry soil is evaluated using the approach proposed by Lu et al. (2007) – see Table 1. Parameters α and β were estimated to be dependent on the contents of sand and clay separates as well as on the soil dry density, that is: (7)α=0.67ϕcl+0.24\alpha = 0.67{\varphi _{cl}} + 0.24(8)β=1.97ϕSa+1.87ρd−1.36ϕSaρd−0.95\beta = 1.97{\varphi _{Sa}} + 1.87{\rho _d} - 1.36{\varphi _{Sa}}{\rho _d} - 0.95 where ϕ_Cl and ϕ_Sa are the contents of, respectively, clay and sand separates. The value of the density ρ_d should be taken in units [g/cm ³ ]. Note, that this model gives the possibility to estimate the thermal conductivity of the same soils, but in different states of compaction.

In 1912, Wiener (Wiener, 1912) proposed an upper and lower bound for the thermal conductivity of a porous medium consisting of three separate phases: solid, water and air (Fig. 1).

Figure 1

The lower (λ^L) and upper (λ^U) bounds of thermal conductivity are obtained, respectively, for the series and parallel arrangement of the individual phases, relative to the direction of the heat flux, that is: (9)λL=(∑iϕiλi)−1{\lambda ^L} = {\left( {\sum\limits_i {\frac{{{\varphi _i}}}{{{\lambda _i}}}} } \right)^{ - 1}}(10)λU=∑iϕiλi{\lambda ^U} = \sum\limits_i {{\varphi _i}{\lambda _i}} where ϕ_i and λ_i are respectively the volume fraction and thermal conductivity of phase i. It is commonly known (e.g., Różański, 2018) that these bounds can only be treated as a rough estimate of soil thermal conductivity. It is due to the fact that for large contrast in conductivities, the bounds are ‘very wide’. In the case of unsaturated soil, the contrast between skeleton, water and air is of the two orders of magnitude.

In the work of Mickley (1951), it was assumed that elementary volume of a soil is a unit cube with separated ‘channels’, filled with air and/or water (Fig. 2). In the case of dry soil (S_r=0), after assuming the direction of heat flow (Fig. 2a), the layer system proposed by Mickley implies that heat flows simultaneously through four separate columns: (1) the air-filled space between the outer walls of a cube of sectional area a² and unit length, (2) the soil skeleton between the outer walls of a cube of sectional area (1−a²) and unit length, and (3) and (4) two columns (air and soil skeleton) in series, each with sectional area a(1−a²) and length a and (1−a), respectively, for air and soil skeleton.

Figure 2

The thermal conductivities at different states of saturation are evaluated using the following relations (Mickley, 1951): –

dry soil(11)λdry=λaa2+λs(1−a)2+λsλa(2a−2a2)λsa+λa(1−a),{\lambda ^{{dry}}} = {\lambda _a}{a^2} + {\lambda _s}{\left( {1 - a} \right)^2} + \frac{{{\lambda _s}{\lambda _a}\left( {2a - 2{a^2}} \right)}}{{{\lambda _s}a + {\lambda _a}(1 - a)}},where λ_a is the thermal conductivity of air,

The geometric variables a and c are to be determined from the soil porosity n and the degree of saturation S_r by solving the equations below: (14)3a2−2a3=n3{a^2} - 2{a^3} = n and (15)3c2−2c3=n(1−Sr).3{c^2} - 2{c^3} = n(1 - {S_r}).

Gemant (1952) proposed a model based on an idealized grain shape – a cube with a side a, the three walls of which form ‘pyramids’ with a square base. It was assumed that grains are in contact with each other only by the top of ‘pyramids’ and the water accumulates around these contact points forming the so-called thermal bridges between the particles (Fig. 3).

Figure 3

Gemant (1952) proposed the following formula for the conductivity of the soil (16)1λ=[(1−a)a]43arctan[λs−λwλw]12(h2)13[λw(λs−λw)]12+1−zwλsaf{b2a}\frac{1}{\lambda } = \frac{{{{\left[ {\frac{{(1 - a)}}{a}} \right]}^{\frac{4}{3}}}{arctan}{{\left[ {\frac{{{\lambda _s} - {\lambda _w}}}{{{\lambda _w}}}} \right]}^{\frac{1}{2}}}}}{{{{\left( {\frac{h}{2}} \right)}^{\frac{1}{3}}}{{\left[ {{\lambda _w}({\lambda _s} - {\lambda _w})} \right]}^{\frac{1}{2}}}}} + \frac{{1 - {z_w}}}{{{\lambda _s}a}}f\left\{ {\frac{{{b^2}}}{a}} \right\} where λ_s and λ_w [Wm⁻¹K⁻¹ ] are the thermal conductivities of the soil skeleton and water, respectively. Remaining parameters are expressed by following relations: (17)a=4.869ρda = 4.869\sqrt {{\rho _d}} (18)h=9.988⋅10−3ρdw−h0h = 9.988 \cdot {10^{ - 3}}{\rho _d} - {h_0}(19)zw=(1−aa)23(h2)13{z_w} = {\left( {\frac{{1 - a}}{a}} \right)^{\frac{2}{3}}}{\left( {\frac{h}{2}} \right)^{\frac{1}{3}}}(20)b2=(a1−a)23(h2)23{b^2} = {\left( {\frac{a}{{1 - a}}} \right)^{\frac{2}{3}}}{\left( {\frac{h}{2}} \right)^{\frac{2}{3}}}

In equations above, h[−] represents the amount of water accumulated around contact points and h₀ is the water adsorbed as a film on the surfaces of solid particles. Gemant assumed that h₀ decreases as the temperature increases following the theory of water adsorption of glass (Farouki, 1981), for example: for T=20°C, h₀=0.01. Analogically, namely using a separate nomogram, the values of f{b²/a} are estimated. In the work of Różański (2018), the following approximation of this function was introduced: (21)f{b2a}≈(0.976+0.2355 ln(b2a))−1f\left\{ {\frac{{{b^2}}}{a}} \right\} \approx {\left( {0.976 + 0.2355\,{\text{ln}}\left( {\frac{{{b^2}}}{a}} \right)} \right)^{ - 1}}

De Vries model (de Vries, 1963) is based on the Maxwell equation defining the effective electrical conductivity of a mixture of spherical particles dispersed in a continuous fluid (Eucken, 1932). De Vries assumed that the soil medium could be modelled as a mixture of N ellipsoidal particles dispersed in a continuous matrix (being water or air, depending on the degree of saturation). The thermal conductivity of the soil medium is expressed as: (22)λ=ϕ0λ0+∑i=1NFiϕiλiϕ0+∑i=1NFiϕi\lambda = \frac{{{\varphi _0}{\lambda _0} + \sum\nolimits_{i = 1}^N {{F_i}{\varphi _i}{\lambda _i}} }}{{{\varphi _0} + \sum\nolimits_{i = 1}^N {{F_i}{\varphi _i}} }} where the subscripts 0 and i indicate the parameters characterizing, respectively, a continuous matrix or ellipsoidal inclusions dispersed within it. Moreover, the weighting parameters F_i are calculated using following formula: (23)Fi=23[1+(λi−λ0λ0)ga]−1+13[1+(λi−λ0λ0)(1−2ga)]−1{F_i} = \frac{2}{3}{\left[ {1 + \left( {\frac{{{\lambda _i} - {\lambda _0}}}{{{\lambda _0}}}} \right){g_a}} \right]^{ - 1}} + \frac{1}{3}{\left[ {1 + \left( {\frac{{{\lambda _i} - {\lambda _0}}}{{{\lambda _0}}}} \right)\left( {1 - 2{g_a}} \right)} \right]^{ - 1}} where g_a is the parameter defining the shape of the ellipsoidal inclusion (de Vries proposed g_a=0.125).

Gori (1983) proposed a theoretical model to determine the thermal conductivity of soils in different saturation states. Gori assumed that the soil (as a multi-phase medium) consists of periodic cubic cells with a ‘centrally’ located skeleton (solid) surrounded by air (dry soil) or water (saturated soil) or both air and water (partially saturated soil). A geometry of the representative unit cell, at various states of saturation, is presented in Fig. 4.

Figure 4

Due to the fact that the calculation procedure using this model, compared to most other proposals, is relatively difficult, the use of this model in engineering practice is quite limited. In addition, the disadvantage of the model is that while the formulation itself results directly from theoretical considerations, the water content values determining the geometry of the periodicity cell (Fig. 4) were determined empirically depending on the type of soil. In the computational sense, determination of the thermal conductivity of dry soils seems to be relatively simple (Fig. 4a). This prediction depends only on the porosity of the soil and the thermal conductivities of the air and soil skeleton and is given by the following equation: (24)1λdry=ψ−1λaψ+ψλa(ψ2−1)+λs\frac{1}{{{\lambda ^{{dry}}}}} = \frac{{\psi - 1}}{{{\lambda _a}\psi }} + \frac{\psi }{{{\lambda _a}\left( {{\psi ^2} - 1} \right) + {\lambda _s}}} where (25)ψ=n1−n3\psi = \sqrt[3]{{\frac{n}{{1 - n}}}}

In the work of Tong et al. (2009), the thermal conductivity of a multicomponent porous medium was considered. A new model based on Wiener as well as Hashin-Shtrikmann

For a macroscopically isotropic medium the bounds of Hashin-Shtrikman provide the narrowest range of estimates, which are expressed only in terms of volume fractions of constituents (e.g., Łydżba, 2011).

Haigh (2012) model is used to determine the thermal conductivity of sandy soils in different water saturation states. The equation describing thermal conductivity was determined on the basis of the assumed simplified geometry of the soil microstructure. A problem of 1D heat flow between two identical grains of spherical shape and radius R was considered (Fig. 5). It was assumed that grains are placed in a cylindrical cell of radius R and length R(1 + ξ) – the contact between grains is only when ξ = 0.

Figure 5

The equation describing the thermal conductivity of the soil is expressed by the following formula: (27)λλs=2(1+ξ)2{αw(1−αw)2ln[(1+ξ)+(αw−1)xξ+αw]+αa(1−αa)2ln[(1+ξ)(1+ξ)+(αa−1)x]}++2(1+ξ)(1−αw)(1−αa)[(αw−αa)x−(1−αa)αw]\begin{array}{*{20}{c}} {\frac{\lambda }{{{\lambda _s}}} = 2{{\left( {1 + \xi } \right)}^2}\left\{ {\frac{{{\alpha _w}}}{{{{\left( {1 - {\alpha _w}} \right)}^2}}}{\text{ln}}\left[ {\frac{{(1 + \xi ) + ({\alpha _w} - 1)x}}{{\xi + {\alpha _w}}}} \right] + \frac{{{\alpha _a}}}{{{{\left( {1 - {\alpha _a}} \right)}^2}}}{\text{ln}}\left[ {\frac{{(1 + \xi )}}{{(1 + \xi ) + ({\alpha _a} - 1)x}}} \right]} \right\} + } \\ { + \frac{{2(1 + \xi )}}{{(1 - {\alpha _w})(1 - {\alpha _a})}}\left[ {({\alpha _w} - {\alpha _a})x - (1 - {\alpha _a}){\alpha _w}} \right]} \end{array} where α_w=λ_w/λ_s and α_a=λ_a/λ_s. In addition, ξ and x are the values depending on porosity and saturation degree and they should be computed using the formulas below: (28)ξ=3n−13(1−n)\xi = \frac{{3n - 1}}{{3(1 - n)}}(29)x=(1+ξ2)(1+cosϑ−3sinϑ)x = \left( {\frac{{1 + \xi }}{2}} \right)(1 + {\text{cos}}\vartheta - \sqrt 3 {\text{sin}}\vartheta )(30)cosϑ=2(1+3ξ)(1−Sr)−(1+ξ)3(1+ξ)3{\text{cos}}\vartheta = \frac{{2(1 + 3\xi )(1 - {S_r}) - {{\left( {1 + \xi } \right)}^3}}}{{{{(1 + \xi )}^3}}}