Numerical Analysis of the Influence of Foundation Replacement Materials on the Hydrothermal Variation and Deformation Process of Highway Subgrades in Permafrost Regions

Abstract

Affected by global warming, permafrost thawing in Northeast China promotes issues including highway subgrade instability and settlement. The traditional design concept based on protecting permafrost is unsuitable for regional highway construction. Based on the design concept of allowing permafrost thawing and the thermodynamic characteristics of a block–stone layer structure, a new subgrade structure using a large block–stone layer to replace the permafrost layer in a foundation is proposed and has successfully been practiced in the Walagan–Xilinji section of the Beijing–Mohe Highway to reduce subgrade settlement. To compare and study the improvement in the new structure on the subgrade stability, a coupling model of liquid water, vapor, heat and deformation is proposed to simulate the hydrothermal variation and deformation mechanism of different structures within 20 years of highway completion. The results show that the proposed block–stone structure can effectively reduce the permafrost degradation rate and liquid water content in the active layer to improve subgrade deformation. During the freezing period, when the water in the active layer under the subgrade slope and natural ground surface refreezes, two types of freezing forms, scattered ice crystals and continuous ice lenses, will form, which have different retardation coefficients for hydrothermal migration. These forms are discussed separately, and the subgrade deformation is corrected. From 2019 to 2039, the maximum cumulative settlement and the maximum transverse deformation of the replacement block–stone, breccia and gravel subgrades are –0.211 cm and +0.111 cm, –23.467 cm and –1.209 cm, and –33.793 cm and –2.207 cm, respectively. The replacement block–stone subgrade structure can not only reduce the cumulative settlement and frost heave but also reduce the transverse deformation and longitudinal cracks to effectively improve subgrade stability. However, both the vertical deformation and transverse deformation of the other two subgrades are too large, and the embankment fill layer will undergo transverse deformation in the opposite direction, which will cause sliding failure to the subgrades. Therefore, these two subgrade structures cannot be used in permafrost regions. The research results provide a reference for solving the settlement and deformation problems of subgrades in degraded permafrost regions and contribute to the development and application of complex numerical models related to water, heat and deformation in cold regions.

1. Introduction

From 2012 to 2020, the area of permafrost under the Earth’s surface decreased from 18 × 10⁶ km² to 16 × 10⁶ km² [1,2], and approximately 11% of permafrost disappeared permanently, leading to more subgrade thawing and settlement deformation [3,4,5]. The survey results of subgrade diseases in permafrost regions such as Siberia, Alaska and the Qinghai–Tibet Plateau show that the settlement caused by the thawing of shallow frozen soil on the ground surface has become the main form of subgrade damage, and the amount of deformation has increased with the increasing thawing depth of frozen soil [6,7,8]. The sixth IPCC climate change assessment report points out [9] that the mean annual air temperature (MAAT) between 2011 and 2020 was 1.09 °C higher than that from 1850–1900, making it the fastest warming decade in human history. In 2019, the highest air temperature in Siberia reached 38 °C, and the degradation of permafrost is accelerating at an unprecedented rate. As a result of climate and environmental change, permafrost is extremely sensitive to temperature change. On the one hand, the strength of the foundation soil decreases after permafrost thaws, and the uneven distribution of water in the subgrade will lead to uneven settlement deformation. Once the deformation exceeds the material limit, longitudinal cracks will occur [10,11,12]. On the other hand, the freezing and thawing action in the active layer is essentially a dynamic process of mutual coupling of water, heat and deformation. Even if there is no external load, the hydrothermal variation in frozen soil can cause obvious stress and deformation. The combined actions of these two aspects may eventually reduce subgrade stability [13,14,15,16]. Considering the background on the accelerated degradation of permafrost [9,17], it is of great practical significance to explore the interaction mechanism between water and heat and the deformation of the permafrost subgrade and establish a road structure system in cold regions to adapt to permafrost degradation.

At present, the prevention and control measures for permafrost subgrade deformation are mainly divided into two categories [18,19,20,21,22]. One category consists of passive protection measures, which mainly reduce heat exchange through measures such as blocking and sealing to maintain the original condition of ground temperature or slow down the degradation of permafrost. Conventional methods include setting thermal insulation materials in the embankment, reasonably increasing the embankment height or setting sunshades on the embankment slope to reduce solar radiation. However, it is difficult to prevent the degradation of permafrost under embankments by relying solely on passive protection measures [23,24,25]. The other category is positive cooling measures, mainly through engineering measures to actively reduce the temperature of permafrost. Conventional methods include setting ventilation ducts, thermosyphons or air–cooled crushed rocks in the embankment [13,14,15,26,27,28]. Among them, setting a crushed–rock layer in the embankment and using strong air convection heat transfer brought by high wind speed can significantly reduce the embankment temperature [28,29], which has a good application effect in low–temperature permafrost regions, such as Alaska and the Qinghai–Tibet Plateau, with more wind in cold seasons and less wind in warm seasons [30,31,32]. The above two measures are mainly used to reduce subgrade deformation by protecting permafrost. However, the field monitoring results [33,34,35] show that the application effect of these measures in high–temperature permafrost regions is not ideal. To date, the problem of subgrade deformation in degraded permafrost regions has not been effectively solved [36,37].

To study the prevention and treatment effect of the above measures on subgrade deformation, researchers have proposed heat transfer, water–heat coupling and water–heat–mechanical coupling models of frozen soil subgrades [38,39,40], which have been widely used in the calculation of subgrade stability in cold regions. However, most studies did not consider the influence of changes in liquid water content, water migration and the freezing forms of water in the active layer on subgrade deformation. These investigations can provide a reference for subgrade thermal–mechanical stability analysis to a certain extent but cannot predict the hydrothermal condition and deformation process of subgrades.

In view of the settlement and deformation of subgrades in degraded permafrost regions, based on the engineering practice of the Walakan–Xilinji section of the Beijing–Mohe Highway in the Da Xing’an Mountains and the design concept of allowing permafrost to thaw, a new subgrade structure, which uses large–size block stones to replace the permafrost layer under the natural ground surface, is proposed to reduce subgrade deformation. To analyze the working mechanism of this new structure and verify the prevention and control effect on the subgrade deformation, the replacement breccia subgrade structure and the replacement gravel subgrade structure are taken as controls, and a coupling model of liquid water, vapor, heat and deformation of unsaturated frozen soil subgrade is established to simulate the hydrothermal variation and deformation process of the subgrades within 20 years of the completion of the highway. In addition, the influence of two freezing forms of water in the active layer in winter (dispersed ice crystals and continuous ice lenses) on the subgrade hydrothermal condition are discussed, and the subgrade deformation is corrected so that the numerical simulation results can better reflect the actual subgrade deformation.

2. Study Area

This paper takes the Walagan–Xilinji section of the Beijing–Mohe Highway in China as the research object. The project location (52°38′42.88″ N, 124°27′04.71″ E to 52°56′16.90″ N, 122°28′33.41″ E) is shown in Figure 1. This highway is an ordinary class II highway with asphalt concrete pavement, and the width of the embankment is 12 m. Construction started in April 2016 and was completed in August 2019.

