# Analysis of heat transfer and water flow with phase change in saturated
porous media by bond-based peridynamics

Petr Nikolaev  Corresponding author<EMAIL_ADDRESS>Department of Mechanical, Aerospace and Civil Engineering, School of
Engineering, The University of Manchester, United Kingdom Majid Sedighi
Corresponding author<EMAIL_ADDRESS>Department of Mechanical,
Aerospace and Civil Engineering, School of Engineering, The University of
Manchester, United Kingdom Andrey P Jivkov Department of Mechanical, Aerospace
and Civil Engineering, School of Engineering, The University of Manchester,
United Kingdom Lee Margetts  Department of Mechanical, Aerospace and Civil
Engineering, School of Engineering, The University of Manchester, United
Kingdom

###### Abstract

A wide range of natural and industrial processes involve heat and mass
transport in porous media. In some important cases the transported substance
may undergo phase change, e.g. from liquid to solid and vice versa in the case
of freezing and thawing of soils. The predictive modelling of such phenomena
faces physical (multiple physical processes taking place) and mathematical
(evolving interface with step change of properties) challenges. In this work,
we develop and test a non-local approach based on bond-based peridynamics
which addresses the challenges successfully. Our formulation allows for
predicting the location of the interface between phases, and for calculating
the temperature and pressure distributions within the saturated porous medium
under the conditions of pressure driven water flow. The formulation is
verified against existing analytical solutions for 1D problems, as well as
finite element transient solutions for 2D problems. The agreement found by the
verification exercise demonstrates the accuracy of the proposed methodology.
The detailed coupled description of heat and hydraulic processes can be
considered as a critical step towards a thermo-hydro-mechanical model, which
will allow, for example, description of the hydrological behaviour of
permafrost soils and the frost heave phenomenon.

Keywords:  heat conduction; coupled problem; phase change; peridynamics;
freezing; melting; porous medium

##  1 Introduction

Heat and mass transfer with phase change occur in many natural and engineered
solids. A non-exhaustive list of examples include seasonal and artificial
freezing of soils, casting of metals, polymers and ceramics, and latent heat
thermal energy storage.

Our focus is on porous media, with specific application to soils. In such
systems, the phenomenon emerges from the coupling the multi-phase heat
transfer with fluid flow through a porous substance. This involves strong
physical non-linearities (rapid change of physical and mechanical properties
between phases) and geometric discontinuities (inter-phase boundaries,
cracks), making the mathematical description a challenging task. The lack of
general analytical solutions [1] calls for appropriate numerical methods with
specific attention to the strong non-linearities. To date, this problem has
been approached by numerical methods based on the local (differential)
formulation of the processes involved, e.g. by the finite element method
(FEM). Such formulations are used in many engineering applications, for
example in the analysis of: underground water flow in permafrost regions [2,
3]; seasonally frozen soils [4, 5]; evolution of methane hydrates under the
seabed [6]; nuclear waste disposal design [7, 8]; and artificial ground
freezing for mining and civil engineering purposes [9, 10, 11, 12].

However, accounting for the strong non-linearities and discontinuities is
challenging for methods that are based on local mathematical formulations. The
continued growth of computational power opens the opportunity to address these
challenges by using methods based on non-local formulations, which are
considered more computationally demanding. One non-local approach with
increasing popularity is the Peridynamics (PD). In PD, the partial
differential equations of any classical local theory are replaced by a set of
integral-differential equations [13], resulting in a mathematically consistent
formulation, even in the presence of strong non-linearities and
discontinuities.

Peridynamics was originally proposed for describing mechanical behavior of
solids [14, 15], and subsequently extended to a variety of diffusion problems
[16, 13, 17]. It was first developed as a bond-based PD, and later generalized
to a state-based PD in two versions - ordinary and non-ordinary [18, 13].
Recently, an element-based PD formulation was proposed [19]. A correspondence
between PD and continuum formulations is made by the concept of the
peridynamics differential operator (PDDO) [20, 21].

In this paper, we construct a bond-based PD model for transient heat transfer
coupled with water flow and phase change in saturated porous media. The first
PD formulation of heat conduction, which also considered electromigration, was
made in [22]. Bond-based formulations of transient heat conduction were
originally developed in [23, 16], see also [24, 25]. Further developments of
the bond-based PD heat transfer models were discussed in [26, 27, 28].
However, to date, the phase change phenomenon has not been included in the
developments by this approach. In [29], the one-dimensional phase change
problem was solved using the PDDO framework, which prompted us to develop a
bond-based PD model for obtaining such solutions when coupled with fluid flow.

A PD formulation of pressure driven water flow in porous media was originally
discussed in [30, 31]. This approach was recently used to develop a wide range
of coupled problems that consider chemical diffusion in partially and fully
saturated porous media e.g. [32, 33, 34, 35]. These demonstrated the
simplicity and universality of the approach for dealing with such problems.
However, the process of thermal diffusion has not yet been coupled with water
flow in the analysis of heat and mass transfer in saturated porous media.

Our development solves the outstanding issues of coupling between different
processes, following ideas from [32, 34]. The present paper is structured as
follows. In Section 2 we consider the classical local differential equations
of heat transfer with phase change under pressure driven water flow in
saturated porous media. In Section 3, we describe the basic concepts of bond-
based peridynamics, and extend them to the case of heat transfer with phase
change. Here, the formulation is also coupled with the existing model for
pressure driven water flow. Section 4 presents our numerical implementation of
the developed PD model for solution of 1D and 2D problems on regular square
meshes. In Section 5, we consider several test problems and compare the
solutions based on our model with existing analytical (1D) and numerical (2D
by finite elements) solutions. The agreement between results demonstrates the
accuracy of the proposed methodology. Importantly, our model is also
successfully tested for 2D heat transfer with phase change with water flow due
to high pressure gradients, a condition that is challenging for other methods
including finite elements. Conclusions are drawn in Section 6.

The developed PD approach can be coupled with mechanical models that describe,
for example, ocean ice cracking [36, 37, 38], cracking of concrete structures
and rocks [35, 39, 40], as well as corrosion pitting [41].

##  2 Thermo-hydraulic model of porous medium with phase change

In this study we consider a fully saturated porous medium with water as the
pore-filling liquid. However, the considerations are applicable to any
Newtonian fluid. The porous medium consists of three phases - solid matrix
particles, liquid water, and solid water (ice). The physical parameters of
these phases will be presented with lower indices $s$, $w$, and $i$,
respectively. The hydraulic behavior of the liquid phase is governed by the
mass conservation law, and the thermal behavior of the three-component
structure is governed by the energy conservation law. The equations describing
these laws can be constructed in the framework of the volume average method,
see, for example, [9]. The equations in our presentation are written with
respect to a Cartesian coordinate system.

The conservation of water mass in the two phases can be written as

$\frac{\partial}{\partial
t}\left(n\rho_{w}S_{w}+n\rho_{i}S_{i}\right)=-\nabla\cdot\left(\rho_{w}\mathbf{u}\right),$
(1)

where $S(\mathbf{x},t)$ is the saturation degree, for which $S_{w}+S_{i}=1$;
$n$ is the porosity; $\rho$ is the density; $\mathbf{u}(\mathbf{x},t)$ is the
flux of liquid water; $\mathbf{x}$ is the coordinate vector of a material
point; and $t$ is the time.

In this work we assume a non-deformable pore system, so that the porosity $n$
is constant. The compressibility of ice is much lower than the compressibility
of water, allowing us to consider the ice density to be approximately
constant. With these assumptions, and taking into account that
$S_{w}+S_{i}=1$, Eq. (1) takes the form:

$nS_{w}\frac{\partial\rho_{w}}{\partial
t}+n\left(\rho_{w}-\rho_{i}\right)\frac{\partial S_{w}}{\partial
t}=-\nabla\cdot\left(\rho_{w}\mathbf{u}\right).$ (2)

The water flow is described by Darcy’s law:

$\mathbf{u}(\mathbf{x},t)=-\frac{k_{r}k_{int}}{\mu}\nabla\Phi,$ (3)

where $k_{r}(\mathbf{x},t)$ is the liquid water relative permeability;
$k_{int}$ is the porous medium intrinsic permeability; $\mu$ is the liquid
water dynamic viscosity; and $\Phi(\mathbf{x},t)$ is the flow potential, which
in terms of fluid pressure is given by:

$\Phi(\mathbf{x},t)=p+g\rho_{w}z,$ (4)

where $p(\mathbf{x},t)$ is the pressure; $g$ is the gravitational acceleration
and $z(\mathbf{x},t)$ is the height of a liquid column. A number of relations
have been proposed between the relative permeability, $k_{r}$ and the liquid
water content. Most of the common relations for the frozen soils are presented
in [42, 43]. In this paper, the relative permeability $k_{r}$ is defined
according to [44] by:

