# Effect of soil saturation and grain size on coupled hydrothermal flow in fine sands based on X-ray µ CT imaging

Kui Liu, Fleur A Loveridge, Richard Boardman, William Powrie

► To cite this version:

Kui Liu, Fleur A Loveridge, Richard Boardman, William Powrie. Effect of soil saturation and grain size on coupled hydrothermal flow in fine sands based on X-ray µ CT imaging. Geomechanics for Energy and the Environment, 2022, 32, pp.100380. 10.1016/j.gete.2022.100380. hal-03882556

HAL Id: hal-03882556

https://hal.science/hal-03882556

Submitted on 2 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effect of soil saturation and grain size on coupled hydrothermal flow in fine sands based on X-ray µCT imaging

Kui Liu<sup>a,*,1</sup>, Fleur A. Loveridge<sup>b</sup>, Richard Boardman<sup>c</sup>, William Powrie<sup>d</sup>

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

**Article history:**

Received 19 November 2021

Received in revised form 29 May 2022

Accepted 7 July 2022

Available online 12 July 2022

**Editors-in-Chief:**

Professor Lyesse Laloui and Professor Tomasz

Hueckel

Dataset link: https://doi.org/10.5518/1200

**Keywords:**

Hydrothermal flow

Unsaturated soils

X-ray µCT techniques

Soil saturation

Grain size

## ABSTRACT

Coupled hydrothermal flow can occur in soils, for example in applications such as ground heat storage and nuclear waste disposal. Therefore, approaches to quantitative analysis of water transfer in response to imposed thermal gradients are required, especially in unsaturated conditions. Analysis methods also require validation by laboratory and field data, which can be hard to obtain. This paper explores the possibility of using X-ray µCT techniques to observe and quantify water content changes in soils under thermal gradients. Specimens of a fine sand and a silty fine sand were prepared at degrees of saturation between 20% and 50%, before being subjected to heating from their base. Repeated scans, set up to balance image quality and scan duration, were carried out during the heating process, and Gaussian decomposition techniques were used to determine the changing soil phase proportions throughout the experiments. Based on these results and the accompanying numerical simulation of the experiments, it is shown that rapid vapour diffusion plays a more significant role than liquid flow in all cases. The rate of water content and hence degree of saturation change was more rapid in the less saturated specimens, especially for the fine sand. In practical terms, these moisture changes would result in reduction in thermal conductivity, especially in the soils of lower saturation. As well as providing insight into the dominant water transfer processes, the experiments show the feasibility of applying X-ray µCT techniques to thermal problems in soil mechanics.

© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Non-isothermal conditions are becoming increasingly important in geotechnical engineering. Buried cables, nuclear waste, or underground thermal energy storage facilities can release heat into the surrounding soils. In saturated soils, conduction through solids or forced advection by groundwater in more permeable strata are usually the key heat transfer mechanisms. These processes are well understood and the application of convection-diffusion equation can give a close match to field data (e.g. Refs. 1-3).

Early work suggested that heat and mass transfer were only likely to be coupled when the water was moving as a liquid.⁴ However, later research contradicts this. Radiation has been shown to be relevant in very coarse materials.⁵ In fine soils (fine sand or smaller grain sizes) with a degree of saturation <60%,

thermal movement of liquid moisture and diffusion of vapour are now believed to be of key importance for heat transfer.⁶ This is consistent with the observation of Hutcheon⁷ that a pressure gradient can be induced by net condensation and evaporation within the liquid phase, thus creating a liquid flow in opposition to the direction of vapour diffusion. These coupled processes are less studied in geotechnical engineering, although more commonly seen to be of importance in soil science.⁸

The classical theory for coupled diffusion of vapour and flow of liquid water in soils was set out by Philip & De Vries,⁹ based on adaptions of Darcy's Law and Fick's Law. However, in some cases the theory does not give a good fit to experimental data (e.g. Ref. 10), resulting in the introduction of empirical adjustment factors. Limitations of the theory for applications in geotechnical engineering include the absence of (i) convection of the gas and liquid phases, (ii) evaporation and condensation, owing to the assumption of water-vapour equilibrium, (iii) hysteresis, and (iv) a deformable solid matrix.¹¹,¹² The inclusion of non-equilibrium mass transfer is especially critical.¹³ The importance of including this in analysis has been shown numerous times in different fields, including in materials engineering,¹⁴ in soil science, e.g. Ref. 12, and for ground heat storage applications.¹⁵,¹⁶

<footnote><sup>a</sup> Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSE-ENV/SEDRE/LETIS, 92260, Fontenay-aux-Roses, France</footnote>

<footnote><sup>b</sup> School of Civil Engineering, Faculty of Engineering and Physical Sciences, University of Leeds, Leeds, UK</footnote>

<footnote><sup>c</sup> µ-VIS X-ray Imaging Centre, Faculty of Engineering and Physical Sciences, University of Southampton, Southampton, UK</footnote>

<footnote><sup>d</sup> School of Engineering, Faculty of Engineering and Physical Sciences, University of Southampton, Southampton, UK</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: kl5g14@soton.ac.uk (K. Liu).</footnote>

<footnote n="1">¹ Formerly Faculty of Engineering and Physical Sciences, University of Southampton, Southampton, UK.</footnote>
Nonetheless, the complexity of fully coupling heat and mass transfer including liquid water and vapour flow means that relatively little work has been carried out on the subject within the geotechnical sphere.⁸ There are also relatively few recent datasets available for testing theory, with notable recent exceptions including the laboratory scale physical models of Moradi et al.¹⁵ and Smits et al.¹². Experimental data at the micro scale can be particularly challenging to obtain.

In recent years, developments in X-ray CT scanning, a non-destructive method that allows imaging of the internal features of objects,¹⁷ have allowed higher resolution images to be achieved in a shorter timeframe. This means that the experimental technique has expanded its application in geotechnical engineering to include quantification of transient processes in finer grained soils. For example, iterative-algorithm reconstruction approaches allow larger datasets to be reconstructed from fewer frames,¹⁸ and synchrotron X-ray CT enables scans to be completed in seconds.¹⁹ There is also an increasing use of CT scanners in real-time experiments to track movement of soil particles (e.g. Ref. 20), changes in moisture conditions,²¹⁻²³ and soil response to mechanical loading.²⁴,²⁵

Taking advantage of these developments, this paper aims to (i) explore the potential for using X-ray CT techniques to observe and quantify the processes of liquid water movement in soils under thermal gradients; and (ii) make observations about the relative importance of liquid and vapour flows for soils of different grain sizes and degrees of saturation. These observations can then be used to draw conclusions about the processes that should be considered in applications in geotechnical engineering in different geological conditions.

## 2. Experimental method

### 2.1. Overall approach

Soil heating experiments were carried out via X-ray micro-CT (µCT) scans (e.g. Refs. 26, 27) in which small specimens were heated at their lower boundary and the consequent changes in water content and degree of saturation interpreted from the resulting scans. The experiments were carried out on two types of soils, a fine sand and a silty fine sand, at four different initial degrees of saturation. The following sections describe the materials used (Section 2.2), the experimental CT setups (Section 2.3) and the CT image reconstruction and segmentation processes (Section 2.4).

It proved difficult to maintain a fully closed system within the soil during the experiments, and additional physical experiments (Section 2.5) demonstrated that it was not possible to achieve conservation of mass within the scanned region. Consequently, the experimental results were also used to validate a numerical simulation. This numerical simulation (described in Section 3) allowed fuller insights to be derived from the experimental results.

### 2.2. Soil materials and specimen preparation

The fine sand used was fraction E Leighton Buzzard sand²⁸ and the silty fine sand used was a mixture of 80% fraction E Leighton Buzzard sand and 20% HIQ5 silica flour. The particle size distributions of both materials are given in Fig. 1. For each soil type, four different degrees of saturation were used in the experiments, with initial saturation of 20%, 30%, 40% and 50%. These types of soils and saturation ranges were selected to facilitate investigation into the conditions in which heat transfer mechanisms other than diffusion dominate. The mechanism proposed⁶ suggest that vapour diffusion could be important in fine sands and silty fine

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Particle size distribution (PSD) for the two types of soil used in the experiments.</figcaption></figure>

sands, but this domain has not been significantly investigated in previous research.

Cylindrical soil specimens with the dimension of 5 mm in diameter and 10 mm in height were prepared in special purpose specimen containers with acrylic side walls and an aluminum base to support one dimensional heat transfer. The small specimen sizes reflected the need for high spatial and temporal resolutions of the soil scanning due to the expected rapid heating.

