Received: 28 January 2024

Accepted: 5 May 2024

DOI: 10.1002/vzj2.20350
Vadose Zone Journal

ORIGINAL ARTICLE
S p e c i a l S e c t i o n : Tr i b u t e t o R i e n v a n G e n u c h t e n , R e c i p i e n t o f t h e 2 0 2 3 Wo l f P r i z e f o r A g r i c u l t u r e

Pore unit cell network modeling of the thermal conductivity
dynamics in unsaturated sandy soils: Unveiling the role ofspanning-wetting phase cluster
Rasoul Mirghafari1
Amir Hossein Helforoosh2
Ehsan Nikooee2
Ghassem Habibagahi2
Amir Raoof3
Martinus Theodorus van Genuchten3,4
1 School of Engineering, Computing andMathematics, Oxford Brookes University,
Oxford, UK
2 Department of Civil and Environmental

Engineering, Shiraz University, Shiraz, Iran3 Department of Earth Sciences, Utrecht

University, Utrecht, The Netherlands
4 Department of Nuclear Engineering,

Abstract
As the world struggles with climate change and energy crises, understanding the roleof soil in the food–water–energy nexus becomes increasingly critical. Accurately
estimating the soil thermal conductivity drying curve is essential for assessing theimpacts of temperature on soil biota and crop growth, environmental changes due toforest fires and global warming, and for designing geo-energy extraction techniques

Federal University of Rio de Janeiro, UFRJ,
Rio de Janeiro, Brazil

such as geothermal energy piles. Existing empirical models often fail to accurately
estimate the soil thermal conductivity (TC), particularly in pendular soil moisture

Correspondence
Ehsan Nikooee, Department of Civil andEnvironmental Engineering, Shiraz
University, P.O. Box: 71348-51156, Shiraz,
Iran.
Email: enikooee@shirazu.ac.ir

regimes where they do not capture sharp changes in TC. This study introduces anovel approach using a pore unit cell network model to more accurately describe
the dynamics of TC in variably saturated soils. A quadratic parallel scheme within

Assigned to Associate Editor Jan Hopmans.

each soil pore unit cell links the TCs of solid, water, and air to the overall effective
conductivity. By modeling air invasion in the pore network model and employing
the proposed equation, we determined the unsaturated soil TC based on varying
local conductivities. The model effectively captures the significant decrease in conductivity in the pendular saturation regime, associated with the shrinkage of thespanning-wetting cluster. Quantitative analyses showed a substantial improvementin prediction accuracy compared to existing models, especially under varying moisture conditions. Our findings have significant implications for better characterizing
soil thermal and hydraulic properties, which are crucial for resource management ina changing climate and advancing geo-energy technologies.

Abbreviations: GA, genetic algorithm; MAPE, mean absolute percentage error; NMR, nuclear magnetic resonance; PNM, pore network model; PUCNM,
pore unit cell network model; RMSE, root mean squared error; SWRC, soil–water retention curve; TC, thermal conductivity; TCDC, soil thermal conductivity
drying curve.
This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided
the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
© 2024 The Authors. Vadose Zone Journal published by Wiley Periodicals LLC on behalf of Soil Science Society of America.
Vadose Zone J. 2024;e20350.
https://doi.org/10.1002/vzj2.20350

wileyonlinelibrary.com/journal/vzj2

1 of 19

1

MIRGHAFARI ET AL.

Vadose Zone Journal

INTRODUCTION

The global population is projected to exceed 9 billion by 2050
and approach 10.4 billion by the mid-2080s (McKersie, 2015;
UN, 2023). With this growth, the demand for energy, water,
and food will intensify, underscoring the need for more efficient management of our natural resources. Precise modelingof the water–food–energy nexus is vital for optimal resource
utilization and adapting to rapidly changing climate conditions. A critical component of applicable models is accurate
simulation of soil heat and water transport to allow more
precise studies of the interplay between the anthroposphere,hydrosphere, and atmosphere.
Soil thermal properties are essential to accurately modelingheat transfer processes in soils. They are important for devising appropriate irrigation models for different soils (Zhang
et al., 2019), managing soil moisture in general (Sanuade
et al., 2020), estimating soil surface evaporation, assessingorganic matter decomposition (Khaledi et al., 2023), andenhancing crop growth (Abu-Hamdeh, 2000; Abu-Hamdeh& Reeder, 2000). They are increasingly important also fordesigning energy-efficient building systems and geothermal
piles (Akrouch et al., 2016; Bayat et al., 2021; Moradshahi
et al., 2021). As climate change influences the frequency offorest fires (Flannigan et al., 2000), accurate models of soil
heat transfer are critical to predict soil temperature variations,temperature-induced changes in soil quality (e.g., organic
matter content), and carbon transfer to the atmosphere (Caonet al., 2014; Conant et al., 2011; Qi et al., 2016).
The importance of accurately estimating soil thermal
properties is highlighted in numerous applications in theenvironmental and soil sciences, geoengineering, and atmospheric sciences. Examples include soil temperature variations and their effects on soil biota and plant growth (Heinze
et al., 2017), understanding frost progression in cold regions(Hayashi et al., 2007), designing buried pipelines and pipeline
thermal backfills (Jia et al., 2019; Newson & Brunning, 2004;
Patwa et al., 2022), backfilling buried high-voltage power
cables (Ocłoń, 2021), and developing bridge deck de-icing
systems (Lei et al., 2020; Yu et al., 2016), among otherapplications. In coarse-grained soils, heat conduction is the
predominant mechanism of soil heat transfer, making the estimation of soil thermal conductivity (TC) essential for thesetypes of applications. Soils are multiphase systems involving
the air, water, and solid phases, each having significantly different TCs: ranging from 1.0 to 9.0 W m−1 K−1 for soil grains
depending on mineral type, 0.6 W m−1 K−1 for water, and0.025 W m−1 K−1 for air (Côté & Konrad, 2005; Farouki,
1981; Johansen, 1977, among others). The multiphase natureof soils deems it necessary to obtain estimates of the overall effective soil TC, denoted hereafter by λ. The TC hence
depends on many factors such as soil mineralogy, pore structure, bulk density, and the inherently transient nature of theCore Ideas∙ Spanning-wetting cluster shrinkage in the pendular regime significantly decreases soil thermal
conductivity.
∙ Soil thermal conductivity is influenced by both
the distribution and specific values of local poresaturation.
∙ Pore unit cell network modeling leads to moreaccurate estimates of the unsaturated thermal conductivity.

soil moisture regime (Gamage et al., 2019; He et al., 2021;
Lehmann et al., 2003; Mustamo et al., 2019; Zhao & Si,2019). Furthermore, recent climate changes have led to new
atmospheric and hydrological conditions, characterized byshort bursts of intense rainfall and longer drought periods(Hirabayashi et al., 2008). These changing conditions underscore the importance of accurately characterizing soil thermaland hydraulic properties and their interplay, crucial for applications such as managing earth dams in extreme weather(Vahedifard et al., 2015). In this context, different effective
TCs are expected due to changing moisture contents for thesame soil type as a consequence of varying conditions ofintense rainfall and drought.
The relationship between the soil thermal conductivity,
hereafter denoted by λ, and the soil moisture content is often
encapsulated in what is known as the soil thermal conductivity
drying curve (TCDC). Various empirical or quasi-empiricalequations for TCDC have been proposed to date (Campbell
et al., 1994; Côté & Konrad, 2005; Hu et al., 2001; Likos,
2014a; Lu et al., 2007; Somerton et al., 1974). Recently, shifts
have occurred toward the use of pore-scale models in attempts
to obtain more precise predictions of the soil hydraulic andthermal properties. For example, Likos (2014b) utilized aunit pore model to determine the thermal conductivity for
different pores based on the bundle of capillary tubes concept, along with a homogenization technique for upscaling,
to find the TCDC. Although this model marks a significant
advancement in physically based modeling of the TCDC, it
has yet to fully address different features pertinent to soil pore
structures made up of an assemblage of pore bodies and porethroats (Daneshian et al., 2021; Khaksar et al., 2013; Nikooee
et al., 2014). For a more realistic representation of pore-scalesoil moisture distributions, and consequently, accurate fluid
saturation estimates at each matric suction, a model integrating both pore body and pore throat elements is essential.
While drying processes are governed mostly by pore throat
radii, they generally have minimal volume compared to porebodies.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

2 of 19

Vadose Zone Journal

In this study, a pore network model (PNM) that accounts
for various topological and geometrical characteristics of thesoil pore structure is employed to determine the TCDC. The
model effectively estimates soil moisture contents and distributions, thus providing a direct link between the soil water
retention curve (SWRC) and the TCDC. It also elucidates
how the evolution of continuous water clusters (or alternatively discontinuous water paths, and isolated pores) leads
to three distinct saturation regimes with concomitant variations in the soil TC with saturation. We first present a brief
review of previously proposed equations for the TCDC andtheir limitations. The proposed pore scale algorithm is thendescribed in detail. Calculated TCDCs and SWRCs are subsequently compared against experimental data as well as thevarious TCDC equations. The discussion section explores theapplication scope of the method, while the study concludes
with key insights, findings, and its future outlook.