$k_{r}=10^{-n\Omega\left(1-S_{w}\right)},$ (5)

where $\Omega$ is an empirical parameter.

Substitution of Eq. (3) into Eq. (2) yields:

$nS_{w}\dfrac{\partial\rho_{w}}{\partial
t}+n\left(\rho_{w}-\rho_{i}\right)\frac{\partial S_{w}}{\partial
t}=\frac{k_{int}}{\mu}\nabla\cdot\left(\rho_{w}k_{r}\nabla\Phi\right),$ (6)

which is the complete equation for the conservation of mass for water flow
given by Eq. (4) and considering phase change.

The energy conservation involves heat conduction and heat convection by a
moving medium with phase change. Following [9], this is presented by

$\left(\rho C\right)_{eq}\frac{\partial T}{\partial
t}=-\nabla\cdot\mathbf{q}-\mathbf{u}\cdot\nabla\left(\rho_{w}C_{w}nS_{w}T\right),$
(7)

where $\mathbf{q}(\mathbf{x},t)$ is the heat flux due to heat conduction;
$T(\mathbf{x},t)$ is the temperature; $C$ is the specific heat capacity at a
constant pressure; and $\left(\rho C\right)_{eq}(\mathbf{x},t)$ is the
equivalent heat capacity of the three-phase medium defined by

$\left(\rho
C\right)_{eq}=nS_{w}\rho_{w}C_{w}+n\left(1-S_{w}\right)\rho_{i}C_{i}+\left(1-n\right)\rho_{s}C_{s}+n\rho_{i}L\frac{\partial
S_{w}}{\partial T},$ (8)

where $L$ is the latent heat of water solidification.

The dependence of the liquid water content on temperature below the freezing
point is different for different materials. Specifically for soils,
approximate relations for the soil freezing characteristic curves are
discussed in [45, 46]. In the present study we consider a Weibull-type
relation, used in [44], according to which

$S_{w}=\left\\{\begin{array}[]{c}{1}\\\
{(1-S_{w_{res}})e^{-\left(\frac{T-T_{f}}{\Delta
T}\right)^{2}}+S_{w_{res}}}\end{array}\quad\begin{array}[]{c}{T>T_{f}}\\\
{T\leq T_{f}}\end{array}\right.$ (9)

where $T_{f}$ is temperature at the onset of water solidification;
$S_{w_{res}}$ is the residual saturation; $\Delta T$ is the Weibull scale
parameter governing the smoothness of the $S_{w}$ function (the Weibul shape
parameter is 2). This relation ensures fast and effective calculation. The
derivative of $S_{w}$ with respect to temperature below freezing point is
required for calculating the equivalent heat capacity in Eq. (8). This is
given by

$\frac{\partial S_{w}}{\partial T}=-2(1-S_{w_{res}})\frac{T-T_{f}}{\Delta
T^{2}}e^{-\left(\frac{T-T_{f}}{\Delta T}\right)^{2}}.$ (10)

The heat flux due to the heat conduction is governed by Fourier’s law:

$\mathbf{q}(x,t)=-\lambda\nabla T,$ (11)

where $\lambda(\textbf{x},t)$ is the material thermal conductivity. The
thermal conductivity of the three-phase medium is given by

$\lambda=\lambda_{w}^{nS_{w}}\lambda_{i}^{n\left(1-S_{w}\right)}\lambda_{s}^{\left(1-n\right)}.$
(12)

Substitution of Eqns. (3), (8) and (11) into Eq. (7) gives:

$\left[nS_{w}\rho_{w}C_{w}+n\left(1-S_{w}\right)\rho_{i}C_{i}+\left(1-n\right)\rho_{s}C_{s}+n\rho_{i}L\frac{\partial
S_{w}}{\partial T}\right]\frac{\partial T}{\partial t}=\\\
-\nabla\left(\lambda\nabla T\right)+\rho_{w}C_{w}\nabla
T\cdot\frac{k_{r}k_{int}}{\mu}\nabla\Phi,$ (13)

which is the complete equation for the conservation of energy under flow
potential given by Eq. (4) and considering phase change.

Equations (6) and (13) describe the processes of heat transfer and water flow
with phase change in a local (differential) form. These will be reformulated
in a non-local form in the next section.

##  3 Bond-based peridynamics formulations for heat transfer with phase change
and water flow

In bond-based Peridynamics, a body occupying a region $\mathbf{R}$ is
considered as a collection of an arbitrary number of individual particles with
associated volumes and masses. Each peridynamic particle is labeled by a
position vector $\mathbf{x}$ with respect to a fixed (background) Cartesian
coordinate system. A particle at position x interacts with (and is connected
to) all particles at positions $\mathbf{x}^{\prime}$ within a certain spatial
region $\mathbf{H_{x}}$ referred to as the horizon of the particle at position
x. The horizon may have any shape and size [25], and these may vary with the
position $\mathbf{x}$, but in the majority of applications, as in this work,
it is assumed to be of the form
$\mathbf{H_{x}}=\\{\textbf{x}^{\prime}:||\textbf{x}^{\prime}-\textbf{x}||\leq\delta\\}$,
where the constant $\delta$ is referred to as the horizon radius.

With this setting, the term ’bond’ refers to the interactions between two
particles at spatial positions $\mathbf{x}$ and
$\mathbf{x}^{\prime}\in\mathbf{H_{x}}$. Considering the nature of the
processes to be modelled, we will refer to the bonds as transport or
’t-bonds’. The peridynamic heat and fluid flux densities (fluxes per unit
volume) along a ’t-bond’ depend on the distance between points x and x’,
represented by the distance vector
$\mathbf{\xi}=\left(\mathbf{x^{\prime}}-\mathbf{x}\right)$.

With respect to the water phases, the region $\mathbf{R}$ can be divided into
two regions – solid and liquid – as illustrated in Fig. 1. The solid region
contains solid water (ice), and it may also contain the liquid water, which
amount is defined by the soil freezing characteristic curve (9). The liquid
region does not contain the solid water. The regions are separated by the
surface with the temperature $T_{f}$.

Figure 1: Illustration of phase regions, particle horizons and transport bonds
(‘t-bonds’)

Depending on the region, where connected particles are located, we define
three different PD t-bonds: solid t-bonds, interfacial t-bonds connecting
different regions, and liquid t-bonds, see Fig. 1. If it is necessary, their
transport properties can be considered differently. Thus, for a porous
material, in which the solid region does not contain liquid water or its
transfer can be neglected, due to the very low relative permeability
associated with it (in the present study defined by relation (5)), it is
suitable to assume that the solid t-bonds transport only heat, and do not
participate in liquid transfer. Such bonds can be called ’impermeable’. The
liquid and interfacial t-bonds transfer both heat and water. The conditions
that attribute t-bonds to the ’impermeable’ type can be chosen based on the
analysis of considered problem. The t-bonds classification may be used to
improve efficiency of the developed computational codes by not considering
water flow through impermeable t-bonds.

###  3.1 PD formulation of water flow

The PD fluid flux density (flux per unit volume) along a t-bond between
particles at positions $\mathbf{x}$ and $\mathbf{x}^{\prime}$ per unit volume
of $\mathbf{x}^{\prime}$ is described by [30]:

$\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\mu}\frac{\Phi\left(\mathbf{x}^{\prime},t\right)-\Phi\left(\mathbf{x},t\right)}{\left\|\xi\right\|}\frac{\xi}{\left\|\xi\right\|},$
(14)

where $\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ is the
volumetric flux of liquid water; $\Phi\left(\mathbf{x},t\right)$ and
$\Phi\left(\mathbf{x}^{\prime},t\right)$ are flow potentials at points
$\mathbf{x}$ and $\mathbf{x}^{\prime}$, respectively;
$\mathbf{\xi}/\left\|\mathbf{\xi}\right\|$ is the unit vector along the
t-bond; and $\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ is
the relative permeability of the t-bond between the points $\mathbf{x}$ and
$\mathbf{x}^{\prime}$.

Since the t-bond relative permeability depends on the locations of its end
points within the two different regions, which properties can vary
significantly, we define it in the following way:

$\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{2k_{r}\left(\mathbf{x}^{\prime},t\right)k_{r}\left(\mathbf{x},t\right)}{k_{r}\left(\mathbf{x}^{\prime},t\right)+k_{r}\left(\mathbf{x},t\right)}$
(15)

This definition ensures that, the solid and interfacial t-bonds would have
$\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\approx 0$,
which cannot be achieved with a simple average.

The mass conservation equation for a t-bond between points $\mathbf{x}$ and
$\mathbf{x}^{\prime}$ is given by [30, 34]

$nS_{w}\frac{\partial\rho_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\partial
t}+R_{b}\left(\mathbf{x}^{\prime},\mathbf{x}\right)=\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\mu}\rho_{w}\frac{\Phi\left(\mathbf{x^{\prime}},t\right)-\Phi\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}},$
(16)

