Damage Statistical Empirical Model for Fractured Rock under Freezing-Thawing Cycle and Loading

The natural rock mass prevailingly exists in the form of a fractured rock mass, and freezing-thawing failure of the fractured rock mass is also frequently encountered during geotechnical projects in cold regions. The previous researches and reports in freezing-thawing field principally focused on intact rocks, while rock joints and fractures were rarely considered, which causes great inconvenience to the safety design and stability assessment of engineering. In response to the special climatic conditions of cold regions, the freezing-thawing damage and degradation mechanism of fractured rock were studied in this paper based on existing laboratory experiments and damage mechanics theory. Primarily, a brief review of the progressive damage process of rock in the conventional triaxial compression experiment was given, as well as the determination methods of four characteristic stresses in the prepeak curve. Then, from the microcosmic perspective, the maximum tensile strain yield criterion was used to reflect the microunit strength which was assumed to statistically satisfy the Weibull distribution, deriving the damage evolution equation of fractured rock under the freezing-thawing cycle and load conditions and quantificationally describing the damage evolution law. Consequently, the statistical empirical constitutive relation of fractured rock considering freezing-thawing and loading damages was established. Ultimately, by combining the existing conventional triaxial compression experimental data of freezing-thawing single fractured rocks with the determination methods of characteristic stresses, the relevant constitutive parameters were solved, and the theoretical constitutive relation curves of the fractured rock after freezing-thawing cycles were obtained, which were compared with the experimental results to verify the validity of the established empirical constitutive relation. The study findings can provide a theoretical basis for revealing the freezing-thawing failure mechanism of the fractured rock mass to some extent.


Introduction
China is one of the countries with the largest cold region distribution, accounting for about 75% of the total land area of the country, where an increasing number of rock mass projects have been launched in cold regions driven by national policies in recent years [1][2][3][4]. In cold regions, the macroscopic damage, frost heave failure, and instability caused by the freezing-thawing cycle of fissure water under alternating temperature changes are considered as the main weathering process of rock mass [5][6][7], which has a decisive influence on the stability of rock mass [8][9][10][11]. Therefore, the researches on the mechanism of damage and degradation of rock mass under the freezing-thawing cycle condition are of great significance for engineering construction in cold regions [12][13][14], which have attracted the attention of numerous experts and scholars [15][16][17]. For example, Mutlutürk et al. [18] proposed an attenuation function model to describe the integrity loss of the freezing-thawing rock materials. Tan et al. [19] conducted research on the strength change of granite after freezing-thawing cycles. Besides, Wang et al. [20], Liu et al. [21], Luo et al. [22], and Zhou et al. [23] have studied the dynamic mechanical responses of freezing-thawing rock from different aspects, and the achievements on other rock properties have also been made [24,25]. In respect of the constitutive model of rock, Zhang established the freezing-thawing and loading coupling damage evolution equation and damage constitutive model based on the failure and deformation characteristics of red sandstone [26,27]. Fang et al. [28] proposed a method to solve the constitutive model parameters of random freeze-thaw cycles. According to the internal state variables theory, Wang et al. [29] established a general thermomechanical water migration-coupled plastic constitutive model of rock subjected to freezing-thawing.
These researches have detailedly discussed the damage theory, microstructure, and macromechanical properties of freezing-thawing rocks, as well as the influencing factors, through theoretical derivation and laboratory experiments, which undoubtedly deepen the recognitions and understandings of the freezing-thawing mechanism of rock. Unfortunately, all of those were carried out on intact rocks, with less relevance with fractured rocks, which is the majority of natural rock mass and widely distributed in cold regions. Additionally, the current focus of freezing-thawing fractured rock mostly locates in the crack extension mechanism [30][31][32][33][34], while the freezing-thawing constitutive relation is seldom investigated.
The task of this study is to establish the empirical freezing-thawing constitutive relation of fractured rock from the perspective of damage statistics. In this paper, the progressive failure process of rock in the conventional triaxial compression experiment was firstly reviewed, and the determination methods of four characteristic stress values were listed, which can be applied to accurately calculate the elastic modulus of rock. Secondly, based on the damage statistics theory and Lemaitre strain equivalence principle, the evolution equation and empirical constitutive relation of loading and freezing-thawing damages for fractured rock were deduced by the macroscopic phenomenological damage mechanics method, as well as the constitutive parameter expressions. Finally, on the basis of the existing conventional triaxial experimental results of single fractured rocks under freezing-thawing cycle and loading conditions, the damage statistical empirical constitutive curves of fractured rocks after a specific freezing-thawing cycle were obtained, which presents a great consistency with experi-mental results, verifying the validity of the established empirical constitutive relation.

