Abstract
The in situ frost-heaving (FH) process and characteristics of the shallow layer of a residual soil slope in the intermittently frozen zone were simulated by a modified three-dimensional particle flow code (PFC3D) model, of which the mesoscopic parameters of soil and ice particles were calibrated through the indoor experiments. In this model, the in situ FH process was gradually achieved by expanding the volume of ice particles (divided into 24 times of expansion), and the process was terminated when the monitored porosity was stable. These countermeasures avoided the stress accumulation and effectively realized the simulation of the in situ FH process. The results found that the displacement occurred firstly and got the largest final value at the surface angle (SA) under the in situ FH effect, followed by that at the foot, and it gradually extended to the interior based on these two regions. The vertical tension was present at the SA, and the major force type in the lateral interlayer was pressure. In addition, the FH effect seemed to be strongly related to the frozen depth, and a sliding surface was found in a steeper slope. Finally, the smaller stone appeared to be favorable to the slope stability, but it was reduced by the larger stone to some extent.
1. Introduction
Frozen regions are usually classified into three categories: permafrost, seasonally, and intermittently frozen regions [1]. The intermittently frozen region is defined as an area where the near-surface soil is frozen annually from one to fifteen days per year [1, 2]. This kind of area is widely distributed in the world, accounting for about 6.6% of the exposed lands in the coldest month of the year in the Northern Hemisphere [1, 3, 4].
Taking China, for example, the intermittently frozen region is mostly distributed in the south of the Yangtze River and the north of the Pearl River Basin, covering an area of about 1.9 million square kilometers and accounting for about 20% of the total national area [5]. It should be noted that the short-term frozen behaviors of this area mainly resulted from the Siberian cold current in the cold seasons, rather than the latitude-induced factors in the other two frozen regions. In addition, accompanied by the cold current, convective rainfall will occur, leading to significantly increased soil water content, and at the same time, the soil will be frozen if the duration of the cold air is long enough. Moreover, due to the uncertain characteristics of the cold current, the temperature will change unpredictably.
Therefore, compared with the soil in the permafrost and seasonally frozen regions, there are five unique characteristics of the intermittently frozen soil [6]: (1) the total formation period of frozen soil is short, commonly no more than 3 months annually, and is acutely affected by local short-term climate change. (2) The frequency of the freeze-thaw cycles is high, which resulted from the frequently fluctuated ambient temperature (within ±10°C) around the freezing point [6–9]. It should be noted that Que et al. [5] observed 16 freeze-thaw times a day, 1 in two days, and 1 in six days in December 2013 to mid-February 2014 in Wuyi Mountain, Fujian Province, China. (3) The soil moisture content is relatively high, usually above 15% based on the field testing. (4) The depth of the frozen soil layer is shallow, as is due to the short duration of the low temperature. (5) The in situ frost-heaving (FH) effects are much more prominent than that induced by water migration, which is mainly ascribed to the high-frequency FH process, making it too short for water migration to occur.
Hence, widely spread shallow surface spalling of soil slope (shown in Figure 1, taken at a slope in Naping, Fujian, China) was promoted under these prominent climatic characteristics, which would further develop into the instability of slopes, causing potential safety issues to the vehicles and pedestrians along with the projects [5,10]. Therefore, there is a great need to carry out in-depth research on the in situ FH behavior of soil slope in the intermittently frozen region.
(a)
(b)
As reported by [11, 12], the instability of the shallow soil slope in winter was closely related to the short-term repeated freezing and thawing. However, the underlying mechanism has not been thoroughly studied, especially for the soil behaviors under such high-frequency FH processes in a short time frozen area.
Discrete element method, like particle flow code (PFC), has been rapidly developed in recent decades and was found to be a useful software to investigate the geological problems [13–17]. Recently, some research has been done in the simulation of soil freezing by PFC, which was reviewed below.
Zhou and Lai [18] analyzed the mechanical properties of frozen sand soil by two-dimensional PFC (PFC2D), where the cementation of ice in soil was simulated by the contact bonding model. Yin et al. [19] considered the cementation of ice occurring in a limited range through three-dimensional PFC (PFC3D) when investigating the cementation behavior of frozen loess under different confining pressures. Huang [20] studied the microscopic characteristics of frozen soil of cement-modified silty clay under uniaxial compression by two-dimensional PFC. In addition, Bao [21] and Wang [22] had also attempted to simulate the FH process by PFC, but some issues were not fully considered. Specifically, the number of ice particles generated was extremely less, and the stress releasing process was somewhat inconsistent with facts, which resulted in the inaccurate modeling of the whole FH process. The prominent problem of the latter one lies in the assumption of the ice encapsulated rock particles, where the expansion of ice particles was replaced by increasing the rock size. Although this method avoided the complex contact problems, it neglected the skeleton effect, resulting in the inconsistency of the actual FH effect.
As aforementioned, on the one hand, it showed that PFC could be an efficient software for simulating the in situ FH behavior of frozen soil in some way [23–25]. However, on the other hand, due to the defects of the assumptions proposed in the former models, the accuracy of mechanical behaviors of the in situ FH process needs to be further discussed. Technically, with the increase of ice particle volume instantaneously, the considerable instantaneous stress would be initiated due to the suddenly increased overlap in the model, which is likely to cause the collapse of the whole soil structure in the PFC modeling. Therefore, to achieve a more realistic simulation of the in situ FH process through instantaneously enlarging the radius of ice particles, there are the following three specific difficulties to be solved urgently: (1) the selection of parameters such as ice particle strength and contact, (2) the release of instantaneous stress caused by frozen expansion, and (3) the boundary stopping condition.
In view of this, the present paper addresses two questions: how to realize the in situ FH simulation process by PFC3D? What are the in situ FH characteristics of frozen soil slope? The objectives of the present article are threefold: (1) to calibrate the mesoparameters of models based on indoor FH tests, (2) to improve the modeling process by optimizing the approaches of stress releasing and the terminal conditions, and (3) to investigate the mesoscopic in situ FH characteristics under various influencing factors based on the proposed models.
2. Methodology
2.1. Experimental Design
The experimental design is illustrated in Figure 2. First, in order to acquire the FH expansion rate, an indoor FH experiment will be conducted. And then, the numerical model is going to be built and described in detail, including the geometry and contact modes, improved algorithm, simulation process, and layout of the measuring balls. By adjusting the simulated results to that of the indoor FH test, the calibrated mesoscopic parameters of the simulation model were obtained. Finally, in order to investigate the influence of the in situ FH effect on the slope deformation, various factors, including FH depth, slope gradient, and block rocks, will be discussed.
2.2. Indoor Test
2.2.1. Equipment
As shown in Figure 3, the main instruments were an automatic low-temperature freeze-thaw machine (CLD-1, Tianjin Gangyuan test instrument factory), a dial gauge (Harbin Measuring & Cutting Tool Group Co. Ltd.), electric demolding machine (DIM-150, Nanjing Soil Instrument Co., Ltd.), and a sample mold (∅39.1 × 120.0 mm, self-made).
2.2.2. Sample Preparation
The dry density and water content of the soil sample were 1.49 g/cm3 and 21.5%, respectively. The gradation of soil sample was shown in Table 1. Based on the JTG E40-2007 [26], the sample was prepared by the pressure sampling method by indoor tests, and the final size was ∅39.1 × 80.0 mm.
2.2.3. Frost-Heaving Testing Procedures
First, the soil sample prepared was placed in the mold, and the equipment was installed as shown in Figure 3. It is noted that, in order to reduce the boundary friction, Vaseline was applied evenly around the interior of the mold, and at the same time, plastic film was wrapped on the specimen. Then, the specimen was put under the cooling bath (1°C) for 10 h to stabilize the displacement and the temperature. Finally, the frost-heaving test was conducted under a cooling rate of 2°C and maintained at −7°C till the end of the experiment. The displacement was recorded for each hour.
It is noted that, as found in the former research [6], the freezing temperature and time of the experimental soil with a water content of 21.5% were −1.11°C and 4.5 h, respectively. Therefore, to ensure the frost effect, the minimum testing temperature and the lasting time were set as −7°C and 12 h, respectively.
2.2.4. Testing Results
Figure 4 shows the relationship between FH displacement and cooling time. It can be seen that, with increasing cooling time, the displacement increases gradually and freezes completely in about 9 hours. By calculation, the expansion rate was 2.75%.
2.3. Calibration of Material Mesoscopic Parameters
2.3.1. Fundamental Assumptions for PFC Simulation
Frost-heaving deformation happens when water, in and outside the soil, migrates and freezes, which results in space-enlarging and redistribution of soil particles. However, due to the complexity of the FH phenomenon and the high-frequency FH cycles in intermittently frozen area, the following assumptions were made in simulation:(1)The frozen water in the soil sample was treated as rigid ice particles, and their volume was considered as a part of the soil pore.(2)The volume expansion of ice particles in soil was realized by gradually increasing their radius.(3)The in situ FH was only considered in this research.
2.3.2. Geometry, Preparation of the Specimen, and Contact Model
(1) Geometry. The radius of soil particles ranged from 0.65 to 1.3 mm. The size of the numerical model was ∅29.25 × 60.36 mm. To be noted, the lateral model boundary was a rigid column wall, with enough height for volume expansion. Before the frost heaving, the upper and bottom walls were set in order to generate the soil and ice particle. However, after the stress balance, the upper wall was deleted, and then the model was ready for simulation.
(2) Preparation of Specimen. There are two things that should be considered in specimen generation. The first one was the porosity of the specimen. The initial porosity and the mass water content of the soil were 0.4322 and 21.5%, respectively. Because we assume that the ice particles were rigid, the actual porosity in the simulation was 0.1119 by conversion (the dry density of soil was 1.49 g/cm3; 0.4322–21.5% × 1.49 = 0.1119), which was too small, leaving the model hard to stabilize resulting from the significant contact force inside the model. Hence, we generated a model with a porosity of 0.4322 (except for the ice particle) and a mass water content of 9.5%, and the actual porosity was 0.2907 (including ice particles, 0.4322–9.5% × 1.49 = 0.2907).
The second consideration was the radius of the ice particle. Restricted by the computing capacity, the simulation efficiency should be balanced. That is, when the particle radius is too small (about 0.3 mm) and the number is too large (about 40,000), the computation time will be enlarged significantly. On the contrary, the accuracy will not be satisfied. To balance the time and accuracy and based on the corresponding water content, the radius of ice particles was determined as 0.34 mm (close to the result of Bao, [21]); the particle numbers of ice and soil are 5678 and 25,594, respectively.
(3) Contact Model. There are three kinds of contacts among soil samples, named soil-soil, soil-ice, and ice-ice, shown in Figure 5. The parallel contact model was utilized to simulate the cementation induced by ice [20–22].
2.3.3. Algorithm Modification
As mentioned in the introduction part, some researchers have tried to realize the FH process by PFC. However, there were still some limitations that made the results seem to be not satisfactory. In general, two problems should be considered in the FH process modeling by enlarging the radius of ice particles: the generation and release of instantaneous stress caused by expansion and the terminal condition of the calculation process. Based on in situ FH simulation experiences through the PFC, major modifications in this paper are as follows.
Firstly, simulate the in situ FH process by expanding ice particles for 24 times, which helps to reduce the instantaneous stress due to the small overlap in each expansion. Secondly, reduce the kinetic energy by every 50 cycles after each expansion. Finally, stop the operation processes when the monitored porosity of each layer is slightly changed and stable. The following shows details of each step.
(1) Determining the Expansion Coefficient of Radius. The overall radius expansion coefficient, λ, of ice particles is defined as follows:where and are the final and initial radii of ice particles after and before expansion, respectively. Due to the fact that the radius expansion at one time (i.e., enlarging the radius of the ice particles from initial to final values in one expansion cycle) would probably result in excessive stress, leading to slope instability and simulation failure, the in situ FH process is realized by multiple expansions. The expansion coefficient in each time is determined bywhere N is the expansion times, taken as 24 in this research. The expansion coefficient was addressed by adjusting the simulation according to indoor testing results, as described in detail in the following part.
The frost-heaving was achieved by expanding the radii of the ice particles, as was calculated bywhere V1 and V2 are the volumes of ice particles before and after expansion; d and h are the radius and height of the model; and ρd are the water content and the dry density of soil. Here, they are 9.5% and 1.49 g/cm3, respectively. In addition, η (2.75%, shown before) and m present the expansion rate of the indoor soil specimen and the expansion coefficient of the radius of the simulated ice particle. Restricted by the computation capacity of the computer, the 21.5% water content was not achieved. Instead, it was replaced by 9.5%. However, the expansion effect was still maintained as that of the original one, of which the expansion rate was 2.75%. When adjusting the expansion rate of the simulating model to the value, the expansion coefficient of the ice particle radius was 1.06.
(2) Determining the Mesoscopic Parameters. The mesoscopic parameters were obtained by adjusting the simulation results to the indoor FH test (mentioned in 2.2.4). Through many attempts, when sample expansion rate and radius expansion coefficient of ice particles reach 2.75% and 1.058 (addressed previously), respectively, mesoscopic parameters of soil particle, ice particle, and contact model are obtained, shown in Table 2.
2.4. Numerical Model
2.4.1. Geometry
The detailed parameters of the geometry model are shown in Table 3.
2.4.2. Layered Slope and Frost-Heaving
(1)As shown in Figure 6(a), the balanced model was divided into four layers from top to bottom, with the thickness of 10%, 10%, 20%, and 60% of the model height (29.25 mm, shown in Table 3) and named T1, T2, T3, and inner, respectively. At the same time, slope layers were made parallel to the slope surface, labeled, respectively, Q1, Q2, and Q3, and the ice and soil particles in each layer were identified and marked with different colors.(2)In the actual frost-heaving process, the shallow surface of the slope would cause first frost heave, and then with the continuation of low temperature, the depth of frost heaving deepens. For simplifying the process, a layered FH method was utilized. For example, the second layer started to expand after the first layer ended the 8th expansion, and the third and fourth layers did in the same way, as shown in Figure 6(b). Taking four-time expansion of ice particle as one freezing cycle, 10%, 20%, 40%, and 100% FH height need 6, 8, 10, and 12 freezing cycles, respectively.
2.4.3. Layout of Measurement Spheres
In order to obtain stress and displacement values of the surface and top, a series of measurement balls with a radius of 0.003 mm, as shown in Figure 7, are placed in the center of the shallow layer and numbered separately.
(a)
(b)
3. Results and Discussion
3.1. FH Behavior of Shallow Surface of Soil Slopes
3.1.1. Displacement Field
Figures 8 and 9 show the horizontal and vertical displacement fields at different freezing times, respectively. It can be seen that (1) before the second freezing period, there is no obvious displacement, which was mainly due to the fact that the initial porosity was filled by the particle volume expansion. (2) The displacement occurs firstly at the surface angle (SA), followed by the foot. And then, it gradually extends to the interior of the slope taking these two positions as the center. (3) The SA is the most influenced region under the FH effect, both in horizontal and vertical displacement.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.1.2. Contacts Analysis
Figure 10 shows the distribution of the contact force chain, pressure, and tension of slope under 20% frozen depth. It is noted that the force chain shown in black and red segments presents, respectively, compression and tension stress, and the larger the segment, the stronger the stress. It can be found from the figure that (1) the contact force in the frozen area is obviously the largest. (2) Pressure takes up the major proportion of force in the model under the FH and gravitation effects. (3) The tension in the slope surface signifies that the particles are suffering expansive force and presents a trend of outward movement.
(a)
(b)
(c)
3.1.3. Local Stress Analysis
The horizontal and vertical stresses of the shallow layer were detected by numbers 4, 5, 6, 7, 8, and 9 measurement balls and plotted in Figure 11. Some conclusions can be drawn from this figure. Firstly, it presents pressure in the top layer, and the force increases gradually during the FH process developing. Secondly, stress declines as the position are closer to the SA, which was mainly due to the smaller constraints when closing to the SA. Therefore, stress can be released easily. This explained why, as shown above, the particles in the SA experienced the obvious displacement firstly and the largest final displacement. Moreover, the stress level in the vertical direction is one order of magnitude smaller than that of the horizontal, which indicates that the interlayer stress is the major force type under the FH effect. Finally, the soil on the slope surface is expanded to some extent under the FH effect, which can be seen from the tension in Figure 11.
(a)
(b)
3.2. Influence of Different Parameters on Slope Deformation
3.2.1. FH Depth
(1) Displacement Field. The FH depth was set at 10%, 20%, 40%, and 100% of the slope, respectively, to investigate its effect on the displacement field. As shown in Figure 12, under the FH effect, there is no obvious deformation and shape-changing. But the relatively larger displacement all exhibits on the SA and at the foot. The average slope displacement of various layers along the slope surface under different frozen depths is shown in Table 4. As noted, the displacement increases gradually with the increasing frozen depth but decreases with the layer deepening. The results show that the FH effect mainly influenced the surface layer and had little effect on the slope body.
(a)
(b)
(c)
(d)
(2) Contacts Analysis. Figure 13 shows the force contact chains of the slope. It can be found that the FH effects on the slopes with different depths were similar. But slopes with a larger frozen depth enhanced their impact on the unfrozen region. As noted in Table 5, the contact number grows with the FH deepening. However, the ratio of tension to pressure declines to some extent. This is mainly due to the fact that the gap between particles was shortened or even became overlapped during the FH process, leading to the contact force change from tension to pressure.
(a)
(b)
3.2.2. Grade of Side Slope
Compared with the 45° slope under 20% frozen depth, the concentrated displacement area for a steeper slope (60°, shown in Figure 14) only exists in the SA, but the vertical displacement behavior was similar. For the 60° slope, the displacement along the slope surface appears firstly, followed by the influencing area toward the interior and top surface, and finally, a sliding surface region generated. In addition, from Table 6, it is found that the average displacement of the 60° slope is much larger.
(a)
(b)
(c)
(d)
3.2.3. Blocks in Shallow Layer
(1) Generation and Layout of Blocks. Rocks laid on or embedded in the soil slope surface layer are commonly found, which usually leads to rockfall, slide, and other unfavorable disasters in the freezing-thawing season. In the following program, Clump was utilized to simulate blocks in slopes. Clump is one of the primary geometry elements in PFC3D and is a rigid collection of several rigid spherical pebbles. Clumps can translate and rotate, and their motion obeys the equations of motion [27]. Clump is widely used to simulate rigid objects in geomechanical simulation, such as rocks or gravels. [15, 28].
Here, 20% height of the slope was assumed to be frozen, and the consequences of the blocks on slopes under the FH effect were analyzed. The detailed layout is shown in Figure 15.
(a)
(b)
(c)
(2) Displacement Field. Figure 16 presents the displacement field of slopes under various conditions. At the same time, the average displacement of the soil particles in the frozen area and the blocks and the contact force between particle and blocks are shown in Table 7. As noted, (1) the blocks show a slight influence on the displacement behavior of the shallow layer but vary in the interior (Figure 16). (2) The existence of the small block in the slope enhanced the influencing area (shown in warmer color in Figure 16), of which the effect is similar to that being embedded in the slope (Figure 16(b)). But the large and middle blocks (Figures 16(c), 16(d), and 16(f)) have less impact on the slope compared to the smaller ones. (3) The influencing area was remarkably decreased by the large block both in 45° and 60° slopes (Figures 16(d) and 16(f)). Not only its displacement but also that of particles nearby is obviously lower than before. It signifies that large blocks in some way decreased the displacement of the slopes, which can be mainly due to the fact that the buried depth of the massive rock was slightly deeper than the depth of the freezing line, and its bond and friction were also greater for its larger surface area and mass. However, the smaller blocks were relatively easily affected by the FH force compared with massive rocks, thereby causing the movement of the nearby particles. This indicates that the FH effect was more harmful to the slope with smaller rocks but can be reduced by the large rock to some extent. (4) As shown in Table 6, the average contact force between particles and rocks decreases with the increasing volume of rocks, but it was greater than that of the interparticles, indicating the concentration of stress around the rocks.
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusion
Based on the fact that the improved PFC3D modeling, observation, and discussion were carried out on the in situ FH characteristics of soil slope, the following shows the major findings:(1)For the developed PFC model, the instantaneous stress induced by the in situ FH effect was reduced by 24-expansion cycle times. In addition, the kinetic energy was released at the early stage to stable the models. Moreover, the computing terminal conditions could help the model remain stable after each expansion. These countermeasures avoided stress accumulation and effectively realized the simulation of the FH process.(2)The FH process has a significant effect on the slope. Firstly, the displacement occurred firstly and got the largest final displacement at the SA, followed by that of the foot. The displacement gradually extended to the interior of the slope, taking the aforementioned two positions as the center. The SA was the most influenced region, in both horizontal and vertical displacements. Moreover, the vertical tension only occurred in the SA region. Finally, the FH effect had a strong relationship to the frozen depth.(3)The analysis of various factors on the slope showed that the displacement towards the opposite direction of the slope increased with the deepening frozen depth. The displacement influencing area of a steeper slope extended from the SA and eventually formed a sliding surface. In addition, the FH effects caused damage to the slope with smaller stones and were reduced by the larger stone. In the meantime, it led to stress concentration between the rock and soil particles.
Data Availability
The data used to support the findings of this study are included in the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project was financially supported by the National Natural Science Foundation of China (41302215 and 41772297) and the National Natural Science Foundation of Fujian Province, China (2018J01771). In addition, the teammates who make considerable efforts in various ways are also acknowledged, especially Rui Huang, Jianbo Yi, Mingshuai Sun, Cong Shi, and others.