where $R_{b}\left(\mathbf{x}^{\prime},\mathbf{x}^{\prime}\right)$ is a
volumetric mass source term, that can be represented by

$R_{b}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=n\left(\rho_{w}-\rho_{i}\right)\frac{\partial
S_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\partial t}.$ (17)

The mass conservation at point $\mathbf{x}$ is obtained by integrating Eq.
(16) over its horizon $\mathbf{H_{x}}$

$\int_{\mathbf{H_{x}}}nS_{w}\frac{\partial\rho_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\partial
t}{\rm
d}V_{\mathbf{x}^{\prime}}+\int_{\mathbf{H_{x}}}n\left(\rho_{w}-\rho_{i}\right)\frac{\partial
S_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\partial t}{\rm
d}V_{\mathbf{x}^{\prime}}=\\\
\int_{\mathbf{H_{x}}}\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\mu}\rho_{w}\frac{\Phi\left(\mathbf{x}^{\prime},t\right)-\Phi\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}{\rm
d}V_{\mathbf{x}^{\prime}}$ (18)

The relation between the water density at point $\mathbf{x}$ and time $t$ and
the water density in all adjacent t-bonds is given by [34]

$\int_{\mathbf{H_{x}}}nS_{w}\frac{\partial\rho_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\partial
t}{\rm
d}V_{\mathbf{x}^{\prime}}=nS_{w}\frac{\partial\rho_{w}\left(\mathbf{x},t\right)}{\partial
t}V_{\mathbf{x}}$ (19)

Similarly, according to [30], the mass source term at point $\mathbf{x}$ can
be obtained as the average in all adjacent t-bonds:

$\int_{\mathbf{H_{x}}}n\left(\rho_{w}-\rho_{i}\right)\frac{\partial
S_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\partial t}{\rm
d}V_{\mathbf{x}^{\prime}}=n\left(\rho_{w}-\rho_{i}\right)\frac{\partial
S_{w}\left(\mathbf{x},t\right)}{\partial t}V_{\mathbf{x}}$ (20)

Combining Eqns. (16) – (20) leads to the mass conservation equation for liquid
water at point $\mathbf{x}$ in terms of density and flow potential

$nS_{w}\frac{\partial}{\partial
t}\rho_{w}\left(\mathbf{x},t\right)+n\left(\rho_{w}-\rho_{i}\right)\frac{\partial}{\partial
t}S_{w}\left(\mathbf{x},t\right)=\\\
\int_{\mathbf{H_{x}}}\frac{\rho_{w}}{\mu}\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{V_{\mathbf{x}}}\frac{\Phi\left(\mathbf{x}^{\prime},t\right)-\Phi\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}{\rm
d}V_{\mathbf{x}^{\prime}}$ (21)

Under the assumption that $\rho_{w}\left(\mathbf{x},t\right)=const$, the last
equation becomes

$n\left(\rho_{w}-\rho_{i}\right)\frac{\partial}{\partial
t}S_{w}\left(\mathbf{x},t\right)=\int_{\mathbf{H_{x}}}\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{V_{\mathbf{x}}}\frac{\rho_{w}}{\mu}\frac{\Phi\left(\mathbf{x}^{\prime},t\right)-\Phi\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}{\rm
d}V_{\mathbf{x}^{\prime}}$ (22)

The PD micro-permeability of a t-bond can be defined as:

$K_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\mu
V_{\mathbf{x}}}$ (23)

Using the approach of uniform constant influence functions [32] the relations
between PD micro-hydraulic conductivity
$K_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ and the macroscopic
properties for 1D and 2D cases is derived as [34]

$K_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\mu\delta}$
(24)
$K_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{4k_{int}\overline{k_{r}}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\mu\pi\delta^{2}}$
(25)

For the known pressure field, that can be found by solution of the equation
(21), following [30], we can define the water flux at any point x by the
relation:

$\mathbf{u}(\mathbf{x},t)=\int_{\mathbf{H_{x}}}K_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\rho_{w}\frac{\Phi\left(\mathbf{x}^{\prime},t\right)-\Phi\left(\mathbf{x},t\right)}{\left\|\xi\right\|}\frac{\mathbf{\xi}}{\left\|\mathbf{\xi}\right\|}{\rm
d}V_{\mathbf{x}^{\prime}}$ (26)

For the 2-D case, the vector $\mathbf{u}(\mathbf{x},t)$ can be considered as
two components: $u_{x}(\mathbf{x},t)$ and $u_{y}(\mathbf{x},t)$ along the $x$
and $y$ axis, respectively.

###  3.2 PD formulation of heat transfer with phase change

The heat flux along a t-bond,
$J\left(T,\mathbf{x}^{\prime},\mathbf{x},t\right)$, can be defined according
to [32, 34] as:

$J=\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\mathbf{\xi}\right\|}\cdot\frac{\mathbf{\xi}}{\left\|\mathbf{\xi}\right\|}-\\\
\rho_{w}C_{w}\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\left[T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)\right]\cdot\frac{\mathbf{\xi}}{\left\|\mathbf{\xi}\right\|}$
(27)

where $T\left(\mathbf{x}^{\prime},t\right)$ and $T\left(\mathbf{x},t\right)$
are the temperature at points $\mathbf{x}^{\prime}$ and $\mathbf{x}$ at time
$t$, respectively;
$\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ is the
thermal conductivity of the t-bond; and
$\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ is the average
water flux along the t-bond.

The thermal conductivity of a t-bond depends on the locations of the two
connected particles in the three considered regions, and we define it with a
simple interpolation by

$\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{\lambda(\mathbf{x},t)+\lambda(\mathbf{x}^{\prime},t)}{2},$
(28)

where $\lambda(\mathbf{x},t)$ and $\lambda(\mathbf{x}^{\prime},t)$ are the
thermal conductivities at the peridynamic particles located at $\mathbf{x}$
and $\mathbf{x}^{\prime}$, respectively.

Similarly, and following [34, 47], the water velocity along a t-bond is
defined by

$\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{\mathbf{u}\left(\mathbf{x},t\right)+\mathbf{u}\left(\mathbf{x}^{\prime},t\right)}{2},$
(29)

where the liquid velocity vectors $\mathbf{u}\left(\mathbf{x},t\right)$ and
$\mathbf{u}\left(\mathbf{x}^{\prime},t\right)$ at particles at positions
$\mathbf{x}$ and $\mathbf{x}^{\prime}$, respectively, are found by solving the
equation (26) for a given pressure field. This relation ensures that both
liquid and inetrfacial types of t-bonds participate in convective heat
transfer.

The peridynamic version of the energy conservation equation is written for
each t-bond in the following form [16, 23, 24]

$\overline{\left(\rho
C\right)_{eq}}(\mathbf{x}^{\prime},\mathbf{x},t)\frac{\partial}{\partial
t}\overline{T}(\mathbf{x}^{\prime},\mathbf{x},t)=\frac{J\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\left\|\xi\right\|}$
(30)

where $\overline{T}(\mathbf{x}^{\prime},\mathbf{x},t)$ is the t-bond
temperature taken as the average of the temperatures at points $\mathbf{x}$
and $\mathbf{x}^{\prime}$, and $\overline{\left(\rho
C\right)_{eq}}(\mathbf{x}^{\prime},\mathbf{x},t)$ is the t-bond equivalent
heat capacity taken as the average of the equivalent heat capacities at points
$\mathbf{x}$ and $\mathbf{x}^{\prime}$.

Inserting Eq. (27) into Eq. (30) provides the energy conservation in a t-bond
in terms of temperature and water velocity

$\overline{\left(\rho
C\right)_{eq}}(\mathbf{x}^{\prime},\mathbf{x},t)\frac{\partial\overline{T}(\mathbf{x}^{\prime},\mathbf{x},t)}{\partial
t}=\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}\cdot\frac{\xi}{\left\|\xi\right\|}-\\\
\rho_{w}C_{w}\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|}\cdot\frac{\xi}{\left\|\xi\right\|}$
(31)

The energy conservation for a particle at position $\mathbf{x}$ involves the
fluxes in all adjacent t-bonds (bonds to particles within the horizon
$\mathbf{H_{x}}$) and is obtained by integrating Eq. (30) over the horizon:

$\int_{\mathbf{H_{x}}}\overline{\left(\rho
C\right)_{eq}}(\mathbf{x}^{\prime},\mathbf{x},t)\frac{\partial\overline{T}(\mathbf{x}^{\prime},\mathbf{x},t)}{\partial
t}{\rm d}V_{\mathbf{x}^{\prime}}=\\\
\int_{\mathbf{H_{x}}}\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}\cdot\frac{\xi}{\left\|\xi\right\|}-\\\
\rho_{w}C_{w}\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|}\cdot\frac{\xi}{\left\|\xi\right\|}]{\rm
d}V_{\mathbf{x}^{\prime}}$ (32)

