#### Abstract

In this paper, the law of ice-rich permafrost embankment thaw consolidation is studied based on three-dimensional nonlinear large strain thaw consolidation theory. To avoid problems associated with numerical simulation efficiency and stability when a nonlinear stress-strain relationship is employed, a segment interpolation function is used to implement the nonlinear relationship between the compression modulus and the void ratio, and the corresponding simulation strategy is proposed. Through a comparison of the monitoring and calculated results, it is indicated that the calculation accuracy on ice-rich embankment thaw settlement can be notably improved after nonlinear theory is implemented with the proposed numerical simulation method. A further analysis of the calculated results indicates that the interactive effects between the thermal and mechanical fields can be more reasonably described by nonlinear theory than by linear theory. It is also determined that the postthaw pore water in the shallow embankment dissipates in the early operation period, while in the following long operation period, the development of the permafrost embankment thaw settlement is mainly due to the dissipation of newly postthaw pore water at the thaw depth or the permafrost table. This is one of the main differences in the law of permafrost embankment thaw settlement compared with that of unfrozen embankments.

#### 1. Introduction

In cold regions, a large number of engineering investigations show that structures built on frozen soil will greatly change the original thermal balance of the frozen soil, leading to significant thawing of permafrost. Along with the drainage of pore water, remarkably uneven settlement of the underneath postthaw layers can often be observed [1–3]. Intensive investigations have indicated that the thaw settlement of permafrost is one of the main causes of infrastructure damage, especially for ice-rich permafrost embankments [4, 5]. Therefore, studying the thaw consolidation law of permafrost and accurately predicting the process are practically significant.

Thaw settlement of frozen soil is a complex physical and mechanical process with continuous thawing of frozen soil and the drainage of postthaw pore water [6]. Thaw consolidation theories based on thermal conduction equations and consolidation theories are usually used to describe the thaw consolidation process. Based on the Neumann analytical solution, Morgenstern and Nixon [7] obtained the semiempirical formula of the thawing front of frozen soil, and a one-dimensional thaw consolidation theory was established by using the Terzaghi consolidation theory in the postthaw area. On the basis of Morgensternʼs work, Nixon further applied the theory to problems with changing thermal boundaries [8]. However, due to the small strain assumption, the prediction errors of the abovementioned theories become very large when large soil deformation is encountered [9]. To overcome the defects of the small strain theories, Foriero and Ladanyi [10] proposed a one-dimensional large strain thaw consolidation theory based on convection coordinates. The theories proposed above are based on linear elastic stress-strain relationships. The prediction errors of these theories still cannot be ignored when ice-rich permafrost is involved. The stress-strain relationship of ice-rich frozen soil shows a strong nonlinearity when it is postthaw. Obviously, the linear elastic stress-strain relationship cannot reasonably describe the nonlinear mechanical behavior [11–13]. Correspondingly, the calculation accuracy is greatly reduced by using elastic theory. Considering this, Dumais and Konrad [12] proposed a nonlinear large strain theory by taking the effective stress-void ratio-hydraulic conductivity nonlinear relationship into account. Based on this theory, Dumais and Konrad [13] studied the law of thaw consolidation of the Inuvik experimental warm-oil pipeline built on ice-rich permafrost foundations. In the abovementioned one-dimensional theories, the void ratio is used as an independent variable, which leads to large calculation errors when describing practical problems with complex boundaries. For the convenience of engineering applications, a three-dimensional thaw consolidation theory for large strain is established based on Biotʼs three-dimensional consolidation theory and the thermal conduction equation considering ice-water phase changes, and the applicability of the theory is verified via laboratory tests [9]. Yao et al. [14] further incorporated a nonlinear stress-strain relationship between the compression modulus and the void ratio into the three-dimensional large strain thaw consolidation theory. However, when the nonlinear stress-strain relationship was directly employed in the numerical analysis on complex boundary problems, the authors found that the extremely low initial compression modulus after the frozen soil thawed usually led to a grid collapse in the numerical calculation. Subsequently, the calculation efficiency and stability are considerably reduced. To avoid this problem, the corresponding numerical implementation strategies must be employed.