Specimens were initially prepared in saturated conditions using wet pluviation.²⁹ The dry fine sand and silty fine sand were drained through a long and narrow funnel into the specimen container filled with de-aired water. The specimens were then dried in controlled conditions (20 °C and relative humidity of 50%) to achieve the pre-calculated gravimetric weight for every target degree of saturation. Repeatability of specimen preparation was checked in terms of soil porosity. For the fine sand this was 40 ± 1% and for the silty fine sand it was 33 ± 2%. A double-layered Parafilm³⁰ seal was then applied to the top of the completed specimens.

### 2.3. Experimental setup

Experimental work was carried out by X-ray µCT imaging using the Hutch 225 kVp³¹ with a horizontally-mounted molybdenum reflection target. This ameliorates any potential influence from the heel effect artefact, which can be pronounced with vertically mounted reflection targets.³² The soil specimen was mounted on the top of a Peltier heating unit from Deben, normal to the manipulator of the scanners. Temperature (T) and relative humidity (H), monitored using a hygrothermograph, were consistently at 20 ± 1 °C and 50 ± 1% respectively during scanning. No contrasting agent was needed thanks to the use of a molybdenum X-ray target.

Scan settings were chosen to achieve a high contrast between the different phases in the soil specimen, balancing the competing demands of a small-grain-size resolution, and to minimise scanning time. The major changes in water content in a soil specimen of this size occur within one hour. For each specimen, therefore, the experiment ran for nine scans, with each scan taking 7 min using the scan settings shown in Table 1. This was the best achievable compromise between (high) image quality and (short) scan duration.

At the start of the experiment the temperature of the Peltier unit was set to 20 °C, consistent with the ambient temperature
<table><caption>**Table 1**
Scan settings used in the heating experiments.</caption><thead><tr><td>Scanner</td><td>Energy</td><td>Power</td><td>Exposure time</td><td>Frame (per projection)</td><td>Projection count</td><td>Resolution (voxel size)</td></tr></thead><tbody><tr><td>Hutch</td><td>70 kV</td><td>6.9 W</td><td>354 ms</td><td>1</td><td>1201</td><td>~6 µm</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Reconstructed image slices for fine sand (left) and silty fine sand (right).</figcaption></figure>

so that no thermal gradient was imposed prior to scanning. The first scan was used to capture these initial conditions. The temperature of the Peltier unit was maintained at 50 °C from the start of the second scan until completion of the experiment after the ninth and final scan. Temperature control was via a remote workstation, and the experiment and nine scans ran continuously.

## 2.4. Image reconstruction and segmentation

Raw 2D projection data from image acquisition were reconstructed into 3D volume data using the cone-beam³³ filtered back-projection algorithm¹⁷ implemented in Nikon's CT Pro v2.2 software.³⁴ This was chosen over algebraic reconstruction techniques³⁵, following Liu et al.³⁶ because use of Gaussian decomposition for segmentation of the data (see below), meant there was no need to suppress the noise from the undersampled data. From visual observation of the reconstructed image slices for the two types of soil specimens (see Fig. 2), soil solid, liquid water and gas phases could be recognised. Examination of the point spread function³⁷ for the image data confirmed the image quality, returning values between 5 and 10.

Data pre-processing was implemented before quantitative segmentation into soil phases. A portion from the top and bottom regions of the reconstructed data was cropped to reduce the potential for errors caused by image artefacts.³⁸ 1500 image slices extracted from the entire 2000-slice image were retained for analysis, with the 250 slices nearest the top and base discarded. Additionally, a representative circular region of interest (ROI) in a horizontal plane was chosen from the middle of the specimen, to avoid the influence of the specimen boundary.

Gaussian decomposition was used for segmentation and post-processing, which was implemented on the raw 32-bit data.³⁸,³⁹ This approach has benefits over established thresholding methods when the scan resolution dimension is not sufficiently small compared with the soil grain sizes, then the partial volume effect could be pronounced.

The density of a material is correlated to the greyscale in X-ray CT images, with the histogram of the grey value (GV) data

from any scan usually following a Gaussian or normal distribution for each material present. In a single Gaussian formula there are three unknown variables: mean grey value ($\mu$), standard deviation ($\sigma$) and voxel count ($a$). The general form of a single Gaussian function can therefore be written:

$$f(x) = \frac{a}{\sqrt{2\pi}\sigma} e^{\frac{-(x-\mu)^2}{2\sigma^2}} \quad (1)$$

For a three-phase soil material, the overall GV histogram is equal to the sum of the individual Gaussian functions for the individual materials present. Hence Gaussian decomposition can be used to back calculate the proportions of the three phases.

Gaussian decomposition was first applied to the global GV histogram from additional scans of a dry and a fully saturated fine sand specimen carried out under the conditions in Table 1. This allowed calculation of the $\mu$ coefficients for each material ($\mu_g$ for gas, $\mu_w$ for liquid water and $\mu_s$ for solid), which are constant for each specimen. The voxel count of soil solid ($a_s$) was then determined from the initial scan data in each heating experiment by validating against the weight of the solid particles. $a_s$, $\mu_g$, $\mu_w$ and $\mu_s$ were then used as fixed inputs to the Gaussian decomposition of the subsequent scans in each experiment. The standard deviation was assumed to be the same for all phases, as all scans were conducted under the same scan settings.

The input coefficients and variables are presented in Table 2 for the fine sand and in Table 3 for the silty fine sand. Different mean GVs were calculated for the specimen corresponding to each degree of saturation, according to the calibration process for result optimisation. This also considered the potential variations in each specimen, for example the existence of dense sand materials, and scan conditions (e.g. subtle variations in GV from the X-ray source instability) along with the variance of water vapour in soils with different degrees of saturation which can be classified as either gas or liquid water according to their densities. The remaining coefficients in Tables 2 and 3 can then be determined by fitting with the GV histograms.

An example of Gaussian decomposition is presented in Fig. 3. The three phases have individual curves, which when summed
<table><caption>**Table 2**
Gaussian decomposition input coefficients for fine sand.</caption><thead><tr><th>Initial saturation</th><th>a<sub>g</sub></th><th>μ<sub>g</sub></th><th>a<sub>w</sub></th><th>μ<sub>w</sub></th><th>a<sub>s</sub></th><th>μ<sub>s</sub></th><th>σ</th></tr></thead><tbody><tr><td>20%</td><td>-</td><td>82</td><td>-</td><td>180</td><td>3.10 × 10<sup>7</sup></td><td>335</td><td>-</td></tr><tr><td>30%</td><td>-</td><td>68</td><td>-</td><td>150</td><td>3.13 × 10<sup>7</sup></td><td>350</td><td>-</td></tr><tr><td>40%</td><td>-</td><td>63</td><td>-</td><td>150</td><td>3.78 × 10<sup>7</sup></td><td>315</td><td>-</td></tr><tr><td>50%</td><td>-</td><td>68</td><td>-</td><td>120</td><td>3.62 × 10<sup>7</sup></td><td>315</td><td>-</td></tr></tbody></table>

<table><caption>**Table 3**
Gaussian decomposition input coefficients for silty fine sand.</caption><thead><tr><th>Initial saturation</th><th>a<sub>g</sub></th><th>μ<sub>g</sub></th><th>a<sub>w</sub></th><th>μ<sub>w</sub></th><th>a<sub>s</sub></th><th>μ<sub>s</sub></th><th>σ</th></tr></thead><tbody><tr><td>20%</td><td>-</td><td>63</td><td>-</td><td>177</td><td>2.82 × 10<sup>7</sup></td><td>334</td><td>-</td></tr><tr><td>30%</td><td>-</td><td>87</td><td>-</td><td>175</td><td>2.33 × 10<sup>7</sup></td><td>296</td><td>-</td></tr><tr><td>40%</td><td>-</td><td>80</td><td>-</td><td>176</td><td>3.13 × 10<sup>7</sup></td><td>330</td><td>-</td></tr><tr><td>50%</td><td>-</td><td>70</td><td>-</td><td>182</td><td>2.49 × 10<sup>7</sup></td><td>290</td><td>-</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Example of Gaussian decomposition for silty fine sand at 30% saturation. The thick black line is the raw CT data; the grey line is the fitted superposition of the three individual Gaussian distributions.</figcaption></figure>

give the total GV curve. The gas phase (with the lowest density) is at the left of the GV axis and the solid phase (with the highest density) at the right, with liquid water in middle. The area underneath each curve represents the voxel count corresponding to that phase portion.

## 2.5. Parafilm tests

