Understanding of the thermal, hydraulic and mechanical behavior of soils under freezing and thawing is important in solving problems of permafrost engineering and artificial ground freezing. Accurate mathematical descriptions of heat and moisture transport and the associated mechanical deformations is critical for developing such understanding. Phase change and transition during freezing and thawing are characterized by a strong non-linearity due to the rapid variations of thermo-physical properties of the material and the release of the latent heat of solidification. Such processes cannot be easily and accurately modelled via the use of local (differential) formulations [1,2]. A non-local formulation, developed within the framework of peridynamics [3–5] is proposed here for ground freezing and thawing. It expands the domain of applicability of the first work on non-local modelling of heat transfer with water flow in saturated porous media [6] into modelling the thermo-hydro-mechanical behavior of unsaturated soils.
 The mathematical description of unsaturated soil behavior is based on three governing equations of conservations of water mass, energy and linear momentum, which are complemented by constitutive relationships that represent the change of soil properties.
 Using the local mathematical model of unsaturated soils [7,8], a peridynamics mathematical framework is formulated. Following [6], the bond-based peridynamics considers a soil body as a collection of particles, see Fig. 1 for the 1D case. The term `bond' refers to the interactions between two points, located at depths and . The point interacts with all other points within a certain finite region , called the horizon of z.
 The peridynamic form of the water mass conservation is:
 where is the water pressure head; is time; is the volumetric content of a soil component; is the density of a soil component; and is the vertical coordinate. The micro-diffusivity function can be defined following [9]. Here, the subscripts i, l, a and s refer to ice, liquid water, air and soil particles respectively.
 The peridynamic form of the energy conservation is:
 where is temperature; is the latent heat of water solidification; is the average thermal conductivity of soil; is the vector of water flux; is the average heat capacity of unsaturated soils and is the specific heat capacity of liquid water. Here, is a micro-conductivity function and is a micro-velocity function, for more details see [6,10].
 The peridynamic form of the linear momentum conservation is:
 where is the acceleration of the peridynamics particle ; is a pairwise force density vector between particles and that depends on its initial and deformed positions Z and Z’; is the horizon volume associated with the particle ; is the body force on the particle . The elastic stress-strain relation for unsaturated soils can be stated as:
 where is the stress increment tensor; is soil compressive modulus for unsaturated conditions, and is the strain increment tensor due to the elastic deformation of the soil matrix, that can be find from the relation for the total strain increment tensor of unsaturated soils . Here is the strain increment tensor corresponding to the change of soil volume due to the frost heave; and is the strain increment tensor corresponding to the plastic deformations of soils.
 The numerical implementation of the peridynamic formulation is verified by comparing the results of simulations with analytical solutions and the formulation is validated by comparing these results with experimental data under the condition of one-dimensional heat and water transfer obtained by laboratory experiments. In the experiments, soil cylinders were cooled from one surface by a source of negative temperature. As a result, the freezing front propagated along the cylinder. The temperature profile in the column was changing, and water migration was observed. The developed model is shown to predict the results of these experiments well.
 The proposed peridynamic model is the first to combine the flow of liquid water with heat transfer and mechanical deformations in the presence of a change of the water phase state in unsaturated soils. It is demonstrated that for the soils with relatively high initial water content, the developed approach describes effectively the change of liquid water content and predict soil deformation. The set of integral-differential equations are solved by a simple numerical scheme.