In this paper, a segment interpolation function is used to implement the nonlinear relationship between the compression modulus and the void ratio in numerical calculations, and a corresponding simulation strategy is employed to avoid the efficiency and stability problems of numerical calculations. The applicability of the proposed nonlinear numerical simulation method is verified through a comparison with in situ monitored data and the calculation results of linear elastic theory.

#### 2. Nonlinear Thaw Consolidation Theory and Numerical Implementation

##### 2.1. Nonlinear Large-Strain Thaw Consolidation Theory

It is assumed that the postthaw regime is saturated. In the postthaw regime, the hydraulic permeability and mechanical properties are considered to be independent of the temperature, and the medium volume changes related to temperature are not considered. Displacement and drainage are assumed to be zero in the frozen domain. Subsequently, consolidation only occurs in the postthaw domain.

Based on the abovementioned assumption, the thaw consolidation theory for saturated soils includes two parts: (I) the thermal conduction equation with the ice-water phase change is used to determine the postthaw area [1, 9]:where *T* is the temperature (°C), *t* is the time (sec), (W/m^{3}) is the volumetric heat source intensity, and *c* (J/kg·°C) and *λ* (W/m·°C) are the specific heat and thermal conductivity, respectively. For the considerable effect of the ice-water phase change, *c* and *λ* are temperature-dependent, where details can be found in the literature [1, 15].

(II) Nonlinear large strain consolidation theory for saturated soils is employed in the postthaw area to calculate the coupling process of soil skeleton compression and pore water dissipation. Consolidation theory for saturated soils includes [1, 9]where is the corrected total stress rate, is the instantaneous velocity of the material point, (m/s^{2}) is the gravitational acceleration, is the mediumʼs density, *E* is Youngʼs modulus, is Poissonʼs ratio, *u* is the pore pressure, *δ*_{ij} is the Kronecker symbol, *k* (m/s) is the hydraulic permeability of soil, (kg/m^{3}) is the fluid density, (1/s) is the volumetric fluid source intensity, (m/s) is the superficial velocity of the fluid relative to the soil skeleton, and is the symmetric deformation tensor, which is defined as

In consolidation theory, the constitutive equation is critical for determining the calculation accuracy of thaw consolidation settlement and pore water pressure. Most recent thaw consolidation theories are based on the linear elastic stress-strain relationship, as shown in equation (3). When ice-rich frozen soil is involved, the linear elastic stress-strain relationship usually causes large errors. Correspondingly, the linear stress-strain relationship (equation (3)) was replaced with a nonlinear stress-strain relationship as follows.

In equation (3), Youngʼs modulus can be expressed by the compression modulus aswhere is the effective pressure under *K*_{0} compression conditions and is the volumetric strain. The relationship between the effective pressure and the void ratio iswhere

By substituting equations (10) and (11) into (8), we have

Equation (12) is the nonlinear relationship between the compression modulus and the void ratio derived from the relationship obtained from the *K*_{0} compression tests. When equation (12) is incorporated into equation (3), the nonlinear stress-strain relationship of postthaw ice-rich frozen soil can be described [14].

##### 2.2. Numerical Implementation

###### 2.2.1. Nonlinear Stress-Strain Relationship

The nonlinear relationship between the void ratio and the compression modulus is used to describe the nonlinear stress-strain relationship of postthaw ice-rich frozen soil. However, direct application of equation (12) to the numerical calculations will considerably reduce the calculation efficiency and stability. The compression modulus of soil with high ice content will increase gradually from a very low initial value after thawing. In the early stage of compression, the compression modulus is extremely small, and a small increase in the effective stress will lead to a sharp decrease in the void ratio, which often leads to a rapid decrease or disappearance of the mesh scale in practical calculations. The rapid reduction or disappearance of the mesh size will further reduce the calculation efficiency and stability [16, 17]. Therefore, a piecewise linear interpolation function is used to numerically implement the nonlinear relationship between the compression and void ratio to guarantee the calculation accuracy and efficiency.

According to equation (12), for a certain void ratio , the corresponding compression modulus iswhere *m* = 0, 1, 2, 3, …, *n*. The curve of the compression modulus and the void ratio is divided into *n* sections. Between two different points ( and ), the compression modulus can be calculated by a linear interpolation function, i.e.,where the void ratio can be determined by the volume strain (), which can be calculated with equation (11) as