2
PREVIOUSLY PROPOSEDEQUATIONS FOR THE TCDC
Most existing models for λ are based on the concept of a normalized thermal conductivity, as first introduced by Johansen(1977). These models interpolate λ between fully saturated
and fully dry states of the soil medium, hereafter designatedby 𝜆sat and 𝜆dry , respectively. A general expression for λ is
given as:𝜆(𝑆) = 𝑓 (𝑆)(𝜆sat − 𝜆dry ) + 𝜆dry ,(1)where 𝑓 (𝑆) is a function of fluid saturation, S. The simplest
forms for 𝑓 (𝑆) include saturation or its square root (Somerton
et al., 1974) as follows:𝜆(𝑆) =√

𝑆(𝜆sat − 𝜆dry ) + 𝜆dry ,(2)

Johansen (1975, 1977) suggested a logarithmic function for𝑓 (𝑆):
𝑓 (𝑆) = 0.7 log(𝑆) + 1 medium and fine sands,(3)𝑓 (𝑆) = log(𝑆) + 1 fine - textured soils.

(4)

Côté and Konrad (2005) proposed another relationshipbased on the Kersten (1949) number, 𝐾𝑒 , after analyzinga large database consisting of various soil types, including
unfrozen and frozen gravels, sands, silts, clays, peat, andcrushed rocks. The Kersten number is defined as:
𝑓 (𝑆) = 𝐾𝑒 =𝜅𝑆,
1 + (𝜅 − 1)𝑆(5)

3 of 19

where 𝜅 is a dimensionless parameter that varies with material
type, with a recommended value of 3.55 for fine sands.
Another form of 𝑓 (𝑆) based on heat pulse measurements
involving a dozen sandy/silty soils was developed by Lu et al.
(2007):
𝑓 (𝑆) = exp[𝛼(1 − 𝑆 𝛼−1.33 )],(6)

where 𝛼 is a material parameter dependent upon soil type,with distinct values for coarse- and fine-grained soils, that is,
0.96 and 0.27, respectively (Dong et al., 2015; Likos, 2014b).
Hu et al. (2001) furthermore used a log-linear function givenby:𝑓 (𝑆) = [0.9879 + 0.181 ln(𝑆)].

(7)

Another commonly used physically based model for thethermal conductivity that accommodates various pore geometries and connections was proposed by de Vries (1963). His
model considers the effective thermal conductivity of a partially saturated porous medium as the weighted sum of thethermal conductivities of water, air, and solid phases (Ghanbarian & Daigle, 2016; Hopmans & Dane, 1986a; Tarnawski
et al., 2021). Subsequently, Campbell et al. (1994) modified
the de Vries (1963) model by formulating the weight function
for the conductivities of the endmembers (i.e., solid, water,
and air) as continuous functions of saturation.
All of the above models attempt to relate bulk saturation
to λ. However, as illustrated in Figure 1, the TCDC generally exhibits three distinct regions. A notable feature is the
sharp decrease in λ that begins at around 20% saturation,
marking the onset of the pendular regime where the liquid
phase tends to become discontinuous (Fu et al., 2024; Likos,
2014b). This results in a two-part trend in the TCDC within
this region. Unfortunately, the above equations do not adequately capture this phenomenon, which is likely influenced
by pore-scale physics and the distribution of water between
soil grains. We hypothesize, examined with the PNM, that
increasing soil suctions and the resulting air invasion reduce
the extent and number of water flow paths connecting the air
and water reservoirs in the PNM. We anticipate that primarily
mobile water and water within connected paths contribute toλ. In contrast, trapped water in pores surrounded by air-filled
pores is presumed to conduct heat close to dry pores since
their contribution is obscured by surrounding air-filled poresor pore clusters.
To the best of our knowledge, the thermal conductivity ofunsaturated soils has not yet been explored from a pore scale
network modeling perspective (i.e., incorporating a network
composed of both pore body and pore throat elements). This
study aims to correlate pore-scale saturation with the thermal conductivity of a single pore, then to upscale results inorder to obtain the effective medium λ using the PNM. The

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

MIRGHAFARI ET AL.

Vadose Zone Journal(2014). Table 1 presents the index soil properties and physicalcharacteristics of the selected soils.

3.2
PNM construction and computational
methodA PNM was used next to simulate the SWRCs and TCDCs.
The PNM along with its structural parameters was calibratedusing a genetic algorithm (GA). We note that the calibration
process only used measured SWRC data, with which we subsequently simulated the TCDCs. Our aim during calibration
hence was to achieve the best possible match between the
PNM simulations and the measured drying SWRC data. The
simulations involved the following steps:
1. Generation of the PNM using initial estimates of thestructural parameters.
2. Optimization of the structural parameters using SWRC
experimental data in conjunction with GA.
3. Simulation of the drying processes.
4. Calculation of the unsaturated TCDC.

F I G U R E 1 Schematics of the TCDC and the SWRC for a typicalsand, depicting the relationships between soil saturation (S), matricsuction, and the thermal conductivity, λ. The top schematic, based on
Likos (2014b), showcases the significant changes in λ as the soiltransitions through different saturation regimes, emphasizing the sharp
decrease when entering the pendular regime.

approach involves pore-scale modeling of air invasion into awater-saturated network and determining λ with the PNM ateach saturation level. With a pore skeleton that accurately captures the drying process of the PNM, this method is expected
to yield more accurate TCDC estimations. Additionally, weinvestigate the impact of pore-scale characteristics such as
coordination number and pore body and pore throat structural/mathematical parameters on the TCDCs. The databaseof sandy soils used in this study and the PNM approach aredetailed in the following sections.

3

MATERIALS AND METHODS

Database of hydraulic and thermal
3.1
properties of sandy soils
For this study, we compiled a database of measured values ofthe drying branch of the SWRCs and TCDCs of eight sandy
soils. These soils are primarily classified as silty sand (SM)and poorly-graded sand (SP) according to the Unified Soil
Classification System (ASTM, 1985). The data were sourced
from studies by Likos (2014a, 2014b), Yao (2014), and Oh

3.2.1

The PNM generation algorithm