where $V_{\mathbf{x}^{\prime}}$ is the horizon volume of the particle at
position $\mathbf{x}^{\prime}$.

The relationship between temperature at point $\mathbf{x}$ and time $t$ and
the temperature in all adjacent bonds can be written as [23]:

$\int_{\mathbf{H_{x}}}\overline{\left(\rho
C\right)_{eq}}(\mathbf{x}^{\prime},\mathbf{x},t)\frac{\partial\overline{T}(\mathbf{x}^{\prime},\mathbf{x},t)}{\partial
t}{\rm d}V_{\mathbf{x}^{\prime}}=\overline{\left(\rho
C\right)_{eq}}(\mathbf{x},t)\frac{\partial T(\mathbf{x},t)}{\partial
t}V_{\mathbf{x}}$ (33)

where $V_{\mathbf{x}}$ is the horizon volume of particle $\mathbf{x}$. Note,
that under the assumption of constant horizon radius
$V_{\mathbf{x}^{\prime}}=V_{\mathbf{x}}$, and specifically for 1D problems
$V_{\mathbf{x}}=2\delta$, and for 2D problems $V_{\mathbf{x}}=\pi\delta^{2}$.

Combining Eqns. (27) – (33) leads to the following equation for the evolution
of temperature at particle $\mathbf{x}$

$\overline{\left(\rho C\right)_{eq}}(\mathbf{x},t)\frac{\partial
T(\mathbf{x},t)}{\partial
t}=\int_{\mathbf{H_{x}}}\frac{\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{V_{\mathbf{x}}}\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}\cdot\frac{\xi}{\left\|\xi\right\|}-\\\
\rho_{w}C_{w}\frac{\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{V_{\mathbf{x}}}\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|}\cdot\frac{\xi}{\left\|\xi\right\|}{\rm
d}V_{\mathbf{x}^{\prime}}$ (34)

The peridynamics microscopic heat conductivity and the peridynamics
microscopic liquid velocity at particle $\mathbf{x}$ are defined by [32]:

$\Lambda\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{V_{\mathbf{x}}}$
(35)
$U\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{\overline{u}\left(\mathbf{x^{\prime}},\mathbf{x},t\right)}{V_{\mathbf{x}}},$
(36)

Therefore, Eq. (34) is written as:

$\overline{\left(\rho C\right)_{eq}}(\mathbf{x},t)\frac{\partial
T(\mathbf{x},t)}{\partial
t}=\int_{\mathbf{H_{x}}}\Lambda\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|^{2}}-\\\
\rho_{w}C_{w}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)U\left(\mathbf{x}^{\prime},\mathbf{x},t\right)\frac{T\left(\mathbf{x}^{\prime},t\right)-T\left(\mathbf{x},t\right)}{\left\|\xi\right\|}{\rm
d}V_{\mathbf{x}^{\prime}}$ (37)

Using the approach of uniform constant influence functions proposed in [23,
24], the relationship between PD microscopic heat conductivity
$\Lambda\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ and macroscopic (which
is experimentally measured) heat conductivity $\overline{\lambda}$ for the 1D
and 2D cases can be derived by equating PD and classical local solutions. The
results are

$\Lambda\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\delta}$
(38)
$\Lambda\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{4\overline{\lambda}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\pi\delta^{2}}$
(39)

Similarly, the relationship between PD microscopic water velocity
$U\left(\mathbf{x}^{\prime},\mathbf{x},t\right)$ and the macroscopic velocity
$\overline{u}$ for the 1D and 2D cases are [32]:

$U\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{2\delta}$
(40)
$U\left(\mathbf{x}^{\prime},\mathbf{x},t\right)=\frac{2\overline{u}\left(\mathbf{x}^{\prime},\mathbf{x},t\right)}{\pi\delta^{2}}$
(41)

##  4 Numerical implementation

Each peridynamic particle has an associated volume, which in the cases
considered here is represented by length in 1D and area in 2D. In our
implementation, we consider that all particles have equal volumes. In such
case, the spatial discretisation of the domain is by a uniform grid with step
equal to the particle size, $\Delta x$, which is the length in 1D and the
square root of the area in 2D. The temporal discretisation is implemented by
considering equally spaced time instances, $t^{n}$ ($n\in\mathbb{N}$), with
time interval between the instances $\Delta t=t^{n+1}-t^{n}$.

Since the pressure gradients in the majority of engineering problems discussed
in Section 1 are insignificant, we consider
$\rho_{w}\left(\mathbf{x},t\right)=const$ and discretise Eq. (22). The spatial
discretisation of the conservation of mass in particle at position
$\mathbf{x}_{\alpha}$ ($\alpha\in\mathbb{N}$) is written as

$n\frac{\left(\rho_{w}-\rho_{i}\right)}{\rho_{w}}\frac{\partial
S_{w}\left(\mathbf{x}_{\alpha},t\right)}{\partial
t}=\sum_{\mathbf{x}_{\beta}}K_{w}\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t\right)\frac{\Phi\left(\mathbf{x}_{\beta},t\right)-\Phi\left(\mathbf{x}_{\alpha},t\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}$
(42)

where $\mathbf{x}_{\beta}$ ($\beta\in\mathbb{N}$) are the positions of
particles in the horizon of $\textbf{x}_{\alpha}$, $V_{\alpha\beta}$ is the
portion of the volume associated with $\mathbf{x}_{\beta}$ within the horizon
of $\mathbf{x}_{\alpha}$, and
$K_{w}\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t\right)$ is calculated by
Eq. (15). We consider a horizontal water table, so that
$\Phi(\mathbf{x},t)=p(\mathbf{x},t)$, see Eq. (4), for which Eq. (42) becomes:

$n\frac{\left(\rho_{w}-\rho_{i}\right)}{\rho_{w}}\frac{\partial
S_{w}\left(\mathbf{x}_{\alpha},t\right)}{\partial
t}=\sum_{\mathbf{x}_{\beta}}K_{w}\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t\right)\frac{p\left(\mathbf{x}_{\beta},t\right)-p\left(\mathbf{x}_{\alpha},t\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}$
(43)

The time derivative of the water saturation in Eq. (43) is calculated using
the value from the previous time step, so that the fully discretised
conservation of mass becomes

$\frac{S_{w}\left(\mathbf{x}_{\alpha},t^{n}\right)-S_{w}\left(\mathbf{x}_{\alpha},t^{n-1}\right)}{\Delta
t}=\\\
\sum_{p}\frac{K_{w}\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t^{n}\right)\rho_{w}}{n\left(\rho_{w}-\rho_{i}\right)}\frac{p\left(\mathbf{x}_{\beta},t^{n}\right)-p\left(\mathbf{x}_{\alpha},t^{n}\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}$
(44)

For the known pressure field, that is defined by the solution of a system of
linear equations (44), we can define the vector of water flux
$\mathbf{u}(\mathbf{x},t)$ by the numerical solution of the relation (26). It
is suitable to find the components of this vector separately, for that, the
following system can be solved:

$\mathbf{u}(\mathbf{x},t)=\begin{bmatrix}u_{x}\\\
u_{y}\end{bmatrix}=\begin{bmatrix}\sum_{p}\frac{K_{w}\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t^{n}\right)\rho_{w}}{n\left(\rho_{w}-\rho_{i}\right)}\frac{p\left(\mathbf{x}_{\beta},t^{n}\right)-p\left(\mathbf{x}_{\alpha},t^{n}\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}\cdot
cos(\gamma_{\alpha,\beta})\\\
\sum_{p}\frac{K_{w}\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t^{n}\right)\rho_{w}}{n\left(\rho_{w}-\rho_{i}\right)}\frac{p\left(\mathbf{x}_{\beta},t^{n}\right)-p\left(\mathbf{x}_{\alpha},t^{n}\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}\cdot
sin(\gamma_{\alpha,\beta})\end{bmatrix}$ (45)

where $\gamma_{\alpha,\beta}$ is the angle between the unit vector along a
considered t-bond and the unit vector j that corresponds to horizontal axes
$x$.

After the water flux vectors $\mathbf{u}(\mathbf{x},t)$ are defined for all
particle within the considered domain $\mathbf{R}$, the average water velocity
along every t-bond, that is necessary to consider the conservation of energy,
can be calculated by the relation (29). It finally lets us define the
microscopic liquid velocity according to (40) or (41).

The spatial discretisation of the conservation of energy, Eq. (37), for a
particle at position $\mathbf{x}_{\alpha}$ is

$\overline{\left(\rho C\right)_{eq}}(\mathbf{x}_{\alpha},t)\frac{\partial
T(\mathbf{x}_{\alpha},t)}{\partial
t}=\sum_{\mathbf{x}_{\beta}}\left[\Lambda\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t\right)\frac{T\left(\mathbf{x}_{\beta},t\right)-T\left(\mathbf{x}_{\alpha},t\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}\right.-\\\
\left.\rho_{w}C_{w}U\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t\right)\frac{T\left(\textbf{x}_{\beta},t\right)-T\left(\mathbf{x}_{\alpha},t\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|}V_{\alpha\beta}\right]$
(46)

The time derivative of temperature in Eq. (46) is approximated by the forward
Euler method. Thus, the fully discretised conservation of energy becomes

$\frac{T\left(\mathbf{x}_{\alpha},t^{n+1}\right)-T\left(\mathbf{x}_{\alpha},t^{n}\right)}{\Delta
t}=\\\
\sum_{p}\left[\frac{\Lambda\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t^{n}\right)}{\left(\rho
C\right)_{eq}\left(\mathbf{x}_{\alpha},t^{n}\right)}\frac{T\left(\mathbf{x}_{\beta},t^{n}\right)-T\left(\mathbf{x}_{\alpha},t^{n}\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|^{2}}V_{\alpha\beta}-\right.\\\
\left.\frac{\rho_{w}C_{w}U\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t^{n}\right)}{\left(\rho
C\right)_{eq}\left(\mathbf{x}_{\alpha},t^{n}\right)}\frac{T\left(\mathbf{x}_{\beta},t^{n}\right)-T\left(\mathbf{x}_{\alpha},t^{n}\right)}{\left\|\mathbf{x}_{\beta}-\mathbf{x}_{\alpha}\right\|}V_{\alpha\beta}\right]$
(47)