Equation (15) is transformed into

In this paper, the incremental form is used to calculate the stress and strain. Therefore, after each time step Δ*t*, the incremental expression of the void ratio corresponding to the volume strain is

The corresponding incremental form of the compression modulus is

In the thaw consolidation calculation, equation (18) is used to calculate the compression modulus according to the actual void ratio of each calculating zone.

###### 2.2.2. Calculation Process

① The temperature field of a time step Δ*t* is calculated, and the postthaw area is determined. ② Consolidation calculation in the postthaw area. According to the volumetric strain of the soil in the thawed area, the compression modulus of each calculation zone is determined using equation (18). On this basis, the consolidation calculation at the same time step Δ*t* is performed.

In each calculation cycle, ① and ② are performed until the target calculation time is reached.

When equation (18) is used to calculate the thaw consolidation problem of frozen soil under a concentrated load, the low compression modulus in the initial stage of consolidation still results in a decrease in the computational efficiency or a poor computational stability [18]. According to previous studies, the rapid consolidation of soils near the concentrated loads will be completed in the early stage of thawing consolidation [19]. To ensure both the calculation efficiency and stability, the compression modulus of the grid near the concentrated load can be taken as the secant modulus of the soil under the corresponding pressure conditions. According to equation (19), this secant modulus can be expressed as

In the case of an underlying high compressibility soil, the compression modulus of the corresponding area can be directly calculated using equation (18) because the volume of the ice-rich layer does not change rapidly after thawing.

#### 3. Case Study

##### 3.1. Site Description

An embankment section in the Kaixinling hilly area on the Qinghai-Tibet Plateau is taken as the study object. The soil strata are listed in Table 1. According to the real size of the embankment, the dimensions of the numerical model are shown in Figure 1. The distribution of the temperature monitoring boreholes and displacement monitoring points is listed in Figure 1. The thermistors are located 0.5 m from the surface with a density of 1/0.5 m. The displacement monitoring points are located in the center of the embankment.

The Qinghai-Tibet Highway, which was originally built in the 1950s, has been rebuilt several times [15]. The surface thermal conditions have changed greatly, and the lack of early embankment ground temperature data has made it difficult to simulate embankment thaw consolidation in real time. Therefore, the natural ground temperature of the embankment section is used as the initial temperature condition, and the measured surface temperature is taken as the boundary condition to simulate the thaw consolidation process of the selected permafrost embankment section over 50 years. The thermal and mechanical parameters of each layer of the embankment are shown in Table 1. The temperature parameters are obtained according to the Interim Provisions for engineering design of the Qinghai-Tibet Plateau permafrost regions [20]. The permeability is obtained according to the literature [1, 21]. The compression modulus of the embankment fill and mudstone in the embankment section is obtained according to the work of Wang et al. [15]. To further analyze the difference between the linear elastic consolidation theory and nonlinear consolidation theory, the thawing consolidation of silty clay and clay layers with a high ice content in Table 1 was calculated with different stress-strain relationships. Only an elastic stress-strain relationship is adopted in other low ice content soil layers (embankment fill and mudstone). In the elastic mode, when the average effective stress of the subsoil is 100 kPa, the compression modulus is obtained when the average effective stress is 100 kPa, which can be obtained from equation (19). The calculation results are listed in Table 1. In the nonlinear mode, the initial compression modulus of this part is calculated using equation (12). The initial void ratio of the silty clay is 1.08, and *λ* is 0.052. The initial void ratio of clay is 1.35, and *λ* is 0.080 [14].

##### 3.2. Boundary Conditions and Accuracy Validity

