Abstract As climate change intensifies, soil water flow, heat transfer, and solute transport in the active, unfrozen zones within permafrost and seasonally frozen ground exhibit progressively more complex interactions that are difficult to elucidate with measurements alone. For example, frozen conditions impede water flow and solute transport in soil, while heat and mass transfer are significantly affected by high thermal inertia generated from water‐ice phase change during the freeze‐thaw cycle. To assist in understanding these subsurface processes, the current study presents a coupled two‐dimensional model, which examines heat conduction‐convection with water‐ice phase change, soil water (liquid water and vapour) and groundwater flow, ad v ective‐dispersive solute transport with sorption, and soil deformation (frost heave and thaw settlement) in variably saturated soils subjected to freeze‐thaw actions. This coupled multiphysics problem is numerically solved using the finite element method. The model's performance is first verified by comparison to a well‐documented freezing test on unsaturated soil in a laboratory environment obtained from the literature. Then based on the proposed model, we quantify the impacts of freeze‐thaw cycles on the distribution of temperature, water content, displacement history, and solute concentration in three distinct soil types, including sand, silt and clay textures. The influence of fluctuations in the air temperature, groundwater level, hydraulic conductivity, and solute transport parameters was also comparatively studied.The results show that (i) there is a significant bidirectional exchange between groundwater in the saturated zone and soil water in the vadose zone during freeze‐thaw periods, and its magnitude increases with the combined influence of higher hydraulic conductivity and higher capillarity; (ii) the rapid dewatering ahead of the freezing front causes local volume shrinkage within the non‐frozen region when the freezing front propagates downward during the freezing stage and this volume shrinkage reduces the impact of frost heave due to ice formation. This gradually recovers when the thawed water replenishes the water loss zone during the thawing stage; and (iii) the profiles of soil moisture, temperature, displacement, and solute concentration during freeze‐thaw cycles are sensitive to the changes in amplitude and freeze‐thaw period of the sinusoidal varying air temperature near the ground surface, hydraulic conductivity of soil texture, and the initial groundwater levels. Our modelling framework and simulation results highlight the need to account for coupled thermal‐hydraulic‐mechanical‐chemical behaviours to better understand soil water and groundwater dynamics during freeze‐thaw cycles and further help explain the observed changes in water cycles and landscape evolution in cold regions. This article is protected by copyright. All rights reserved.

Subsurface pipe failures in cold regions are generally believed to be exacerbated by differential strain in shallow soils induced by seasonal freeze and thaw cycles. The transient stress–strain fields resulting from soil water phase change may influence the occurrence of local buried pipe breaks including those related to urban water mains. This work proposes that freezing-induced frost loading results in uneven stress–strain distributions along the buried water mains placing them at risk of bending, breaking, and/or leaking. A coupled thermal-hydraulic-mechanical (THM) model was developed to illustrate the interactions among moisture, temperature, and stress–strain fields within variably saturated freezing soils. Several typical cases involving highly frost-susceptible and lower frost-susceptible soils underlying roadbeds were examined. Results show that the magnitude of frost-induced compressive stress and strain changes between different frost-susceptible soils can vary significantly. Such substantial differences in stress–strain fields would increase the breakage risk of water mains buried within different types of soils. Furthermore, even water mains buried within soils with low frost-susceptibility are at risk when additional sources of soil water exist and are available to migrate to the freezing front. To reduce the risk of damage to buried pipe-like infrastructure, such as municipal water mains, from soil freezing phenomena, the selected backfill material should have fairly consistent frost susceptibility or a broad zone of transition should be considered between materials with significantly different frost susceptibility. In addition, buried pipes should be kept as far away from external sources of subsurface water as possible considering the potential for the water source to exacerbate the level of risk to the pipe.

