Conduction, radiation, and convection are the three main modes of heat transfer. Heat conduction is the main mode of heat transfer in salinized frozen soil in nature, and the effects of convection and radiation can be neglected (Luo et al., 2023; He et al., 2021). Assuming that the salinized frozen soil is an isotropic continuous medium, unfrozen water, crystalline salt, soil particles, and ice crystals are incompressible. Ignoring convective heat due to moisture seepage, according to Fourier’s law of heat conduction and the law of energy conservation, the crystallization of moisture salts will release heat (Wang et al., 2016; Zhang et al., 2020). Treating the phase changes of water and salt crystallization as heat sources, the differential equation of heat conduction in salinized frozen soil is Equation 1 (Deng et al., 2021): ρ C θ ∂ T ∂ t = λ θ ∇ 2 T + L c ∂ θ c ∂ t + L i ρ i ∂ θ i ∂ t , (1) where ∇ is the Hamiltonian operator for two-dimensional ∂ / ∂ x , ∂ / ∂ y ; T is the temperature of the soil; t is the time; L c and L i are the latent heat of phase change of crystalline salt and ice; ρ c and ρ i are the densities of the crystalline salt and ice, respectively; θ c and θ i are the volume contents of the crystalline salt and the pore ice, respectively; λ is the thermal conductivity unit of W/(m-K), which is expressed by the following Equations 2, 3: λ θ = λ s θ s λ i θ i λ c θ c λ u θ u , (2) where λ s is the thermal conductivity of soil; λ i is the thermal conductivity of ice; λ c is the thermal conductivity of crystalline salt; λ u is the thermal conductivity of water; θ s , θ i , θ c , θ u , are the proportions of the corresponding substances; C is the specific heat capacity, which can be expressed by the sum of the components of the volume-specific heat capacity: ρ C θ = ρ s C s θ s + ρ i C i θ i + ρ c C c θ c + ρ w C u θ u , (3) where C s is the specific heat capacity of the earth; C i is the specific heat capacity of ice; C c is the specific heat capacity of crystalline salt; C u is the specific heat capacity of water.

Unfrozen water always exists in the frozen soil. The migration of water follows Darcy’s law (Ming et al., 2022; Andersland et al., 1996). According to Richard’s Equation 4 (Chai et al., 2018; Burt and Williams, 1976), the total water content of saline soil comprises three components: ice, unfrozen water, and mannite crystalline water: θ T = θ u + ρ i ρ w θ i + ρ c M 180 ρ w M c θ c , (4) where θ T the total water content, M c is the molecular mass of the crystalline salt; M 180 is the molecular mass of the bound water in the crystalline salt.

The ratio of ice to unfrozen water volume content as a function of temperature with temperature reduction is as follows (Bai et al., 2015; Watanabe and Osada, 2016) (Equation 5): B i = θ i θ u 1.1 T T i B − 1 T < T i 0   T ≥ T i , (5) where B is a constant related to the type of soil and salt content. The values of 0.61, 0.47, and 0.56 are taken for sandy, pulverized, and clay soils, respectively, and T i is the freezing temperature.

The differential equation for water infiltration in salted frozen soil is Equation 6 ∂ θ u ∂ t + ρ i ρ w ∂ θ i ∂ t + ρ c M 180 ρ w M c ∂ θ c ∂ t = ∇ D θ u ∇ θ + k g θ u , (6) where θ u is the volume content of unfrozen water in the frozen soil; k g is the permeability coefficient of the unsaturated soil in the direction of the gravitational acceleration; D θ u is the moisture diffusion coefficient of the soil, which can be calculated by the following Equations 7–9 (Chen and Zhang, 2020): D θ u = k u i C θ u , (7) where C θ u is the specific water capacity (1/m); k u i is the permeability coefficient of salinized frozen soil. This paper adopts the method of calculating the specific water capacity from the Van Genuchten (VG) model: S e = θ − θ r θ s − θ r = 1 1 + α h m n m , (8) C θ u = a 0 m / 1 − m S e 1 m 1 − S e 1 m m , (9) where S e is the relative saturation of frozen soil, θ s is the saturated water content of the soil, and θ r is the residual water content of the soil.

The permeability coefficient of salinized frozen soil k u i is closely related to the pore diameter in the soil, with the pores in the soil being simplified as capillary tubes with different diameters.

As the temperature decreases, the solubility of salt in solution decreases, leading to the crystallization and precipitation of salt. The relationship between the volume fraction of crystallized salt and the volume fraction of unfrozen water, as well as the solution concentration, is established (Huang et al., 2008) as follows (Equation 10): θ c = C − n T θ u M c ρ c   C ≥ n T 0   C < n T , (10) where n(T) is the solubility curve of easily soluble salts in the pores of the soil, which is a function of temperature.