Additional tests were conducted to investigate whether the Parafilm lid, which was applied to the top of each specimen for sealing purposes, really acted to prevent mass transfer from the top of the specimen and hence ensured a closed soil system. The experiments were prompted by observations of liquid water droplets on the underside of the Parafilm at the end of the heating experiments, which suggested that although the specimen-Parafilm system may have been sealed, water was able to leave the soil specimen itself. This would mean that there was no longer mass conservation within a soil specimen.

Additional soil specimens were prepared using the procedure outlined in Section 2.2 above. These specimens were then weighed, sealed with Parafilm, and subjected to heating in the same way as described in Section 2.3, using a 50 °C basal temperature. However, these tests were not carried out inside the scanner. In this case, the specimens were heated for one hour as described in Section 2.3, the Parafilm lid was removed and the specimens immediately re-weighed to determine the mass of water that was now outside the soil specimen, having re-condensed on the Parafilm lid.

## 3. Numerical model

A numerical simulation was developed using COMSOL Multiphysics,⁴⁰ based on the non-equilibrium approach.¹²,¹⁵,⁴¹ The approach includes liquid/gas phase change and has been validated against physical experiments of evaporation,¹² and heat transfer around ground heat exchangers.¹⁵

Fig. 4 illustrates the key physical processes included in the model, which is based on the principles of mass conservation for the flow of liquid and gas (air and water vapour) and energy conservation for heat transfer. The governing equations are set out below.

### 3.1. Governing equations

The governing equations were obtained by substituting the relevant flow law(s) into each conservation equation. Darcy's law was used for liquid water and gas flow, Fick's law was used for vapour flow, and Fourier's law was used for heat flow.

For liquid water and gas, the governing equations in terms of mass balances¹³,⁴² are respectively:

$$ \rho_w \frac{\partial \theta_w}{\partial p_c} \frac{\partial p_c}{\partial t} + \theta_w \frac{\partial \rho_w}{\partial t} + \nabla \left( -\frac{\rho_w k_{rw} k_s}{\mu_w} (\nabla p_w + \rho_w g) \right) = -R_{gw} \quad (2) $$

$$ \rho_g \frac{\partial \theta_g}{\partial p_c} \frac{\partial p_c}{\partial t} + \theta_g \frac{\partial \rho_g}{\partial t} + \nabla \left( -\frac{\rho_g k_{rg} k_s}{\mu_g} (\nabla p_g + \rho_g g) \right) = R_{gw} \quad (3) $$

where $P_c$ is the capillary pressure ($P_c = P_g - P_w$) (kPa); $P_w$ is the liquid water pressure and $P_g$ is the gas pressure; $\rho_w$ and $\rho_g$ are the densities of liquid water and gas, respectively (kg/m³); $\mu_w$ and $\mu_g$ are the dynamic viscosities of liquid water and gas, respectively (Pa·s); $\theta_w$ and $\theta_g$ are the volumetric water content and gas content, respectively; $R_{gw}$ is the phase change rate between water and vapour due to evaporation or condensation (kg/m³ s); $k_s$ is the intrinsic permeability of the granular matrix (m²); $k_{rw}$ and $k_{rg}$ are the relative permeabilities to liquid water and gas, respectively; g is the gravitational acceleration (m/s²); t is the time (s).

Fixed parameters (dynamic viscosity, intrinsic permeability) are given in Table 4.

The relationship between volumetric water content ($\theta_w$) and capillary pressure ($P_c$) is as proposed by van Genuchten⁴³ and Mualem⁴⁴:

$$ \theta_w = \theta_r + (\theta_s - \theta_r) [1 + (\alpha P_c)^N]^m \quad (4) $$

where $\alpha$ (kPa⁻¹) is the air-entry value; N and m are the model coefficients related to the pore size distribution and the geometry of the solid matrix ($m = 1 - 1/N$, Ref. 43).

Values of the relevant parameters are given in Table 5.

The relative permeabilities of water and gas can be obtained from the effective volumetric water content using Eqs. (5) and (6) for the liquid phase gas phase, respectively.⁴³

$$ k_{rw} = \theta_e^\gamma [1 - (1 - \theta_e^{1/m})^m]^2 \quad (5) $$

$$ k_{rg} = (1 - \theta_e)^\gamma (1 - \theta_e^{1/m})^{2m} \quad (6) $$

where an effective water content term ($\theta_e$) is expressed as $\theta_e = \frac{\theta_w - \theta_r}{\theta_s - \theta_r}$; $\theta_s$ (equivalent to the porosity, n = 0.4) and $\theta_r$ (0.01, Ref. 44) are the saturated and residual volumetric water content, respectively.; $\gamma$ is a model fitting parameter related to pore connectivity or tortuosity.

The governing equation for water vapour can be expressed as:

$$ \frac{\partial \rho_g \theta_g w_v}{\partial t} + \nabla \left( \rho_g \theta_g w_v u_g - D_v \rho_g \nabla w_v \right) = R_{gw} \quad (7) $$
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Schematic diagram summarising the physical mechanisms underpinning the numerical model.</figcaption></figure>

<table><caption>**Table 4**

Material Parameters used in numerical simulations.</caption><thead><tr><th>Material parameter</th><th>Density, ρ (kg/m³)</th><th>Intrinsic permeability, k<sub>s</sub> (m²)</th><th>Dynamic viscosity, μ (Pa s)</th><th>Specific heat capacity, c<sub>p</sub> (J/(kg K))</th><th>Thermal conductivity, k<sub>T</sub> (W/(m K))</th></tr></thead><tbody><tr><td>Soil solid</td><td>2650</td><td>–</td><td>–</td><td>800</td><td>7.7</td></tr><tr><td>Liquid water</td><td>See Table 6</td><td>2 × 10<sup>-12</sup></td><td>See Table 6</td><td>4181</td><td>0.6</td></tr><tr><td>Water vapour</td><td>ρ<sub>ve</sub> (= ρ<sub>v,s</sub> H)<br>See Table 6</td><td>2 × 10<sup>-12</sup></td><td>–</td><td>1864</td><td>–</td></tr><tr><td>Gas</td><td>Table 6</td><td>2 × 10<sup>-12</sup></td><td>18.6 × 10<sup>-6</sup></td><td>From Eq. (11)</td><td>0.024</td></tr></tbody></table>

<table><caption>**Table 5**

Empirical coefficients used in numerical simulations.</caption><thead><tr><th>Empirical coefficient</th><th>b</th><th>B (s/m²)</th><th>N</th><th>m</th><th>α (1/kPa)</th><th>γ</th></tr></thead><tbody><tr><td>Liquid flow</td><td>–</td><td>000021</td><td>6.77</td><td>m = 1–1/N</td><td>0.15</td><td>0.5</td></tr><tr><td>Vapour flow</td><td>9</td><td>000021</td><td>6.77</td><td>m = 1–1/N</td><td>0.15</td><td>0.5</td></tr><tr><td>Gas flow</td><td>–</td><td>000021</td><td>6.77</td><td>m = 1–1/N</td><td>0.15</td><td>0.5</td></tr><tr><td>Heat flow</td><td>–</td><td>000021</td><td>6.77</td><td>m = 1–1/N</td><td>0.15</td><td>0.5</td></tr></tbody></table>

where $w_v$ is the vapour mass fraction in gas phase; $u_g$ is the mean pore velocity of gas (m/s); $D_v$ is the vapour diffusion coefficient in a soil (m²/s).

The vapour diffusion coefficient in a soil can be expressed as $D_v = D_{va}n^2(S_g/n)^{2+3/b}$,⁴⁵,⁴⁶ where $D_{va}$ is the vapour diffusion coefficient in air (Table 6) and $S_g$ is the saturation of gas phase. The term $2 + 3/b$ is an analogous to the Burdine⁴⁷ and Campbell⁴⁸ tortuosity model for describing unsaturated hydraulic conductivity,⁴⁹ where b is the Campbell pore-size distribution parameter (Table 5).

The governing equation for heat flow is expressed as:

$$C_v \frac{\partial T}{\partial t} + \nabla \left( (\rho_w c_{pw}) u_w \theta_w T + (\rho_g c_{pg}) u_g \theta_g T - (k_T \nabla T) \right) = -L R_{gw} \quad (8)$$