• A novel analytical-numerical scheme for calculating temperature profiles in porous media with temperature-dependent thermal properties during the freezing process; • The hybrid analytical-numerical method can deal with different types of nonlinear soil freezing functions; • Neumann's two-layer solution underestimates the penetration rate and depth of the freezing front; • The profiles of temperature, equivalent thermal conductivity and diffusivity, conductive heat flux, and dynamics of the freezing front were significantly impacted by the shape of the unfrozen water content curve and the magnitude of soil grain thermal conductivity. The freeze-thaw cycle associated with climatic seasonality is a common phenomenon in cold regions affecting a wide range of subsurface processes. Due to the complex and highly nonlinear nature of the associated hydrologic processes, transient freeze-thaw dynamics are conventionally quantified in a numerical way. Here we present a hybrid analytical-numerical scheme for solving one-dimensional soil (or porous media) temperature profiles when the soil profile is subjected to unidirectional freezing (or thawing) conditions. This scheme divides the partially-frozen soil into multi-layers, each with constant thermal parameters and fixed-temperature boundaries. Temperature profiles within each layer were obtained by solving multiple moving-boundary problems. The proposed hybrid analytical-numerical scheme was tested into a freezing test of silty clay in a permafrost region on the Qinghai-Tibetan Plateau, and its solution was in good agreement with the finite element numerical solution. Results show that the proposed multi-layer method adapted well to the changes in unfrozen water content and thermal properties of soil over a wide range of subzero temperatures. By contrast, the freezing front's migration rate and penetration depth calculated by Neumann's classical solution, which only considers two zones (frozen and unfrozen), was found to be underestimated. As for our proposed multi-layer solution, by dividing the subsurface domain into many layers with smaller proportion ratios (thinner layers close to the freezing front), there was a slower penetration rate of the freezing front resulting in shallower penetration depth. The predicted profiles of temperature, thermal conductivity and diffusivity, heat flux, and dynamics of the freezing front were significantly impacted by the shape of the soil freezing curves and the magnitude of soil grain thermal conductivity, especially for the accuracy of long-term predictions.

• Coupled modelling of water flow, heat transfer, water-ice phase change and ice lens formation in deformable, variably saturated freezing soils. • Moisture, vapour, temperature and stress-strain fields significantly interact with each other and should be fully accounted for within the modeling platform. • The large increases in effective stress ahead of the freezing front causes substantial compaction in the unfrozen zone. Although many frost heave and freezing soil models have been developed in the past decades, saturated conditions are commonly assumed and/or the behavior of pore ice rather than ice lenses are conventionally predicted. This study presents a fully coupled thermal-hydraulic-mechanical (THM) model for variably saturated freezing soil, which examines a number of processes. These include heat conduction and convection, phase change, water (moisture) movement through cryosuction, and the development of independent ice lenses. Instead of directly solving for the pore pressure distributions, the void ratio is considered as a dependent variable related to the degree of water saturation. Both the stress-deformation and ice lens segregation are inextricably linked to the evolution of the void ratio as well. The coupled mechanism and performance of the model is first verified by comparison with laboratory freezing experiment observations obtained from literature and then is further evaluated by a series of parametric analyses. The results show that the calculated profiles of temperature, water content and frost heave are in good agreement with literature experimental data, demonstrating that the proposed THM coupling model appropriately represents the mechanisms of heat-moisture-deformation in variably saturated freezing soil. In addition, the sensitivity analysis illustrates that in the test cases considered, thermally-induced cryosuction due to phase change is the main driving force for water migrating towards the freezing front. Also, ahead of the freezing front, a significant increase in effective stress developed due to the elevated negative pore pressure and expansion of ice lenses causing substantial consolidation and reduction in porosity in the unfrozen zone. As the freezing front penetrated with time, the temperature, moisture, vapour and stress-strain fields interact with each other. The distribution of water vapour was mainly controlled by the temperature gradient and location of the freezing front. Both the initial degree of saturation and hydraulic conductivity affected the distribution of pore pressure and displacements. Higher compression moduli and lower overburden load led to greater frost heave but exerted little influence on the temperature field. Finally, the two-sided freezing scenario for soils underlain by permafrost made the middle ice-poor zone highly compacted with ice lenses accumulating near both freezing boundaries.