If it is assumed that the salt solution is contained in a capillary tube with a radius of Ri (i = 1 − n), the solution is always saturated after the salt precipitation. If the supersaturation of the solution is not taken into account, the following Equation 11 is obtained: C ′ = n T   C > n T   C   C ≤ n T , (11) where C is the initial mass concentration of solution, C′ is the subsequent mass concentration of solution after cooling, Ri is the initial pore radius, and R_i′ is the subsequent pore radius after the salt precipitation due to cooling. Assuming that the salt crystals are evenly deposited on the pore walls after salt precipitation, the quantity of salt crystallization equals the amount of variation in the pore (Equation 12): M c M ρ c π R i u 2 C 100 − C − π R i 2 C ′ 100 − C ′ = π R i u 2 − π R i 2 , (12) where M c is the molecular weight of salt after crystallization, M is the molecular weight of salt, and ρ c is the specific gravity of the crystalline salt. For the same salt, M c M ρ c is a constant value. Let M c M ρ c = α , then Equations 12, 13 can be used to determine the pore radius R i after salt precipitation (Bai et al., 2015): R i = R i u α C − 100 − C 100 − C ′ α C ′ − 100 − C ′ 100 − C . (13)

The following equation can be used to determine the connection between temperature T and pore radius R i (Peng et al., 2021; Wang et al., 2019) as Equation 14: R i = σ l s ρ i L f T i T i − T , (14) where Lf is the latent heat of ice melting, Lf = 3.34 × 105 J/Kg; Ti is the freezing temperature of saline soil; T is the temperature of saline soil; σ l s is the ice water interface energy, where σ _ls = 2.90×10⁻² J/m².

Using the unfrozen water characteristic curve, assuming the temperature range (T0, Tn) is divided into n equal parts, for every decrease of △θu in the unfrozen water content, the corresponding k increases by 1. When k = 0, ks represents the permeability coefficient of saturated thawed soil under the same conditions as Equation 15. k s = π ξ 8 μ ∞ ∑ i = 1 M n i R i 4 , (15) where μ∞ is the viscosity of volumetric water, which is a function of temperature (Thomas and Sansom, 1995) and can be represented as μ ∞ T = 0.6612 T − 229 − 1.562 ; ξ is the ratio of the soil column length H to the actual path of water seepage L( ξ =L/H), ξ i = 1 + 0.41 ⁡ ln 1 ∕ θ i (Ishizaki et al., 1996); n i is the number of pores with the same diameter; R i is the diameter of unfrozen pores.

Therefore, the permeability coefficient of saturated salinized frozen soil is expressed as (Zhang and Guan, 2022) Equation 16 k i = k s μ ∞ μ d ∑ i = k + 1 M n i R i 4 + 4 R i 3 l s , t / ∑ i = 1 M n i R i 4 , (16) where μ d is the effective viscosity. The unfrozen water content in frozen soil equals the sum of liquid water in each pore, which can be expressed as Equation 17 θ u = π ∑ i = k + 1 M n i R i 2 . (17)

The number of pores with diameter R i , denoted as n i , can be calculated from the unfrozen water characteristic curve. k u is the permeability (m/s) of unsaturated soil. In this paper, the calculation method of the permeability coefficient of unsaturated soil is adopted by the VG model in Equation 18: k u = k s S e 0.5 1 − 1 − S e 1 m m 2 = k s A , (18) where m and n are the VG model fitting parameters, which vary for different types of soils.

Therefore, the permeability coefficient of unsaturated salinized frozen soil is expressed as Equation 19 k u i = k u A μ ∞ μ d ∑ i = k + 1 M n i R i 4 + 4 R i 3 l s , t / ∑ i = 1 M n i R i 4 . (19)

Convection, molecular diffusion, and mechanical dispersion are the main modes of salt transport in solutes. Amongst these modes, molecular diffusion and mechanical dispersion are collectively known as hydrodynamic dispersion (Jian et al., 2021; Ma et al., 2016).

Convection, which is the migration of salts along with water in the soil, has a significant effect on soil solute fluxes. The solute flux due to convection is related to the soil water flux and the concentration of the solution and can be described by the following Equation 20 (Ishizaki et al., 1996): J c = q l C , (20) where J c refers to the flux of solute due to convection through a unit area of soil in a unit of time, q l is the flux of liquid water in the soil, and C is the solute concentration.