For this study, we used a regular PNM representing the pore
structure of the sands. The network was modeled as a threedimensional cubic lattice, where cubic pore bodies at thenodes are interconnected by cylindrical throats. The dimension of the PNM was determined by evaluating errors in
predictions of the measured λ values (Figure 2). Formulationsof the mean absolute percentage error (MAPE) and root meansquared error (RMSE) performance criteria provided in thisfigure are introduced in Table 2. A network comprising 40 ×
40 × 40 pores proved sufficient for capturing the experimental
λ data and providing a representative elementary volume.
The simulated PNM was considered connected to a wetting
phase (i.e., water) reservoir at the bottom and a nonwetting fluid (i.e., air phase) reservoir at the top. No-flow
conditions were imposed along the lateral boundaries. A lognormal distribution function, consistent with previous studies(Daneshian et al., 2021; Raoof & Hassanizadeh, 2012; Rostami et al., 2015), was employed to generate pore body radii
described by the following equation:𝑓 (𝑥|𝜇, 𝜎) =)((ln(𝑥) − 𝜇)2
,
exp −
2𝜎 2
𝑥𝜎 2𝜋
1
√

(8)

where x denotes pore radii and μ and σ represent the mean andstandard deviation of the lognormal distribution, respectively.
To determine the relationship between neighboring poresand their throat radius, the Acharya et al. (2004) formulas were used as illustrated in Figure 3 and expressed by

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

4 of 19

TA B L E 1

5 of 19

Vadose Zone Journal

Summary of index (physical) and thermal properties of the eight sandy soils considered in this study.

Soil namenp

𝑮𝒔

𝑫𝟓𝟎 (mm)𝑪𝒖

𝑪𝒄

𝝀𝐝𝐫𝐲 (𝑾 ∕𝒎 × 𝑲)

𝝀𝐬𝐚𝐭 (𝑾 ∕𝒎 × 𝑲)

Soil type

Reference

F-75 Sand

0.41

2.65

0.20

2.32

1.08

0.23

2.42

SP

Likos (2014b)

12/20 Sand

0.31

2.65

1.04

1.33

1.00

0.31

2.97

SP

Likos (2014b)

River sand

0.38

2.64

0.33

1.40

0.96

0.32

2.29

SP

Likos (2014a)

SM1

0.39

2.68

0.11

2.55

1.12

0.48

2.20

SM

Yao (2014)

SM2

0.40

2.75

0.13

3.20

1.01

0.38

2.17

SM

Yao (2014)

SP

0.37

2.66

0.21

1.60

1.03

0.40

2.43

SP

Oh (2014)

SP2

0.39

2.68

0.29

1.83

0.94

0.50

2.43

SP

Yao (2014)

SP5

0.37

2.66

0.21

1.59

1.03

0.40

2.43

SP

Yao (2014)

Abbreviations: np : soil porosity, Gs : specific gravity, D50 : median grain size (Dx , grain size at x % passing), Cu : uniformity coefficient (D60 /D10 ), Cc : coefficient ofcurvature (D30 2 /D60 × D10 ), SM: silty sand, SP: poorly-graded sand.

TA B L E 2

Measures of performance used in this study.

Error functions
Root mean squared
error (RMSE)

Equations√ ∑(𝜆exp − 𝜆com )2
RMSE = 𝑁1 𝑁
𝑖=1Mean absolute
percentage error(MAPE)

MAPE = 𝑁1

Coefficient ofdetermination (𝑅2 )

𝑅2 = 1 −∑𝑁|𝜆exp −𝜆com |

𝑖=1

𝜆exp

∑

(𝜆exp −𝜆com )2∑
̄2(𝜆exp −𝜆)Note: λexp and λcom denote experimental and computed values of thermal conductivity data, respectively, and 𝜆̄ specifies the average value of the observed
datapoints.

F I G U R E 2 Variations in the (a) mean absolute percentage error(MAPE) and (b) root mean squared error (RMSE) performancemeasures of the unsaturated soil thermal conductivity drying curve(TCDC) of River sand versus size of the network (i.e., the number ofpores in the pore network model [PNM]). The red line depicts 64,000
pores (40 × 40 × 40 pores).

Equation (9):
(
)
(1∕𝑛)
(1∕𝑛) (−𝑛)𝑟′ 𝑖𝑗 = 𝑝𝑖 𝑝𝑗 𝑝𝑖+ 𝑝𝑗
𝑅′ sin(𝜋∕4)𝑝𝑖 = (1−𝑅′𝑖 cos(𝜋∕4))𝑛𝑖

𝑅′ 𝑗 sin(𝜋∕4)𝑝𝑗 = (1−𝑅′ cos(𝜋∕4))𝑛𝑗𝑛 > 0, 𝑅′′′
𝑖 = 𝑅𝑖 ∕𝑑𝑖𝑗 , 𝑅 𝑗 = 𝑅𝑗 ∕𝑑𝑖𝑗 , 𝑟 𝑖𝑗 = 𝑟𝑖𝑗 ∕𝑑𝑖𝑗

(9)

F I G U R E 3 Schematic of a pore-to-pore connection within thepore network model (PNM), illustrating the geometric and topologicalsetup of the cubic pore bodies and cylindrical pore throats.

where 𝑅𝑖 and 𝑅𝑗 are the radii of two adjacent pore bodies
and 𝑟𝑖𝑗 is the radius of the cylindrical throat connecting them.
The 𝑝𝑖 and 𝑝𝑗 values denote the distance between the centerof pores i and j and the respective end of the arc tangent to
the pore body and intersecting the pore throat (see Acharya
et al. [2004] for further details). The distance between the pore
body centers is denoted by 𝑑𝑖𝑗 , and 𝑛 is a non-dimensional
aspect ratio parameter governing the pore throat radius.
The coordination number Cn , representing the maximum
and minimum number of neighboring pores, to which a porecan connect, varied between 6 and 2. In this study, thecoordination number for various pores is controlled by thecoordination number parameter, 𝑐, which ranges from 0 to 1.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

MIRGHAFARI ET AL.

Vadose Zone JournalThis parameter acts as a threshold: a random number assigned
to each throat determines if the throat is retained or removed
based on its comparison to the parameter 𝑐. Thus, 𝑐 influences the coordination number distribution. Consequently, thestructural PNM parameters include the aspect ratio parameter(𝑛), the coordination number parameter (𝑐), the mean (𝜇) andstandard deviation (𝜎) of the pore radii distribution, and theminimum (𝑅min ) and maximum (𝑅max ) pore body radii.

𝑆𝑖 = −)(
2𝑇𝑠1.
ln 1 −
6.83
𝜓𝑅𝑖

(11)

4. Comparison of simulated and experimental SWRCs fora set of PNM structural parameters. The performance ofthe network is evaluated using the MAPE performance
criterion:
MAPE = 𝑓 (𝑛, 𝑐, 𝜇, 𝜎, 𝑅min , 𝑅max )3.2.2

PNM simulation processes

Simulation of SWRC and calibration of network structuralparameters
The simulation process begins with an initial guess of thestructural parameters, leading to the generation of a preliminary PNM. To accurately describe the soil’s pore structure,a calibration and optimization technique is applied. The
approach involves using the main drying branch of the SWRC
experimental data for calibration. The structural parameters
are iteratively modified until the error in predictions of the
SWRC data falls below a predefined threshold. The GA is thenemployed to optimize the structural characteristics, ensuringthat the simulated PNM closely aligns with the SWRC data.
The SWRC main drying branch is simulated as follows:
1. Generation of an initial PNM using assumed structural parameters and setting the top, bottom, and lateralboundary conditions.
2. Invasion of the non-wetting (air) phase into pores adjacentto the top boundary connected to the air reservoir. An airpressure is imposed along the top of the boundary of thepore network, while the bottom boundary is connected to
the wetting phase reservoir. By stepwise increasing the airpressure, incremental matric suctions are applied to facilitate air invasion while maintaining a continuous drainage
path. A pore is drained when air invades its connecting
throat. The criterion for air invasion into a throat is definedby the Young–Laplace equation, which calculates the air
entrance pressure for each throat during drainage:
𝑃𝑐 throat =2𝑇𝑠 cos 𝜃
,
𝑟𝑖𝑗

(10)

where 𝑟𝑖𝑗 denotes the throat radius linking pores 𝑖 and 𝑗,
𝑇𝑠 represents surface tension, and 𝜃 specifies the receding
contact angle.
3. Calculation of the degree of PNM saturation at each suction level (ψ) by identifying drained pores. Local saturation
Si of each cubic pore body i at different suction values
is determined using a function proposed by Joekar-Niasar
et al. (2010) in which Ri stands for the pore body radius:|𝑁 |1 ∑ ||𝑆exp − 𝑆pnm ||
,
=
𝑁 𝑖=1
𝑆exp

(12)

where 𝑆exp and 𝑆𝑝𝑛𝑚 are the experimental and simulateddegrees of saturation, respectively, and 𝑁 is the number ofdata points (suction levels).
5. Adjustment of structural parameters via GA, continuing until the error between experimental and simulated
SWRCs is minimized, or a maximum number of algorithm generations is reached. The optimized structuralparameters are then used to create the presumably mostrepresentative PNM.

Determination of TCDC
The thermal conductivity (λ) of each pore unit within the
PNM is calculated using a quadratic parallel expression thatcombines the λ values of the solid, air, and water phases.
The equation for calculating the λ of each soil unit embracing a pore is as follows (detailed calculations are provided inAppendix A):[
𝜆(𝑆) =]
𝜅𝑆(𝜆sat − 𝜆dry ) + 𝜆dry ,1 + (𝜅 − 1)𝑆𝜅 = 13,(13)
where S denotes the local saturation of the pore, while theparameter κ is a function of soil structure, and the value
obtained (κ = 13) is therefore suggested for sandy soils. The
average value of the local pore thermal conductivity is thenderived from the entire network and reported at each saturation level. It is important to note that the saturation ofpores (and pore clusters) not belonging to the spanning cluster of the wetting phase (which connects the air and water
reservoirs) is considered non-contributory to the thermal conduction. These isolated pores are surrounded by air-filled
pores, effectively masking their heat conduction. For suchpore unit cells, we assign the λ value of a dry pore unitcell.
To identify the continuous path of the wetting phase,including various isolated wetting phase clusters and themain spanning cluster, a search algorithm commonly usedin computer vision and robotics is employed (Haralick &15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

6 of 19

TA B L E 3

Best-matched pore network model (PNM) structural parameters.

Soil

n

c

μ

σ

Rmin (µm)

Rmax (µm)

F-75 Sand

0.727

0.165

4.45

2.696

14.701112.798

12/20 Sand

2.565

0.986

6.946

1.127

0.059

227.31

River sand

0.015

0.116

5.505

2.696

13.929112.798

SM1

2.022

0.405

3.771

1.57

5.985

154.668

SM2

0.919

0.039

3.648

1.005

4.352

275.071

SP

0.015

0.663

3.84

2.696

18.657112.798

SP2

0.538

0.924

4.034

2.189

13.701103.941

SP5

0.015

0.663

4.034

2.59

4.748

111.798

Shapiro, 1992; Sedgewick, 1998). The algorithm makes it
possible to detect connected objects in an image.
Since the model calculates the local thermal conductivityof a soil unit cell embracing a pore, the term pore unit cell
network model (PUCNM) is used hereafter to better convey
the concept behind the calculations. To summarize, at first,
saturations of different pores within the PNM are calculated,
and then local λ values of the soil units are obtained basedon the local saturations and upscaled by averaging throughthe network. The efficacy and appropriateness of the proposed
approach for determining TCDC are assessed using various
measures of performance as outlined in Table 2.

4

7 of 19

Vadose Zone Journal

RESULTS

The structural parameters of the best-matching PNMs forour eight soils are detailed in Table 3. Figure 4 shows thesimulated SWRCs along with the experimental data. The optimized PNMs accurately replicated the SWRCs of the eight
sandy soils, demonstrating the effectiveness of the selectedstructural parameters.
Figures 5 and 6 present the eight soils-simulated TCDCs
as obtained with the pore unit cell network model (PUCNM)
against those derived with the literature equations and experimental data. The figures show that the PUCNM estimatesclosely align with the experimental TCDC data, especiallyin the pendular regime. The other models within this regime
deviate significantly from the experimental values, often
underestimating them. The performance of the PUCNM is
further contrasted in Figure 5a,b with the pore-scale model of
Likos (2014b). The PUCNM approach provided much better
TCDC estimates compared to the bundle of capillary tubes
model used by Likos (2014b). In the capillary tube model,
pores of the same size drain simultaneously at a specificmatric suction, without considering the presence of both porebodies and pore throats, thereby neglecting connectivity in
the PNM. The PUCNM, on the other hand, incorporates boththe pore body and pore throat elements and thus accounts
for the connectivity between different pore unit cells. This

allows for more precise estimates of local saturation andthe thermal conductivity in especially the pendular moistureregime where λ has its lowest values (Figure 6). The PUCNM
formulation, with its more comprehensive account of the natural soil geometrical and topological features, hence offers amore accurate representation than the simpler capillary tubes
model.
To better understand the superior performance of the
PUCNM pore-scale network model across all saturation
regimes, especially in the pendular regime, an in-depth analysis of the wetting phase saturation evolution within the PNM
was conducted. This included examining the local λ valuesof pores located within the spanning-wetting cluster of thesamples (i.e., the continuous wetting phase cluster connecting the air and water reservoirs), and the ratio of the number ofpores in the spanning-wetting phase cluster to the total number of pores (denoted as RPSC). These aspects are illustrated
in Figure 7.
First and on top in Figure 7 are the simulated and experimental data of the TCDC of River sand, as well as the portionsof all pores of the network that belong to the spanning-wetting
cluster are depicted (denoted by RPSC). In the lower partof the figure, at each saturation, are two illustrations. One
is the soil PNM at a certain soil saturation, with the corresponding local pore saturations presented in different colors.
The second network shows the local λ values of pores in the
spanning-wetting cluster and the extent of the spanning cluster in size. This illustration reveals a significant decrease inboth the extent of the spanning-wetting phase cluster and thepercentage of pores within this cluster at around 20% saturation, which marks the transition to the pendular regime,
approximately corresponding to the residual water content in
the classical SWRC models of Brooks and Corey (1964) andvan Genuchten (1980). Consequently, the effective medium
value of λ, calculated as the average of the local λ values in the PUCNM, is markedly reduced. Eventually, thecontinuous and spanning cluster of the wetting phase connecting the two reservoirs diminishes, causing the PUCNMderived thermal conductivities to approach those of a dry
soil.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

MIRGHAFARI ET AL.

Vadose Zone Journal

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

F I G U R E 4 Pore unit cell network model (PUCNM) simulated and experimental data points of soil water retention curves (SWRCs) of various
sandy soils: (a) F-75 sand, (b) 12/20 sand, (c) River sand, (d) SM1, (e) SM2, (f) SP, (g) SP2, and (h) SP5. PUCNM, pore unit cell network model.

Figure 8 presents, for the different approaches, the performance metrics in terms of the R2 , RMSE, and the MAPE.
The analysis reveals that none of the models provided perfectly accurate TCDC estimates for all soils. The model by
Hu et al. (2001) tended to overestimate experimental data,
while the other models by Somerton et al. (1974), Johansen(1975), Côté and Konrad (2005), and Lu et al. (2007) generally underestimated the thermal conductivity. However, the
PUCNM demonstrated a notable ability to capture the general trend of the TCDC experimental data across the different
flow regimes, exhibiting the highest R2 and lowest RMSE and
MAPE values compared to the other approaches.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

8 of 19

Vadose Zone Journal

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

9 of 19F I G U R E 5 Calculated thermal conductivity drying curves (TCDCs) versus experimental data for soils (a) F-75 sand, (b) 12/20 sand, (c) River
sand, (d) SM1, (e) SM2, (f) SP, (g) SP2, and (h) SP5 and comparisons with other models. PUCNM, pore unit cell network model. Likos (2014)
stands for pore scale simulations of Likos (2014b).

5

DISCUSSION OF THE RESULTS

Figure 7 illustrates that in coarse-grained soils, such as sandy
soils, distinct regimes of the thermal conductivity varying
with saturation are observable. The sharp decrease in theunsaturated thermal conductivity is notably marked by thetransition from the funicular to the pendular regime. The
PUCNM simulations show that the connected spanningcluster of the pore network begins to shrink at this transition.
Specifically, the proportion of total pores in the network

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

Vadose Zone Journal

F I G U R E 6 Predicted λ values of the thermal conductivity drying
curves (TCDCs) versus experimental data for the eight sandy soils,
estimations through other formulas and the pore-scale model of Likos(2014b).

that belong to this cluster, denoted by RPSC in Figure 7,
starts to decrease. Once there is no connected wetting phasecluster, the thermal conductivity values remain unchanged.
This direct correlation between heat flow and water flowis also evident in evaporation phenomena (Fu et al., 2024).
Coussot (2000) utilized the concept of water connectivity
loss in the pendular regime to define the second regime (often
called the falling-rate regime) in evaporation curves. As air
becomes trapped in the pores, water persists in the form ofwater bridges, causing the heat-conducting water pathways
to narrow and eventually disconnect. Cho and Santamarina(2001) also identified the boundary between the pendular
and funicular regimes as a critical point that defines different
regimes for electrical conduction and chemical diffusion in
unsaturated porous media.
In our study, we employed Equation (13) (determined froma quadratic parallel scheme as presented in Appendix A)
to correlate the local saturation of pores with the thermal
conductivity of a unit pore. The equation is somewhat similar to the function used by Côté and Konrad (2005). This
raises the question what the outcome would be if the samefunction was applied to the global saturation of the network,
instead of accounting for the local saturation of each pore andthen upscaling. Figure 9 addresses this question by comparing PUCNM results with the predictions of Côté and Konrad(2005) for different values of the 𝜅 parameter in Equation (5),applied to the bulk network saturation. Results clearly showthat the PUCNM calculations align more closely with theexperimental data across all values of the 𝜅 parameter. This
underscores the inability of traditional models to capture thecritical role of the spanning cluster and its influence on the

MIRGHAFARI ET AL.

network thermal conductivity, thus highlighting the necessityof PNM. It also emphasizes the importance of considering the
local saturation values of individual pores and their impact on
a macroscopic parameter like the soil thermal conductivity.
In scenarios with varying grain size distributions, local poresaturations can differ significantly from the overall network
saturation. Consequently, the average λ value derived fromvariable local saturations can substantially deviate from thethermal conductivity estimated directly from a single (bulk)
network saturation value. Mathematically, for a nonlinear
function f (e.g., for mapping saturation to λ), the expectedvalue of the function for a set of input values (like saturation)
does not equate to the function’s value when applied to theaverage of the inputs, that tis, 𝐸(𝑓 (𝑠𝑖 )) ≠ 𝑓 (𝐸(𝑠𝑖 )). This discrepancy highlights the critical role of pore-scale models in
upscaling and accurately estimating macroscopic unsaturated
properties such as the hydraulic and thermal conductivities.
Figure 10 explores the impact of the PNM structural parameters, notably the coordination number parameter (c) and the
aspect ratio parameter (n), on thermal conductivity. Results
suggest that a decreased threshold probability (indicating
lower connectivity) leads to a higher saturation threshold forthe onset of the decline in λ. In other words, less connected
networks experience a reduction in the spanning cluster size
at higher saturations. Similarly, increasing the aspect ratio toreflect smaller throat sizes implies a more gradual change innetwork saturation with suction until some high suction isreached. At this point, more extensive pore drainage occurs,
causing the spanning cluster to retract and λ to decrease at
higher saturations.
As can be seen, the aspect ratio parameter (n) significantly
affects the unsaturated soil thermal conductivity. It should be
noted that in the proposed approach, the aspect ratio parameter
does not directly enter the calculations; instead, it affects the
distribution of saturation within the soil PNM, thereby influencing the local thermal conductivities of soil pore unit cells
and, ultimately, the bulk soil value of λ.

6
MODEL LIMITATIONS, SCOPE OF
APPLICATIONS, AND FUTURE OUTLOOKThe limitations of the proposed model, its cope of application, and potential for further developments can be discussed
in light of different physical phenomena as well as soil
characteristics, grouped into the following four items:
1. Temperature effects: The proposed model was formulated
for the thermal conductivity of unsaturated sandy soils
at ambient temperature. However, the thermal conductivity is known to vary with temperature (Campbell et al.,
1994; de Vries, 1963; Hopmans & Dane, 1986a; Sepaskhah & Boersma, 1979; Xu et al., 2019). Particularly, in

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

10 of 19

Vadose Zone Journal

11 of 19F I G U R E 7 Evolution of local saturation (refer to the PNMs illustrated in the first and third columns) and spanning-wetting phase cluster local
thermal conductivities in the PUCNM (see the networks presented in the second and fourth columns) during drying of River sand at different
saturations: (a) S = 1.0, (b) S = 0.88, (c) S = 0.54, (d) S = 0.28, (e) S = 0.14, (f) S = 0.1, (g) S = 0.08, and (h) S = 0.03. RPSC gives the ratio of thenumber of pores in the spanning-wetting phase cluster to the total number of pores. PUCNM, pore unit cell network model.

medium- to fine-textured soils, vapor transfer has beenreported to occur at high temperatures (Sepaskhah &Boersma, 1979), while latent heat transfer may becomea significant portion of total heat transfer at relatively
low water contents and high temperatures (Tarnawski &Leong, 2000). To more precisely estimate the thermal conductivity at higher temperatures, these phenomena mustbe accounted for (i.e., at high temperatures, the vapor
conductivity should be substituted for air conductivity
when determining the unit cell thermal conductivity in the
PUCNM). Therefore, the development of the model for

predicting the unsaturated thermal conductivity at higher
temperatures, especially for fine-grained soils, requires
careful incorporation of relevant physics (such as vapor
transfer and temperature-dependent changes in soil structure, such as the double layer thickness, which indirectly
affects the SWRC parameters as discussed by Nimmo
and Miller (1986), Hopmans and Dane (1986b), Grant
and Salehzadeh (1996), François and Ettahiri (2012), andKhaksar et al. (2015), among others).
2. Hysteresis effects: The PUCNM technique was evaluated for determining the unsaturated thermal conductivity

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

Vadose Zone Journal

F I G U R E 8 Comparison of the R2 , root mean squared error(RMSE), and mean absolute percentage error (MAPE) performancemeasures of the various models in describing the thermal conductivity
drying curves (TCDCs) of eight sandy soils. NA, not available. Likos(2014) stands for pore scale simulations of Likos (2014b).

MIRGHAFARI ET AL.

during drying. The model effectively revealed how wetting phase connectivity and the extent of the connected
wetting phase cluster influence thermal conductivity
changes across different capillary regimes, especially during the transition from a funicular to a pendular regime. We
found that the formation of entrapped air and a reduction
in network connectivity lead to a decrease in unsaturated
thermal conductivity. However, the amount of entrapped
air varies from a soil undergoing drying to the same
soil at the same saturation reached through wetting or
scanning processes. This variation results in different air
and water saturations across the soil pore network, a
key aspect of the hysteresis phenomenon (Chen et al.,2015; Kaluarachchi & Parker, 1987; Lenhard & Parker,
1987; Poulovassilis, 1970; Stonestrom & Rubin, 1989).
Sharma and Mohamed (2003) noted that the amount ofentrapped air is determined by the reversal degree of saturation during a drying-wetting cycle, which specifies thesaturation level at the cycle’s turning point. Variations
in air entrapment significantly impact water flow (Fayer& Hillel, 1986a, 1986b; Stauffer & Dracos, 1984) and,due to a strong correlation with air-filled porosity, also
affect the unsaturated thermal conductivity (Ochsner et al.,2001). Experimental studies have shown that the degreeof hysteresis in unsaturated thermal conductivity curves
heavily depends on soil structure and soil type (Ali et al.,2014; Cardoso et al., 2018; Farouki, 1981; Nguyen et al.,2017, 2018; Smits et al., 2010; Tang et al., 2008; Vu
et al., 2023; Xu et al., 2020; Zhou et al., 2019). Soils
with different grain size distributions and types, including those prepared in undisturbed or compacted states,
exhibit varying pore structures, leading to different levelsof air entrapment during drying-wetting cycles and, consequently, varying degrees of hysteresis in the unsaturated
thermal conductivity (λ). Given that PUCNM incorporates
elements of pore bodies and pore throats, which govern
the processes of drying and wetting, this model appears
to be a promising tool for investigating unsaturated thermal conductivity in future studies on hysteresis effects.
However, the model must incorporate fluid phase displacement mechanisms that govern wetting and scanning
processes, such as cooperative pore filling and snap-off
during imbibition (Mirghafari et al., 2024; Rostami et al.,2015).
3. The role of soil structure, connectivity and tortuosity: We validated the proposed model using data fromcoarse-grained soils, specifically sandy soils. Fine-grainedsoils, especially those that are structured and aggregated,
generally possess a more complex soil pore network
comprising two distinct pore categories: intra- and interaggregate pores. Studies have demonstrated that heattransfer decreases as aggregate size increases (Usowicz
et al., 2013; Zhou et al., 2019). Consequently, the path

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

12 of 19

(b)

(c)

(d)

(e)

(f)

(g)

(h)

F I G U R E 9 Performance of pore unit cell network model (PUCNM) compared to direct estimation of the thermal conductivity using the Côté
and Konrad (2005) model with different κ values for soils (a) F-75 sand, (b) 12/20 sand, (c) River sand, (d) SM1, (e) SM2, (f) SP, (g) SP2, and (h)SP5.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License(a)13 of 19
Vadose Zone Journal

MIRGHAFARI ET AL.

MIRGHAFARI ET AL.

Vadose Zone Journal

(a)

(b)

F I G U R E 1 0 Sensitivity analysis of simulated thermal conductivity drying curves (TCDCs) for different values of pore unit cell networkmodel (PUCNM) structural parameters (a) the coordination number parameter (c) and (b) the aspect ratio parameter (n).

F I G U R E 1 1 Comparison of the average of 𝑅2 and root mean squared error (RMSE) of all equations and the pore unit cell network model(PUCNM) model used in this study.

of heat transfer and its tortuosity are influenced by thesize of aggregates, or more precisely, by the distribution of intra- and interaggregate pore sizes. To extend
the proposed PUCNM technique to other soil types, it isessential to account for variations in soil pore structure
and to consider a wider range of pore sizes, extending
down to nanometer-scale pores. The development of multiscale pore-scale models, a practice already underway infields such as petroleum engineering for studying twophase flow in shales (Alafnan & Yucel Akkutlu, 2020;
Cui et al., 2019), offers valuable insights for investigating the physical properties of fine-grained structured
soils.
4. Surface water adsorption: To adapt the PUCNM for finegrained soils, it is crucial to account for the presence ofadsorbed water. The current model, designed primarily
for sandy soils, considers only mobile and capillary water
within a single pore class. In fine-grained soils, distinguishing between adsorbed (including hydration water)
and capillary (mobile) water is essential (Dong et al., 2015;
Lu, 2016). Hydration water does not enhance connectiv-

ity among grains, which explains the pronounced stasis
behavior or the extended tail in the unsaturated thermal
conductivity profiles of fine-grained soils. As the finegrained soil content, especially clay, increases, so does thelength of this tail and the corresponding threshold watercontent. Exceeding this threshold leads to an increase in
the mobile water content, thereby establishing connectivity between grains and significantly enhancing the overallsoil thermal conductivity. The soil’s plasticity, the tendency of minerals to adsorb water, and the specific surface
area of fine soil particles govern this behavior. Recent
advancements in nuclear magnetic resonance (NMR) haveenabled the quantification of adsorbed water in unsaturated porous media such as wood and soil (Cocusse
et al., 2022; Eizaguirre et al., 2024). These advancements
in dynamic NMR relaxometry, combined with dynamic
X-ray (or neutron) tomography together with thermal
imaging, provide insights into the local heat paths and theirdevelopment within the soil porous medium. Until theseadvanced experimental techniques become more accessible, advancements in pore-scale models such as PUCNM

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

14 of 19are vital for understanding the dynamics of the thermal
conductivity in unsaturated porous media.

7

CONCLUSIONSIn this study, we utilized a PUCNM to estimate the TCDC
of coarse-textured soils. The model demonstrated promising
accuracy in estimating the unsaturated TCDC of sands as indicated by the average R2 and RMSE values across all studied
soils (Figure 11). The PUCNM results showed an impressive R2 of 95%, indicating highly accurate TCDC estimations
across the full range of saturation. The model excelled especially in the pendular regime (saturations below 20%), whereother models failed to accurately describe the experimental
data. The enhanced performance of the proposed model isdue to its ability to differentiate between disconnected wetting clusters and pores situated on the spanning (continuous)
water cluster within the soil PNM. Since the thermal conductivity (λ) of isolated and disconnected wetting clusters is
obscured by surrounding air-filled pores, the extent of thespanning cluster and the percentage of pores belonging tothis cluster significantly control the overall unsaturated soil
λ value. A marked shrinkage of the spanning cluster in thependular regime leads here then to a substantial decrease in
the unsaturated λ.
The PUCNM results also underscored the importance ofpore-level saturation and the resulting local λ. The modelleads to more accurate global λ estimates compared to models using bulk network saturation. The relevance of pore-scalemodeling, as demonstrated by PUCNM, in upscaling macroscopic characteristics such as the TCDC is evident from ourfindings.
Future research should focus on incorporating various porescale displacement processes to more accurately capture fluid
saturation and thermal conductivity changes throughout theentire hydraulic cycle (i.e., to enable simulating unsaturated
thermal properties in not only drying but also in wetting
curves). Using advanced imaging techniques, such as thermal imaging combined with neutron tomography and X-ray
computed tomography, together with dynamic NMR relaxometry, will provide more insights into heat conduction
processes in the unsaturated zone across different saturation
regimes. With the advent of 4D tomography and thermal
imaging, unique opportunities arise for observing these processes in real time. Until such advancements are widespread,the insights gained from PNM can greatly enhance ourunderstanding of unsaturated soil pore-scale characteristics
and their impact on macroscopic properties. The PUCNM
approach in this study was applied to very coarse-textured(sandy) soils. Its broader applicability, such as to fine-grainedsoils, requires the incorporation of corner and film flow inthe pore throats and the use of multi-structural soil PNMs

Vadose Zone Journal

15 of 19that consider a broader range of pore sizes (e.g., interand intra-aggregate pores typical of aggregated or structured
media).
AU T H O R C O N T R I B U T I O N S
Rasoul Mirghafari: Investigation; methodology; software;
validation; visualization; writing—original draft; writing—review and editing. Amir Hossein Helforoosh: Investigation;
software; validation. Ehsan Nikooee: Conceptualization;
methodology; software; supervision; writing—original draft;writing—review and editing. Ghassem Habibagahi: Supervision; writing—review and editing. Amir Raoof: Methodology; writing—review and editing. Martinus Theodorus van
Genuchten: Writing—review and editing.
C O N F L I C T O F I N T E R E ST STAT E M E N T
The authors declare no conflicts of interest.
ORCID
Ehsan Nikooee https://orcid.org/0000-0003-2365-1132
Martinus Theodorus van Genuchten https://orcid.org/
0000-0003-1654-8858
REFERENCES
Abu-Hamdeh, N. H. (2000). Effect of tillage treatments on soil thermal
conductivity for some Jordanian clay loam and loam soils. Soil and
Tillage Research, 56(3–4), 145–151. https://doi.org/10.1016/S01671987(00)00129-X
Abu-Hamdeh, N. H., & Reeder, R. C. (2000). Soil thermal conductivity
effects of density, moisture, salt concentration, and organic matter.
Soil Science Society of America Journal, 64(4), 1285–1290. https://doi.org/10.2136/sssaj2000.6441285x
Acharya, R. C., van der Zee, S. E. A. T. M., & Leijnse, A. (2004).
Porosity–permeability properties generated with a new 2-parameter
3D hydraulic pore-network model for consolidated and unconsolidated porous media. Advances in Water Resources, 27, 707–723.
https://doi.org/10.1016/j.advwatres.2004.05.002
Akrouch, G. A., Sánchez, M., & Briaud, J. L. (2016). An experimental,
analytical and numerical study on the thermal efficiency of energy
piles in unsaturated soils. Computers and Geotechnics, 71, 207–220.
https://doi.org/10.1016/j.compgeo.2015.08.009
Alafnan, S., & Yucel Akkutlu, I. (2020). The transport mechanisms anddynamics of shale via multiscale multiphysics pore network modeling. Journal of Energy Resources Technology, 142(9), 093003. https://doi.org/10.1115/1.4046522
Ali, A., Mohamed, M., Aal, M. A., Schellart, A., & Tait, S. (2014). Thermal and hydraulic properties of sandy soils during drying and wetting
cycles. New frontiers in geotechnical engineering (pp. 129–138).
ASCE.
ASTM. (1985). Classification of soils for engineering purposes: Annual
book of ASTM standards (D 2487-83, Vol. 04). American Society for
Testing and Materials.
Bayat, H., Ebrahimzadeh, G., & Mohanty, B. P. (2021). Investigating thecapability of estimating soil thermal conductivity using topographical attributes for the Southern Great Plains, USA. Soil and Tillage
Research, 206, 104811. https://doi.org/10.1016/j.still.2020.104811

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

Vadose Zone Journal

Brooks, R., & Corey, T. (1964). Hydraulic properties of porous media(Hydrology Papers No. 3). Colorado State University Fort Collins
Colorado.
Campbell, G. S., Jungbauer, J. D., Jr, Bidlake, W. R., & Hungerford,
R. D. (1994). Predicting the effect of temperature on soil thermal
conductivity. Soil Science, 158(5), 307–313. https://doi.org/10.1097/
00010694-199411000-00001
Caon, L., Vallejo, V. R., Ritsema, C. J., & Geissen, V. (2014). Effects ofwildfire on soil nutrients in Mediterranean ecosystems. Earth-Science
Reviews, 139, 47–58. https://doi.org/10.1016/j.earscirev.2014.09.001
Cardoso, R., Sousa, M., & Vieira, A. (2018). Suction effect on
the thermal properties of compacted kaolin. In 7th International conference on unsaturated soils. International
Society for Soil Mechanics and Geotechnical Engineering.
https://www.issmge.org/publications/publication/suction-effecton-the-thermal-properties-of-compacted-kaolin
Chen, P., Wei, C., & Ma, T. (2015). Analytical model of soil-water
characteristics considering the effect of air entrapment. International
Journal of Geomechanics, 15(6), 04014102. https://doi.org/10.1061/(ASCE)GM.1943-5622.0000462
Cho, G. C., & Santamarina, J. C. (2001). Unsaturated particulate materials—Particle-level studies. Journal of Geotechnical andGeoenvironmental Engineering, 127(1), 84–96. https://doi.org/10.
1061/(asce)1090-0241(2001)127:1(84)Cocusse, M., Rosales, M., Maillet, B., Sidi-Boulenouar, R., Julien,E., Caré, S., & Coussot, P. (2022). Two-step diffusion in cellularhygroscopic (vascular plant-like) materials. Science Advances, 8(19),eabm7830. https://doi.org/10.1126/sciadv.abm7830
Conant, R. T., Ryan, M. G., Ågren, G. I., Birge, H. E., Davidson, E. A.,
Eliasson, P. E., Evans, S. E., Frey, S. D., Giardina, C. P., Hopkins,
F. M., Hyvönen, R., Kirschbaum, M. U. F., Lavallee, J. M., Leifeld,
J., Parton, W. J., Megan Steinweg, J., Wallenstein, M. D., Martin
Wetterstedt, J. Å., & Bradford, M. A. (2011). Temperature and soilorganic matter decomposition rates–synthesis of current knowledge
and a way forward. Global Change Biology, 17(11), 3392–3404.
https://doi.org/10.1111/j.1365-2486.2011.02496.x
Cosenza, P., Guérin, R., & Tabbagh, A. (2003). Relationship between
thermal conductivity and water content of soils using numerical modelling. European Journal of Soil Science, 54(3), 581–588. https://doi.
org/10.1046/j.1365-2389.2003.00539.x
Côté, J., & Konrad, J. M. (2005). A generalized thermal conductivity
model for soils and construction materials. Canadian Geotechnical
Journal, 42(2), 443–458. https://doi.org/10.1139/t04-106
Coussot, P. (2000). Scaling approach of the convective drying of aporous medium. The European Physical Journal B—Condensed Matter and Complex Systems, 15(3), 557–566. https://doi.org/10.1007/
s100510051160
Cui, R., Feng, Q., Chen, H., Zhang, W., & Wang, S. (2019). Multiscalerandom pore network modeling of oil-water two-phase slip flow inshale matrix. Journal of Petroleum Science and Engineering, 175, 46–
59. https://doi.org/10.1016/j.petrol.2018.12.026
Daneshian, B., Habibagahi, G., & Nikooee, E. (2021). Determination of
unsaturated hydraulic conductivity of sandy soils: A new pore network approach. Acta Geotechnica, 16, 449–466. https://doi.org/10.
1007/s11440-020-01088-3
de Vries, D. A. (1963). Thermal properties of soils. In W. R. Van Wijk
(Ed.), Physics of plant environment (pp. 210–235). North-Holland.
Dong, Y., McCartney, J. S., & Lu, N. (2015). Critical review of thermal
conductivity models for unsaturated soils. Geotechnical and Geolog-

MIRGHAFARI ET AL.

ical Engineering, 33, 207–221. https://doi.org/10.1007/s10706-0159843-2
Eizaguirre, P., Tang, A. M., Maillet, B., Sidi-Boulenouar, R., Talandier,J., Pereira, J.-M., Vu, M. N., Chabot, B., Dangla, P., Bornert, M.,& Aimedieu, P. (2024). Exploring two regimes of water mobility inunsaturated expansive clay using NMR relaxometry. Applied Clay
Science, 251, 107324. https://doi.org/10.1016/j.clay.2024.107324
Farouki, O. (1981). Thermal properties of soils. US Army Corps ofEngineers, Cold Regions Research and Engineering Laboratory.
Fayer, M. J., & Hillel, D. (1986a). Air encapsulation: I. Measurement ina field soil. Soil Science Society of America Journal, 50(3), 568–572.
https://doi.org/10.2136/sssaj1986.03615995005000030005x
Fayer, M. J., & Hillel, D. (1986b). Air encapsulation: II. Profile water
storage and shallow water table fluctuations. Soil Science Society ofAmerica Journal, 50(3), 572–577. https://doi.org/10.2136/sssaj1986.
03615995005000030006x
Fedaoui, K., Madani, S., & KANIT, T. (2016). Prediction of effective
thermal conductivity of heterogeneous random multi-phase composites. UPB Scientific Bulletin Series D: Mechanical Engineering,78(3), 91–100.
Flannigan, M. D., Stocks, B. J., & Wotton, B. M. (2000). Climate change
and forest fires. Science of the Total Environment, 262(3), 221–229.
https://doi.org/10.1016/S0048-9697(00)00524-6
François, B., & Ettahiri, S. (2012). Role of the soil mineralogy on thetemperature dependence of the water retention curve. Unsaturated
soils: Research and applications (Vol. 1, pp. 173–178). Springer.
Fu, Y., Ghanbarian, B., Horton, R., & Heitman, J. (2024). New insights
into the correlation between soil thermal conductivity and water
retention in unsaturated soils. Vadose Zone Journal, 23(1), e20297.
https://doi.org/10.1002/vzj2.20297
Gamage, D. N. V., Biswas, A., & Strachan, I. B. (2019). Spatial variability of soil thermal properties and their relationships with physicalproperties at field scale. Soil and Tillage Research, 193, 50–58.
https://doi.org/10.1016/j.still.2019.05.012
Ghanbarian, B., & Daigle, H. (2016). Thermal conductivity inporous media: Percolation-based effective-medium approximation.
Water Resources Research, 52, 295–314. https://doi.org/10.1002/
2015WR017236
Grant, S. A., & Salehzadeh, A. (1996). Calculation of temperature effects
on wetting coefficients of porous solids and their capillary pressurefunctions. Water Resources Research, 32(2), 261–270. https://doi.org/
10.1029/95WR02915
Haralick, R. M., & Shapiro, L. G. (1992). Computer and robot vision.
Addison-Wesley.
Hayashi, M., Goeller, N., Quinton, W. L., & Wright, N. (2007). A simple
heat-conduction method for simulating the frost-table depth in hydrological models. Hydrological Processes: An International Journal,
21(19), 2610–2622. https://doi.org/10.1002/hyp.6792
He, H., Liu, L., Dyck, M., Si, B., & Lv, J. (2021). Modelling dry soil
thermal conductivity. Soil and Tillage Research, 213, 105093. https://doi.org/10.1016/j.still.2021.105093
Heinze, J., Gensch, S., Weber, E., & Joshi, J. (2017). Soil temperature modifies effects of soil biota on plant growth. Journal of Plant
Ecology, 10(5), 808–821.
Hirabayashi, Y., Kanae, S., Emori, S., Oki, T., & Kimoto, M. (2008).
Global projections of changing risks of floods and droughts in achanging climate. Hydrological Sciences Journal, 53(4), 754–772.
https://doi.org/10.1623/hysj.53.4.754

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

16 of 19

Hopmans, J. W., & Dane, J. H. (1986a). Thermal conductivity of two
porous media as a function of water content, temperature, and density.
Soil Science, 142(4), 187–195. https://doi.org/10.1097/00010694198610000-00001
Hopmans, J. W., & Dane, J. H. (1986b). Temperature dependence of soil water retention curves. Soil Science Society ofAmerica Journal, 50(3), 562–567. https://doi.org/10.2136/sssaj1986.
03615995005000030004x
Hu, X. J., Du, J. H., Lei, S. Y., & Wang, B. X. (2001). A model forthe thermal conductivity of unconsolidated porous media based on
capillary pressure-saturation relation. International Journal of Heatand Mass Transfer, 44(1), 247–251. https://doi.org/10.1016/S00179310(00)00079-X
Jia, G. S., Tao, Z. Y., Meng, X. Z., Ma, C. F., Chai, J. C., & Jin, L. W.
(2019). Review of effective thermal conductivity models of rock-soil
for geothermal energy applications. Geothermics, 77, 1–11. https://doi.org/10.1016/j.geothermics.2018.08.001
Joekar-Niasar, V., Hassanizadeh, S. M., & Dahle, H. K. (2010). Nonequilibrium effects in capillarity and interfacial area in two-phase
flow: Dynamic pore-network modelling. Journal of Fluid Mechanics,
655, 38–71. https://doi.org/10.1017/S0022112010000704
Johansen, O. (1975). Varmeledningsevne av jordarter [Doctoral dissertation, Norwegian Institute of Technology].
Johansen, O. (1977). Thermal conductivity of soils. Cold Regions
Research and Engineering Laboratory.
Kaluarachchi, J. J., & Parker, J. C. (1987). Effects of hysteresis
with air entrapment on water flow in the unsaturated zone. Water
Resources Research, 23(10), 1967–1976. https://doi.org/10.1029/
WR023i010p01967
Kersten, M. S. (1949). Thermal properties of soils (Bulletin No. 28).
Engineering Experiment Station, Institute of Technology, Universityof Minnesota.
Khaksar, H., Habibagahi, G., & Nikooee, E. (2013). SWRC modeling in unsaturated soils: A pore network approach. In Poromechanics V: Proceedings of the Fifth Biot Conference on Poromechanics (pp. 1570–1579). ASCE. https://ascelibrary.org/doi/10.1061/
9780784412992.186
Khaksar, H., Nikooee, E., Habibagahi, G., & Hassanizadeh, S. M. (2015).
Modeling temperature effects on soil water retention curve: A pore
network approach. Geotechnical engineering for infrastructure anddevelopment (pp. 3267–3272). ICE. https://www.icevirtuallibrary.
com/doi/abs/10.1680/ecsmge.60678.vol6.511
Khaledi, S., Delbari, M., Galavi, H., Bagheri, H., & Chari, M. M. (2023).
Effects of biochar particle size, biochar application rate, and moisture content on thermal properties of an unsaturated sandy loam soil.
Soil and Tillage Research, 226, 105579. https://doi.org/10.1016/j.still.
2022.105579
Lehmann, P., Stahli, M., Papritz, A., Gygi, A., & Fluhler, H. (2003).
A fractal approach to model soil structure and to calculate thermal conductivity of soils. Transport in Porous Media, 52, 313–332.
https://doi.org/10.1023/A:1023595527999
Lei, G., Yu, X., Li, T., Habibzadeh-Bigdarvish, O., Wang, X., Mrinal, M.,& Luo, C. (2020). Feasibility study of a new attached multi-loop CO2
heat pipe for bridge deck de-icing using geothermal energy. Journal of
Cleaner Production, 275, 123160. https://doi.org/10.1016/j.jclepro.
2020.123160
Lenhard, R. J., & Parker, J. C. (1987). A model for hysteretic constitutive relations governing multiphase flow: 2. Permeability-saturation

Vadose Zone Journal

17 of 19

relations. Water Resources Research, 23(12), 2197–2206. https://doi.
org/10.1029/WR023i012p02197
Likos, W. J. (2014a). Modeling thermal conductivity dryout curves fromsoil-water characteristic curves. Journal of Geotechnical and Geoenvironmental Engineering, 140(5), 04013056. https://doi.org/10.1061/(ASCE)GT.1943-5606.0001078
Likos, W. J. (2014b). Pore-scale model for thermal conductivity ofunsaturated sand. Geotechnical and Geological Engineering, 33(2),179–192. https://doi.org/10.1007/s10706-014-9744-9
Lu, N. (2016). Generalized soil water retention equation for adsorption and capillarity. Journal of Geotechnical and Geoenvironmental Engineering, 142(10), 04016051. https://doi.org/10.1061/(ASCE)
GT.1943-5606.0001524
Lu, S., Ren, T., Gong, Y., & Horton, R. (2007). An improved model
for predicting soil thermal conductivity from water content at room
temperature. Soil Science Society of America Journal, 71(1), 8–14.
https://doi.org/10.2136/sssaj2006.0041
McKersie, B. (2015). Planning for food security in a changing climate.
Journal of Experimental Botany, 66(12), 3435–3450. https://doi.org/
10.1093/jxb/eru547
Mirghafari, R., Nikooee, E., Raoof, A., & Habibagahi, G. (2024). Determination of a pedotransfer function for specific air water interfacial
area in sandy soils: A pore network-informed multigene genetic programming approach. Vadose Zone Journal, e20352. https://doi.org/
10.1002/vzj2.20352
Moradshahi, A., Faizal, M., Bouazza, A., & McCartney, J. S. (2021).
Effect of nearby piles and soil properties on thermal behaviour of afield-scale energy pile. Canadian Geotechnical Journal, 58(9), 1351–
1364. https://doi.org/10.1139/cgj-2020-0353
Mustamo, P., Ronkanen, A. K., Berglund, Ö., Berglund, K., & Kløve, B.
(2019). Thermal conductivity of unfrozen and partially frozen managed peat soils. Soil and Tillage Research, 191, 245–255. https://doi.
org/10.1016/j.still.2019.02.017
Newson, T. A., & Brunning, P. (2004). Thermal conductivity of deepwater offshore sediments. International Journal of Offshore and Polar
Engineering, 14(4), 310–314.
Nguyen, V. T., Heindl, H., Frost, J. D., Pereira, J. M., & Tang, A.
M. (2018). Suction and thermal conductivity of unsaturated loessfrom Northern France. UNSAT. https://www.issmge.org/uploads/
publications/62/63/4-17.UNSAT2018_085.pdf
Nguyen, V. T., Heindl, H., Pereira, J. M., Tang, A. M., & Frost, J. D.
(2017). Water retention and thermal conductivity of a natural unsaturated loess. Géotechnique Letters, 7(4), 286–291. https://doi.org/10.
1680/jgele.17.00037
Nikooee, E., Habibagahi, G., Khaksar, H., Hassanizadeh, S., & Raoof,
A. (2014). Pore network modeling of unsaturated soils: Fundamentals,
recent advancements and future perspectives. Numerical Methods inGeotechnical Engineering-NUMGE, 2014, 1007–1012. https://doi.
org/10.1201/b17017-179
Nimmo, J. R., & Miller, E. E. (1986). The temperature dependence ofisothermal moisture vs. potential characteristics of soils. Soil Science Society of America Journal, 50(5), 1105–1113. https://doi.org/
10.2136/sssaj1986.03615995005000050004x
Ochsner, T. E., Horton, R., & Ren, T. (2001). A new perspective on soil
thermal properties. Soil Science Society of America Journal, 65(6),1641–1647. https://doi.org/10.2136/sssaj2001.1641
Ocłoń, P. (2021). The effect of soil thermal conductivity and cableampacity on the thermal performance and material costs of under-

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.

Vadose Zone Journal

ground transmission line. Energy, 231, 120803. https://doi.org/10.
1016/j.energy.2021.120803
Oh, H. (2014). Thermal resistivity dry-out curves for thirteen sandy soils[Master’s thesis, University of Wisconsin-Madison].
Patwa, D., Bordoloi, U., Dubey, A. A., Ravi, K., Sekharan, S., & Kalita,
P. (2022). Energy-efficient biochar production for thermal backfill
applications. Science of the Total Environment, 833, 155253. https://doi.org/10.1016/j.scitotenv.2022.155253
Population. United Nations. (2023). https://www.un.org/en/globalissues/population
Poulovassilis, A. (1970). The effect of the entrapped air on the hysteresis curves of a porous body and on its hydraulic conductivity.
Soil Science, 109(3), 154–162. https://doi.org/10.1097/00010694197003000-00003
Qi, R., Li, J., Lin, Z., Li, Z., Li, Y., Yang, X., Zhang, J., & Zhao, B.
(2016). Temperature effects on soil organic carbon, soil labile organiccarbon fractions, and soil enzyme activities under long-term fertilization regimes. Applied Soil Ecology, 102, 36–45. https://doi.org/10.
1016/j.apsoil.2016.02.004
Rao, S. E., Ray, L., Khan, T., & Ravi, G. (2022). Thermal conductivity,
density and porosity of sedimentary and metamorphic rocks from theLower and Higher Himalaya, Western Himalaya, India. GeophysicalJournal International, 231(1), 459–473. https://doi.org/10.1093/gji/ggac176
Raoof, A., & Hassanizadeh, S. M. (2012). A new formulation for porenetwork modeling of two-phase flow. Water Resources Research,48(1), W01514. https://doi.org/10.1029/2010WR010180
Rostami, A., Habibagahi, G., Ajdari, M., & Nikooee, E. (2015). Pore
network investigation on hysteresis phenomena and influence of stressstate on the SWRC. International Journal of Geomechanics, 15(5),04014072. https://doi.org/10.1061/(ASCE)GM.1943-5622.0000315
Sanuade, O. A., Hassan, A. M., Akanji, A. O., Olaojo, A. A., Oladunjoye,
M. A., & Abdulraheem, A. (2020). New empirical equation to estimate the soil moisture content based on thermal properties using
machine learning techniques. Arabian Journal of Geosciences, 13,
1–14. https://doi.org/10.1007/s12517-020-05375-x
Sedgewick, R. (1998). Algorithms in C (3rd ed.). Addison-Wesley.
Sepaskhah, A. R., & Boersma, L. (1979). Thermal conductivity of soils
as a function of temperature and water content. Soil Science Societyof America Journal, 43, 439–444. https://doi.org/10.2136/sssaj1979.
03615995004300030003x
Sharma, R. S., & Mohamed, M. H. (2003). An experimental investigationof LNAPL migration in an unsaturated/saturated sand. Engineering
Geology, 70(3–4), 305–313. https://doi.org/10.1016/S0013-7952(03)00098-X
Smits, K. M., Sakaki, T., Limsuwat, A., & Illangasekare, T. H. (2010).
Thermal conductivity of sands under varying moisture and porosity in drainage–wetting cycles. Vadose Zone Journal, 9(1), 172–180.
https://doi.org/10.2136/vzj2009.0095
Somerton, W. H., El-Shaarani, A. H., & Mobarak, S. M. (1974, April
4–5). High temperature behavior of rocks associated with geothermal type reservoirs [Paper presentation]. SPE California Regional
Meeting, San Francisco, CA, USA. https://doi.org/10.2118/4897-MS
Stauffer, F., & Dracos, T. (1984). Local infiltration into layered soil andresponse of the water table, experiment and simulation. Frontiers inhydrology (pp. 228–242). Water Resources Publication.

MIRGHAFARI ET AL.

Stonestrom, D. A., & Rubin, J. (1989). Water content dependence oftrapped air in two soils. Water Resources Research, 25(9), 1947–1958.
https://doi.org/10.1029/WR025i009p01947
Tang, A. M., Cui, Y. J., & Le, T. T. (2008). A study on the thermal conductivity of compacted bentonites. Applied Clay Science, 41(3–4),
181–189. https://doi.org/10.1016/j.clay.2007.11.001
Tarnawski, V. R., & Leong, W. H. (2000). Thermal conductivityof soils at very low moisture content and moderate temperatures.
Transport in Porous Media, 41(2), 137–147. https://doi.org/10.1023/
A:1006738727206
Tarnawski, V. R., Wagner, B., Leong, W. H., McCombie, M., Coppa,
P., & Bovesecchi, G. (2021). Soil thermal conductivity model by deVries: Re-examination and validation analysis. European Journal of
Soil Science, 72(5), 1940–1953. https://doi.org/10.1111/ejss.13117
Usowicz, B., Lipiec, J., Usowicz, J. B., & Marczewski, W. (2013).
Effects of aggregate size on soil thermal conductivity: Comparison of measured and model-predicted data. International Journal of
Heat and Mass Transfer, 57(2), 536–541. https://doi.org/10.1016/j.
ijheatmasstransfer.2012.10.067
Vahedifard, F., AghaKouchak, A., & Robinson, J. D. (2015). Drought
threatens California’s levees. Science, 349(6250), 799–799. https://doi.org/10.1126/science.349.6250.799-a
Van Genuchten, M. T. (1980). A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.
03615995004400050002x
Vu, Q. H., Pereira, J. M., & Tang, A. M. (2023). Effect of clay content on
the thermal conductivity of unfrozen and frozen sandy soils. International Journal of Heat and Mass Transfer, 206, 123923. https://doi.
org/10.1016/j.ijheatmasstransfer.2023.123923
Xu, Y., Zeng, Z., & Lv, H. (2019). Effect of temperature on thermal
conductivity of lateritic clays over a wide temperature range. International Journal of Heat and Mass Transfer, 138, 562–570. https://doi.org/10.1016/j.ijheatmasstransfer.2019.04.077
Xu, Y., Zeng, Z., & Lv, H. (2020). Comparative study on thermal properties of undisturbed and compacted lateritic soils subjected to drying
and wetting. Engineering Geology, 277, 105800. https://doi.org/10.
1016/j.enggeo.2020.105800
Yao, J. (2014). Measurement and modeling of thermal propertiesof coarse-grained soils [Master’s thesis, University of WisconsinMadison].
Yu, X., Zhang, N., Pradhan, A., & Puppala, A. J. (2016). Geothermalenergy for bridge deck and pavement deicing—A brief review. GeoChicago, 2016, 598–609.
Zhang, J., Hou, J., Zhang, H., Meng, C., Zhang, X., & Wei, C. (2019).
Low soil temperature inhibits yield of rice under drip irrigation. Journal of Soil Science and Plant Nutrition, 19, 228–236. https://doi.org/
10.1007/s42729-019-0008-x
Zhao, Y., & Si, B. (2019). Thermal properties of sandy and peat soils
under unfrozen and frozen conditions. Soil and Tillage Research, 189,
64–72. https://doi.org/10.1016/j.still.2018.12.026
Zhou, Y., Yan, C., Tang, A. M., Duan, C., & Dong, S. (2019). Mesoscopic prediction on the effective thermal conductivity of unsaturated
clayey soils with double porosity system. International Journal of
Heat and Mass Transfer, 130, 747–756. https://doi.org/10.1016/j.
ijheatmasstransfer.2018.11.001

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

18 of 19

Vadose Zone Journal

19 of 19

How to cite this article: Mirghafari, R., Helforoosh,A. H., Nikooee, E., Habibagahi, G., Raoof, A., & van
Genuchten, M. T. (2024). Pore unit cell networkmodeling of the thermal conductivity dynamics inunsaturated sandy soils: Unveiling the role ofspanning-wetting phase cluster. Vadose Zone Journal,
e20350. https://doi.org/10.1002/vzj2.20350

APPENDIX A
To determine the thermal conductivity (λ) of each soil unit
cell embracing a pore, the λ values of the solid, air, and waterphases were combined. Given the unknown value of λ of different sand types, an approximate mean value for sandstones(λs = 4.23 W/m K as reported by Rao et al. [2022]) wasused. The thermal conductivities of water and air at room temperature were taken as λw = 0.56 and λa = 0.026 W/m K,
respectively.
According to Cosenza et al. (2003), a quadratic parallel
expression is recommended for complex multiphase materials to link the effective thermal conductivity of a unit cell
containing solid, air, and water phases to their respective bulk
values, through their volume fractions (φ). This approach
has later been validated against 3D numerical heat flow
simulations, showing promising results (Fedaoui et al.,2016). In this approach, the λ values of the solid, air, andwater phases in a unit cell are combined using the following
quadratic parallel expression:

𝜆unit cell =

( 3
∑
𝑖=1

)2
(√ )
𝜑𝑖
𝜆𝑖
.

(A1)

In our study, to reduce computational costs of running
the PUCNM and to directly use soil thermal properties,F I G U R E A 1 Comparison of effective unit cell thermal
conductivities (λ) values (calculated using a quadratic parallel scheme)with values obtained from the proposed semi-empirical equation foreight sands and at different saturations (S = 0.1–1).

this formulation was applied to eight sand types across different saturation levels (0.1–1). We found that the valuesof λ (unit cell) at these saturation levels can be empirically related to λdry and λsat of the soils as follows(Figure A1):]
𝜅𝑆
(𝜆sat − 𝜆dry ) + 𝜆dry ,𝜆(𝑆) =
1 + (𝜅 − 1)𝑆[𝜅 = 13.

(A2)
This equation is similar in form to the model suggested by
Côté and Konrad (2005). However, in this study, we derived
the equation based on a quadratic parallel conceptualizationof solid, water, and air in a single pore unit cell. The equation
connects local pore saturation (S) to the local conductivity (λ)
of a pore unit cell, which could otherwise be obtained directly
through quadratic parallel conceptualization.

15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20350 by Martinus van Genuchten - CAPES , Wiley Online Library on [09/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

MIRGHAFARI ET AL.