Frost-Heaving Cracking Sensitivity of Single-Flaw Rock Mass Based on a Numerical Experimental Method

School of Resources and Civil Engineering, Northeastern University, Shenyang 110819, China State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China Hongda Blasting Engineering Group Co., Ltd., Guangzhou 510623, China


Introduction
Rocks are heterogeneous composite materials, which can develop defects, such as cracks, holes, joints, weak surfaces, and interlayers as a result of long-term geological tectonic activity. Infrastructure construction projects at high altitudes or in cold regions often face the problem of freezing damage to fractured rock mass, which is a critical engineering challenge that cannot be ignored. There are many factors affecting the engineering stability of rock masses in cold regions [1][2][3], including microcracks (microscopic defects) and joint surfaces (macroscopic defects) [4,5]. Under high rainfall and subsurface flow conditions, water typically accumulates in fissures. During the cooling process, the freezing of water induces a volume expansion of approximately 9% [6,7]. When cracked ice is restricted by the surrounding rock medium, it generates volumetric dilatancy pressure, also known as, frost-heaving pressure. The frost-heaving pressure tends to concentrate the stress at the crack tip, which is enough to induce the formation of a secondary crack that grows and breaks through; this may cause damage or even the destruction of the rock mass, thus hindering the construction and implementation of rock mass engineering in cold regions [8].
Currently, the frost-heaving and freezing-thawing engineering disasters involving fractured rock masses in cold regions are mainly studied through theoretical calculations and experimental measurements of the frost-heaving pressure, as well as through analysis of the mechanical properties of the rock material after freeze-thawing.
There are two research ideas that are typically applied in theoretical calculations of the frost-heaving pressure. One involves determining the frost-heaving pressure based on the circular tunnel model [9,10]. For example, Lai et al. [11] proposed a viscoelastic analytical computation of tunnel frost-heaving pressure in cold regions by using a Laplace transform method. Gao et al. [12] derived the analytical elastic-plastic solution of the frost-heaving pressure of a tunnel in a cold region using a continuity equation, and they proposed a semianalytical method to calculate the stress in the plastic zone of the surrounding rock. Feng et al. [13] developed a theoretical model to obtain the elastic-plastic solutions from the stress and deformation analysis around the tunnel in cold regions where the entire surrounding rock comprised a nonfrozen elastic zone, a frozen elastic zone, a frozen plastic zone, and a support zone. Zhao et al. [14] established an elastic-plastic analytical determination of the frost-heaving pressure by considering the transversely isotropic frost heave of the surrounding rock. The results indicated that the transversely isotropic frost heave of the surrounding rock had a significant impact on the frostheaving pressure.
The second approach to such theoretical calculations involves computing the frost-heaving pressure based on a single flaw [15][16][17]. For example, Kang et al. [18] proposed a theoretical method to predict crack propagation under freezing pressure and far-field stress, and the authors determined the direction and length of crack propagation based on Griffith's theory of brittle fracture mechanics. Considering the influence of thermal stress on flaw deformation, Huang et al. [19] derived a mathematical model for calculating the frost-heaving pressure of elliptical flaws in lowpermeability rock. Tan et al. [20] expanded on previous research and evaluated the influence of external loading conditions and other parameters.
Aiming to empirically determine the frost-heaving pressure, Davidson and Nye [21] exploited the photoelastic effect to measure the stress generated in transparent materials following ice expansion, and they also discussed the influence of the frost-heaving pressure on the freezing process of real rock materials. Matsuoka and Murton [22] studied the effects of porosity and lithology on the frost-heaving deformation of saturated rock in open and closed systems and discussed the law of flaw propagation during the melting and freezing stages. A new experimental device was designed by Amanuma et al. [23] to simultaneously measure the frostheaving pressure in the temperature gradient direction and perpendicular to the temperature gradient direction. Xia et al. [24] conducted transverse isotropic frost-heaving experiments involving saturated rock under unidirectional freezing conditions and deduced the formula for tunnel frost-heaving pressure. Based on experimental research investigating the frost-heaving pressure characteristics of fractured red sandstone, Shan et al. [25] analyzed the factors affecting the frost-heaving pressure and cracking mode and discussed the frost-heaving expansion mechanism of a single-flaw rock mass.
Several research groups have also studied the mechanical properties of rock materials after freezing-thawing [26,27]. For example, Jia et al. [28] studied the changes in the compressive tensile strength of sandstone based on freezingthawing tests, and they also monitored the strain changes of specimens during a single freeze-thaw cycle to evaluate the damage mechanism. Huang et al. [29] analyzed the strength and failure characteristics of rock materials containing single cracks under freeze-thaw and uniaxial compression conditions, and they concluded that the new frost-heaving crack grew from the crack tip along the coplanar direction because of the boundary effect; then, most of it expanded toward the short side of the specimen. Bai et al. [30] carried out experimental research regarding the strength, deformation, and crack evolution characteristics of red sandstone samples containing double flaws and subjected to triaxial compression conditions at various temperatures and confining pressures.
The mechanical behavior of brittle rock is affected by the initiation, propagation, and coalescence of preexisting rock flaws. The development cracking mode and microcrack location of a secondary crack of fractured rock mass caused by frost-heaving pressure has important guiding significance for construction projects in cold regions. Considering the heterogeneity of rock materials, this study establishes a thermo-mechanical model based on elastic mechanics, thermodynamics, and fracture mechanics and proposes a novel numerical experimental method for evaluating frostheaving cracking via RFPA 2D -thermal analysis. The numerical experimental system reproduces the physical frostheaving cracking process, and the feasibility of this approach is verified (Section 2). Additionally, the failure mode of the single preexisting water-saturated flaw in the rock mass during frost-heaving, and the effects of the inclination angle, flaw width, crack length, and cooling rate are discussed (Sections 3.1-3.4). Ultimately, a prediction regression model of the secondary crack is derived by combining multiple stepwise regressions (Section 4). Finally, the few limitations of this study are discussed, and several significant conclusions are drawn (Section 5).