Molecular diffusion is the process of solute transport by ions or molecules due to Brownian motion (Wu et al., 2017; Xiao et al., 2015). This process is irreversible. It makes a solution in the presence of a concentration gradient, and the solution tends to be homogeneous. Molecular diffusion in soil pore solutions can also be described by Fick’s first law as Equation 21: J d = − D 0 a e b θ u ∂ C ∂ z , (21) where J d is the molecular diffusion flux of the solute, ∂ C / ∂ z is the concentration gradient, D 0 is the molecular diffusion coefficient of solute in free water, e is the soil pore ratio, a and b are empirical parameters, a is generally taken between 0.005 and 0.001, and b is generally taken as 10.

Mechanical dispersion is due to the non-homogeneity of the soil. Therefore, the water flow in the soil pores has different magnitudes and directions, and the solute spreads to other ranges. The solute flux due to mechanical dispersion can be approximated by Fick’s law as Equation 22 (Deng, 2006): J h = − α v ∂ C ∂ z , (22) where J h is the solute flux due to mechanical dispersion, α is an empirical coefficient and is related to the nature of the soil itself, taking the value of 0.2–0.55 cm; v = q / θ , q can be obtained from Darcy’s law, q = D d θ / d y .

The solute hydrodynamic dispersion effect can be described by the following Equation 23: J D = J d + J h = − D 0 a e b θ u + α v ∂ C ∂ z , (23) where the total salt content comprises the readily soluble salt content in unfrozen water, the crystalline salt precipitated by ice crystallization, and the crystalline salt after solubility reduction as Equation 24: θ w c = θ u c + ρ i ρ w θ i c + ρ c M c θ c . (24)

Therefore, the salt field control (Equation 25) is (Liu et al., 2021; Zhao and Luo, 2019) ∂ θ u c + ρ i ρ w θ i c + ρ c M c θ c ∂ t = ∇ D 0 a e b θ u + α v ∂ c ∂ z + q c . (25)

The results of the model calculations of sulfate and chloride saline soils under open conditions were compared with the experimental results to verify the validity and feasibility of the above numerical coupled hydrothermal salt model. During the temperature drop process of the soil sample, the unfrozen water content is determined by Formula 17, the freezing temperature T_i of the salinized soil, and the crystallization of salt is obtained from the water–salt (Na₂SO₄, NaCl) binary phase diagram. As shown in Figure 1, the freezing temperature T_i curves of saline soil are AB(Na₂SO₄) and AB′(NaCl), and the solubility curves of salt n(T) are BC(Na₂SO₄) and BC′(NaCl). These two curves intersect at the eutectic point B (Niu and Cheng, 2002). When the soil temperature drops to T_X, the corresponding amount of unfrozen water content, the freezing temperature, and salt crystallization are calculated from the relationship.