The study area in this paper is located in the hinterland of the Da Xing’an Mountains and the southern edge of the permafrost region in Eurasia. The climate is alternately affected by the Siberian high pressure and continental monsoon climate. The MAAT is low, approximately –2.3 °C. Affected by regional environmental changes, a wide range of inversion layers are distributed in this area, and permafrost is distributed in discontinuous island shapes [41,42,43]. Permafrost varies in thickness from several meters to tens of meters and is commonly found in swamps with thick peat accumulations. Since 2004, due to the severe degradation of permafrost, the surface frost number in the study area has dropped below 0.5 (Figure 2), and a distorted ground temperature curve has appeared [43]. Due to the high forest coverage along the highway in the study area [44,45], the surface roughness coefficient is large, and the wind speed is low. In addition, the slope of the embankment is covered by snow in cold seasons [46,47], which makes it difficult for heat to be discharged from the embankment, thus causing great difficulty in permafrost protection. Therefore, traditional prevention and control measures based on the design concept of protecting permafrost under subgrade conditions are unsuitable for regional highway construction.

3. Model Description

3.1. Geometric Model

Because the subgrade has linear characteristics along the length direction, the longitudinal influence is ignored. Based on the K308 + 200 section parameters of the Beijing–Mohe Highway (52°39′41.39″ N, 124°19′36.16″ E), a two–dimensional geometric model is adopted and divided into five calculation units (Figure 3). Part I is the subgrade fill layer above the ground surface, with a height of 3.5 m and a width of 6 m. The slope ratio of the subgrade is 1:1.5. Part II is the waterproof and heat insulation layer of peat clay below the ground surface, with a thickness of 1.5 m. Part III is an undisturbed breccia layer. Part IV is a moderately weathered andesite layer, which is a less frozen soil layer. Part V is the replacement layer, which belongs to the near ground surface rich frozen soil layer, with a thickness of 4.5 m. To reduce the influence of the boundary effect, the model horizontally extends 30 m from the right boundary of the replacement layer to the right and vertically extends 30 m downward from the bottom of the replacement layer.

To compare the influence of different replacement materials on subgrade stability, the replacement materials in part V are divided into three groups for comparison. Model 1: the replacement block–stone subgrade structure. The material in part V is block stones (diameter range from 20–40 cm), which do not contain fine particles. Model 2: considers the replacement breccia subgrade structure. The material in part V is breccia, which contains a large amount of fine–grained soil. Model 3: the replacement gravel subgrade structure. The material in part V is gravel, which contains a small amount of fine–grained soil.

To facilitate an understanding of the structure and composition of the subgrade soil layers, two site photos are provided (Figure 4).

3.2. Mathematical Model

Based on the nonisothermal water–heat–vapor coupled migration model, the water migration equation of the unsaturated frozen soil subgrade is established as follows [48]:

$$\frac{\partial {\theta }_w}{\partial t}+\frac{{\rho }_v}{{\rho }_w}\frac{\partial {\theta }_v}{\partial t}+\frac{{\rho }_i}{{\rho }_w}\frac{\partial {\theta }_i}{\partial t}=\mathit{\nabla }\left[K_{\mathrm{wh}}\left(\mathit{\nabla }h+1\right)+K_{\mathrm{wT}}\mathit{\nabla }T+K_{\mathrm{vh}}\mathit{\nabla }h+K_{\mathrm{vT}}\mathit{\nabla }T\right]$$

(1)

where θ_w is the volumetric liquid water content, expressed as ${\theta }_w=a{\left|T\right|}^b$; θ_v and θ_i are the volumetric contents of vapor and ice, respectively; ρ_w, ρ_v and ρ_i are the densities of liquid water, vapor and ice, respectively; h is the height of the water head; K_wh and K_wT are the isothermal and nonisothermal hydraulic conductivities caused by the matrix potential and temperature gradient, respectively; K_vh and K_vT are the isothermal and nonisothermal vapor conductivities, respectively.

According to the van Genuchten model, the expression of the soil–water characteristic curve of unsaturated soil can be expressed as follows:

$$\mathsf{\Theta }={\left[1+{(-\alpha h)}^n\right]}^{-m}$$

(2)

where Θ is the effective saturation and α, n and m are the model fitting parameters. Under isothermal conditions, the permeability coefficient of unsaturated soil can be expressed as follows:

$$K_{\mathrm{wh}}={10}^{-\mathsf{\Omega }{\theta }_i}{K}_{\mathrm{s}}{\mathsf{\Theta }}^l{\left[1-{\left(1-{\mathsf{\Theta }}^{1/m}\right)}^m\right]}^2$$

(3)

where Ω is an empirical parameter with dimensions of 1, l is a Mualem model parameter, and K_s is a saturated permeability coefficient.

The nonisothermal liquid water permeability coefficient K_wT, the vapor permeability coefficient K_vT caused by the temperature gradient, and the vapor permeability coefficient K_vh caused by the matrix potential gradient are listed as follows:

$$\left\{\begin{array}{c}{K}_{\mathrm{wT}}=K_{\mathrm{wh}}h{G}_{\mathrm{wT}}\frac{1}{{\gamma }_0}\frac{d\gamma }{dT}\\ K_{\mathrm{vh}}=\frac{D}{{\rho }_{\mathrm{w}}}{\rho }_{\mathrm{v}}\frac{Mg}{RT}{H}_{\mathrm{r}}\\ K_{\mathrm{vT}}=\frac{D}{{\rho }_{\mathrm{w}}}\eta H_{\mathrm{r}}\frac{d{\rho }_{\mathrm{v}}}{dT}\end{array}\right.$$

(4)

where G_wT is the empirical parameter for evaluating the influence of temperature on the soil–water characteristic curve and γ is the surface tension at 25 °C. M and g are the molar mass of water and the acceleration of gravity, respectively. R is the gas constant, H_r is the relative humidity, η is the water vapor diffusion enhancement factor, and D is water vapor diffusion in the soil.

According to energy conservation, the heat transfer equation of the unsaturated frozen soil subgrade can be expressed as follows [48]:

$$C\frac{\partial T}{\partial t}+L_{\mathrm{i}}{\rho }_{\mathrm{i}}\frac{\partial {\theta }_{\mathrm{i}}}{\partial t}+L_{\mathrm{v}}{\rho }_{\mathrm{v}}\frac{\partial {\theta }_{\mathrm{v}}}{\partial t}=\mathit{\nabla }\left(\lambda \mathit{\nabla }T\right)-C_{\mathrm{w}}\left[\mathit{\nabla }\left(q_{\mathrm{w}}T\right)\right]-C_{\mathrm{v}}\left[\mathit{\nabla }\left(q_{\mathrm{v}}T\right)\right]-L_{\mathrm{v}}{\rho }_{\mathrm{v}}\mathit{\nabla }{q}_{\mathrm{v}}$$

(5)

where C is the equivalent specific heat of the soil, λ is the equivalent thermal conductivity of the soil, and L_v is the latent heat of evaporation of water. The expressions of other hydrothermal–related parameters are shown in

Table 1. Expressions of other hydrothermal–related parameters [21,48].

Parameters	Expression
Surface tension	$\mathit{\gamma }=75.6-0.1425T-2.38\times {10}^{-4}{T}^2$
Saturated vapor density	${\rho }_{\mathrm{v}}=\exp \left(31.37-6014.79/T-7.92\times {10}^{-3}T\right)\times {10}^{-3}/T$
Enhancement factor	$\eta =9.5+3{\theta }_{\mathrm{w}}{\theta }_{\mathrm{sat}}{}^{-1}-8.5\left(1/\exp {\left(1+2.6{\left(f_{\mathrm{c}}\right)}^{-0.5}{\theta }_{\mathrm{w}}{\theta }_{\mathrm{s}}{}^{-1}\right)}^4\right)$
Relative humidity	$H_{\mathrm{r}}=\exp \left(hM\mathrm{g}/RT\right)$
Water head	$h=L_{\mathrm{i}}\left(T-T_{\mathrm{f}}\right)/\mathrm{g}{T}_{\mathrm{f}}$
Equivalent specific heat of the soil	$C=C_{\mathrm{s}}{\theta }_{\mathrm{s}}+C_{\mathrm{w}}{\theta }_{\mathrm{w}}+C_{\mathrm{i}}{\theta }_{\mathrm{i}}+C_{\mathrm{v}}{\theta }_{\mathrm{v}}$
Latent heat of water evaporation	$L_{\mathrm{v}}=2.501\times {10}^6-2369.2T$
Effective thermal conductivity Volumetric vapor content Vapor diffusion	$$\begin{gathered} \lambda ={\lambda }_{\mathrm{s}}{}^{{\theta }_{\mathrm{s}}}{\lambda }_{\mathrm{w}}{}^{{\theta }_{\mathrm{w}}}{\lambda }_{\mathrm{i}}{}^{{\theta }_{\mathrm{i}}}{\lambda }_{\mathrm{v}}{}^{{\theta }_{\mathrm{v}}} \\ {\theta }_{\mathrm{v}}={\rho }_{\mathrm{v}}{H}_{\mathrm{r}}\left({\theta }_{\mathrm{sat}}-{\theta }_{\mathrm{w}}\right)/{\rho }_{\mathrm{w}} \\ D=2.12\times {10}^{-5}{\left(T/273.15\right)}^2{\theta }_{\mathrm{v}}^{10/3}/{\theta }_{\mathrm{sat}}^2 \end{gathered}$$
Unfrozen water flux	$q_{\mathrm{w}}=-K_{\mathrm{wh}}\left(\mathit{\nabla }h+1\right)-K_{\mathrm{wT}}\mathit{\nabla }T$
Vapor flux	$q_{\mathrm{v}}=-K_{\mathrm{vh}}\mathit{\nabla }h-K_{\mathrm{vT}}\mathit{\nabla }T$

Note: In Table 1, λ_s, λ_w, λ_i and λ_v are the thermal conductivity of soil particles, unfrozen water, ice and vapor, respectively. f_c is the mass fraction of clay in the soil, and L_i is the latent heat of ice–water phase change. θ_s is the volume content of soil particles, θ_sat is the saturated water content of the soil, and T_f is the freezing temperature of soil (273.15 K). C_s and C_i are the specific heat of soil particles and ice, respectively.

.

According to the theory of viscoplasticity, the stress–strain relationship of frozen soil subgrades is as follows [49,50]:

$$\left\{\mathrm{\Delta }\sigma \right\}=\left[D_{\mathrm{T}}\right]\left(\left\{\mathrm{\Delta }\epsilon \right\}-\left\{\mathrm{\Delta }{\epsilon }_{\mathrm{vp}}\right\}-\left\{\mathrm{\Delta }{\epsilon }_{\mathrm{v}}\right\}\right)$$

(6)

$$\left[D_{\mathrm{T}}\right]=\frac{E_{\mathrm{T}}\left(1-v_{\mathrm{T}}\right)}{\left(1+v_{\mathrm{T}}\right)\left(1-2v_{\mathrm{T}}\right)}\left[\begin{array}{ccc}1& \frac{v_{\mathrm{T}}}{1-v_{\mathrm{T}}}& 0\\ \frac{v_{\mathrm{T}}}{1-v_{\mathrm{T}}}& 1& 0\\ 0& 0& \frac{1-2v_{\mathrm{T}}}{2\left(1-v_{\mathrm{T}}\right)}\end{array}\right]$$

(7)

where [D_T] is the elasticity matrix related to temperature and water content, E_T is the elastic modulus and v_T is Poisson’s ratio. {Δε}, {Δε_vp} and {Δε_v} are the total strain increment vector, the viscoplastic strain increment vector, and the frost heave strain increment vector, respectively.

In the two–dimensional stress state, the viscoplastic strain rate is as follows [49]:

$$\dot{{\epsilon }_{\mathrm{vp}}}={\gamma }_T⟨\mathrm{\Phi }\left(F\right)⟩\frac{\partial Q}{\partial \left\{\sigma \right\}}$$

(8)

where γ_T is the viscosity parameter, Q is the plastic potential function, and Φ(F) is a scalar function, which can be written as [49]

$$\mathrm{\Phi }\left(F\right)=\frac{F-F_0}{F_0}$$

(9)

$$⟨\mathrm{\Phi }\left(F\right)⟩=\left\{\begin{array}{cc}\mathrm{\Phi }\left(F\right)& \mathrm{\Phi }\left(F\right)>0\\ 0& \mathrm{\Phi }\left(F\right)\le 0\end{array}\right.$$

(10)

where F is the yield function and F₀ is the yield stress. Assuming that the plasticity criterion applies to frozen soil, according to the D–P yield criterion, Q is equivalent to F.

The strain caused by the water phase transition can be described as

$${\epsilon }^v=\delta \left[\left(1+\zeta \right)\left({\theta }_0+\mathrm{\Delta }{\theta }_{\mathrm{L}}-{\theta }_{\mathrm{w}}\right)-\left(n-{\theta }_{\mathrm{w}}\right)\right]$$

(11)

where ζ is the in situ frost heave rate, n is the porosity of the porous media, θ₀ is the initial volumetric water content, Δθ_L is the volumetric content of migrated water, and δ is defined as follows:

$$\delta =\left\{\begin{array}{cc}1& \left(1+\zeta \right)\left({\theta }_0+\mathrm{\Delta }{\theta }_{\mathrm{L}}-{\theta }_{\mathrm{w}}\right)-\left(n-{\theta }_{\mathrm{w}}\right)>0\\ 0& \left(1+\zeta \right)\left({\theta }_0+\mathrm{\Delta }{\theta }_{\mathrm{L}}-{\theta }_{\mathrm{w}}\right)-\left(n-{\theta }_{\mathrm{w}}\right)\le 0\end{array}\right.$$

(12)

Under plane strain conditions, {ε_v} can be expressed as follows [49]:

$$\left\{{\epsilon }_{\mathrm{v}}\right\}=\frac{{\epsilon }^v}{3}{\left\{\begin{array}{ccc}1+v_{\mathrm{T}}& 1+v_{\mathrm{T}}& 0\end{array}\right\}}^T$$

(13)

Equations (1)–(13) constitute the numerical model of liquid water, vapor, heat and deformation coupling of the unsaturated frozen soil subgrade.

3.3. Parameters of the Soil Layers

The geometric shape of the subgrade and the division of the calculation domain are shown in Figure 3, the physical parameters of the soil layer are shown in

Table 2. Physical parameters of the soil layers.

Soil Layers	λs (W/m·K)	Cs (J/kg·K)	Ks (m/s)	ρ (kg/m3)	θ0 (%)	n0	a	b
Embankment fill	2.116	1028.6	8.522 × 10−6	1940	14.2	0.233	0.082	−0.28
Peat clay	1.215	1404.7	3.043 × 10−6	1300	18.6	0.438	0.123	−0.20
Breccia	1.741	1226.6	7.496 × 10−6	1710	21.2	0.481	0.133	−0.25
Moderately weathered andesite	2.022	1161.4	5.137 × 10−6	1800	15.2	0.352	0.107	−0.18
Gravel	2.253	981.2	4.525 × 10−5	2080	14.7	0.313	0.085	−0.29
Block–stone layer	2.642	923.0	8.404 × 10-4	2700	10.2	0.350	0.067	−0.21

[13,18,19,51,52,53], and the mechanical parameters are shown in

Table 3. Mechanical parameters of the soil layers.

Soil Layers	a1 (MPa)	b1	a2	b2	a3 (o)	b3	a4 (MPa)	b4	c1 (MPa)	d1 (MPa)	e1 (MPa)	c2	c3 (o)	d2 (o)	e2 (o)	c4	c5 (MPa)	d3 (MPa)	e3 (MPa)
Embankment fill	48.43	27.26	0.35	−0.007	15	0.75	0.014	0.032	−31	−5.0	48.43	1.1	−3	−2	15	1.1	−0.013	−0.008	0.014
Peat clay	4.52	2.16	0.40	−0.008	12	0.67	0.068	0.108	−4	−0.4	4.52	1.2	−6	−4	12	1.2	−0.025	−0.016	0.068
Breccia	6.27	2.84	0.38	−0.006	18	0.55	0.046	0.065	−6	−0.6	6.27	1.2	−7	−5	18	1.2	−0.031	−0.028	0.046
Moderately weathered andesite	40.11	22.41	0.25	−0.004	26	0.84	0.140	0.149	−21	−4.1	40.11	1.2	−4	−3	26	1.2	−0.047	−0.035	0.140
Gravel	5.38	2.56	0.42	−0.004	8	0.08	0.001	0.046	−96	−13.2	5.38	1.1	−15	−8	8	1.1	−0.098	−0.079	0.001
Block–stone layer	86	24	0.30	−0.001	30	0.86	0	0.005	−16	−0.2	86	1.1	−1	−0.4	30	1.1	−0.001	−0.001	0

[22,49,50,54,55,56,57,58,59] The expressions of other mechanical parameters are as follows [39]:

$$E_{\mathrm{T}}=\left\{\begin{array}{cc}{a}_1+b_1{\left|T\right|}^{0.6}& T\le T_{\mathrm{f}}\\ c_1{\theta }_{\mathrm{w}}^2+d_1{\theta }_{\mathrm{w}}+e_1& T>T_{\mathrm{f}}\end{array}\right.$$

(14)

$$v_{\mathrm{T}}=\left\{\begin{array}{cc}{a}_2+b_2\left|T\right|& T\le T_{\mathrm{f}}\\ a_2& T>T_{\mathrm{f}}\end{array}\right.$$

(15)

$${\phi }_{\mathrm{T}}=\left\{\begin{array}{cc}{a}_3+b_3{\left|T\right|}^{c_2}& T\le T_{\mathrm{f}}\\ c_3{\theta }_{\mathrm{w}}^2+d_2{\theta }_{\mathrm{w}}+e_2& T>T_{\mathrm{f}}\end{array}\right.$$

(16)

$$c_{\mathrm{T}}=\left\{\begin{array}{cc}{a}_4+b_4{\left|T\right|}^{c_4}& T\le T_{\mathrm{f}}\\ c_5{\theta }_{\mathrm{w}}^2+d_3{\theta }_{\mathrm{w}}+e_3& T>T_{\mathrm{f}}\end{array}\right.$$

(17)

where a₁–a₄, b₁–b₄, c₁–c₅, d₁–d₃ and e₁–e₃ are all test fitting parameters.

3.4. Boundary Conditions

The temperature boundaries for the top surface of the subgrade, the side slope of the subgrade and the natural ground surface can be expressed as follows:

$$T=T_0+A\sin \left(2\pi t+{\alpha }_0\right)+\mathrm{\Delta }Tt$$

(18)

where α₀ is the initial phase angle, determined by the completion time of the subgrade. ΔT is the climate warming rate, taken as 0.052 °C/a [17]. The mean annual temperature T₀ and annual temperature amplitude A for each boundary are shown in

Table 4. Parameters of the temperature boundary conditions.

Variables	T0 (°C)	A (°C)
Natural ground surface: IJ	0.29	16.24
Side slope of the subgrade: JK	0.52	18.11
Top surface of the subgrade: AK	1.13	19.62

. The side ABCDE is a symmetrical boundary, the temperature boundary FGHI is adiabatic, and the geothermal flux at the bottom boundary EF is 0.03 W/m² [43]. The water boundary IJK is permeable, AK and FGHI are impermeable [50,60,61], and the bottom boundary EF has a 15.20% liquid water supply. The horizontal displacement for the lateral FGHI and the vertical displacement of the bottom boundary EF are restrained. The top surface of the subgrade, the side slope of the subgrade and the natural ground surface (boundaries AK, JK and IJ) are all free boundaries, and the displacement is not restricted.

3.5. Initial Conditions and Model Verification

To obtain the initial conditions of the model calculation, the model is solved according to the initial water, heat and displacement boundary conditions without considering climate warming. When the ground temperature outside the influence range of the subgrade under the natural ground surface tends to be stable, the calculated results of borehole ground temperature match the measured results, and then, the hydrothermal state of the subgrade will be used as the initial conditions of numerical simulation. Figure 5 shows that the measured values of ground temperature at the 1 m borehole outside the right slope toe of the subgrade from July 2018 to October 2019 are in good agreement with the calculated values. The water and heat state of the subgrade on 15 October 2019, are used as the initial conditions for the numerical simulation, and the water, heat and deformation of the subgrade within 20 years after the completion of the subgrade are simulated and calculated under the conditions of climate warming.

4. Results

4.1. Verification of the Temperature Field after Subgrade Completion

To further verify the accuracy of the numerical simulation, the temperature field after subgrade completion is verified. Figure 6 shows the comparison between the measured and simulated ground temperature values at the 1 m borehole outside the right slope toe from January 2020 to October 2021. The matching results are good and verify the accuracy of the model prediction.

4.2. Laboratory Test Verification

The numerical model was verified by one–sided freezing tests of unsaturated silt [54]. The height of the soil column is 12 cm, and the initial water content is 0.126. The top temperature of the soil column was –3 °C, the bottom temperature was 1 °C, and heat insulation treatment was carried out. During the test, the bottom boundary was supplied with pressureless water. Three parallel samples were designed for testing, and the total water content of the soil columns was measured by slicing after freezing for 5 h, 10 h and 40 h. The test and simulation results show that (Figure 7) the measured values match the simulated values, indicating that the model can accurately predict the hydrothermal change and deformation process of freezing soil.

4.3. Analysis of the Variation in Thermal Conditions in Different Subgrade Structures

To compare the effects of different replacement materials on the thermal condition of the subgrades, the changes in the temperature field of the three subgrade structures in cold and warm seasons were numerically analyzed.

The maximum air temperature in warm seasons in the Da Xing’an Mountains is generally reached in July every year, and the maximum seasonal thawing depth is generally reached in October. By comparing and analyzing the temperature contour maps on October 15 of the three subgrade structures in different years (Figure 8), the results show that the permafrost tables under the road centerlines of the three subgrade structures decrease to varying degrees with increasing time. From 2019 to 2029, the permafrost tables under the replacement block–stone subgrade, the replacement breccia subgrade and the replacement gravel subgrade move down from –1.421 m to –4.976 m, –5.906 m and –6.780 m, respectively. The average degradation rates of permafrost are 0.355 m/a, 0.449 m/a and 0.536 m/a. From 2029 to 2039, the permafrost tables under the road centerlines of the three subgrade structures drop to –5.695 m, –8.284 m and –9.789 m. The average degradation rates of permafrost are 0.072 m/a, 0.238 m/a and 0.301 m/a. The permafrost table of the replacement block–stone subgrade is higher than that of the other two subgrades, and the permafrost degradation rate is smaller.

This phenomenon is related to the heat transfer process in the subgrade and the heat transfer mechanism of the replacement materials. The block–stone layer is a porous medium with a large number of interconnected pores and a good thermal insulation effect [30,31,32,33,34]. In warm seasons, when the heat from the ground surface from solar radiation is transferred to the permafrost foundation, the low thermal conductivity of the block–stone layer (part V) makes less heat enter part V. Even if this heat enters part V, it is difficult for heat to diffuse to the deep foundation (part IV) [13,19]. Heat accumulates on the upper part of part V and causes the local temperature to rise. In cold seasons, after the ground surface temperature decreases, the temperature gradient between the upper part of part V and the ground surface is larger and the heat transfer efficiency is higher, so that part of the heat accumulated in the upper part of part V can be effectively discharged from the subgrade. This thermal regulation mechanism of the block–stone layer can greatly reduce the total heat entering the deep foundation [26,28], which improves the permafrost table and subgrade stability. The thermal regulation effect of the block–stone layer is directly proportional to its thickness in warm seasons [62], while in cold seasons, it is directly proportional to the amount of heat accumulated in the upper part of the block–stone layer. The more heat that accumulates, the higher the thermal regulation efficiency and the higher the permafrost table. Therefore, with increasing time, the decline rate of the permafrost table is smaller than that under the other two subgrades.

The replacement breccia and replacement gravel have relatively high thermal conductivities due to their relatively low porosity [63,64]. In warm seasons, when the heat from the ground surface is transferred to the permafrost foundation, the two materials cannot effectively reduce the total heat entering the replacement layer, resulting in most of the heat finally entering the deep foundation (part IV). Therefore, less heat accumulates in the replacement layer, and the temperature decreases. In cold seasons, after the ground surface temperature drops, the heat transfer rate between the replacement layer and the ground surface is relatively low, resulting in the heat in the foundation not being effectively discharged from the subgrade. Therefore, the permafrost tables under the road centerlines of these two subgrade structures are lower than those of the replacement block–stone subgrade.

The maximum seasonal freezing depth in the Da Xing’an Mountains is generally reached in March every year. By comparing and analyzing the temperature contour maps of the three subgrade structures on March 15 in different years (Figure 9), the results show that in 2039, the thicknesses of the thawed interlayer under the replacement block–stone subgrade, the replacement breccia subgrade and the replacement gravel subgrade are 5.016 m, 7.535 m and 8.918 m, respectively. The thawed interlayer thickness of the replacement block–stone subgrade is thinner than that of the other two subgrades, which shows that the block–stone layer can also play a significant role in thermal regulation in cold seasons. In 2039, the bottom of the thawed interlayer is still located in the block–stone layer, which shows that the block–stone layer can significantly reduce the degradation rate of permafrost. However, the permafrost in parts III and V in the replacement breccia subgrade and the replacement gravel subgrade is completely thawed. Due to the high ice content in these two soil layers, the thawing of the permafrost increases the subgrade settlement [55,65]. Therefore, these two subgrade structures can neither effectively reduce the temperature of the subgrade nor improve the deformation.

4.4. Analysis of the Variation in Liquid Water Content in Different Subgrade Structures

The degradation of permafrost can lead to an increase in the liquid water content in the subgrade. The strength damage caused by the increase in water content in the soil is the main reason for subgrade instability [65,66].

By comparing and analyzing the vertical distribution curves of the volumetric liquid water content under the road centerlines in different subgrade structures (Figure 10), in the early stage (2019–2021), the water contents in the replacement layer (part V) and the fill layer (part I) of the three subgrade structures are differentiated, which shows that the water content of the replacement block–stone subgrade in the above zones is significantly lower than that of the other two subgrades. From the long–term change process (2019–2039), this differentiation phenomenon becomes increasing obvious; that is, the water content of the replacement block–stone subgrade increases slightly in the above zones, while the other two subgrades demonstrate a significant water content increase. This phenomenon eventually leads to an uneven distribution of water content in part IV, which further reduces the foundation stability [66]. Among them, in part V, the water content of the replacement breccia subgrade is the largest, while in part I, the water content of the replacement gravel subgrade is the largest, which will lead to road fragmentation and foundation instability after repeated freezing and thawing [49,55].

The block–stone layer has the worst water–holding capacity [66], and the breccia layer has the strongest water–holding capacity because it contains a large amount of fine–grained soil, while the gravel layer has a slightly poor water–holding capacity because it contains a small amount of fine–grained soil [63,64]. The water–holding capacity of a material is proportional to its water content. Therefore, in part V, the liquid water content in the block–stone layer is always less than that in the gravel layer and the breccia layer.

In part I, the water content of the embankment fill layer changes greatly due to seasonal variations in the ground surface temperature and water supply in the thawed interlayer. However, by comparing the water contents in cold seasons and warm seasons, the water content on the top surface in part I of the replacement block–stone subgrade is significantly lower than that in the other two subgrades. This is related to the process of water migration and the size of the driving force [55,65,66]. In cold seasons, the temperature in the active layer decreases, and the water content in the thawed interlayer mainly migrates to the active layer under the action of a temperature gradient. In warm seasons, the temperature in the active layer increases, and the water in the thawed interlayer mainly migrates to the active layer under the action of matrix potential. The lack of fine particles in the block–stone layer results in an extremely low matrix potential. In addition, the temperature gradient between the block–stone layer and the active layer in cold seasons is relatively small. As a result, the water content supply in part I of the block–stone subgrade is lower. The water content in part I of the replacement gravel subgrade is greater than that of the replacement breccia subgrade because the larger temperature gradient between the gravel layer and active layer leads to a greater water supply in the active layer.

4.5. Deformation Analysis of Different Subgrade Structures

The thawing of permafrost leads to a reduction in soil strength. Under the action of its weight, the subgrade undergoes vertical deformation (settlement), which is the main source of deformation of the permafrost subgrade [56]. Generally, the deeper the thawing depth of permafrost is, the greater the settlement [49,50].

The permafrost table in the Da Xing’an Mountains is the lowest in October every year, and the subgrade settlement is the largest at this time. Comparing and analyzing the contour maps of vertical deformation of the three subgrade structures on October 15 in different years (Figure 11), the results show that the settlement of the replacement block–stone subgrade is mainly located at the shoulder, and the maximum settlement is –0.235 cm (the negative sign represents only the deformation direction, the same as below). The uplift deformation is mainly located at the slope toe, and the deformation first increases and then decreases with time. Additionally, the deformation reaches a maximum value of +0.695 cm in 2024, then begins to decrease, and finally drops to +0.346 cm in 2039. The maximum deformations located at the road centerline of the replacement breccia subgrade and the replacement gravel subgrade are mainly settlement, and the maximum values are –22.114 cm and –30.908 cm, respectively. There will also be uplift deformation at the slope toe of the two subgrades, and the deformation trends are the same as those of the replacement block–stone subgrade. This is the result of the compaction of embankment fill soil at the foot of the slope [49]. Unlike settlement, the influence of uplift deformation on subgrade stability is quite small and can be almost ignored.

Generally, the influence of vertical deformation on subgrade stability is often of greater concern, while the influence of transverse deformation is usually ignored.

The maximum transverse deformation of the replacement block–stone subgrade is mainly concentrated near the top of the right side of the upper boundary of the block–stone layer (Figure 12). From 2019 to 2024, the soil near the top moved +1.319 cm to the outside of the subgrade and then remained stable without increasing. This is because the interlocking structure formed between the block–stone particles make the block–stone layer have a certain rigidity and can effectively resist the local uneven deformation of the foundation [55,65]. Therefore, the maximum transverse deformation occurs at the top of the upper boundary of the block–stone layer. The transverse deformation trends of the replacement breccia subgrade and the replacement gravel subgrade are similar, and the deformation is mainly divided into two parts. One part is the transverse deformation moving toward the road centerline, which is located within 2.1 m below the top surface of the subgrade, with maximum deformations of –1.271 cm and –1.835 cm. The other part is the transverse deformation moving to the outside of the subgrade, which is mainly located inside the replacement layer (part V), with maximum deformations of +3.126 cm and +3.784 cm. This transverse deformation in the opposite direction is likely to lead to sliding failure of the subgrade. In addition, after ice crystals melt, breccia and gravel have great strength damage and cannot effectively resist uneven foundation deformation [63,64]. Therefore, the maximum transverse deformation does not appear at the boundary of the replacement layer.

5. Discussion

When the water in soil freezes, ice crystals discontinuously distributed on the surface of soil particles are generated [66]. Ice crystals have a retardation effect on water migration, which is measured by the retardation coefficient $\left(I={10}^{-\mathsf{\Omega }{\theta }_i}\right)$ [49,67]. However, the water in the active layer forms a relatively continuous distribution of ice lenses when it freezes in cold seasons. The retardation effect of ice lenses on the migration of water is much greater than that of discontinuously distributed ice crystals [68,69,70], but there is currently no specific calculation method.

To analyze the difference in the influence of the two freezing forms of water in the active layer in cold seasons on the effect of water migration, the discussion is divided into two cases. In the first case (Case 1), it is assumed that only discontinuous ice crystals form in the active layer in cold seasons, and then, the retardation effect of water migration is measured by the retardation coefficient. At this time, the convective flux of the convection boundaries (embankment slope and natural ground surface) is not zero. In the second case (Case 2), it is assumed that only completely closed and continuously distributed ice lenses form in the active layer in cold seasons. At this time, ice lenses can block water migration, and the convection flux at the convective boundaries is zero. However, in general, ice lenses do not completely block water migration, so the retardation effect of ice in the active layer on water migration should be between Case 1 and Case 2 [65,66].

5.1. Comparison of Permafrost Table Changes in Three Subgrade Structures in Case 1 and Case 2

First, the change trends of the permafrost table under the road centerlines in the three subgrade structures in Case 1 are compared and analyzed (Figure 13). The results show that the permafrost table under the replacement block–stone subgrade is higher than that of the other two subgrades, and the difference increases with time. In 2039, the permafrost table under the replacement block–stone subgrade, replacement breccia subgrade and replacement gravel subgrade will be reduced from –1.421 m to –5.695 m, –8.284 m and –9.789 m, respectively. The average degradation rates will be 0.214 m/a, 0.343 m/a and 0.418 m/a. The replacement block–stone subgrade structure can effectively reduce the degradation rate of permafrost.

Second, the change trends of the permafrost table under the three subgrade structures in Case 2 are compared and analyzed (Figure 13). The results show that in 2039, the permafrost table under the replacement block–stone subgrade, replacement breccia subgrade and replacement gravel subgrade will be reduced to –5.807 m, –8.839 m and –10.539 m, respectively. The average degradation rates will be 0.219 m/a, 0.371 m/a and 0.456 m/a. The permafrost table under the replacement block–stone subgrade is higher, and the degradation rate of permafrost is smaller. However, compared with Case 1, the permafrost tables under the three subgrades in Case 2 have decreased to varying degrees. This is because in Case 2, a part of the heat in the subgrades cannot be discharged by convection in winter, thus lowering the permafrost tables. However, only the thermal regulation efficiency of the block–stone layer is affected, and the essence of slowing the decline of the permafrost table does not change.

5.2. Comparison of the Volumetric Liquid Water Content of the Three Subgrade Structures in Case 1 and Case 2

In Case 1, the vertical distribution of the volumetric liquid water content under the road centerlines of the three subgrade structures (Figure 10) is analyzed in Section 4.4 and is not repeated here.

The relationship between the liquid water contents in parts I and V of the three subgrade structures from 2019 to 2039 in Case 2 is the same as that in Case 1 (Figure 14). Compared with Case 1, the distribution law of the liquid water content does not change in Case 2, but the absolute quantity changes. Due to the large increase in water contents in parts I and V of the replacement breccia subgrade and replacement gravel subgrade, a water supply in part IV is needed, resulting in an uneven distribution of water content in part IV, which further reduces the foundation stability [55,66,67].

The vertical distributions of the volumetric liquid water contents under the road centerlines of the three subgrade structures in cold seasons in different years are compared and analyzed (Figure 15). The results show that the liquid water content in part I of the three subgrade structures in Case 2 is greater than that in Case 1. This indicates that a considerable amount of water cannot be discharged from the subgrade in cold seasons when a closed and continuous ice lens forms in the active layer under the convection boundaries. This water remains in the subgrade, which increases the degradation rate of permafrost and then increases the subgrade settlement [19,35].

The permafrost degradation rate, liquid water content and settlement interact and promote each other. Therefore, the greater the liquid water content is, the greater the subgrade settlement [49,68].

5.3. Comparison of the Deformation at the Road Centerline and the Shoulder of the Three Subgrade Structures in Case 1 and Case 2

The difference between the two cases is that the convective fluxes at the convective boundaries are different, which includes both water convection and heat convection. The change in water and heat content in the subgrade affect the soil strength and then, the deformation [3,4,70].

In Case 1, the comparison of vertical deformations at the road centerline and the shoulder of the three subgrade structures (Figure 16) shows that from 2019 to 2039, the cumulative settlements at the road centerline and the shoulder (the following is replaced by the above two vertical deformations) for the replacement block–stone subgrade are –0.195 cm and –0.235 cm, respectively. The above two vertical deformations for the replacement breccia subgrade are –22.114 cm and –15.079 cm. The above two vertical deformations for the replacement gravel subgrade are –30.909 cm and –20.299 cm. In Case 2 (Figure 17), the above two vertical deformations for the replacement block–stone subgrade are –0.109 cm and –0.186 cm. The above two vertical deformations for the replacement breccia subgrade are –24.820 cm and –16.863 cm. The above two vertical deformations for the replacement gravel subgrade are –36.677 cm and –23.838 cm. The above two vertical deformations in Case 1 and Case 2 are different. The differences in cumulative settlement (the deformation in Case 2 minus the deformation in Case 1, and the following is replaced by the above two vertical deformation differences in Case 1 and Case 2) at the road centerline and the shoulder in Case 1 and Case 2 for the replacement block–stone subgrade are –0.086 cm and –0.049 cm, respectively. The above two vertical deformation differences in Case 1 and Case 2 for the replacement breccia subgrade are 2.706 cm and 1.784 cm, respectively. The above two vertical deformation differences in Case 1 and Case 2 for the replacement gravel subgrade are 5.768 cm and 3.539 cm, respectively.

The above two vertical deformations in Case 2 are smaller than those in Case 1 for the replacement block–stone subgrade. There are two reasons. First, the content of liquid water in the active layer in Case 2 is greater than that in Case 1. More liquid water participates in the ice–water phase transition, resulting in a greater frost heave in Case 2 [71,72,73,74]. Second, the block–stone layer can reduce the uneven settlement caused by water migration and the ice–water phase transition [19,65], so the increased amount of liquid water in Case 2 contributes little to the subgrade settlement deformation. The combined action of these two reasons finally leads to a decrease in cumulative settlement on the top surface of the subgrade and an increase in frost heave in Case 2. However, the cumulative settlement and frost heave of the replacement block–stone subgrade are significantly smaller than those of the other two subgrades, so the subgrade stability is better. The above two vertical deformations in Case 2 are larger than those in Case 1 for the replacement breccia subgrade and the replacement gravel subgrade. This is because a considerable amount of liquid water in Case 2 cannot be discharged from the subgrade in cold seasons. This liquid water will remain in the subgrades, which will reduce the foundation strength [55], resulting in an increase in subgrade settlement. Therefore, the vertical deformation stability of these two subgrade structures is worse.

In addition, the above two vertical deformations in Case 1 and Case 2 have the following common points. The maximum cumulative settlement of the replacement block–stone subgrade is located at the shoulder, while the maximum cumulative settlement of the replacement breccia subgrade and that of the replacement gravel subgrade are located at the road centerline. This is because there is a large internal friction angle between block–stone particles, which can form a firm interlocking structure and effectively resist the uneven settlement of the foundation caused by the thawing of permafrost [75]. At the same time, as the subgrade slope has no lateral earth pressure support, the cumulative settlement at the shoulder is slightly larger than that at the road centerline. In contrast, after the permafrost under the replacement breccia subgrade and the replacement gravel subgrade thaws, the internal friction angle, cohesion and elastic modulus of the breccia and gravel decrease rapidly and show obvious plasticity [63,64]. Under the self–weight load of the subgrade, the settlement deformation at the road centerline increases rapidly. These two replacement materials cannot form a strong interlocking structure, and it is difficult for them to resist the uneven deformation of the foundation. Therefore, the cumulative settlement at the road centerline is greater than that at the shoulder. Note that the gravel has liquefaction characteristics [76,77,78], and the strength of gravel decreases sharply after the ice crystals thaw, resulting in the cumulative settlement of the replacement gravel subgrade being greater than that of the replacement breccia subgrade.

Vertical deformation is the main factor affecting subgrade stability, while transverse deformation may cause longitudinal cracks in the subgrade, which will affect its stability and durability and even damage the subgrade structure [29,79].

In Case 1, the comparison of the transverse deformation at the road centerline and the shoulder of the three subgrade structures (Figure 18) shows that from 2019 to 2039, the maximum transverse deformations at the road centerline and the shoulder (the following is replaced by the above two transverse deformations) for the replacement block–stone subgrade are –0.111 cm and +0.128 cm, respectively. The above two transverse deformations for the replacement breccia subgrade are +0.749 cm and –1.271 cm. The above two transverse deformations for the replacement gravel subgrade are −0.509 cm and –2.043 cm. In Case 2 (Figure 19), the above two transverse deformations for the replacement block–stone subgrade are –0.056 cm and –0.146 cm. The above two transverse deformations for the replacement breccia subgrade are +0.278 cm and –1.113 cm. The above two transverse deformations for the replacement gravel subgrade are +0.821 cm and –2.516 cm. The above two transverse deformations in Case 1 and Case 2 are different. The maximum differences in transverse deformation in Case 1 and Case 2 (the deformation in Case 2 minus the deformation in Case 1, and the following is replaced by the above two maximum differences in transverse deformation in Case 1 and Case 2) at the road centerline and the shoulder for the replacement block–stone subgrade are –0.055 cm and –0.274 cm, respectively. The above two maximum differences in transverse deformation in Case 1 and Case 2 for the replacement breccia subgrade are +0.471 cm and –0.158 cm, respectively. The above two maximum differences in transverse deformation in Case 1 and Case 2 for the replacement gravel subgrade are +1.330 cm and –0.473 cm, respectively.

Figure 15, Figure 16, Figure 17 and Figure 18 show that the vertical and transverse deformations of the three subgrades are relatively large from 2019 to 2026, and the deformation rate and size are also different. From 2026 to 2039, the vertical and transverse deformations will fluctuate slightly and gradually become stable [49,68]. Although the final deformations are different, the change trends are the same, which shows that the influence of Case 1 and Case 2 on subgrade deformation is mainly reflected in the rapid deformation stage. This result is because from 2019 to 2026, due to the rapid decline in the permafrost table, the water and heat exchange between the thawed interlayer and the active layer is in the drastic change stage (Figure 7, Figure 8 and Figure 9). The rapid discharge of water from the thawed interlayer will lead to a rapid increase in subgrade settlement [49,55]. This will cause uneven water distribution (Figure 10 and Figure 14) and further lead to uneven subgrade settlement (Figure 11). The more uneven the water distribution is, the more obvious the uneven settlement, which eventually leads to changes in the size and direction of transverse deformation. Due to the lower liquid water contents in parts I and V of the replacement block–stone subgrade structure, water migration and the phase transition have less impact on the subgrade deformation. In addition, the block–stone layer can resist the uneven deformation of the foundation well. Therefore, the vertical and transverse deformations are smaller, and the subgrade stability is better. Due to the poor water stability of gravel [76,77,78], the uneven distribution of water causes the replacement gravel subgrade to exhibit the worst stability after the thawing of permafrost. Although the liquid water content in part V of the replacement breccia subgrade is the largest, the mechanical stability of breccia is better than that of gravel [63,64]. In addition, the liquid water content in part I is lower, so the vertical and transverse deformations are smaller than those of the replacement gravel subgrade. From 2026 to 2039, the decline rate of the permafrost table of the three subgrades gradually decreases, and the water and heat exchange between parts I and V gradually slows down. At this time, the three subgrades enter the slow deformation stage, and the vertical deformation and transverse deformations fluctuate little and tend to become stable gradually [49,55,65].

In addition, the following two points need to be emphasized. One is that excessive transverse deformation of the subgrade may cause longitudinal cracks [29]. Generally, when the transverse deformation exceeds the material limit deformation, the width of the longitudinal crack increases with increasing transverse deformation [79,80,81]. The transverse deformation of the top surface of the replacement gravel subgrade is the largest, and the subgrade structure may be the first to reach the limit deformation. The transverse deformation of the top surface of the replacement block–stone subgrade is the smallest, and it is the least likely to reach the limit deformation. The transverse deformation of the top surface of the replacement breccia subgrade lies between the two. Therefore, it is speculated that the width of the longitudinal cracks on the top surface of the replacement gravel subgrade may be the largest, followed by the replacement breccia subgrade, and the width of the longitudinal cracks on the top surface of the replacement block–stone subgrade is the smallest (and in some cases, no cracks exist). Second, the settlement at the road centerline of the replacement breccia subgrade and the replacement gravel subgrade is the largest, resulting in tensile stress on the surrounding soil. As a result, the transverse deformation within 2.1 m below the top surface of the subgrade fill layer is negative. Therefore, the longitudinal cracks for these two subgrades will expand to the road centerline. However, the transverse deformation at the shoulder for the replacement block–stone subgrade is the largest, so the longitudinal cracks will expand to the shoulder [80,81].

Since the actual retardation effect of ice in the active layer on water migration is between Case 1 and Case 2, the average values (Figure 20 and Figure 21) of vertical and transverse deformations at the road centerline and the shoulder in Case 1 and Case 2 are used as the actual deformations of the subgrades.

Summarizing the actual maximum deformations of the top surface of the three subgrade structures (

Table 5. Actual maximum deformations of the top surface of the three subgrade structures from 2019 to 2039.

Subgrade Structures	Cumulative Settlement at the Road Centerline (cm)	Cumulative Settlement at the Shoulder (cm)	Maximum Transverse Deformation at the Road Centerline (cm)	Maximum Transverse Deformation at the Shoulder (cm)
Replacement block–stone subgrade structure	−0.153	−0.211	−0.077	+0.111
Replacement breccia subgrade structure	−23.467	−15.971	+0.495	−1.209
Replacement gravel subgrade structure	−33.793	−22.068	+0.395	−2.207

), the results show that the maximum cumulative settlement and maximum transverse deformation of the replacement block–stone subgrade are –0.211 cm and 0.111 cm, respectively. The maximum cumulative settlement and maximum transverse deformation of the replacement breccia subgrade are –23.467 cm and –1.209 cm, respectively. The maximum cumulative settlement and maximum transverse deformation of the replacement gravel subgrade are –33.793 cm and –2.207 cm, respectively. The replacement block–stone subgrade structure can not only reduce the cumulative settlement and frost heave but also reduce the transverse deformation and longitudinal cracks to improve the overall stability of the subgrade. However, the replacement breccia subgrade structure and the replacement gravel subgrade structure can neither reduce the cumulative settlement and frost heave nor reduce the transverse deformation and longitudinal cracks. Therefore, these two subgrade structures cannot be used in permafrost regions.

6. Conclusions

In view of the thawing settlement deformation of the permafrost subgrade caused by climate warming and construction disturbance, based on the engineering practice of the Walagan–Xilinji section of the Beijing–Mohe Highway and the design concept of allowing permafrost thawing, a new subgrade structure using a large block stone to replace the permafrost layer of the foundation is proposed to reduce subgrade settlement. To compare and explore the influences of different foundation replacement materials on the stability of highway subgrades in permafrost regions, the water, heat and deformation of the replacement block–stone subgrade, replacement breccia subgrade and replacement gravel subgrade from 2019 to 2039 are numerically analyzed. In addition, the following main conclusions are obtained:

(1)

The replacement block–stone subgrade structure has a higher permafrost table in warm seasons and a thinner thawed interlayer in cold seasons. Compared with the replacement breccia layer and the replacement gravel layer, the replacement block–stone layer can effectively reduce the total heat entering the deep foundation in warm seasons. Its good thermal regulation performance reduces the impact of seasonal temperature changes on subgrade temperature to improve the permafrost table.

(2)

The low thermal conductivity of the block–stone layer reduces the temperature gradient between the block–stone layer and the active layer. In addition, the very small matrix potential of the block–stone layer leads to a very small driving force of water migration. The reduction in the driving force of water migration reduces the water supply from the thawed interlayer to the active layer, which reduces the liquid water content in the active layer and finally, reduces the transverse and vertical deformations of the subgrade.

(3)

The replacement block–stone layer can effectively resist the local uneven deformation of the foundation, so the settlement of the replacement block–stone subgrade structure is the smallest. After the thawing of the permafrost in the breccia layer and the gravel layer, the strength of these two replacement materials decreases rapidly, and it is difficult for them to resist the uneven deformation of the foundation. Therefore, the cumulative settlement of the two subgrades is relatively large. In addition, the liquefaction characteristics of the gravel lead to a greater settlement of the replacement gravel subgrade than that of the replacement breccia subgrade.

(4)

In cold seasons, when the water in the active layer freezes, two forms of dispersed ice crystals and continuous ice lenses form, which have different retardation effects on water migration. We discussed these effects and corrected the subgrade deformation. The results show that from 2019 to 2039, the maximum cumulative settlement and the maximum transverse deformation of the replacement block–stone subgrade are –0.211 cm and +0.111 cm, respectively. The maximum cumulative settlement and the maximum transverse deformation of the replacement breccia subgrade are –23.467 cm and −1.209 cm, respectively. The maximum cumulative settlement and the maximum transverse deformation of the replacement gravel subgrade are –33.793 cm and –2.207 cm, respectively. The replacement block–stone subgrade structure can not only reduce the cumulative settlement and frost heave but also reduce the transverse deformation and longitudinal cracks to improve the overall stability of the subgrade. In contrast, the vertical and transverse deformation of the replacement breccia subgrade and the replacement gravel subgrade are too large, and even the subgrade fill layer will undergo transverse deformation in the opposite direction, which will cause sliding failure. Therefore, these two subgrade structures cannot be used in permafrost regions.

In addition, since the study area is located in permafrost regions with high temperatures and high ice contents, the deformation process of the subgrade is greatly affected by changes in water and heat. The numerical model fully considers the influence of water migration, the phase transition process and the phase transition forms on the stability of the subgrade. It is of great significance to understand the interactions among water, heat and deformation in the permafrost subgrade, which can provide a reference for the development and application of hydrothermal–related models in cold regions. In the case of the accelerated degradation of permafrost, there are more subgrade thawing settlement problems. The results of the study can provide a reference for subgrade design and protection in degraded permafrost regions.