Equations (44) and (47) represent the numerical implementation of the
developed bond-based PD model. Based on the known
$S_{w}\left(\mathbf{x},t^{n-1}\right)$ and
$S_{w}\left(\mathbf{x},t^{n}\right)$, the solution of Eq. (44) provides the
pressure field at $t^{n}$, that defines the average water flux along every
t-bond (29). Taking it into account, the vectors of water flux at every
particle can be found by the relation (45). Finally, the solution to Eq. (47)
provides the temperature field at $t^{n+1}$. The detailed algorithm of the
proposed numerical implementation is presented in Fig . 2.

Figure 2: The general flow-chart of the developed PD model

According to [32], we should pay special attention when we calculate the terms
for which $\alpha=\beta$. For the 1D dimensional case, it is proposed to
estimate the average value from the closest neighbors [32]:

$\frac{1}{2}\left[\Lambda\left(\textbf{x}_{\alpha},\textbf{x}_{\alpha+1},t\right)\frac{T\left(\textbf{x}_{\alpha+1},t\right)-T\left(\textbf{x}_{\alpha},t\right)}{\left\|\textbf{x}_{\alpha+1}-\textbf{x}_{\alpha}\right\|^{2}}\right.+\\\
\Lambda\left(\textbf{x}_{\alpha},\textbf{x}_{i-1},t\right)\frac{T\left(\textbf{x}_{i-1},t\right)-T\left(\textbf{x}_{\alpha},t\right)}{\left\|\textbf{x}_{i-1}-\textbf{x}_{\alpha}\right\|^{2}}-\\\
\rho_{w}C_{w}\left(\textbf{x}_{\alpha},t\right)U\left(\textbf{x}_{\alpha},\textbf{x}_{i+1},t\right)\frac{T\left(\textbf{x}_{i+1},t\right)-T\left(\textbf{x}_{\alpha},t\right)}{\left\|\textbf{x}_{i+1}-\textbf{x}_{\alpha}\right\|}-\\\
\left.\rho_{w}C_{w}\left(\textbf{x}_{\alpha},t\right)U\left(\textbf{x}_{\alpha},\textbf{x}_{i-1},t\right)\frac{T\left(\textbf{x}_{i-1},t\right)-T\left(\textbf{x}_{\alpha},t\right)}{\left\|\textbf{x}_{i-1}-\textbf{x}_{\alpha}\right\|}\right]$
(48)

For the 2D case, we consider the eight neighbor nodes of the node
$\textbf{x}_{\alpha}$:

$\frac{1}{8}\sum_{\beta=1}^{8}\Lambda\left(\textbf{x}_{\alpha},\textbf{x}_{\beta},t\right)\frac{T\left(\textbf{x}_{\beta},t\right)-T\left(\textbf{x}_{\alpha},t\right)}{\left\|\textbf{x}_{\beta}-\textbf{x}_{\alpha}\right\|^{2}}-\\\
\frac{1}{8}\rho_{w}C_{w}\sum_{\beta=1}^{8}U\left(\textbf{x}_{\alpha},\textbf{x}_{\beta},t\right)\frac{T\left(\textbf{x}_{\beta},t\right)-T\left(\textbf{x}_{\alpha},t\right)}{\left\|\textbf{x}_{\beta}-\textbf{x}_{\alpha}\right\|}$
(49)

##  5 Verification

To assess the accuracy of the implementation, we consider several problems,
where we calculate the temperature distribution at different time instances
and compare the results with alternative solutions obtained using analytical
expressions and/or the FEM. The first problem was designed to verify the
implementation of heat conduction with phase change only; the second problem
was chosen to test the implementation of convective heat transport only; and
the third problem was included to test the implementation of fully coupled
convective-conductive heat transport with phase change under pressure driven
water flow. In addition we analyse a fourth problem, involving high water
velocities, to demonstrate the application of our implementation for
situations that are challenging for numerical methods that are based on local
(differential) formulations.

In all problems, the medium undergoing phase change is water. The thermo-
physical properties of the three phases - solid, water and ice - are presented
in Table 1, where data was taken from [48]. The PD boundary conditions were
applied using a boundary layer with thickness equal to the horizon radius
$\delta$ that was added at the boundary ends in 1D problems, and along the
domain perimeter in 2D problem (details are given in [32]).

Table 1: Thermo-physical properties of phases in test problems Physical properties | Parameter values
---|---
Porosity, n | 0.3
Thermal conductivity of liquid water, $\lambda_{w}$ , (W m-1 ${}^{\circ}{\rm C}$-1) | 0.6
Thermal conductivity of ice, $\lambda_{i}$ , (W m-1 ${}^{\circ}{\rm C}$-1) | 2.14
Thermal conductivity of solid grains, $\lambda_{s}$ , (W m-1 ${}^{\circ}{\rm C}$-1) | 9
Specific heat capacity of liquid water, $C_{w}$ , (J kg-1 ${}^{\circ}{\rm C}$-1) | 4182
Specific heat capacity of ice, $C_{i}$ , (J kg-1 ${}^{\circ}{\rm C}$-1) | 2060
Specific heat capacity of solid grains, $C_{s}$ , (J kg-1 ${}^{\circ}{\rm C}$-1) | 835
Liquid water density, $\rho_{w}$ , (kg m-3) | 1000
Ice density, $\rho_{i}$ , (kg m-3) | 920
Solid grain density, $\rho_{s}$ , (kg m-3) | 2650
Dynamic viscosity of liquid water, $\mu$ , (kg m-1 s-1) |  1.793$\mathrm{\cdot}$10-3
Latent heat of solidification, $L$ , (J kg-1) | 334000
Residual saturation, $S_{w_{res}}$ | 0
$\Delta T$ , ${}^{\circ}{\rm C}$ | 0.5
Water solidification temperature, $T_{f}$ , ${}^{\circ}{\rm C}$ | 0.0
Intrinsic permeability, $k_{int}$ , m2 |  1.3$\mathrm{\cdot}$10-10
$\Omega$ | 50

###  5.1 Heat transfer with phase change: 1D problem

The problem of heat transfer with phase change has an analytical solution in
1D. The problem is formulated for a semi-infinite domain, $0\leq x\leq\infty$,
with initial condition $T(x,0)=T_{\infty}$ and boundary conditions
$T(0,t)=T_{0}$ (constant source temperature) and $T(\infty,t)=T_{\infty}$. The
solution is [1]:

$\begin{array}[]{l}{T\left(x,t\right)=\left\\{\begin{array}[]{c}{T_{0}+\left(T\left(x_{c},t\right)-T_{0}\right)\frac{{\rm
erf}\left(\frac{x}{\left(4\alpha_{i}t\right)^{1/2}}\right)}{{\rm
erf}\left(\frac{\beta}{\left(4\alpha_{i}\right)^{1/2}}\right)},\quad 0\leq
x\leq x_{c}\left(t\right)}\\\
{T_{\infty}+\left(T\left(x_{c},t\right)-T_{\infty}\right)\frac{{\rm
erfc}\left(\frac{x}{\left(4\alpha_{w}t\right)^{1/2}}\right)}{{\rm
erfc}\left(\frac{\beta}{\left(4\alpha_{w}\right)^{1/2}}\right)},\quad
x_{c}\left(t\right)\leq x\leq\infty}\end{array}\right.}\\\ {}\end{array}$ (50)

where $x_{c}\left(t\right)=\beta\sqrt{\alpha_{i}t}$ is the position of the
freezing front at time $t$;
$\alpha_{i}=\lambda_{i}/\left(C_{i}\rho_{w}\right)$ and
$\alpha_{w}=\lambda_{w}/\left(C_{w}\rho_{w}\right)$ are the thermal
conductivities of ice and water, respectively; ${\rm
erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-p^{2}}dp$ is the error function;
${\rm erfc}(x)=1-{\rm erf}(x)$ is the complementary error function; and
$\beta$ is the root of the transcendental equation:

$\frac{e^{-\beta^{2}}}{{\rm
erf}\left(\beta\right)}+\frac{\lambda_{w}}{\lambda_{i}}\left(\frac{\alpha_{i}}{\alpha_{w}}\right)^{1/2}\frac{T\left(x_{c},t\right)-T\left(\infty,t\right)}{T\left(x_{c},t\right)-T\left(0,t\right)}\frac{e^{-\beta^{2}\left(\alpha_{i}/\alpha_{w}\right)}}{{\rm
erfc}\left(\beta\sqrt{\alpha_{i}/\alpha_{w}}\right)}=\\\ \frac{\beta
L\sqrt{\pi}}{C_{i}\left[T\left(x_{c},t\right)-T\left(0,t\right)\right]}$ (51)

We considered the freezing of water and compared our result with the solution
provided by Eq. (50). The simulated domain had finite length of 0.3 m. The
particle size was $\Delta x=0.002$ m, the horizon radius was $\delta=0.006$ m,
and the time step was $\Delta t=2$ s. In the calculations we neglected the
change of density during the phase change, since the analytical solution did
not take it into account. Hence, the density of liquid and solid water was
taken as $\rho_{w}$. The specific values of the initial and boundary
conditions were $T\left(x,0\right)=T\left(\infty,t\right)=8^{\circ}{\rm C}$
and $T\left(0,t\right)=-20{}^{\circ}{\rm C}$. Equation (42) was solved with
$n=1$ (no solid phase) and
$u\left(\mathbf{x}_{\beta},\mathbf{x}_{\alpha},t\right)=0$ (no water flow).
This model simulated a physical time for the process of 10 hours.

It is noted that in this problem we compare the temperature field in an
infinite domain, given by the analytical solution, with the temperature field
in a finite domain, provided by the PD solution. The results are presented in
Fig. 3, where the left graph shows the displacement of the phase front with
time, and the right graph shows the temperature distributions at time
instances 2, 4, 6, 8 and 10 hours. There is very close agreement between the
two solutions, providing further confidence on the implementation of heat
conduction with phase change in the bond-based PD model. It is noted that the
proposed bond-based solution provides results with the same accuracy as the
PDDO model that was developed in [29].

Figure 3: One-dimensional heat transfer with phase change according to the PD
model and analytical results. The upper figure shows the phase front
propagation. The lower figure shows the profiles of temperature distribution
for different times. The solid lines are the developed PD solution; the dots
are the analytical solution.

###  5.2 Heat transfer with water flow: 1D problem

This problem of convective-conduction thermal diffusion also has an analytical
solution in 1D. The problem is again formulated for a semi-infinite domain,
$0\leq x\leq\infty$, with initial condition $T(x,0)=T_{\infty}$ and boundary
conditions $T(0,t)=T_{0}$ (constant source temperature) and
$T(\infty,t)=T_{\infty}$, supplemented by fluid velocity, $v$. The solution is
given by [1]:

$T\left(x,t\right)-T_{\infty}=\frac{T_{0}}{2}\left({\rm
erfc}\left(\frac{x-vt}{2\sqrt{\alpha_{w}t}}\right)+\exp\left(\frac{vt}{\alpha_{w}}\right){\rm
erfc}\left(\frac{x+vt}{2\sqrt{\alpha_{w}t}}\right)\right)$ (52)

We simulated the process in a domain of finite length 0.4 m, with particle
size $\Delta x=0.005$ m, horizon radius $\delta=0.015$ m, and time step
$\Delta t=10$ s. The specific values of the initial and boundary conditions
were $T\left(x,0\right)=T\left(\infty,t\right)=40^{\circ}{\rm C}$ and
$T\left(0,t\right)=5^{\circ}{\rm C}$. Furthermore, the water was assumed to
have constant velocity $v=0.0025$ cm/s.

The results from our simulations and from the analytical solution are
presented in Fig. 4 for time instances 0.5, 1, 1.5 and 2 hours. There is close
agreement between the two alternative solutions which provides further
confidence that the PD formulation for convective-conductive heat transport
has been correctly implemented. These results are also in close agreement with
similar solutions for chemical diffusion-advection problems presented in [32,
34] .

Figure 4: One-dimensional heat transfer with water flow according to the PD
model and analytical results. The figure shows the profiles of temperature
distribution for different time intervals. The solid lines are the developed
PD solution; the dots are the analytical solution.

###  5.3 Heat transfer with water flow and phase change: 2D problem

In this sub-section we assess the performance of the developed model applied
to heat transfer with phase change in a fully saturated porous medium with
pressure driven water flow. Analytical solutions for such problems do not
exist and we compare our PD results, with results obtained using the finite
element method (Comsol Multiphysics [48]).

We adopted a benchmark problem proposed in [48] to estimate the efficacy and
accuracy of different FEM packages. This problem describes the thawing of an
initially frozen inclusion in a porous medium, subject to pressure driven
water flow under constant positive temperature. The domain and its boundary
conditions are illustrated in Fig. 5 and geometric and physical parameters are
presented in Table 2. We note that these values are selected for our work and
differ from the test problem in [48].

Table 2: Parameters of the frozen inclusion model Parameter | Value
---|---
Domain length, $L_{x}$ , m | 0.75
Domain height, $L_{y}$ , m | 0.5
Horizontal position of inclusion centre, $L_{cx}$ , m | 0.28
Vertical position of inclusion centre, $L_{cy}$ , | 0.25
Inclusion length, $L_{i,x}$ | 0.075
Inclusion height, $L_{i,y}$ | 0.075
Initial temperature of unfrozen domain, $T_{domain}$ , ${}^{\circ}{\rm C}$ | 15.0
Initial temperature of frozen inclusion, $T_{inclusion}$ , ${}^{\circ}{\rm C}$ | -5.0
Temperature of infused liquid, $T_{in}$ , ${}^{\circ}{\rm C}$ | 15.0
Imposed pressure (3th / 4th example) $P_{0}$ ,Pa | 5 / 150
Excess pressure (3th / 4th example), $\Delta P$ , Pa | 370 / 11100
Figure 5: Schematic of frozen inclusion problem: geometry, boundary and
initial conditions

In our simulations, the particle size was $\Delta x=0.005$ m, the horizon
radius was $\delta=0.015$ and the time step was $\Delta t=5$ s. In this
problem, we assumed that the t-bonds in the model can be divided into groups:
impermeable solid, liquid and interfacial; for more details, see Section 3. A
t-bond was classified as impermeable solid if it connects nodes with
temperature $T\leq T_{f}-\Delta T$. We supposed that $\Delta T=0.5^{\circ}{\rm
C}$ by analysing the soil freezing characteristic curve (9) and the relation
for relative permeability (5). It showed that by this temperature the relative
permeability of soils $k_{rel}$ became negligible. Therefore, these
impermeable bonds were not used for fluid transfer, and the nodes with such
temperature were not considered in pressure estimation.

The results of PD and FEM simulations are presented in Figs. 6 – 9. Firstly,
Fig. 6 shows the temperature distribution in the domain at several time
instances. This illustrates how the inclusion is changing size and eventually
disappearing, and how the temperature field is deforming due to the water
flow.

Figure 6: Temperature distribution, ${}^{\circ}{\rm C}$, at different moments
of time calculated with the proposed PD model

Comparison between the PD and FEM results is presented in Fig. 7 (a) for the
X–X cross section, (see Fig. 5) at time instances 0 s, 500 s, 1000 s, 1500 s,
2000 s, 2500 s and 3000 s. On this figure, we also present the distribution of
temperature within the inclusion at the initial time instances 0 s, 50 s, 100
s and 250 s. The results indicate a good agreement between the FEM and PD
simulations. Small differences are observed in the temperature distribution
within the inclusion as its temperature increases to $T_{f}$ (results up to
time 500 s). This may be explained by the fact that the release of the latent
heat of solidification is described differently in the two numerical methods.
The temperature distribution for positive temperatures are in very good
agreement. Notably, the agreement can be improved by decreasing the particle
size and the time step.

Figure 7: Temperature distribution along the X–X cross section at different
moments of time, calculated by PD and FEM simulations. The solid lines are the
developed PD solution; the dots are the analytical solution. (a) is the
temperature distribution within the domain; (b) is the temperature
distribution within the inclusion.

Our model allows us to determine the pressure distribution within the domain
at any moment of time. The results are presented in the plots of Fig. 8 for
time instances 0 s, 250 s, 500 s, 1000 s, 2000 s and 3000 s. All ice is melted
by 1300 s leading to steady-state pressure field. These figures demonstrate
clearly how the pressure field is changing with changing the size of frozen
inclusion. The pressure within the inclusion is not calculated, as the water
saturation within it provides very low relative permeability, so the t-bonds
connected to the inclusion’s nodes can be considered as an impermeable solid.

Figure 8: Pressure distribution, Pa, at different moments of time calculated
by the proposed PD model

Comparison between the PD and FEM results for the pressure distribution along
the X–X cross-section for the time instances 0 s, 500 s, 1000 s and 2000 s is
presented in Fig. 9. These plots provide good evidence of the capability of
the PD approach in describing the behaviour of porous material with an
evolving discontinuity, as during the simulated period of time, the number of
bonds that are participating in water transfer are constantly changing.

The presented agreement between two numerical methods indicates that the non-
local PD formulation represents accurately the water flow and heat transfer
with phase change in a saturated porous medium.

Figure 9: Pressure distribution along the X–X cross-section at different
moments of time, calculated by PD and FEM simulations. The solid lines are the
developed PD solution; the dots are the analytical solution.

It should be noted that the developed PD model can provide a solution with the
same high accuracy for the case when all bonds are participating in water
transfer. In this case, the model defines the pressure distribution both in
the liquid region and within the frozen inclusion for every time instance. As
the permeability of the frozen soils is extremely low, the phase change that
coincides with density change, occurs under conditions of nearly constant
volume that, according to the Clausius–Clapeyron equation, leads to the
appearance of extreme negative pressure, see [49, 42]. This phenomenon is
observed in the obtained PD results; however, the FEM model cannot track it at
all, see Fig. (9). The obtained PD pressure values cannot be directly
described by the existing forms of the Clausius–Clapeyron equation. However,
for the saturated soils with an incompressible solid matrix that are
considered in the present study, taking into account this negative pressure
does not have a significant physical meaning and does not affect the accuracy
of the solution. This effect is out the scope of the present paper. However,
in planned future work considering partially saturated soils with deformable
soil matrices, the appearance of extreme pressure values and gradients within
the frozen soils have to be considered.

Such coupled problems are challenging for the FEM. It is widely accepted that
the mesh size and the time step of the model strongly affect the computations
and the results. Some geometric model parameters lead to instability of
computations and divergent solutions. The FEM solution is particularly
sensitive to the flow velocity. This is the reason for our selection of a
relatively low pressure gradient in this example. Higher values may lead to
computational instability and incorrect results. For example, we could not
obtain a converging solution for the present problem for pressures larger than
$\approx 3000$ Pa by FEM. In contrast, the PD implementation can effectively
deal with such conditions, as will be shown in the next sub-section.

###  5.4 Heat transfer with high velocity water flow and phase change

To illustrate the effectiveness of the developed PD model for the case of
convection-dominated heat transfer with phase change, we consider the
following example. We suppose that a frozen body is divided into two parts by
a narrow unfrozen ’crack’. The initial thickness of this crack is 0.09 m and
its shape is given by a $sin$-function with a period 0.5 m and an amplitude
0.15 m. The domain and its boundary conditions are illustrated in Fig. 5. The
values of the initial and boundary conditions are included in Table 2.

Figure 10: Schematic of an unfrozen ’crack’ in a frozen body: geometry,
boundary and initial conditions

The PD simulations were performed with particle size $\Delta x=0.005$ m,
horizon radius $\delta=0.015$, and time step $\Delta t=2$ s. The temperature
distribution in the domain at several time instances is presented in Fig. 11.
The figure illustrates how the crack changes its shape and expands due to the
heat flux from the high velocity water flow. Initially, the temperature within
the ’crack’ decreases from its initial value due to heat exchange with the
frozen part. The constant inflow of warm liquid increases the temperature
within the ’crack’ over time and after 2000 s the crack has changed fully.

Figure 11: Temperature distribution, ${}^{\circ}{\rm C}$, at different moments
of time, according to the proposed PD model

This example illustrates the ability of our bond-based PD implementation to
describe convection-dominated heat transfer in a porous medium with phase
change. It shows clearly how fast inflow of warm water can melt a significant
amount of frozen water in a very short period of time - a situation that we
found to be very challenging using computational methods based on a local
formulation.

##  6 Conclusions

This paper presented a non-local formulation, based on bond-based
peridynamics, for heat transfer with phase change in saturated porous media
under conditions of pressure driven water flow. It allows for the tracking of
the interface between phases and for calculating the temperature and pressure
distributions within the medium. The model and its numerical implementation
were verified by comparison with results obtained using existing analytical
solutions for 1D problems, and with transient finite element solutions for a
2D problem. The agreement found by the verification exercises demonstrates the
correct implementation of the proposed methodology.

The detailed description of coupled heat and hydraulic processes is a critical
step towards a thermo-hydro-mechanical model, which will allow for example, to
describe the hydrological behaviour of permafrost soils and the frost heave
phenomenon that greatly affects civil and mining construction works.

## Acknowledgements

The financial support received by Petr Nikolaev in the form of a President
doctoral scholarship award (PDS award) by the University of Manchester is
gratefully acknowledged. Jivkov and Margetts are grateful for the financial
support of the Engineering and Physical Sciences Research Council, UK (EPSRC)
via grant EP/N026136/1.

## References

  * [1] David W Hahn and M Necati Özisik.  Heat conduction.  John Wiley & Sons, Inc., 3rd edition, 2012.
  * [2] J.C. C. Rowland, B.J. J. Travis, and C.J. J. Wilson.  The role of advective heat transport in talik development beneath lakes and ponds in discontinuous permafrost.  Geophysical Research Letters, 38(17):n/a–n/a, 9 2011.
  * [3] Laurent Orgogozo, Anatoly S. Prokushkin, Oleg S. Pokrovsky, Christophe Grenier, Michel Quintard, Jérôme Viers, and Stéphane Audry.  Water and energy transfer modeling in a permafrost-dominated, forested catchment of central siberia: The key role of rooting depth.  Permafrost and Periglacial Processes, 30(2):75–89, 2019.
  * [4] Hanli Wan, Jianmin Bian, Han Zhang, and Yihan Li.  Assessment of future climate change impacts on water-heat-salt migration in unsaturated frozen soil using CoupModel.  Frontiers of Environmental Science & Engineering, 15(1):10, feb 2021.
  * [5] Xudong Zhang, Yajun Wu, Encheng Zhai, and Peng Ye.  Coupling analysis of the heat-water dynamics and frozen depth in a seasonally frozen zone.  Journal of Hydrology, 593:125603, 2021.
  * [6] J M Frederick and B A Buffett.  Taliks in relict submarine permafrost and methane hydrate deposits: Pathways for gas escape under present and future conditions.  Journal of Geophysical Research: Earth Surface, 119(2):106–122, 2 2014.
  * [7] Patrik Vidstrand, Sven Follin, Jan-Olof Selroos, Jens-Ove Näslund, and Ingvar Rhén.  Modeling of groundwater flow at depth in crystalline rock beneath a moving ice-sheet margin, exemplified by the Fennoscandian Shield, Sweden.  Hydrogeology Journal, 21(1):239–255, 2 2013.
  * [8] Christophe Grenier, Damien Régnier, Emmanuel Mouche, Hakim Benabderrahmane, François Costard, and Philippe Davy.  Impact of permafrost development on groundwater flow patterns: a numerical study considering freezing cycles on a two-dimensional vertical cut through a generic river-plain system.  Hydrogeology Journal, 21(1):257–270, 2 2013.
  * [9] M. Vitel, A. Rouabhi, M. Tijani, and F. Guérin.  Modeling heat and mass transfer during ground freezing subjected to high seepage velocities.  Computers and Geotechnics, 73:1–15, 2016.
  * [10] M. M. Zhou and G. Meschke.  A three-phase thermo-hydro-mechanical finite element model for freezing soils.  International Journal for Numerical and Analytical Methods in Geomechanics, 37(18):3173–3193, 2013.
  * [11] H. Tounsi, A. Rouabhi, and E. Jahangir.  Thermo-hydro-mechanical modeling of artificial ground freezing taking into account the salinity of the saturating fluid.  Computers and Geotechnics, 119(September 2019):103382, 3 2020.
  * [12] Ahmad Zueter, Aurelien Nie-Rouquette, Mahmoud A. Alzoubi, and Agus P. Sasmito.  Thermal and hydraulic analysis of selective artificial ground freezing using air insulation: Experiment and modeling.  Computers and Geotechnics, 120(July 2019):103416, 4 2020.
  * [13] Erdogan Madenci and Erkan Oterkus.  Peridynamic Theory and Its Applications, volume 9781461484.  Springer New York, New York, NY, 2014.
  * [14] S.A. Silling.  Reformulation of elasticity theory for discontinuities and long-range forces.  Journal of the Mechanics and Physics of Solids, 48(1):175–209, 1 2000.
  * [15] Stewart A Silling, M Epton, O Weckner, Ji Xu, and E23481501120 Askari.  Peridynamic states and constitutive modeling.  Journal of Elasticity, 88(2):151–184, 2007.
  * [16] Florin Bobaru and Monchai Duangpanya.  A peridynamic formulation for transient heat conduction in bodies with evolving discontinuities.  Journal of Computational Physics, 231(7):2764–2785, 2012.
  * [17] Ali Javili, Rico Morasata, Erkan Oterkus, and Selda Oterkus.  Peridynamics review.  Mathematics and Mechanics of Solids, 24(11):3714–3739, 2019.
  * [18] Abigail Agwai, Ibrahim Guven, and Erdogan Madenci.  Predicting crack propagation with peridynamics: a comparative study.  International Journal of Fracture, 171(1):65–78, sep 2011.
  * [19] Shuo Liu, Guodong Fang, Jun Liang, Maoqing Fu, and Bing Wang.  A new type of peridynamics: Element-based peridynamics.  Computer Methods in Applied Mechanics and Engineering, 366:113098, 7 2020.
  * [20] Erdogan Madenci, Atila Barut, and Michael Futch.  Peridynamic differential operator and its applications.  Computer Methods in Applied Mechanics and Engineering, 304:408–451, jun 2016.
  * [21] Erdogan Madenci, Atila Barut, and Mehmet Dorduncu.  Peridynamic Differential Operator for Numerical Analysis.  Springer International Publishing, Cham, 2019.
  * [22] Walter Gerstle, Stewart Silling, David Read, Vinod Tewary, and Richard Lehoucq.  Peridynamic simulation of electromigration.  Comput Mater Continua, 8(2):75–92, 2008.
  * [23] Florin Bobaru and Monchai Duangpanya.  The peridynamic formulation for transient heat conduction.  International Journal of Heat and Mass Transfer, 53(19-20):4047–4059, 2010.
  * [24] Ziguang Chen and Florin Bobaru.  Selecting the kernel in a peridynamic formulation: A study for transient heat diffusion.  Computer Physics Communications, 197:51–60, 2015.
  * [25] Xin Gu, Qing Zhang, and Erdogan Madenci.  Refined bond-based peridynamics for thermal diffusion.  Engineering Computations, 36(8):2557–2587, oct 2019.
  * [26] L. J. Wang, J. F. Xu, and J. X. Wang.  The green’s functions for peridynamic non-local diffusion.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2193):20160185, 2016.
  * [27] Akbar Jafari, Reza Bahaaddini, and Hasan Jahanbakhsh.  Numerical analysis of peridynamic and classical models in transient heat transfer, employing galerkin approach.  Heat Transfer—Asian Research, 47(3):531–555, 2018.
  * [28] Linjuan Wang, Jifeng Xu, and Jianxiang Wang.  A peridynamic framework and simulation of non-Fourier and nonlocal heat conduction.  International Journal of Heat and Mass Transfer, 118:1284–1292, mar 2018.
  * [29] Erdogan Madenci, Mehmet Dorduncu, Atila Barut, and Michael Futch.  Numerical solution of linear and nonlinear partial differential equations using the peridynamic differential operator.  Numerical Methods for Partial Differential Equations, 33(5):1726–1753, sep 2017.
  * [30] Amit Katiyar, John T. Foster, Hisanao Ouchi, and Mukul M. Sharma.  A peridynamic formulation of pressure driven convective fluid transport in porous media.  Journal of Computational Physics, 261:209–229, 2014.
  * [31] Rami Jabakhanji and Rabi H. Mohtar.  A peridynamic model of flow in porous media.  Advances in Water Resources, 78(January):22–35, 4 2015.
  * [32] Jiangming Zhao, Ziguang Chen, Javad Mehrmashhadi, and Florin Bobaru.  Construction of a peridynamic model for transient advection-diffusion problems.  International Journal of Heat and Mass Transfer, 126:1253–1266, 11 2018.
  * [33] Majid Sedighi, Huaxiang Yan, and Andrey P Jivkov.  Peridynamics modelling of clay erosion.  Géotechnique, 0(0):1–12, 12 2020.
  * [34] Huaxiang Yan, Majid Sedighi, and Andrey P. Jivkov.  Peridynamics modelling of coupled water flow and chemical transport in unsaturated porous media.  Journal of Hydrology, 591(October):125648, 2020.
  * [35] Wanhan Chen, Xin Gu, Qing Zhang, and Xiaozhou Xia.  A refined thermo-mechanical fully coupled peridynamics with application to concrete cracking.  Engineering Fracture Mechanics, 242:107463, 2021.
  * [36] Yanzhuo Xue, Renwei Liu, Zheng Li, and Duanfeng Han.  A review for numerical simulation methods of ship–ice interaction.  Ocean Engineering, 215:107853, 2020.
  * [37] Wei Lu, Mingyang Li, Bozo Vazic, Selda Oterkus, Erkan Oterkus, and Qing Wang.  Peridynamic Modelling of Fracture in Polycrystalline Ice.  Journal of Mechanics, 36(2):223–234, 2020.
  * [38] Renwei Liu, Yanzhuo Xue, Duanfeng Han, and Baoyu Ni.  Studies on model-scale ice using micro-potential-based peridynamics.  Ocean Engineering, 221:108504, 2021.
  * [39] Jiangming Zhao, Ziguang Chen, Javad Mehrmashhadi, and Florin Bobaru.  A stochastic multiscale peridynamic model for corrosion-induced fracture in reinforced concrete.  Engineering Fracture Mechanics, 229:106969, 2020.
  * [40] Ziguang Chen, Sina Niazi, and Florin Bobaru.  A peridynamic model for brittle damage and fracture in porous materials.  International Journal of Rock Mechanics and Mining Sciences, 122:104059, 2019.
  * [41] Ziguang Chen, Siavash Jafarzadeh, Jiangming Zhao, and Florin Bobaru.  A coupled mechano-chemical peridynamic model for pit-to-crack transition in stress-corrosion cracking.  Journal of the Mechanics and Physics of Solids, 146:104203, 2021\.
  * [42] Barret L. Kurylyk and Kunio Watanabe.  The mathematical representation of freezing and thawing processes in variably-saturated, non-deformable soils.  Advances in Water Resources, 60:160–177, Oct 2013.
  * [43] Jidong Teng, Han Yan, Sihao Liang, Sheng Zhang, and Daichao Sheng.  Generalising the kozeny-carman equation to frozen soils.  Journal of Hydrology, 594:125885, 2021.
  * [44] Jeffrey M. McKenzie, Clifford I. Voss, and Donald I. Siegel.  Groundwater flow with energy transport and water–ice phase change: Numerical simulations, benchmarks, and application to freezing in peat bogs.  Advances in Water Resources, 30(4):966–983, 4 2007.
  * [45] Zheng Yu, Longcai Yang, and Shunhua Zhou.  Comparative study of relating equations in coupled thermal-hydraulic finite element analyses.  Cold Regions Science and Technology, 161:150–158, 2019.
  * [46] Haiyan Wen, Jun Bi, and Ding Guo.  Evaluation of the calculated unfrozen water contents determined by different measured subzero temperature ranges.  Cold Regions Science and Technology, 170:102927, 2020.
  * [47] Selda Oterkus, Erdogan Madenci, and Erkan Oterkus.  Fully coupled poroelastic peridynamic formulation for fluid-filled fractures.  Engineering Geology, 225:19–28, 7 2017.
  * [48] Christophe Grenier, Hauke Anbergen, Victor Bense, Quentin Chanzy, Ethan Coon, Nathaniel Collier, François Costard, Michel Ferry, Andrew Frampton, Jennifer Frederick, Julio Gonçalvès, Johann Holmén, Anne Jost, Samuel Kokh, Barret Kurylyk, Jeffrey McKenzie, John Molson, Emmanuel Mouche, Laurent Orgogozo, Romain Pannetier, Agnès Rivière, Nicolas Roux, Wolfram Rühaak, Johanna Scheidegger, Jan Olof Selroos, René Therrien, Patrik Vidstrand, and Clifford Voss.  Groundwater flow and heat transport for systems undergoing freeze-thaw: Intercomparison of numerical simulators for 2D test cases.  Advances in Water Resources, 114(January):196–218, 2018.
  * [49] M. Dall’Amico, S. Endrizzi, S. Gruber, and R. Rigon.  A robust and energy-conserving model of freezing variably-saturated soil.  The Cryosphere, 5(2):469–484, Jun 2011.