A numerical simulation of water–salt (Na₂SO₄, NaCl migration in a one-way freezing test with brine recharge is conducted to analyze the water–salt migration and salt crystallization law of saline soil. Through the relevant model parameters, the model is set to be solved with the same initial and boundary conditions as the indoor test, and the numerical simulation results are compared with the measured data of the test to verify model validity.

Both salt crystallization and ice occur during the cooling process of Na₂SO₄ saline soil. However, in the calculation of NaCl saline soil, the concentrations of NaCl solution are preset to be 0.4 mol/kg, 0.6 mol/kg, 0.8 mol/kg, and 1.0 mol/kg, with NaCl contents in the soil of 0.412%, 0.618%, 0.824%, and 1.030%, respectively. The temperature ranges from −8°C to 4°C. As indicated by Figure 1, the NaCl-H₂O phase diagram, during the soil cooling process, NaCl does not crystallize; the NaCl content mainly affects the freezing temperature and the unfrozen water content of the soil. Therefore, neglecting NaCl crystallization, the influence of NaCl content on the soil’s thermal conductivity, latent heat, and the pore diameter of phase transition is disregarded.

These soil parameters (Xu et al., 2021; Zhang et al., 2021b; Tang et al., 2023; Wan et al., 2015) are used for the model calculations (Table 1).

The model assumes the following:

1. The pores are not deformed, and the porosity remains constant;

COMSOL was used to simulate the 2D frozen soil columns using a map illustration of a sub-grid with 1,200 cells, 148 boundary cells, and four vertex cells (Figures 2, 3). The model calculation time is 96 h, the calculation result output step is 1, and the solver step is set to intermediate.

To verify the accuracy of the model, it was configured with the same initial conditions and boundary conditions as the test (Zhang, 2020). Under water replenishment conditions, numerical simulations of unidirectional freezing of Na₂SO₄ saline soil were conducted, and the computed results were compared with test data. The variables of the equation, the boundary conditions, and the initial conditions set are shown in Table 2.

The distribution of salt crystallization amount with different initial Na₂SO₄ concentrations is shown in Figure 10. According to the theoretical calculation model of freezing temperature, the freezing temperature of Na₂SO₄ soil comprises two parts: moisture activity and pore water freezing temperature. When the pore water freezing temperature of the sample is −0.76°C, the freezing temperature of the soil is −2.35°C, −2.5°C, −2.01°C and −2.22°C; at this time, the maximum crystalline amount of the salts is 0.263, 0.428 and 0.483. The maximum crystallization amount of salt is 0.263, 0.428, 0.483, and 0.6, and the salt crystallization amount is large under high initial salt content.

Due to the high initial temperature and large initial salt concentration, the salt crystallized in the upper and lower parts of the soil, and the upper part crystallized more than the lower part of the soil due to the low temperature. As the freezing temperature shows a decreasing, increasing, and then decreasing trend, the location of salt crystallization aggregation displays an increasing, decreasing, and then increasing trend, and the impact on the upper part of the sample is substantial when the location of salt crystallization is high.

The distribution of aqueous salt with different initial water contents (IWC) is shown in Figure 11. When the initial salt concentrations were all 0.6 mol/kg, the initial volumetric water contents were 0.1, 0.21, and 0.3, and the soil salt contents were 0.714%, 1.500%, and 2.143%, respectively. Thus, the freezing temperatures were all approximately −2.5°C. Due to similar initial salt concentrations, the soil freezing temperature is the same under different initial water contents, the location of the soil freezing front is consistent, and the locations of moisture salt migration and aggregation are also the same. The latent heat of moisture phase change has a certain delaying effect on cooling. Thus, the migration time of moisture and salt upward becomes longer when the water content is large, and the location of moisture and salt aggregation is slightly higher.

The distribution of aqueous salt with different initial salt concentrations is shown in Figure 12. When the salt concentration is 0.4 mol/kg, 0.6 mol/kg, 0.8 mol/kg, and 1.0 mol/kg, the initial volumetric water content is 0.21, the salt content of the soil is 1%, 1.5%, 2%, and 2.5%, respectively, and the freezing temperatures of soil are −2.35°C, −2.5°C, −2.01°C, and −2.22°C, respectively. As the salt concentration increases, the soil freezing temperature initially decreases, increases, and then decreases again, and the soil freezing front initially rises, falls, and then rises (2.01°C, −2.22°C). Moreover, the figures indicate that as the soil freezing temperature decreases, increases, and then decreases again with the rise in salt concentration, the soil freezing front rises, decreases, and then rises again, and there is a tendency for the location of the water–salt migration and aggregation to rise, fall, and then rise again. However, the location of water–salt aggregation is the same due to the small difference in freezing temperature.

The distribution of aqueous salt with different initial water contents is shown in Figure 13. When the initial NaCl concentrations were all 0.6 mol/kg, and the initial volumetric water contents were 0.1, 0.21, and 0.3, the NaCl content of the soil was divided into 0.294%, 0.618%, and 0.883%, respectively, and the freezing temperatures were all approximately −2.99°C. Due to the same initial salt concentration, the soil freezing temperature remains consistent under different initial water contents, the location of the soil freezing front is the same, and the location of moisture salt migration and aggregation also remains the same. The latent heat of moisture phase change has a certain delaying effect on the cooling. When the water content is large, the migration time of moisture and salt upward lengthens, and the location of moisture and salt aggregation is slightly elevated.

The distribution of aqueous salt with different initial NaCl concentrations is shown in Figure 14. The initial volumetric water content is 0.21; the NaCl concentrations are 0.4 mol/kg, 0.6 mol/kg, 0.8 mol/kg, and 1.0 mol/kg; the soil salt contents are 0.412%, 0.618%, 0.824%, and 1.030%, respectively; the solution freezing temperatures are taken as −1.46°C, −2.23°C, −2.98°C, and −3.74°C, respectively, and the soil freezing temperatures are −2.22°C, −2.99°C, −3.74°C, and −4.74°C, respectively. With the increase in salt concentration, the soil freezing temperature gradually decreases, the soil freezing front gradually moves upward, and the location of water–salt aggregation steadily rises. Compared with sulfate-salted soil, the permeability coefficient is relatively larger due to the absence of salts blocking the pores, and the upward water–salt migration is high. When the freezing temperature is high, the aggregation of the water–salt position in the upper soil is high, and the water–salt displays an oscillatory distribution; when the freezing temperature is low, the aggregation position of water in the upper soil decreases.