where $C_v$ is the volumetric heat capacity (J/(m³ K)); $c_{pw}$ and $c_{pg}$ are the specific heat capacities of water and gas (Table 4), respectively, (J/(kg K)); $k_T$ is the effective soil thermal conductivity (W/(m K); $u_w$ is the mean pore velocity of water (m/s); L is the latent heat coefficient given in Table 6 (J/kg); T is the temperature (K).

The effective (bulk) soil thermal conductivity ($k_T$, W/(m K)) can be written in terms of the constituent components as Ref. 50:

$$k_T = (1 - n) k_{Ts} + n Sk_{Tw} + n (1 - S) k_{Tg} \quad (9)$$

where the specific thermal conductivities of soil materials are assumed to be constant (Table 4) and independent of temperature effect. $k_{Ts}$ (W/(m K)) is based on Refs. 51–53; and $k_{Tw}$ (W/(m K)) and $k_{Tg}$ (W/(m K)) is based on Ref. 54.

The volumetric heat capacity ($C_v$), is given as Refs. 16, 55:

$$C_v = \rho c_p = (1 - n) \rho_s c_{ps} + nS (\rho_w c_{pw}) + n (1 - S) (\rho_g c_{pg}) \quad (10)$$

where the specific heat capacity values of materials, $c_p$ (J/(kg K)) are assumed to unaffected by temperature. $c_{ps}$ is the specific heat capacity of the soil solid and $c_{pw}$ the specific heat capacity of the water. The specific heat capacity of the gas, $c_{pg}$, is given by Eq. (11):

$$c_{pg} = c_{pv} w_v + c_{pa} (1 - w_v) \quad (11)$$

where $c_{pa}$ (1012 J/(kg K)) is the specific heat capacity of dry air and $c_{pv}$ (J/(kg K)) is the specific heat capacity of the vapour (Table 4).

The governing equations were implemented in COMSOL 5.4 using the standard modules for Darcy's law (for flow of liquid water and gas), the stabilised convection–diffusion equation (for the movement of water vapour in the gas phase), and the general coefficient based partial differential equation module (for heat flow).

## 3.2. Geometry, initial and boundary conditions

A radial plane of the 10 mm high cylinder soil specimen was considered in the axisymmetric model, with the analysis conducted in one dimension in the Z-direction (Fig. 5). The initial temperature of the soil specimen was set to the ambient temperature in the Hutch scanner, 20 °C (293.15 K). To represent the start of heating, the temperature at the base was set to increase linearly to 35 °C (308.15 K) over the time of the first scan, after which it was held constant. Testing of the Peltier unit has shown this to be a good representation of the physical experiments. The initial mass fraction of vapour in the gas phase ($w_v$) was set
<table><caption>**Table 6**
Expressions for parameters used in numerical simulations.</caption><thead><tr><th>Parameter</th><th>Equation</th><th>Reference</th></tr></thead><tbody><tr><td>Relative humidity at equilibrium (H)</td><td>H = exp [P<sub>c</sub>M<sub>w</sub> / (ρ<sub>w</sub>RT)]</td><td>Ref. 9</td></tr><tr><td>Density of liquid water (ρ<sub>w</sub>)</td><td>ρ<sub>w</sub> = 1 - 7.37 × 10<sup>-6</sup>(T - 273.15 - 4)<sup>2</sup> + 3.79 × 10<sup>-8</sup>(T - 273.15 - 4)<sup>3</sup></td><td>Ref. 56</td></tr><tr><td>Density of gas (ρ<sub>g</sub>)</td><td>ρ<sub>g</sub> = P<sub>g</sub>/RT M<sub>g</sub> = P<sub>g</sub>/RT [w<sub>v</sub> (M<sub>a</sub> - M<sub>w</sub> / M<sub>a</sub>M<sub>w</sub>) + 1 / M<sub>a</sub>]<sup>-1</sup></td><td>Ideal gas law</td></tr><tr><td>Viscosity of liquid water (μ<sub>w</sub>)</td><td>μ<sub>w</sub> = 0.6612(T - 229)<sup>-1.562</sup></td><td>Ref. 57</td></tr><tr><td>Latent heat of water vapourisation (L)</td><td>L = 2.501 × 10<sup>6</sup> - 2369.2(T - 273.15)</td><td>Ref. 58</td></tr><tr><td>Vapour diffusion coefficient in air (D<sub>va</sub>)</td><td>D<sub>va</sub> = 2.12 × 10<sup>-5</sup>(T/273.15)<sup>2</sup></td><td>Ref. 59</td></tr><tr><td>Phase change rate (R<sub>gw</sub>)</td><td>R<sub>gw</sub> = (B(θ<sub>w</sub> - θ<sub>r</sub>)RT / M<sub>w</sub>) (ρ<sub>v,s</sub>RH - ρ<sub>g</sub>w<sub>v</sub>)</td><td>Refs. 14, 60</td></tr><tr><td>Density of saturated water vapour (ρ<sub>v,s</sub>)</td><td>ρ<sub>v,s</sub> = exp(31.37 - 6014.79/T - 7.92 × 10<sup>-3</sup>T) / T</td><td>Ref. 59</td></tr></tbody></table>

Note: R is the universal gas constant (8.31 J/mol K); M<sub>w</sub> (0.018 kg/mol), M<sub>a</sub> (0.028 kg/mol) and M<sub>g</sub> (kg/mol) are the molecular weights of water, dry air, and gas phases, respectively; B (s/m<sup>2</sup>) is an empirical fitting parameter given in Table 5.

### 3.3. Key assumptions

The model contains some important general assumptions. Adjacent phases were assumed to have the same local temperature (T<sub>gas</sub> = T<sub>water</sub> = T<sub>soil</sub> = T), although it will vary spatially and temporally in the model. The soil structure was also assumed to be undeformable, so that the porosity of the specimen remained constant with no volume change occurring. Most soil parameters were assumed to be constant and independent of temperature, including the residual water content, viscosity of gas, density and volume of the soil matrix along, and the thermal conductivity for each phase. The P<sub>c</sub> - θ relationship described by the van Genuchten model was also assumed to be temperature invariant.

Temperature dependent parameters are given in Table 6. The temperature dependency of the density and viscosity of liquid water means that buoyancy-driven natural-convection flow was included in the simulation. While the interfacial drag between the fluid phases was neglected, the much more significant drag between the fluid phases and the solid matrix was included in Darcy's law.

### 3.4. Mesh sensitivity

The mesh was formed of triangular elements. A mesh sensitivity study was carried out for the fine sand with an initial degree of saturation of 30%, using total element numbers between approximately 100 and 6500. For all cases within that range, the final degree of saturation at the end of heating was 21.7% compared with 21.1% in the physical experiment. However, cases with fewer than approximately 1000 elements gave inconsistent internal patterns of specimen moisture content, suggesting the element size was too large compared with the scale of the processes taking place within the specimen. After further investigation, a model with 1784 elements (Fig. 5) was chosen to provide more consistent patterns of liquid water distribution, with no further change in distribution occurring with the use of smaller elements.³⁸

### 3.5. Model calibration and validation

The outflow boundary at the top of the numerical model required calibration. This was because the precise amount of water leaving the soil specimen at the top of the sample, and any change in rate of flux with time, could not easily be quantified from

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Schematic of the model domain.</figcaption></figure>

to 0.03, following a sensitivity analysis. P<sub>g</sub> was assumed to be atmospheric pressure (1 atm) and P<sub>w</sub> was determined from the van Genuchten model.⁴³,⁴⁴

The model boundary conditions are summarised in Table 7. The thermal boundary conditions at the outer side and the top of the soil specimen were set to convective. For the acrylic wall at the outer side a heat transfer coefficient of 20 W/m² was used. For the top surface a value of 10 W/m² was applied, corresponding to natural convection.⁶¹ Sensitivity analyses showed that these values do not have a significant impact of the key simulation outputs. Mass transfer was not permitted from the base or the outer edge of the specimen. However, based on observations of liquid water on the base of the Parafilm at the end of the physical experiments, an outflow boundary was assumed for liquid, gas and vapour at the upper surface. Calibration of this boundary is described in Section 3.5.
<table><caption>**Table 7**
Flow boundary conditions.</caption><thead><tr><th>Boundary</th><th>Liquid flow</th><th>Gas flow</th><th>Vapour flow</th><th>Heat flow</th></tr></thead><tbody><tr><td>Top surface</td><td>Outflow</td><td>Outflow</td><td>Outflow</td><td>Convective heat transfer</td></tr><tr><td>Bottom base</td><td>No flux</td><td>No flux</td><td>No flux</td><td>Constant uniform temperature</td></tr><tr><td>Outer side</td><td>No flux</td><td>No flux</td><td>No flux</td><td>Convective heat transfer</td></tr><tr><td>Inner side</td><td>Symmetry</td><td>Symmetry</td><td>Symmetry</td><td>Symmetry</td></tr></tbody></table>

<table><caption>**Table 8**
Outflow boundary conditions used in simulations.</caption><thead><tr><th>Outflow boundary</th><th>Liquid inward mass flux [kg/(m² s)]</th><th>Gas inward mass flux [kg/(m² s)]</th><th>Vapour boundary flux/source [kg/(m² s)]</th><th>Heat flow [W/m²]</th></tr></thead><tbody><tr><td>Case 1</td><td>−ρ<sub>w</sub> × U<sub>(t)</sub></td><td>−ρ<sub>g</sub> × U<sub>(t)</sub></td><td>−ρ<sub>ve</sub> × U<sub>(t)</sub></td><td>10 × 293.15</td></tr><tr><td>Case 2</td><td>0</td><td>−ρ<sub>g</sub> × U<sub>(t)</sub> −ρ<sub>w</sub> × U<sub>(t)</sub></td><td>−ρ<sub>ve</sub> × U<sub>(t)</sub><br>−ρ<sub>w</sub> × U<sub>(t)</sub></td><td>10 × 293.15</td></tr><tr><td>Case 3</td><td>−ρ<sub>w</sub> × U<sub>(t)</sub><br>−ρ<sub>ve</sub> × U<sub>(t)</sub></td><td>0</td><td>0</td><td>10 × 293.15</td></tr><tr><td>Case 4</td><td>−ρ<sub>w</sub> × 0.5 × U<sub>(t)</sub></td><td>−(ρ<sub>g</sub>+0.5 × ρ<sub>w</sub>) × U<sub>(t)</sub></td><td>−(ρ<sub>ve</sub>+0.5 × ρ<sub>w</sub>) × U<sub>(t)</sub></td><td>10 × 293.15</td></tr><tr><td>Case 5</td><td>−ρ<sub>w</sub> × U<sub>(t)</sub><br>−0.5 × ρ<sub>ve</sub> × U<sub>(t)</sub></td><td>−ρ<sub>g</sub> × U<sub>(t)</sub>+<br>0.5 × ρ<sub>ve</sub> × U<sub>(t)</sub></td><td>−0.5 × ρ<sub>ve</sub> × U<sub>(t)</sub></td><td>10 × 293.15</td></tr></tbody></table>

the X-ray CT experiments. However, the Parafilm tests described in Section 2.5 indicated that a mass of water in the order of $10^{-4}$ to $10^{-3}$ grams left each soil specimen over the one-hour heating period. This corresponds to an average rate of evaporation (m³/s per m² specimen area, i.e. m/s) of $10^{-8}-10^{-7}$ m/s at the soil surface, consistent with the order of magnitude of observed steady state soil surface evaporation rates.⁶²

To investigate the nature of the outflow conditions, a constant total flux of $U_{(t)}$, within the range $10^{-8}-10^{-7}$ m/s, was specified as a starting point. Different relative outflows of liquid water, gas and vapour were then trialled based on the five cases given in Table 8. Case 1 considered the same rate (in m/s) for all three phases, Case 2 had no liquid flux, Case 3 had only liquid flux and Cases 4 and 5 had lesser and greater proportions of liquid flow respectively. For each case, the value of $U_{(t)}$ at each scan time in the physical experiment was adjusted manually to achieve the best fit between the simulation results and the global saturation values in the experiments for the four specimens of fine sand. The best match was obtained for Case 2 with no liquid water outflow. Fig. 6 shows the calibrated outflows in terms of $U_{(t)}$ and Fig. 7 shows the match between the simulation and the experimental data. The values of $U_{(t)}$ vary up to a maximum of $4 \times 10^{-7}$ m/s and are therefore consistent with both the results of the Parafilm tests and previous work of soil surface evaporation. Although not shown in the figures, the distribution and changes in distribution of local moisture content were also checked for comparability between the simulation and the experimental data to confirm the plausibility of calibration process.³⁸

Fig. 6 shows that $U_{(t)}$ tended to increase and reaches its maximum value at an early stage when heating was first applied (time > 0 min) and then decreased over time. The specimens of the lower saturation showed a higher initial rate of outflow, which then reduced rapidly. By contrast, the specimens of the higher initial saturation also displayed an initial increasing outflow but became stable at a relatively lower rate that sustained for a longer period. In all cases, equilibrium with no further outflow was reached by the end of the experiments. While Case 2 provided the best fit, it should be noted that this same pattern was observed whilst fitting for all the cases (given in Table 8). Nonetheless, the adoption of Case 2 implies that the role of gas and vapour flow at the top of the sample was much more important than that of liquid water flow, and that these processes occurred much more rapidly upon heating in lower saturation specimens. This topic will be reviewed in the following section, where the experimental results are discussed in more detail.

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Calibrated distributions of $U_{(t)}$ over time for Case 2.</figcaption></figure>

# 4. Results

## 4.1. Experimental results

### 4.1.1. Characteristics of specimens

Both the fine sand and the silty fine sand show a non-uniform distribution of initial moisture over the height as shown in Fig. 8. Typically the GV, which indicates the density of a material, is lower at the top of the specimen. This represents either a higher porosity or a lower water saturation; the latter is more likely, given the drying of specimens used to obtain the target initial saturation values.

### 4.1.2. Variation in saturation with time

Fig. 9 shows how all the eight specimens dried overall with time. (Note that the scale is very different from that used in Fig. 7 and the origin in Fig. 9 is true; but that the experimental data for the fine sand data are the same). For the fine sand specimens, the change in the overall degree of saturation due to heating was approximately 10% regardless of the initial value, although the changes occurred more rapidly when the initial saturation was lower. For the silty fine sand, the overall change in saturation due to heating was less, at about 5%, and it was harder to discern a difference in behaviour depending on initial saturation values.

The local distribution of moisture in the specimens and its change over time in response to heating is investigated with
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison of simulated and experimental results for the change in degree of saturation over time for fine sand with outflow condition Case 2. Initial degree of saturation: (a) 20% (b) 30% (c) 40% (d) 50%.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Example of orthogonal images at initial conditions and GV distributions with depth over time for (a) fine sand of 30% initial saturation and (b) silty fine sand of 30% initial saturation.</figcaption></figure>

reference to isochrones for the degree of saturation with depth (Fig. 10). Each data point corresponding to a certain height at each isochrone (under each scan) was determined from the Gaussian decomposition analysis on the corresponding subsection image data (comprising 300 image slices per subsection), and the middle position of the subsection was taken as its vertical location within the sample height.

Consistent with the specimen overall saturation data (Fig. 9), the coarser specimens and those with a lower initial saturation dry more rapidly and to a greater degree. In all cases, drying is greatest at the base of the specimen closest to the heat source. Here progressive drying is seen with each time step. However, in the upper part of some specimens, particularly the fine sand at lower degrees of saturation, the degree of saturation may increase again later in the test. This is manifest in the crossing of the saturation vs. depth isochrones for different times. This may be a result of the direct movement of liquid water within the specimen, but as discussed in later sections is more likely to

be due to re-condensation of water vapour that had previously evaporated in the lower part of the specimen.

### 4.1.3. Loss of mass

The degree of saturation data from the X-ray CT scans suggests that the fine sand specimens lost approximately 10% saturation and the silty fine sand specimens lost approximately 5% saturation due to heating, regardless of the initial degree of saturation. This can be compared with the changes in saturation calculated based on the mass loss measured during the Parafilm sealing experiments described in Section 2.5. Fig. 11 shows that in these experiments, the change in degree of saturation depends on both the grain size and the initial saturation. For an initial degree of saturation of 50%, the Parafilm test specimens exhibit a change in saturation within a few percent of that in the X-ray CT experiments. However, for the specimens with an initial saturation degree of 20%, the Parafilm tests show a much more reduced degree of drying than the X-ray CT experiments.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Change in overall degree of saturation over time.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Isochrones showing the distribution of degree of saturation with depth within specimens at different times, for the fine sand (grey) and silty fine sand (black) with initial saturation values of (a) 20% (b) 30% (c) 40% (d) 50%.</figcaption></figure>

Both the X-ray CT results and the Parafilm tests data show water leaving the soil specimen but to differing degrees. To understand this difference, we need to remember that the X-ray CT data are only for the ROI in the centre of the specimen, whereas the Parafilm tests cover the entire specimen. The ROI from the X-ray CT data is a cylindrical volume that crops about 1 mm from the top and base of the specimen and approximately 0.5 mm for the sides. The part of the specimen outside the ROI forms a three-dimensional annulus. The ROI and the annulus represent 51% and

49% of the total specimen volume, respectively. Assuming that the specimen is uniform in terms of initial water distribution (since the spatial variability within the annulus is not known), we can use the experimental results to determine a water mass balance for the ROI and the annulus.

Fig. 12 shows the results of the water balance calculation. Initially, prior to heating, 51% of the water in the specimen was within the ROI and 49% within the annulus. After heating, up to 14% of the water mass had left the specimen completely. Based
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Overall change in the degree of saturation for the fine sand and silty fine sand determined from the change in mass during the parafilm tests.</figcaption></figure>

on the mass of liquid water measured on the Parafilm, it is known that up to 6% of this had vaporised and re-condensed. The remaining up to 8% must have escaped as vapour. The amount of water within the ROI decreased, to between 27% and 40% for the fine sand (Fig. 12a), and to between 39% and 45% for the silty fine sand (Fig. 12b). The amount of water in the annulus outside the ROI either increased or decreased depending on the initial degree of saturation. If the initial degree of saturation was high, the water in the annulus reduced. If the initial degree of saturation was low, the amount of water in the annulus increased; these increases corresponded to the cases where the reduction in moisture in the ROI was greatest. This suggests that redistribution of moisture occurred from the ROI to the annulus for specimens possessing a lower initial saturation, and especially for the fine sand (Fig. 12a).

Knowing the specimen porosities of 40% (for the fine sand) and 33% (for the silty fine sand), it is possible to calculate the available void space in the annulus and determine whether it is possible for it to have received this additional water. In all cases the answer is yes. However, the specimens are heated from the base, and given the evidence from Fig. 10 for the greater drying at the base as well as increases in the degree of saturation near the top of the ROI in some cases, this suggests that water is likely to have moved, in particular from the ROI to the upper part of the annulus. If all the transferred liquid water had moved to the upper half of the annulus, the final degree of saturation here would be between 40% and 90%, with the higher values pertaining to cases with a higher initial saturation.

These calculations assume that the specimen is uniform, which it is known not to be. Owing to the nature of specimen preparation, specimens tended to have a greater saturation near the base initially and also may have a higher porosity near the edge. The first may mean that the extent of moisture movement is underestimated, while the second may mean there is more void space at the edge of the specimen to accommodate the moisture movement; consequently the final saturation values for the upper part of the annulus given above could be overestimated. The other place to accommodate moisture would be in any ullage or “headspace” between the top of the specimen and the Parafilm cap. While this should theoretically be zero, very little space would be needed to accommodate the volumes of water involved ($<1 \times 10^{-8} \text{ m}^3$).

## 4.2. Numerical results

Examination of the numerical results allows additional inferences to be made about heat and mass transfer in the different soil phases which cannot be observed directly from the X-ray CT results.

### 4.2.1. Thermally driven vapour flow

Consistent with other studies in different scenarios, e.g. Ref. 63, the results of the numerical analysis conducted for the fine sand case shown that the dominant flow mechanism in the specimens is vapour flow. Fig. 13 illustrates this with arrow plots of liquid and vapour flow directions for the initial saturations of 20% and 50% at 7 min into the heating experiments, close to when the flow would be expected to be the most significant. The flow domain is presented within the full specimen height (z) versus half symmetrical geometry (radius, r), and the colour of the arrows indicates the magnitude of the flow occurring at this time. In Fig. 13 (and for all other analyses), water flow is seen to be small and in a downward direction due to gravity, while vapour flow is much more important and in an upward direction. This is because of the effect of the imposed temperature change, which induces a vapour density difference that drives flow. These results are consistent with the work of Hutcheon,⁷ which showed that the direction of liquid flow is opposite to that of vapour flow, and that the liquid pressure gradient induced by heat is less significant than that of gravity. This effect is greater closer to the heat source, so that the gravity driven liquid water flow is reduced in the lower parts of the specimens (Fig. 13). Liquid flow is also less significant in the less saturated specimen, reflecting the reduced water permeability associated with low saturation levels.

Similarly, the fact that moisture is rapidly moving away from the heat source by vapour diffusion is at least a partially reflection of the high relative vapour permeability. Fig. 14 shows the estimated variations in the relative gas and liquid permeabilities (as Krg and Krw, respectively) over time for two initial degrees of saturation (20% and 50%). The higher initial saturation specimen shows a higher initial liquid permeability with some reduction due to drying. These cases also show a corresponding large increase in vapour permeability along the experiment progress.

In both cases Fig. 14 shows vapour diffusion to be dominant. However, the highest relative gas permeability occurs in the lower saturation specimen. This is consistent with these specimens experiencing the quickest changes in saturation (Fig. 7).

### 4.2.2. Vapour density

Fig. 15 shows the changes in vapour density with time, along with the changes in degree of saturation from the numerical model for the fine sand specimens with the initial saturation of 20% and 50%, respectively. While the change in degree of saturation depends on the initial saturation, being faster for the least saturated specimen, the change in vapour density is rapid in both cases. All the change in vapour density happens within the first 7 min. This is due to the rapid change in specimen temperature, for example the specimen with the saturation of 20% as shown in Fig. 16. Beyond 7 min, there is negligible further change in temperature occurs.

The fact that the degree of saturation continues to change after 7 min, however, suggests that the vapour must be leaving the system (and in reality re-condensing on the Parafilm) and being replaced by a further water evaporation from the specimen, which results in the continued drying of the specimen. This effect is especially important for specimens with a higher initial saturation. It is also consistent with the Parafilm test results, which show that a greater proportion of the moisture leaves the soil and accumulates on the Parafilm in the case of specimens having a higher initial degree of saturation.
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Change in water distribution between the ROI (square symbols), the outer annulus (triangle symbols) and outside of the specimen (circle symbols) specimens at end of heating (open symbols) relative to Initial Conditions (solid symbols); (a) fine sand; (b) silty fine sand.</figcaption></figure>

### 4.2.3. Re-distribution of liquid water

As well as confirming the dominance of vapour diffusion as the means of heat and mass transfer within the specimens, the numerical simulation also allows a greater understanding of the potential for thermal redistribution of moisture due to re-condensation of vapour. Fig. 10 shows that in the experiments on the fine sand with the lowest initial degree of saturation there appeared to be some increase in moisture content in the upper part of the specimen later in the test. This tendency is confirmed by the numerical results as shown in Fig. 17. For the first 14 min the sample is seen to dry significantly. After this time, further drying is seen near the base but in the upper third of the specimen the change in saturation reverses slightly, before drying is observed again throughout the specimen in the later stage of the experiment. This confirms the suggestion that the water that is not escaped as vapour, re-condensed on the Parafilm lid, or retained within the ROI of the X-ray CT image, could be present as liquid water in the top of the specimen (refer also to Fig. 12).

## 5. Discussion

### 5.1. Experimental approach

One of the aims for this paper was to explore the usefulness of the X-ray µCT scanning techniques in energy geotechnics problems in which water may be redistributed in unsaturated fine sands or smaller-grain-size soils. While the results presented do show the ability to observe these processes in fine sands and silty fine sands, there are some limitations. First, the experimental results were dominated by vaporisation and vapour movement. While this cannot be observed directly by X-ray CT techniques, the later re-condensation of moisture can be detected. However, this was only seen in some scenarios during the experiments, owing in part to the challenge of maintaining a closed system. It also reflects the need to crop the X-ray CT data to avoid errors or the influence from image artefacts. These points could be addressed in future work, but require permanent sealing of the specimen, or testing under confined conditions which would add to the complexity of the experimental arrangements.

### 5.2. Implications for applications in energy geotechnics

The results indicate that there is potential for significant rapid vapourisation and vapour flow due to thermal gradients, particularly in soils with a low initial degree of saturation. Resultant drying will increase the soil shearing resistance. Also, there is potential for changes in soil thermal properties which are of importance in problems such as ground heat storage and heat dissipation from buried cables.

The results presented in Section 4 suggest changes in soil saturation of only 5% to 10% is expected to occur, even with heating up to 50 °C. Consideration of the changes in thermal conductivity with the degree of soil saturation measured by other researchers (Fig. 18) shows that the relative significance of the changes in thermal properties depends on the type of soil. Clayed soils will experience rapid changes in thermal conductivity due to small changes in moisture content at degrees of saturation over 40%. On the other hand, silts and sands can experience rapid changes in thermal conductivity at degrees of saturation less than 10%. This suggests that the type of soil materials tested in this study could see a reduction in thermal conductivity by 50% as a result of heating, which would be significant for problems in energy geotechnics.

## 6. Conclusions

This paper has presented the results of heating experiments on soil specimens conducted via X-ray µCT scans and numerical analyses, to investigate the processes of heat and mass transfer problems occurring in unsaturated fine soils. The results indicate that:

* Coarser soil (fine sand) dries more rapidly and to a greater extent than finer soil (silty fine sand). This is due to the relatively larger surface area and the smaller pore size of the finer soil which would take longer to dry at the same saturation under the same other conditions.

* The specimens with the lowest initial saturation experienced the most rapid drying process. This may be due to the fact that vapour diffusion is the dominant flow pattern and a higher relative gas permeability was possessed by the specimen with a lower saturation.
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Flow patterns (height (z) versus half symmetrical geometry (radius, r)); (a) liquid flow (fine sand-20%), (b) vapour flow (fine sand-20%), (c) liquid flow (fine sand-50%) and (d) vapour flow (fine sand-50%).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Evolution of estimated relative gas (Krg) and liquid (Krw) permeabilities over time for different initial degrees of saturation.</figcaption></figure>

* Vaporisation and movement of vapour are the main mechanisms of mass transfer in soil specimens with the grain sizes and degrees of saturation investigated (fine sand and silty fine sand).

* Especially at the lower moisture contents, the lower liquid water permeability and higher vapour permeability are likely to be key factors governing the thermal behaviour of soil specimens.

* Re-distribution of liquid moisture due to evaporation and re-condensation was found to play a lesser role in the cases investigated. However, there was clear evidence for this process occurring in the fine sand, especially at the lowest degrees of saturation. The soil specimens in the tests were not fully sealed, which may or may not be representative of field conditions.

* Overall, the changes of moisture content that occurred within the specimens were relatively small. Nonetheless, this could result in thermal conductivity reductions of up to 50% for silts and sands with low initial saturations.
<figure><img src="image_18.png" /><figcaption>Fig. 15. Changes of overall vapour density with degree of saturation; (a) fine sand-20%, (b) fine sand-50%.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Example of isochrones of temperature distribution through the specimen.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. Distribution of local temporal saturation from numerical simulations for 20% fine sand.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Relationship between thermal conductivity and degree of soil saturation, based on data from Ref. 64.</figcaption></figure>

## CRediT authorship contribution statement

**Kui Liu:** Methodology, Formal analysis, Investigation, Validation, Data curation, Writing – original draft, Visualization. **Fleur A. Loveridge:** Conceptualization, Writing – review & editing, Supervision, Funding acquisition, Visualization. **Richard Boardman:** Writing – review & editing, Methodology, Supervision. **William Powrie:** Writing – review & editing, Supervision, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data access statement

The data associated with this paper are openly available from the University of Leeds Data Repository. https://doi.org/10.5518/1200.

## Acknowledgements

The work reported in this paper forms a part of a project funded by the Royal Academy of Engineering, the Doctoral Training Centre at University of Southampton and EPSRC (EP/G036896/1).

## References

1. Angelotti A, Sterpi D. On the performance of energy walls by monitoring assessment and numerical modelling: a case in Italy. *Environ Geotech.* 2018;7(4):266-273. http://dx.doi.org/10.1680/jenge.18.00037.

2. Farabi-Asl H, Fujii H, Kosukegawa H. Semi-open loop ground source heat pump system: Space-heating tests, numerical modeling and sensitivity analysis. *J Geotherm Res Soc Jpn.* 2017;39(2):81-92. http://dx.doi.org/10.11367/grsj.39.81.

3. Choi JC, Park J, Lee SR. Numerical evaluation of the effects of groundwater flow on borehole heat exchanger arrays. *Renew Energy.* 2013;52:230-240. http://dx.doi.org/10.1016/j.renene.2012.10.028.

4. Joshua WD, Jong ED. Soil moisture movement under temperature gradients. *Can J Soil Sci.* 1973;53(1):49-57. http://dx.doi.org/10.4141/cjss73-006.

5. Fillion M, Côté J, Konrad J. Thermal radiation and conduction properties of materials ranging from sand to rock-fill. *Can Geotech J.* 2011;48(4):532-542. http://dx.doi.org/10.1139/t10-093.

6. Farouki O. *Thermal Properties of Soils*. Clausthal-Zellerfeld, Germany: Trans Tech; 1986:720. In: Series on Rock and Soil Mechanics; Vol. 11.

7. Hutcheon WL. *Moisture Flow Induced by Thermal Gradients Within Unsaturated Soils (Ph.D thesis).* University of Minnesota; 1955.

8. Lekshami KRA, ArnPali DN. A review on coupled heat and water vapour transport in unsaturated soils. In: *Geotechnical Frontiers*. Orlando, Florida: 2017. http://dx.doi.org/10.1061/9780784480472.079.

9. Philip JR, De Vries DA. Moisture movement in porous materials under temperature gradients. *EOS Trans Am Geophys Union*. 1957;38(2):222-232. http://dx.doi.org/10.1029/TR038i002p00222.

10. Parlange MB, Cahill AT, Nielsen DR, Hopmans JW, Wendroth O. Review of heat and water movement in field soils. *Soils Tillage Res.* 1998;47(1-2):5-10. http://dx.doi.org/10.1016/S0167-1987(98)00066-X.

11. De Vries DA. The theory of heat and moisture transfer in porous media revisited. *Int J Heat Mass Transfer*. 1987;30(7):1343-1350. http://dx.doi.org/10.1016/0017-9310(87)90166-9.

12. Smits KM, Cihan A, Sakaki T, Illangasekare TH. Evaporation from soils under thermal boundary conditions: Experimental and modeling investigation to compare equilibrium-and nonequilibrium-based approaches. *Water Resour Res.* 2011;47(5). http://dx.doi.org/10.1029/2010WR009533.

13. Bear J. *Dynamics of Fluids in Porous Media*. New York: American Elsevier Pub; 1972.

14. Zhang J, Datta A. Some considerations in modeling of moisture transport in heating of hygroscopic materials. *Drying Technol.* 2004;22:1983-2008. http://dx.doi.org/10.1081/DRT-200032740.

15. Moradi A, Smits KM, Massey J, Cihan A, McCartney. Impact of coupled heat transfer and water flow on soil borehole thermal energy storage (SBTES) systems: Experimental and modeling investigation. *Geothermics*. 2015;57:56-72. http://dx.doi.org/10.1016/j.geothermics.2015.05.007.

16. Moradi A, Smits KM, Lu N, McCartney JS. Heat transfer in unsaturated soil with application to borehole thermal energy storage. *Vadose Zone J.* 2016;15(10). http://dx.doi.org/10.2136/vzj2016.03.0027.

17. Buzug TM. Computed tomography. In: Springer Handbook of Medical Technology. Berlin, Heidelberg: Springer; 2011.

18. Biguri A, Lindroos R, Bryll R, Towsyfyan H, Deyhle H, El Khalil Harrane I, Boardman R, Mavrogordato M, Dosanjh M, Hancock S, Blumensath T. Arbitrarily large tomography with iterative algorithms on multiple GPUs using the TIGRE toolbox. *J Parallel Distrib Comput.* 2020;146:52-63. http://dx.doi.org/10.1016/j.jpdc.2020.07.004.

19. Atwood RC, Bodey AJ, Price SW, Basham M, Drakopoulos M. A high-throughput system for high-quality tomographic reconstruction of large datasets at Diamond light source. *Phil Trans R Soc A.* 2015;373(2043):20140398. http://dx.doi.org/10.1098/rsta.2014.0398.

20. Hall SA, Bornert M, Desrues J, Pannier Y, Lenoir N, Viggiani G, Bésuelle P. Discrete and continuum analysis of localised deformation in sand using X-ray µCT and volumetric digital image correlation. *Géotechnique.* 2010;60(5):315-322. http://dx.doi.org/10.1680/geot.2010.60.5.315.

21. Cromwell V, Kortum DJ, Bradley DJ. The use of a medical computer tomography (CT) system to observe multiphase flow in porous media. In: *SPE Annual Technical Conference and Exhibition*. OnePetro; 1984. http://dx.doi.org/10.2118/13098-MS.

22. Hirono TMSF, Takahashi M, Nakashima S. Direct imaging of fluid flow in fault-related rocks by X-ray CT. *Geol Soc Lond Spec Publ.* 2003;215(1):107-115. http://dx.doi.org/10.1144/GSL SP.2003.215.01.10.

23. Sato A, Sawada A. Analysis of tracer migration process in the crack by means of X-ray CT. In: Proceedings of the 11th International Congress on Rock Mechanics, Lisbon, 2007, p. 15-8.

24. Andò E, Hall SA, Viggiani G, Desrues J, Bésuelle P. Experimental micromechanics: grain-scale observation of sand deformation. *Geotech Lett.* 2012;2(3):107-112.

25. Beckett CTS, Hall MR, Augarde CE. Macrostructural changes in compacted earthen construction materials under loading. *Acta Geotech.* 2013;8(4):423-438.

26. Taina IA, Heck RJ, Elliot TR. Application of X-ray computed tomography to soil science: A literature review. *Can J Soil Sci.* 2008;88(1):1-19.

27. Van Geet M, Swennen R, Wevers M. Quantitative analysis of reservoir rocks by microfocus X-ray computerised tomography. *Sediment Geol.* 2000;132(1-2):25-36.

28. Naughton PJ, O'Kelly BC. The anisotropy of Leighton Buzzard sand under general stress conditions. In: Proceedings Third International Symposium on Deformation Characteristics of Geomaterials, Vol. 1, 2003, p. 285-91.

29. Raghunandan M, Juneja A, Hsiung B. Preparation of reconstituted sand samples in the laboratory. *Int J Geotech Eng.* 2012;6(1):125-131. http://dx.doi.org/10.3328/IJGE.2012.06.01.125-131.

30. Miller D. Parafilm: Frequently asked questions. 2012 [Online]. Available: https://camblab.info/wp/index.php/parafilm-frequently-asked-questions/ [Accessed 17th September 2018].

31. University of Southampton. µ-VIS: Multidisciplinary, multiscale, micro-tomographic volume imaging. 2017 [Online]. Available: http://www.southampton.ac.uk/muvis/about/index.page [Accessed 26th March 2017].

32. Liu K, Boardman R, Mavrogordato M, Loveridge FA, Powrie W. The importance of the heel effect in X-ray CT imaging of soils. *Environ Geotech.* 2020;40(XXXX):1-16. http://dx.doi.org/10.1680/jenge.20.00048.

33. Feldkamp LA, Davis LC, Kress JW. Practical cone-beam algorithm. *Josa A.* 1984;1(6):612-619. http://dx.doi.org/10.1364/JOSAA.1.000612.
34. Metris. *CT Pro User Manual [Online]*. Hertfordshire: 2008 Available: http://muvis.soton.ac.uk/wiki/images/1/1d/CT_Pro_User_Manual.pdf [Accessed 31st May 2015].

35. Gordon R. A tutorial on ART (algebraic reconstruction techniques). *IEEE Trans Nucl Sci.* 1974;21(3):78-93. http://dx.doi.org/10.1109/TNS.1974.6499238.

36. Liu K, Loveridge FA, Boardman R, Powrie W. Evaluation of reconstruction and segmentation techniques on high temporal resolution µCT scans for geotechnical applications. In: *International Symposium on Energy Geotechnics (SEG-2018)*. Lausanne, Switzerland: 2018.

37. Ketcham RA, Hildebrandt J. Characterizing, measuring, and utilizing the resolution of CT imagery for improved quantification of fine-scale features. *Nucl Instrum Methods Phys Res B.* 2014;324:80-87. http://dx.doi.org/10.1016/j.nimb.2013.08.064.

38. Liu K. *Computed Tomography (CT) Scanning to Quantify the Movement of Moisture in Soils under Thermal Gradients* (EngD thesis). UK: University of Southampton; 2020.

39. Liu K, Loveridge FA, Boardman R, Powrie W. Study of short-term evaporation in sand specimens via micro-focus X-ray computed tomography. In: *2nd International Symposium on Coupled Phenomena in Environmental Geotechnics (CPEG2)*. Leeds, UK: 2017.

40. COMSOL. COMSOL Modeling Software [Online]. Comsol, Inc: 2015 Available: COMSOL.com [Accessed 22nd May 2017].

41. McCartney JS, Baser T. Role of coupled processes in thermal energy storage in the vadose zone. In: *2nd Symposium on Coupled Phenomena in Environmental Geotechnics (CPEG2)*. Leeds, UK: 2017.

42. Richards LA. Capillary conduction of liquids through porous mediums. *Physics*. 1931;1(5):318-333. http://dx.doi.org/10.1063/1.1745010.

43. Van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci Am J.* 1980;44(5):892-898. http://dx.doi.org/10.2136/sssj1980.03615995004400050002x.

44. Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. *Water Resour Res.* 1976;12(3):513-522. http://dx.doi.org/10.1029/WR012i003p00513.

45. Moldrup P, Olesen T, Yoshikawa S, Komatsu T, Rolston DE. Three-porosity model for predicting the gas diffusion coefficient in undisturbed soil. *Soil Sci Am J.* 2004;68(3):750-759. http://dx.doi.org/10.2136/sssj2004.7500.

46. Ungureanu A, Statescu F. Models for predicting the gas diffusion coefficient in undisturbed soil. *Trans Hydrotech.* 2010;55(69):168-172.

47. Burdine N. Relative permeability calculations from pore size distribution data. *J Pet Technol.* 1953;5(03):71-78. http://dx.doi.org/10.2118/225-G.

48. Campbell GS. A simple method for determining unsaturated conductivity from moisture retention data. *Soil Sci.* 1974;117(6):311-314.

49. Moldrup P, Kruse CW, Rolston DE, Yamaguchi T. Modeling diffusion and reaction in soils: III. Predicting gas diffusivity from the Campbell soil-water retention model. *Soil Sci.* 1996;161(6):366-375.

50. Campbell GS, Jungbauer Jr JD, Bidlake WR, Hungerford RD. Predicting the effect of temperature on soil thermal conductivity. *Soil Sci.* 1994;158(5):307-313.

51. Haigh SK. Thermal conductivity of sands. *Géotechnique*. 2012;62(7):617-625. http://dx.doi.org/10.1680/geot.11.P.043.

52. Tarnawski VR, Momose T, Leong WH, Bovesecchi G, Coppa P. Thermal conductivity of standard sands. Part I. Dry-state conditions. *Int J Thermophys.* 2009;30(3):949-968. http://dx.doi.org/10.1007/s10765-009-0596-0.

53. Horai KI. Thermal conductivity of rock - forming minerals. *J Geophys Res.* 1971;76(5):1278-1308. http://dx.doi.org/10.1029/JB076i005p01278.

54. Mitchell J. Conduction phenomena: from theory to geotechnical practice. *Géotechnique*. 1991;41(3):299-340. http://dx.doi.org/10.1680/geot.1991.41.3.299.

55. Suryatriyastuti ME, Mroueh H, Burlon S. Understanding the temperature-induced mechanical behaviour of energy pile foundations. *Renew Sustain Energy Rev.* 2012;16(5):3344-3354. http://dx.doi.org/10.1016/j.rser.2012.02.062.

56. Hillel D. *Fundamentals of Soil Physics*. New York: Academic Press; 1980.

57. Thomas HR, King SD. A non-linear, two-dimensional, potential-based analysis of coupled heat and mass transfer in a porous medium. *Internat J Numer Methods Engrg.* 1994;37(21):3707-3722. http://dx.doi.org/10.1002/nme.1620372108.

58. Monteith J, Unsworth M. *Principles of Environmental Physics*. Burlington, MA: Elsevier Science; 2008.

59. Campbell GS. *Soil Physics with BASIC: Transport Models for Soil-Plant Systems*. New York: Elsevier; 1985.

60. Bixler NE. *NORIA-a Finite Element Computer Program for Analyzing Water, Vapor, Air, and Energy Transport in Porous Media*. NASA STI/Recon Technical Report N, 86; 1985:16531.

61. Rohsenow WM, Hartnett JP, Cho YI. *Handbook of Heat Transfer*. New York: McGraw-Hill; 1998.

62. Marshall TJ, Holmes JW. *Soil Physics*. Cambridge, UK: Cambridge University Press; 1979.

63. Gens A, Sánchez M, Guimarães LDN, Alonso EE, Lloret A, Olivella S, Villar MV, Huertas F. A full-scale in situ heating test for high-level nuclear waste disposal: observations, analysis and interpretation. *Géotechnique*. 2009;59(4):377-399. http://dx.doi.org/10.1680/geot.2009.59.4.377.

64. Dong Y, McCartney JS, Lu N. Critical review of thermal conductivity models for unsaturated soils. *Geotech Geol Eng.* 2015;33(2):207-221.