Based on previous studies [15, 22, 23], the temperature of the pavement, side slope, and natural ground surface follows the equationwhere *T*_{0} is the mean annual temperature, *θ* is the climatic warming velocity (0.02°C/year) [1, 23], *l* is the temperature amplitude, and *m**π* is the initial phase. A regression analysis on the monitoring data is performed to obtain the annual mean temperature and amplitude (as listed in Table 2). The calculation begins at the highest temperature in the summer, then the initial phase can be taken as *π**/*2, and the corresponding ground initial temperature can be specified according to the monitored temperature data in the field in Figure 2. According to Figure 2, the heat flux at the lower boundary is taken as 0.018 W/m^{2} [1].

#### 4. Results and Analysis

##### 4.1. Thaw Depth and Settlement

Figure 3 gives the thaw depth of the calculated results for both the linear and nonlinear modes over 50 years. In Figure 3, the black spot denotes the thaw depth of the monitored data in 2004. The results show that the monitored data in 2004 are approximately the same as the calculated thaw depth in the 30th year with the linear elastic mode. For the nonlinear results, the thaw depth measured in 2004 is approximately the same as the calculated results in the 25th year. Therefore, the calculated settlements after 30 years with the linear mode and 25 years with the nonlinear mode are compared with the monitored settlement.

Figure 4 shows the calculated thaw depth development at the pavement center with both the linear and nonlinear modes. From Figure 4, the thaw depth at the pavement center is approximately linear with time. Figure 5 shows the calculated settlement of different modes and the monitoring data. It can be determined from Figure 5 that the development of thaw consolidation settlements calculated by different stress-strain relationships is linearly related to time. According to previous studies [1, 24–26], the development of thaw settlement of frozen soil is dominantly controlled by the development of the thaw depth. By comparing the results of two different stress-strain relationships in Figure 5, it is determined that the settlement results of the nonlinear mode are larger than those of the linear mode, and the nonlinear results are closer to the monitored results. By comparing the thaw depth and settlement developments, it can be determined that the thaw depth calculated by the nonlinear mode is also larger than that of the linear elastic relationship (Figure 4). This indicates that the development of the soil deformation field will further influence the development of the temperature field. With the development of the thaw settlement, the thermal conduction path between the temperature boundary of the embankment and the subsurface soil will be shortened, and the thaw speed of the soil will also be accelerated [27–31]. Therefore, the rapid thaw settlement rate obtained with the nonlinear mode will inevitably lead to the rapid development of the thaw depth. At the same time, with the rapid development of thaw depth, more postthaw soil will further promote the development of the thaw settlement. Therefore, it can be said that the deformation and temperature fields of frozen soil will influence one another in the thaw settlement process. In addition, the calculation results of different stress-strain relationships show that the thaw consolidation theory based on the nonlinear stress-strain relationship can further improve the calculation accuracy of the thaw settlement for ice-rich permafrost.

What should be noted is that the monitored results are far greater than the predicted results, even though the nonlinear mode has a higher accuracy. There are two main reasons for this. First, this paper only considers embankment thaw consolidation under self-weight conditions. In practice, vehicle loading on the embankment pavement will cause additional local deformation and overall consolidation settlement deformation. Second, according to previous studies [32–34], the creep of high-temperature frozen soil also contributes to part of the embankment settlement. These factors will be taken into consideration in future research work.

##### 4.2. Stress Field

To further investigate the effects of different stress-strain relationships on the thaw consolidation behavior of permafrost embankments, the vertical effective stress fields in postthaw areas are analyzed (Figures 6 and 7). According to a previous study [1], the thaw consolidation process in the embankment fill usually finishes during the initial operation period (less than 5 years), and the subsequent thaw settlement is mainly due to the thaw consolidation of the deeper permafrost layer or the new postthaw layer. In addition, according to the calculated results at 10 and 30 years (Figures 8 and 9), the pore water flows down and toward the outside of the embankment slope side; that is, after the initial operation period, there is no pore water pressure in the embankment fill. Correspondingly, the effective stress calculated in the embankment fill is approximately equal to the total stress, and the effective stress in Figures 6 and 7 can be used to represent the embankment stress state. Therefore, in Figures 6 and 7, the calculated vertical effective stress was used to represent the stress state of the postthaw area.

**(a)**

**(b)**

**(a)**

**(b)**

As shown in Figures 6 and 7, the distribution range and the development of the vertical effective stress are similar at 10 and 30 years, i.e., the stress contour lines are distributed in the shape of a pot, which is due to the effects of the shape of the thaw depth (Figure 3); with the development of time, the range of the contour lines is enlarged, and the distribution form is geometrically consistent. By comparing the stress contour lines with different stress-strain relationships, it can be seen that the area encircled by the maximum contour lines calculated with the nonlinear mode (such as the area encircled by the contour lines at 160 kPa, Figure 7) is larger than that of the linear mode. On one hand, as described above, the thaw depth calculated with the nonlinear mode is larger than that of the linear mode (Figure 3). Correspondingly, the calculated effective stress of the nonlinear theory at the thaw depth is higher than that of linear theory due to the effect of the soil mass self-weight; on the other hand, as has been proven by previous research work [14, 35, 36], under the same overloading (), the consolidation degree (*U*_{t}) (determined with overloading and effective pressure ) obtained via the nonlinear mode is larger than that of the linear mode. According to the definition of consolidation degree at a certain soil layer,

Certainly, a larger consolidation degree means a larger effective pressure. As shown in Figure 10, at a certain subsurface soil layer with a self-weight of , the effective pressure determined by the nonlinear relationship is larger than the linear pressure (). This is why the distribution range of maximum stress obtained by the linear stress-strain relationship is smaller than that obtained by the nonlinear relationship at different time points. At the same time, according to the definition of consolidation degree *U*_{t} (equation (21)), the thaw consolidation degree of ice-rich permafrost embankments would be underestimated using a linear stress-strain relationship. Therefore, the large strain thaw consolidation theory with a nonlinear stress-strain relationship must be used to improve the calculation accuracy of the relevant design parameters in the stability design of permafrost engineering.

##### 4.3. Pore Water Pressure

In the above, it is indicated that the stress-strain relationship has a dominant effect on the distribution of effective stress. The consolidation process is mainly affected by the pore water pressure distribution [37–39]. Figures 8 and 9 show the calculated pore water flow direction and the contour lines at 10 and 30 years with nonlinear stress-strain relationships, where the arrows indicate the pore water flow direction in the postthaw area. As shown in Figures 8 and 9 at different times, the excess pore water in the postthaw area is drained out toward the embankment slope and the natural ground surface. This is due to the effects of the embankment self-weight and boundary conditions (where asphalt pavement is impermeable). With the excess pore water pressure contour lines, it can be seen that most of the excess water pressure is mainly concentrated near the thaw depth area, while in the shallow area, most of the excess water pressure is dissipated at both 10 and 30 years; that is, the postthaw pore water in the shallow embankment dissipates in the early operation period. In the following long operation period, the development of the permafrost embankment thaw settlement is mainly due to the dissipation of the newly postthaw pore water at the thaw depth or permafrost table. This is one of the main differences between the law of permafrost embankment thaw settlement between unfrozen soil embankments.

#### 5. Conclusions

In this paper, a nonlinear large strain thaw consolidation theory was used to study the law of ice-rich permafrost embankments. A segment interpolation function and the corresponding simulation strategy were used to avoid efficiency and stability problems when a nonlinear stress-strain relationship is employed. Based on the measured data of the Qinghai-Tibet Highway, the applicability of the numerical simulation method and the nonlinear theory was verified, and the following conclusions were obtained:

The accuracy of the nonlinear theory after being implemented with the proposed numerical simulation method is obviously higher than that of the linear thaw consolidation theory, which can be applied to calculating the thaw settlement of frozen soil embankments with high ice content.

By analyzing the calculation results of the thaw depth and settlement with both linear and nonlinear theories, it is shown that the calculation accuracies of both the thaw settlement and depth are influenced by the stress-strain relationships due to the interaction between the embankment thermal and mechanical fields.

It is shown that the consolidation degree of soil will be underestimated when the linear stress-strain relationship is used to calculate the thaw settlement of the embankment. Therefore, the thaw consolidation theory based on the nonlinear stress-strain relationship should be employed in the stability design of frozen soil embankments with high ice content to improve the calculation accuracy of the relevant design parameters.

The postthaw pore water in the shallow embankment dissipates in the early operation period. In the following long operation period, the development of the permafrost embankment thaw settlement is mainly due to the dissipation of the new postthaw pore water at the thaw depth or permafrost table.

#### Data Availability

All data, models, and code generated or used during the study appear in the submitted article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors are grateful to the Research Fund of the State Key Laboratory of Eco-hydraulics in the Northwest Arid Region, Xiʼan University of Technology (Grant No. QNZX-2019-07), the National Natural Science Foundation of China (U2034204), and the Fundamental Research Funds for Beijing Universities (X18008).