Progressive Damage Process of Rock in Conventional Triaxial Compression Experiments
During the conventional triaxial compression experiments, the damage process of rock is usually accompanied by crack closure, initiation, extension, and transfixion [35][36][37][38] (see Figure 1), and the development of internal cracks is closely associated with the mechanical properties of the rock itself, which can be divided into the following stages: (1) Stage I: Closure stage-the crack closure stress σ cc The initial defects inside the rock continuously shrink under the action of load, and the axial stress-strain curve (see the blue curve in Figure 1) presents the shape of "convex downward" and "concave upward." The corresponding characteristic stress at the end of this stage is termed as the closure stress σ cc , which is generally determined by the crack volumetric strain method [39]. With respect to this method, however, Peng et al. [35] pointed out that it was strongly dependent on the accurate estimation of elastic parameters (elastic modulus E and Poisson's ratio μ). Figure 1 makes it clear that E is closely related to σ cc , resulting in certain difficulties in the application of volumetric strain method. In this case, the axial strain response method was proposed by considering the nonlinearity of initial stage in axial stress-strain curve, which can not only avoid the dependence on elastic parameters but also eliminate the influence of subjective factors.
(2) Stage II: Linear elasticity stage-the crack initiation stress σ ci Further closure of the initial cracks emerges in this stage, but there is no crack extension or new crack development. The corresponding axial stress-strain curve is shown as an 2 Geofluids inclined straight line, whose slope is defined as the elastic modulus E (see Eq. (1)), and the terminal stress is known as the crack initiation stress σ ci . In most cases, σ ci is considered as approximately 0.4~0.5 times to the peak stress [40].
where ε ci is the crack initiation strain and ε cc is the crack closure strain.
(3) Stage III: Stable extension stage of cracks-the crack damage stress σ cd The initial cracks begin to extend, new cracks gradually develop, and the mechanical performances of rock continuously degrade. On the other hand, obvious deviation from the straight line of the axial stress-strain curve can be observed, and an inflection point appears on the volumetric strain curve (see the red curve in Figure 1), which means the occurrence of rock damage. Therefore, the stress corre-sponding to the inflection point is regarded as the damage stress σ cd .
(4) Stage IV: Unstable extension stage of cracks-the peak stress σ p The microcracks are interconnected to form a network, and the macrofractures begin to appear and interpenetrate, accompanied by energy release. At the end of this stage, both the rock bearing capacity and axial stress-strain curve reach the maximum, i.e., peak stress σ p .

Damage Statistical Empirical Constitutive
Relation of Fractured Rock under Freezing-Thawing Cycle and Loading

Damage Evolution Equation of Intact Rock under
Loading. From the microscopic view, the failure of the rock microunit shows a certain randomness, which can be described by a statistical method [41,42]. Assuming the   Figure 3: Experimental rock specimens [46].         [43], the probability density function Pð εÞ is expressed as follows: where ε is the microstrain and m and ε m are both Weibull parameters, wherein m reflects the brittle properties of rock and ε m represents the macroscopic mean strength of the rock. During loading, the number of damaged microunit N d can be calculated by Eq. (3).
where N is the total number of microunit. If making the loading damage variable of rock D L defined as the ratio of the damaged microunit number to the total microunit number (D L = N d /N), the rock damage evolution equation can be gained as follows:

Damage Evolution Equation of Intact Rock under
Freezing-Thawing Cycle and Loading. In the process of the freezing-thawing cycle, due to the repeated occurrence of water-ice phase transformation and the uneven shrinkage and expansion of mineral particles, cracks rapidly develop and extend, which affects the microstructure of rock and is macroscopically manifested as the degradation response of mechanical properties. Such a progressive failure caused by the accumulation of microdamage can be represented by the freezing-thawing damage D F-T . As previously described, rock damage occurs once the applied load is greater than the value of σ cd . Hence, according to the macroscopic phe-nomenological damage mechanics method, the elastic modulus of rock after freezing-thawing cycles E n can be used as the reference of damage measurement to describe D F-T [26] (see Eq. (5)).
where E 0 is the elastic modulus of fresh rock. From the definitions of damage variables (see Eqs. (4) and (5)), the mechanical properties of rocks are degraded by the freezing-thawing cycle and load in different approaches and resulting in a coupling effect, in which case the coupling damage D can be derived as shown in Eq. (6) underpinned by the generalized damage variable [44].
Substituting Eqs. (4) and (5) into Eq. (6), D can be recast as: Obviously, D = D L before the freezing-thawing cycle treatment to rocksðE n = E 0 Þ ; D = D F−T while no load is applied (ε = 0). According to the Lemaitre strainequivalence principle [45], the constitutive relation of rock damage under conventional triaxial compression condition is: Combining Eq. (7) and Eq. (8), the empirical constitutive relation of rock damage under the freezing-thawing cycle and loading conditions can be obtained: where σ 1 and σ 3 are the maximum and minimum principal stresses, respectively; ε 1 is the maximum principal strain. At the peak of the stress-strain curve, the derivative at the peak equals to 0 (see Eq. (10)), and the expressions of m and ε m can be derived (see Eqs. (11) and (12)). Figure 2 demonstrates that both m and ε m have significant influences on the shape of the theoretical stress-strain curve. However, it is interesting to note that no matter how m or ε m changes, the linear part of curves all coincide, which indicates that the influence of the two is mainly reflected in the other stages except the linear elastic stage. Besides, the

Constitutive Relation Establishment of Fractured Rock under Freezing-Thawing Cycle and Loading.
According to the definition of the coupling damage variable, D = 0 can be deduced for the initial intact rock, which means no macroor microdefects inside the rock. If regarding the macroscopic fractures as a kind of damage, the fractured rocks can be considered as rocks with certain initial damage D 0 , which can also be measured by the elasticity modulus of fractured rocks E C and intact rocks E 0 (see Eq. (13)).
As previously described, the Weibull parameters m and ε m are in close connection with the shape of the theoretical stress-strain curve (see Figure 2), and their values are easily influenced by the elasticity modulus of rocks (see Eqs. (11) and (12)). Therefore, the influence of macroscopic fractures on the curve can be considered by introducing Eqs. (14) and (15).
Thus, Eq. (7) can be rewritten as: The expression of the corresponding empirical constitutive relation, as well as m 0 and ε 0 , can be recast as: In the case of D 0 = 0 (i.e., E C = E 0 ), Eqs. (16) and (17) will be transformed into Eqs. (7) and (9), which illustrates the great applicability of this empirical constitutive relation.

Verification of Empirical Constitutive Relation
To verify the validity and accuracy of the established empirical constitutive relation, the conventional triaxial compression experimental results under the freezing and thawing cycle and loading conditions of sandstone-like  [46] were used, wherein the data with the confining stresses of 2, 4, and 6 MPa (σ 3 = 2, 4, and 6 MPa) under 0, 20, 40, and 60 freezing and thawing cycles (n = 0, 20, 40, and 60) were utilized to obtain the specific theoretical stress-strain curves. The specimens are shown in Figure 3, and the corresponding physical and mechanical parameters are presented in Table 1. All the fracture length is 30 mm, and the dip angles are 30°, 60°, and 90°, denoted as A4, A5, and A6, respectively. The characteristic stresses of intact rock and different fractured rocks are shown in Table 2, wherein σ cd was gained by volumetric strain method, σ cc was determined by axial strain response method, σ p was directly obtained from the axial stress-strain curve, σ ci approximated to 0.4~0.5σ p , and the elasticity modulus E was calculated by Eq. (1). Taking  (1) For a given axial stress-strain curve (see the red curve in Figure 4), σ cd can be gained by the volumetric strain method, and the corresponding ε cd can be located, which point was denoted as P cd (5.11, 50.91) (2) Connecting the origin O with P cd to obtain the expression of the line: ε = 0:1σ (see the green curve in Figure 4) (3) Solving the strain difference Δε between the red curve and green curve, plotting the Δε~ε curve (see the blue curve in Figure 4), and the stress corresponding to the maximum of Δε is just the crack closure stress   Similar to intact rocks, both the σ cd and σ p increase with the increase of σ 3 , but decrease with the growing of cycle number n, while no obvious variation law of σ cc or σ ci can be observed. By substituting the relevant parameters into Eqs. (18) and (19), the m 0 and ε 0 of each curve can be solved (see Table 3), and therefore, the theoretical curves can be obtained (see . In accordance with Figures 5-7, a comprehensive comparison between the experimental and theoretical curves was made. First of all, a general agreement of the two can be observed, the variation trend is basically consistent, the liner parts of curves agree well or keep parallel, and the peaks also coincide, which prove the validity and rationality of   9 Geofluids established empirical constitutive relation. But it is important to recognize that the linear elastic stage of the theoretical curve emerges at the beginning for this empirical constitutive relation cannot describe the nonlinear characteristic of the closure stage. As a result, the experimental data with no obvious closure stage are in greater agreement with the theoretical curves, while those with obvious closure stage greatly deviate from theoretical curves before peak, such as the experimental data of A6 under the conditions of σ 3 = 2 MPa and N = 60 (see Figure 5(d)). In addition, it is assumed that the rock microunits do not bear loads after damage in this paper, that is, the residual strength is not considered. Therefore, the postpeak stage of theoretical curve differs greatly from experimental results.

Conclusions
Since the existing researches on the constitutive relations of freezing-thawing rock materials have rarely considered rock fractures, the evolution equation and empirical constitutive relation of freezing-thawing and loading damages for fractured rock were derived by using the statistical method and damage mechanics theory. And the following conclusions can be drawn: (1) According to the extension law of rock cracks, the prepeak stress-strain curve can be divided into four stages: closure stage, linear elasticity stage, stable extension stage of cracks, and unstable extension stage of cracks, respectively, corresponding to the four characteristic stresses. The determination method of each was briefly reviewed, which is conducive to accurately calculate the elastic modulus of rocks. Besides, both the crack damage stress and peak stress increase with the increase of confining stress, but decrease with the growing of cycle number, while no obvious variation law of the crack closure stress or crack initiation stress can be concluded (2) By combining the macroscopic phenomenological method with microscopic statistical method, the coupling damage evolution equation of the freezingthawing cycle and load was established. On this basis, the equation parameters were modified by considering the rock fractures, and the empirical constitutive relation of freezing-thawing and loading damages for fractured rock was derived in combination with the Lemaitre strain equivalence hypothesis (3) Underpinned by existing conventional triaxial compression experimental data under freezing and thawing cycle and loading conditions of single fractured sandstone-like specimens, the specific empirical constitutive curves were obtained based on the determination methods of characteristic stresses and damage evolution equation parameter expressions, which were then compared with the experimental results. The theoretical curves are generally consistent with the experimental data, demonstrating the validity of the established empirical constitutive rela-tion, which has important practical guiding significance for the construction of geotechnical engineering in cold region Notation σ cc : The crack closure stress, MPa E: Elastic modulus of rock, MPa μ: Poisson's ratio σ ci : The crack initiation stress, MPa ε ci : The crack initiation strain ε cc : The crack closure strain σ cd : The crack damage stress, MPa σ p : The peak stress, MPa ε: The microstrain of rock m: Weibull parameter ε m : Weibull parameter N: Total number of microunit D L : The rock damage caused by load N d : Damaged microunit number D F-T : The freezing-thawing damage E n : The elastic modulus of freezing-thawing rock after n cycles, MPa E 0 : The elastic modulus of fresh intact rock, MPa D: The coupling damage σ 1 : The maximum principal stress, MPa σ 3 : The minimum principal stress, MPa ε 1 : The maximum principal strain ε p : The peak maximum principal strain D 0 : The initial damage of fractured rock E C : The elastic modulus of fractured rock, MPa m 0 : Weibull parameter ε 0 : Weibull parameter ρ: Density, g/cm 3 σ c : Compressive strength, MPa σ t : Tensile strength, MPa c: Cohesion, MPa φ: Internal friction angle Δε: Strain difference

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

Conflicts of Interest
The authors declare no conflicts of interest.