Verifying the Feasibility of Numerical Simulations of Rock Mass Frost-Heaving Cracking
Based on the theory of elastic damage mechanics, RFPA 2D discretizes rock materials into several elements with distinct physical and mechanical properties, while ensuring that their material properties obey the Weibull distribution [31][32][33]. The RFPA 2D -thermal module was developed by Tang and coworkers, and several reports have been published on this topic [34,35]. This module is mainly employed for the simultaneous analysis of the temperature and stress fields for a brittle material. First, the temperature field in the model is computed according to the temperature conditions. Second, the stress field is deduced according to the mechanical conditions. Then, the selected phase transition criterion is used to evaluate the basic elements, and 2 Geofluids the corresponding damage treatment is carried out. The elements are treated as linearly elastic before and after phase transformation. RFPA 2D -thermal can simulate the initiation and propagation of fractured rock mass cracking and reproduce the failure process. When analyzing the frost-heaving cracking of a fractured rock mass, the coupled temperature and stress effects are mainly considered, and the effect of the water-ice phase transition is generally ignored; rather, it is simplified as the damage to the surrounding rock medium induced by the expansion of cracked ice at low temperature.

Theoretical
Model. The first set of heat conduction temperature boundary conditions is used for numerical simulations. The differential equation of heat conduction is described by where k x and k y represent the thermal conductivity of the elements in the X-and Y-directions, respectively; T is the temperature; ρ and c represent the density and specific heat of the material, respectively; t is the time.
For the frost-heaving analysis of a fractured rock mass, the stress-strain field equation is expressed as shown in where σ ij is the stress tensor; F i is the body force; ε ij and u i represent the strain tensor and displacement vector, respectively; β, G, and λ are the thermal stress coefficient, shear modulus, and Lamé coefficient, respectively; δ ij represents the Kronecker function (i = j, In elastic mechanics, the elastic modulus of an element gradually decreases with damage aggravation, and the elastic modulus of a damaged material is defined by During frost-heaving, fractured rock cracking is mainly affected by tensile stress. When the tensile stress in an element reaches its uniaxial tensile strength σ t0 , the relationship in Equation (6) is satisfied.
In this case, tensile damage occurs in the element, the   3 Geofluids damage variable D is defined by the conditions in where E 0 represents the initial elasticity modulus; σ tr is the residual strength of the tensile damage; ε t0 and ε tu are the tensile strain at the elastic limit and the ultimate tensile strain of the element, respectively. The tensile damage constitutive relationship for the frost-heaving cracking of a fractured rock mass can be expressed as shown in To verify the reliability of the rock frost-heaving analysis using the RFPA 2D -thermal module, the simulated results were compared with the experimental study regarding the frost-heaving pressure characteristics of fissured red sandstone conducted by Shan et al. [25]. The test specimen was a standard cylinder (ϕ = 50 mm; l = 100 mm). The width and length of the flaw in the specimen were 1 mm and 20 mm, respectively ( Figure 1). The frost-heaving analysis of the fractured rock was carried out by directly soaking the sample in water and cooling to -28°C over a total freezing time of 25 h. To maintain consistency with the test conditions, a numerical simulation model was established (Figure 1) based on a plane strain model. The temperature boundary conditions were applied to the boundary around the model, and the temperature drop was set to -28°C. The mechanical parameters and thermal properties of the model are summarized in Tables 1 and 2, respectively.

Cracking Mode and Frost-Heaving Pressure
Comparison. Figure 2 shows the cracking morphology of the fractured rock model as the temperature drops to -28°C. Initially, the preexisting closed and water-saturated single-flaw begins to expand and extrude into the surrounding media without causing any damage. As the frost-heaving pressure gradually increases, the extrusion effect of the surrounding medium becomes clearer. Under the applied crack tip stress concentration, microcracks appeared when the stress exceeded the tensile strength of the rock material. As the number of microcracks increases, individual microcracks gradually converge to form macrocracks. The experimentally measured frost-heaving pressure was 0.23 MPa [25]. In the numerical simulation, the minimum principal stress at the center of the cracked ice was defined as the frost-heaving pressure, and the computed value was 0.45 MPa (Figure 3).
In Figures 2(e) and 2(f), the numerical results are similar to the experimental results. Although there are some differences in the cracking mode (i.e., inclination angle of the secondary crack), the numerical simulation results concerning the frost-heaving pressure closely reproduce the test results. The feasibility of the numerical experimental system is thus verified again.

Sensitivity Analysis of a Preexisting Closed and Water-Saturated Single-Flaw during Frost-Heaving
The single flaw is the simplest and most important form of rock mass composition. The mechanical properties of a single-flaw rock mass also comprise the foundation for studying the mechanical properties of complex fractured rock masses. This section discusses and analyzes the effects of the flaw width, length, inclination angle, and cooling rate on the cracking mode. The feasibility of the numerical experimental system for frost-heaving analysis was verified in Section 2.2; to facilitate the computations, the model and material parameters in Tables 1 and 2 were used, and only the BC edge ( Figure 1) was cooled to -28°C.

Influence of the Inclination Angle on the Frost-Heaving
Cracking Process. Considering the inclination angles of 0°, 30°, 45°, 60°, and 90°and assuming a flaw length and width of 20 mm and 1 mm, respectively, the influence of the inclination angle on the frost-heaving cracking of a preexisting water-saturated single-flaw was analyzed. From the real-time monitoring of the rock frost-heaving process at various inclination angles ( Figure 5), it was concluded that the preexisting single-flaws with different incli-nation angles mainly propagated along the flaw-coplanar direction.
With an inclination angle of 0°, the lengths of the left and right sides were approximately 3.475 mm. When the inclination angle was 30°, 45°, or 60°, the lengths of the left and right sides were about 6.243 mm, 7.218 mm, and 8.851 mm, respectively. When the inclination angle was 90°, the secondary crack was about 9.538 mm long; the secondary crack length increases with increasing inclination angle, mainly because the frost-heaving pressure reaches the strength limit of the rock medium sooner at a greater inclination angle. This causes tensile failure, and failure units gradually form a secondary crack. However, a secondary crack must appear at the tip of a preexisting flaw because of the stress concentrations, regardless of the inclination angle.

Influence of the Flaw Width on the Frost-Heaving
Cracking Process. Setting the flaw length as 20 mm and the inclination angle at 60°, the width of the flaw was varied (i.e., 0.6, 1, 1.6, or 2 mm) to analyze the influence of the flaw width on the frost-heaving cracking of a preexisting closed and water-saturated single-flaw. Figure 6 shows that the cracking modes of the secondary cracks produced by frostheaving in a fractured rock mass with various flaw widths were quite diverse. When the flaw width was less than 1 mm, the initiation crack began from the tip and developed gradually along the direction of the flaw, thus forming a secondary crack. When the flaw width was greater than 1 mm, there were two stress concentration points at the end of the preexisting flaw, and secondary cracks developed from these two stress concentration points. It is clear that the flaw width affects the cracking mode of rock frost-heaving. As the flaw width increased, the larger frost-heaving pressure caused by the cooling effect causes a longer secondary crack to develop. When the flaw width was 0.6, 1, or 1.6 mm, the secondary crack length was about 3.855, 8.851, or 18.591 mm, respectively. When the flaw width was 2 mm, the secondary crack ran through the entire model, and the upper covering layer tends to lift; in this case, the secondary crack length is about 22 mm. Figure 6(a), path A-A' was taken as the research object to analyze the change in the minimum principal stress during the frostheaving cracking process as a function of the flaw width. In the module, pressure is positive, and tension is negative. The plots in Figure 7 reveal that the minimum principal stress of the rock medium along path A-A' exhibits a concave trend in all simulated cases. The reason for this type of stress distribution is that the preexisting flaw was located in the middle of the whole model, and the rock medium above the flaw has the greatest influence on the frostheaving pressure. When the flaw width was 0.6, 1, 1.6, or 2 mm, the minimum principal stress value along path A-A' reached -0.08, -0.16, -0.25, and -0.5 MPa, respectively. The larger the flaw width, the greater the frost-heaving pressure. The stress curves fluctuate because the properties of the 5 Geofluids rock materials are not uniform, and low strength elements are destroyed first. The stress acting on a given element disperses around it, thus causing the stress in adjacent elements to vary greatly.

Influence of the Flaw Length on the Frost-Heaving
Cracking Process. Setting the inclination angle to 60°and the flaw width to 1 mm, the flaw length was varied (10, 15, 20, or 25 mm) to determine the influence of the flaw length on the frost-heaving cracking of preexisting closed and water-saturated single-flaws. Figure 8 shows that as the flaw length increased, the secondary cracks gradually lengthen; however, the main trend of secondary crack development still followed the direction of the preexisting cracks. When the flaw length was 10, 15, 20, or 25 mm, the secondary crack length was about 5.6, 8.474, 8.851, or 11.97 mm, respectively. Figure 9 presents the Y-direction (vertical) displacement of preexisting flaws with various lengths along the B-B' path (Figure 8(a)). For example, when the flaw length was 10 mm, the vertical displacement of the upper part of the rock mass exhibited a convex trend (Figure 9(a)). The peak of this curve tends toward the left side, mostly because the preexisting flaw's inclination angle was 60°, and the rock medium in the upper-left part of the specimens is more significantly affected by the frost-heaving pressure, so the peak value appears to be offset. A comparison of the plots in Figures 9(a)-9(d) indicates that as the length of the preexisting flaw increases, a frost-heaving pressure is generated, and the effect on the upper part of the specimen becomes more obvious. As a result, the deformation increases from 1:82 × 10 −3 mm to 5:02 × 10 −3 mm.

Influence of the Cooling Rate on the Rock Frost-Heaving
Cracking Process. The frost-heaving cracking process for preexisting closed and water-saturated single-flaws was analyzed by setting an inclination angle of 45°, a flaw width of 1 mm, a flaw length of 20 mm, and modulating the cooling rate (-0.5, -1, -2, or -4°C/h). Figure 10 shows that the secondary crack length gradually decreased with faster cooling rates, and the development direction still followed the direction of the flaw. This was because the faster the cooling rate, the more failure units generated around the tip of the preexisting flaw under the frost-heaving pressure, which consumes more energy, thus decreasing the length of the 7 Geofluids secondary crack. When the cooling rate was -0.5, -1, -2, or -4°C/h, the secondary crack length was about 9.758, 7.218, 4.903, or 3.583 mm, respectively.

Analysis of the Microcrack Initiation Process.
The graphs in Figure 11 shows the microcrack changes under various cooling rates during the frost-heaving process. It is clear that as the cooling rate increases, the maximum number of accumulated microcracks gradually increases. When the cooling rate was -0.5, -1, -2, or -4°C/h (each step took one hour), the number of microcracks reached about 20, 35, 80, or 130, respectively.

Multiple Stepwise Regression Analysis of the Secondary Crack Length Sensitivity
The variation in the secondary crack length as influenced by the inclination angle, flaw width, crack length, and cooling rate has been studied previously. Herein, a regression model for predicting the secondary crack length has been established by applying a multiple stepwise regression method.

Multiple Stepwise Regression Analytical Model Setup.
The main steps in the multiple stepwise regression analysis are describe here. First, independent variables are introduced, including the inclination angle, flaw width, flaw length, and cooling rate. The independent variable with the most significant effect on the secondary crack length y is introduced each time. Then, other variables are tested one by one to eliminate the variables associated with negligible changes. The multiple stepwise regression equation is shown in where y is the secondary crack length; x 1 , x 2 , x 3 , x 4 represent the independent variables affecting the secondary crack length (i.e., inclination angle, flaw width, crack length, and cooling rate, respectively); β 0 is a constant; β 1 , β 2 , β 3

Geofluids
are the regression coefficients (indicating the degree of influence of the inclination angle, flaw width, crack length, and cooling rate, respectively, on the secondary crack length); ε is the random error.
The described analytical tests indicate that the regression equation of this model is suitable for use in predictions. In the established model, the flaw width is positively correlated with secondary crack length, with a correlation coefficient of 0.817. There is also a positive correlation between the inclination angle and the secondary crack length (correlation coefficient = 0:294). The cooling rate is negatively correlated with the length of the secondary crack, with a correlation coefficient of 0.249. Finally, the flaw length is positively correlated with the secondary crack length, and the correlation coefficient is 0.207. The relative contributions to the secondary crack length follow the order flaw width > inclination angle > cooling rate > flaw length.
The regression prediction equation for the secondary crack length is shown in y = −14:751 + 13:668x 1 + 0:08x 2 + 1:621x 3 + 0:345x 4 : ð10Þ 4.3. Inspection and Application. Equation (10) can be used to predict the length of the secondary cracks in entries 15-17 in Table 3, and the results are compared with the measured values in Table 7 as programs 1-3, respectively. The obtained relative error is within the allowable range, so the established regression prediction equation demonstrates suitable applicability.

Discussion and Conclusions
The simulations presented herein investigated the effects of the inclination angle, flaw width, flaw length, and cooling rate on crack propagation. The regression equation for predicting the secondary crack length was established via multiple stepwise regression analysis. The results confirmed that the frost-heaving pressure in a flaw was enough to drive crack initiation and propagation along the coplanar direction of a preexisting crack under unconfined conditions. Therefore, the frost-heaving process may easily destroy the unconfined rock mass near the surface when crack water is available. Based on the results of this study, the following conclusions can be drawn: (1) During the freeze-heaving process, preexisting closed and water-saturated single-flaws with various inclination angles, flaw widths, flaw lengths, and cooling rates mainly propagate along the flawcoplanar direction. However, the secondary crack   10 Geofluids must appear at the tip of the preexisting flaw because of the stress concentration (2) The cracking mode of a rock mass is influenced by the flaw width of a preexisting single flaw. As the flaw width increases, the secondary crack changes from one main crack to two main cracks, and the minimum principal stress value along the path A-A' changes more dramatically. As the flaw width increases from 0.6 to 2 mm, the maximum minimum principal stress increases from 0.08 to 0.5 MPa (3) The secondary crack length increases with increasing inclination angle, flaw width, or flaw length, and decreases with increasing cooling rates. The entire process of secondary crack formation can be divided into a rupture incubation period, a rupture growth period, and a rupture steady period (4) The influencing factors of flaw width, inclination angle, cooling rate, and flaw length correspond to 0.817, 0.294, 0.249, and 0.207, respectively, leading to an order of contribution to the secondary crack length where flaw width > inclination angle > cooling rate > flaw length. The model for predicting the secondary crack length is thus established as follows: y = -14:751 + 13:668x 1 + 0:08x 2 -1:621x 3 + 0:345x 4 Herein, simulations of the frost-heaving cracking of fractured rock masses are simplified to a model probing the coupled temperature and stress fields. The described work considered the cracking of the surrounding rock medium caused by the expansion of cracked ice under an ambient temperature of 0°C. This study also focused on elucidating the relationship between secondary crack formation and propagation and the preexisting cracked ice; however, it ignored the effect of the water-ice phase transition. Therefore, ongoing studies in our laboratory are aimed at further enriching the established theoretical model and redeveloping the RFPA software to incorporate the water-ice phase transition process.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work; there is no professional or other personal interest of any nature or kind in any product, service and or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.