Statistical Damage Shear Constitutive Model of Rock Joints Under Seepage Pressure

The constitutive relation of rock joints in the percolation environment has always been a frontier research topic in the field of geotechnical engineering, involving nuclear waste disposal, dam foundation construction, tunnel tunneling, and other fields. In this study, a shear statistical damage constitutive model of rock joints considering seepage pressure was proposed, which takes into account the statistical damage, stage deformation characteristics and seepage action in the process of joint shear deformation. By introducing the seepage pressure reduction coefficient and establishing the relationship between the model parameters and the seepage pressure, the prediction of joint shear stress-displacement curves under different seepage pressures was achieved. The comprehensive comparisons between the theoretical results and the experimental results reveals that the model can accurately reflect the shear deformation law of joints and reasonably simulate the shear stress-displacement relationship of joints. And ultimately, the physical significance of the model parameters and their relationship with shear strength characteristics are discussed via the gray correlation method. The model and verification results presented in this study can deepen the understanding of the shear behavior of joints subject to high seepage pressure and contribute to the design and optimization of geotechnical engineering in the percolation environment.


INTRODUCTION
Rock joints, which are ubiquitous in slope, fracture zone, and underground engineering excavation area, have a significantly weaker bearing capacity than the two wall surface rock masses (Yong et al., 2018;Lin et al., 2019a;Wang et al., 2020;Xie et al., 2020a;Zhang et al., 2020a,b). And typically, shear slip failure along joints has dominated the main failure modes of the aforementioned rock engineering under the action of external forces or drastic changes in the geological environment (Meng et al., 2017;Han et al., 2020;Li et al., 2020;Lin et al., 2020a;Yuan et al., 2020). Besides, as the main movement channel of groundwater, the seepage properties of joints also present a powerful influence on the shear mechanical properties of joints. In this case, the investigation on the strength characteristics and seepage characteristics of rock joints under shear load is of great practical significance for guaranteeing the safety and stability of geotechnical engineering and scientific rationality of design.
The salient issues of shear failure of joint under seepage pressure have globally attracted extensive attention from scientists and engineers. Since the establishment of the classical cubic law (Snow, 1969), numerous scholars have carried out explorations on the influence of stress on seepage characteristics of rock joints, including experimental research and numerical model research on the seepage characteristics of rock joint under normal stress, shear stress and composite stress, and obtained very rich results (Raven and Gale, 1985;Lee and Cho, 2002;Yasuhara et al., 2006;Xia et al., 2017;Zhao et al., 2017a,b). For the seepage characteristics of the fracture surfaces with circular convex contact areas, Walsh (1981) applied the effective stress theory to derive the influence coefficient of the contact area. Zimmerman et al. (1992) simplified the fracture surfaces to oval bulges, derived the correction factor for the contact area effect, and extended it to the analysis of irregular contact areas. Assisted with the COMSOL Multiphysics numerical method, Yang et al. (2015) proposed a finite element mathematical model of rock masses considering permeability and elastic anisotropy. Equipped by self-developed apparatus, Liu and Chen (2011) studied the seepage and deformation properties of rock joints subjected to normal stress, as well as lateral stress. Rong et al. (2017) conducted joint shear-seepage coupling tests to three types of granite joints with different roughness, and pointed out that the discharge was sensitive to the changes of normal stress, joint roughness coefficient (JRC), and pressure gradient.
Although these aforesaid valuable contributions lay a good foundation for further understanding of the deformation characteristics of rock joints under hydro-mechanical (HM) coupling conditions, the researches on the shear stressdisplacement relationship of joints under seepage pressure, that is, the constitutive model, is still rarely reported. Olsson and Barton (2001) proposed an improved empirical engineering model that considers the coupling of hydraulic aperture and shear stress-displacement of joints during the shear process. Under the framework of the cubic law and minimum potential energy principle, Zhao et al. (2016) developed a coupled stress-seepage model to reflect the characteristics of the pre-and post-shear stress-displacement curves. Despite enriching the constitutive modeling of rock joints under seepage pressure, however, the theoretical model that directly include the shear stress, normal stress and seepage pressure remain unavailable.
Therefore, in this study, the modified constitutive model proposed by Lin et al. (2019b), as well as the additional introduction of the seepage pressure reduction coefficient β to quantitatively investigate the reduction effect of seepage pressure, was developed to predict the shear stress-displacement curves of joints under high seepage pressure. Particularly, the model considers the stage characteristics of shear behavior of joints under seepage pressure. All the parameters required by the proposed model can be directly determined from the laboratory direct shear tests. The validity of the proposed constitutive model was verified by experimental results.

Constitutive Damage Model of Joint Shear in Seepage Environment
Based on damage statistics theory and joint shear deformation characteristics, Lin et al. (2019b) proposed a shear constitutive model that reflects the whole joint shear deformation characteristics as follows: where k s is the shear stiffness, u y is the yield shear displacement, τ r is the residual shear stress, and m and u 0 are the Weibull distribution parameters. Moreover, the Weibull distribution parameters m and u 0 can be mathematically solved as: where τ p and u p are the peak shear stress and peak shear displacement, respectively. The shear constitutive model of rock joints is sensitive to permeable engineering environment, such as the interaction and superposition of stress field, seepage field, even temperature field and chemical field (Lei et al., 2019;Zhao et al., 2019;Lin et al., 2020b;Xie et al., 2020b). In order to thoroughly understand the problem of multi-field coupling, the method of the system research is needed, including the internal relations of the subsystems and the interactions among the subsystems. However, due to the complexity of this multi-field coupling, it is difficult to clearly reflect the whole coupling process, and it is quite difficulty to establish an accurate mathematical model up to now.
By establishing the relationship between model parameters (k s ,u y , τ p , u p , u 0 , m) and seepage pressure P w , a statistical damage constitutive model considering seepage pressure can be established. Then the joint shear stress-displacement curve under a certain seepage pressure can be predicted by the existing sample shear stress-displacement curve. Based on the above conjecture, on the basis of Equation 1, the relationship between each parameter in Equation 1 and seepage pressure P w is considered where A(P w ), B(P w ), C(P w ), D(P w ), E(P w ), F(P w ) represent the certain correlations between corresponding model parameters and seepage pressure P w . Therefore, as long as the relationship between model parameters and seepage pressure P w are determined, the statistical constitutive model of joint shear damage considering seepage pressure can be established by substituting Equation 3 into Equation 1:

Experimental Results
In this study, experimental data from literature (Qian, 2018;Qian et al., 2018) were selected as examples to establish a statistical damage constitutive model of rock joints considering seepage pressure. The shear-flow coupling tests were performed on three artificial joint samples with different roughness surfaces. Twenty seven samples with the size of 200 × 100 × 100 mm (length × width × height) were produced with different joint morphologies, divided into 3 groups (denoted as L, M, and H) and 9 samples with identical morphology in each group. The joint roughness coefficient (JRC) corresponding to these 3 groups of samples are : L: JRC = 5.657, M: JRC = 10.793, and H: JRC = 13.897, respectively. Note that, these samples were numbered as the rule of X-RJ N, where X stands for the group name (L, M, or H), and N stands for sample number. The values of seepage pressure (p w ) adopted in the tests were 1, 2, and 3 MPa, and the normal stresses of 5, 8, and 10 MPa were applied, respectively. The uniaxial compressive strength (σ c ) and basic friction angle (ϕ b ) of prefabricated joint samples are 77.14 MPa and 32.48 • , respectively.
All direct shear tests strictly follow the latest revised ISRM recommendations (Muralha et al., 2014). The obtained shear test results are presented in Figures 1-3, and the relevant test data are recorded as shown in Table 1.
According to the test results, the shear stress-displacement curves of joints under seepage pressure basically conform to the typical shear stress-displacement curve (Figure 4): elastic deformation stage (1), yield stage (2), post-peak deformation stage (3), residual deformation stage (4), where τ y and u y represent the yield shear stress and yield shear displacement,   τ p and u p stand for the peak shear stress and the peak shear displacement, τ r and u r are the residual shear stress and the residual shear displacement, respectively. Practically, not all characteristics of four stage mentioned above are available in experimental curves due to objective attributes or subjective operations. A typical phenomenon is that a few experimental constitutive relationship curves may not have obvious residual deformation stage (lack of residual shear stress τ r ), that is, the shear stress still keeps a downward trend at the end of the shear process, but this stage still cannot be omitted in the constitutive relationship modeling process.

Determination of the Peak Shear Stress τ p Under Seepage Pressure
It is strongly evident in Figures 1, 3 that with the increase of seepage pressure, the peak shear stress τ p decreases continuously. From the view of influence mechanism, Xia and Sun (2002) pointed out that both the cohesive force c and internal friction angle ϕ of most rocks (except mudstone, shale and other soft rocks) showed no significant changes due to the existence of water, and therefore the reduction of shear strength was mainly caused by the decrease of the effective normal stress on joint surface, i.e., (σ n −P w ), which can be explained by the effective stress principle of geotechnical materials (Xia et al., 2020), namely: where σ n is normal stress, P w is seepage pressure, ϕ and c are internal friction angle and cohesive force, respectively. By means of wavelet transform and 3D laser scanning technology, Qian (2018) established a new peak shear strength criterion considering the influence of joint match coefficient (JMC) and the first order undulations and second order roughness of joints, namely: where JMC is the joint match coefficients, JCS is the joint wall compressive strength, ϕ b is the basic friction angle of the joint, z + 2w and z + 2r are the slope root mean square of first-order waviness and second-order unevenness, respectively.
Actually, in the process of shearing, due to the alternant occurrence of shear shrinkage and shear dilatation on joint surface, the contact interaction emerged between part of the joint surface will trigger the instantaneous variation of the action of seepage pressure. As a consequence, the effective stresses acting on rough joints are dynamic. Referring to the modification of the effective stress principle by Biot (1963) (i.e., Biot coefficient), a seepage pressure reduction coefficient symbolically expressed as β, representing the average effect of seepage pressure on shear strength reduction, was introduced in this section. And therefore, Equation 6 can be recast into: The relevant test data in Table 1 were substituted into Equation 7, and the seepage pressure reduction coefficient β was calculated as 0.825. In order to find a good model to predict the peak shear stress, by comparing Equation 7 with JRC-JMC model, (Xia et al., 2014;Yang et al., 2016), the prediction errors of each model are demonstrated in Table 2: To facilitate the comparison of the error between the experimental peak shear stress and the predicted peak shear stress of different models, the average estimation error e ave was adopted as the precision index, as shown in Equation 9 (Liu et al., 2017): where τ peak,exp is the experimental peak shear stress,τ peak,est is the estimated peak shear stress, and n is the number of joint samples. The comparison results in Table 2 show that the deviations with the measures values of the four models are 8.805% (Equation 7), 11.326% (JRC-JMC), 27.353% (Yang et al., 2016), 29.338% (Xia et al., 2014), respectively. Obviously, the prediction error of Equation 7 is significantly smaller than that of other models. Therefore, Equation 7 is selected as the method to establish the peak shear stress τ p in this study.

Determination of the Peak Shear Displacement u p Under Seepage Pressure
In literature (Qian, 2018), the existing models commonly used to calculate the peak shear displacement u p were compared: where L 0 is the laboratory-scale joint sample length (in meters) and Ln is the field-scale length (i.e., the spacing of cross-joints that gives the effective block size) (in meters). In this modeling attempt, Ln was set equal to L 0 .
The average deviation equation (See Equation 9) is still used to analyze the prediction effects of each model in Equation 10. The comparison results are tabulated in Table 3, it can be seen from Table 3 that Qian's model has the best effect (the average deviation is only 7.49%). Therefore, Qian's model is selected as the method to predict the peak shear displacement u p in this study.

Determination of the Yield Shear Displacement u y Under Seepage Pressure
The yield shear displacement u y is the turning point of the shear stress-displacement curve from the elastic deformation stage to the pre-peak softening stage. At present, there is no general theoretical analytic solution for the determination of the position of this point. With reference to the method of Xia et al. (2014) and Qian (2018), the shear displacement of the joint sample at the moment where the initial dilatancy occurred was taken as the yield shear displacement. The ratio of shear displacement at the time of dilatation to peak shear displacement of each sample can be obtained from reference (Qian, 2018). Finally, u y = 0.38u p was selected as the yield point in this study.

Determination of the Shear Stiffness k s Under Seepage Pressure
The determination of shear stiffness k s (i.e., the slope of elastic deformation stage) strongly depends on the accurate value of yield shear stress τ y . Goodman (1976) conducted a large number of joint direct shear tests, and the results showed that the yield shear stress was approximately equal to 70-90% of the peak shear stress. Amadei et al. (1998) analyzed the direct shear test characteristics of natural joint sample and artificially created joint sample under constant load, and found that 75% of the peak shear stress was suitable for yield shear stress. Hungr and Coates (1978) found that 0.9τ p as the yield shear stress is suitable for engineering applications. Although the above literatures have evaluated the yield stress determination of joints based on experimental data, there is no general method that can accurately capture the yield stress of joints in complex stress states. According to the test data of different types of joints, it is a common method to obtain a summary of the test data with yield stress characteristics from a limited number of strength tests. In section Determination of the Shear Stiffness k s Under Seepage Pressure, the method for determining the yield point has been discussed. The shear stress corresponding to the yield displacement selected in section Determination of the Shear Stiffness k s Under Seepage Pressure and the corresponding peak shear stress were plotted in the Y-X coordinate system to fit the curve. When τ y = 0.7213τ p , the optimal fitting effect can be observed. The specific process is described in detail in reference (Qian, 2018), which is not shown in this section due to space limitations.
In this case, the shear stiffness k s can be determined by the following equation: Therefore, the shear stiffness ks determined in this verification is k s = 1.898 τ p u p .

Determination of the Residual Shear Stress τ r Under Seepage Pressure
As described in Figure 4, the constitutive relationship is ideal, with obvious residual deformation stage, as well as exact residual shear stress (τ r ), which is the key point when using the Equation 1 to simulate the shear stress-displacement curve. However, in the actual stress-displacement curve (Figures 1-3), there is no constant residual stage in the curve. Therefore, the suggested methods by International Society for Rock Mechanics (ISRM) were adopted for determining the residual shear stress τ r . A value of τ corresponding to a variation lower than 5% during 1.0 cm shear displacements is selected. The residual shear stress τ r is generally determined by the following equation: where ϕ r is the residual internal friction angle.
Frontiers in Earth Science | www.frontiersin.org In order to simplify the modeling process, the common practice is to use basic friction angle ϕ b instead of residual friction angle ϕ r . However, on the basis of the test results of a large number of laboratory direct shear tests on natural tensile fractured joints (gneiss joints, limestone joints, granite joints) and corresponding tensile joint replicas, Park et al. (2013) found that the basic internal friction angle (ϕ b ) slightly differs with the residual internal friction angle (ϕ r ). Moreover, an empirical estimation method for the residual internal friction angle (ϕ r ) in the percolation environment is recommended in literature (Park et al., 2013;Qian, 2018), as shown in Equation 14: where ϕ p is the peak dilatancy angle, a and b are dimensionless parameters used to fit test data, respectively.
The peak dilatancy angle ϕ p under higher normal stress can be determined by Equation 15 (Bandis et al., 1983): Due to the limited space, the relevant experimental regression results are not presented here, and the final residual shear stress data used for verification are indicated in Table 4.

Model Verification
In order to evaluate the capability and rationality of the constitutive model proposed in this paper (Equation 4), the experimental example in reference (Qian, 2018) under high seepage pressure was selected as the verification. The experimental parameters (i.e., the parameters in the proposed model with clear shear physics meanings) ks, τ r , u p , τ p in Equation 4   Specifically, M-RJ2 test data was adopted here as an illustrative example for detailed validation. Before the shear-flow coupling test, the relevant roughness parameters and mechanical properties parameters of the M-RJ2 sample can be determined through 3D morphology scanning and relevant basic mechanical properties testing, as shown in Equation 16. Also note that these parameters in Equation 16 have been shown in the previous sections (see Tables 1-4), and the re-display here is only for the specific test (M-RJ2 sample) to detailedly elaborate on the steps that these parameters in Equation 16 are used to solve the specific model parameters (τ p , u p , m, u 0 , etc.).

M-RJ2 sample can be computed by substituting the relevant test data in Equation 16 and Equation 17
: By combining Equations 17 and 18 with Equation 2, it can be calculated that the Weibull distribution parameters u 0 and m of M-RJ2 sample are, respectively: In such cases, all these constitutive model parameters in Equation 4 have been solved, so the predicted specific constitutive model expression of corresponding M-RJ2 sample can be derived as following: The comparative plots for the experimental data (M-RJ2 sample) and predicted model curve (Equation 20) are plotted in Figure 5 below. As shown, there is a great agreement between the experimental data of M-RJ2 sample and predicted model curve.
And the obtained correlation coefficient R 2 value of 0.978 reveals that the proposed model sufficiently simulates the shear characteristics of rock joints under seepage pressure. It is worth noting that the peak values of the predicted curve and the experimental results are not completely consistent because the predicted τ p and u p have certain errors with the experimental results. Of course, this error is relatively acceptable. In sections Determination of the Peak Shear Stress τ p Under Seepage Pressure and Determination of the Peak Shear Displacement u p Under Seepage Pressure, the prediction error of the method used in this study has been compared with these widely used methods such as the JRC-JMC model and Barton model.
In the same way, the proposed model is also used to simulate the remaining experimental results in Figures 1-3, and the predicted results are presented in Figures 6-8. Good agreement

Parametric Analysis
Although the Weibull distribution parameters m and u 0 can be calculated from Equation 2, their physical significance in engineering practice is difficult to be seen from the expression itself. Equation 20 was selected as the data source to analyze the influence of parameters m and u 0 on the constitutive relation curve. The relationship between shear stress-displacement curves under the action of different parameters m and u 0 is schematically shown in Figure 9. It can be clearly seen that the two parameters have no effect on the elastic deformation stage.
Specifically, with the increase of the value of the parameter u 0 , the peak shear stress increases, while the relationship between the parameter m and the peak shear stress shows a trend opposite to that of u 0 . When u 0 increased to 1.7 times of the original value (i.e., the u 0 values varied from 0.49338 to 0.83875), the peak shear stress increased by 18.12%, while m increased to 1.7 times of the original value (i.e., the m values varied from 0.39913 to 0.67852), the peak shear stress decreased by 26.525%, that is, m was more sensitive to the peak shear stress. The shear stress vs. shear displacement curves corresponding to different u 0 are approximately parallel after the yield point, and with the increase of m, the post-peak softening characteristic of joints becomes more obvious, that is, the stress drop from peak shear stress τ p to residual shear stress τ r occurs more rapidly.
The variation trend of model parameters under different normal stress and seepage pressure is extracted, as illustrated in Figures 10, 11. The results show that seepage pressure and normal stress have the opposite effect on model parameters. With the increase of normal stress, both m and u 0 decrease, while with the increase of seepage pressure, both m and u 0 increase.

Gray Relational Analysis
The qualitative research in section Parametric Analysis has found that the Weibull distribution parameters m and u 0 have different influence mechanisms on the peak shear stress τ p . In order to further quantitatively explore the influence degree of m and u 0 on the peak shear stress, the gray relational analysis is introduced. For the analysis of the relationship between objects with incomplete information of behavior mechanism, scarce behavior data and unclear inherent connotation, the gray relational analysis (Deng, 1982), as a mathematical method, can accurately find out the relevancy between influencing factors (comparative sequences) and system characteristic factors (reference sequence).
In this paper, the correlation between the model parameters m and u 0 and the peak shear stress τ p is calculated by using the gray correlation analysis method, and the main influencing factors are sought by comparing the correlation degree. The gray relational analysis steps are as follows (Ning et al., 2019): (1) Build the original data matrix The peak shear stress τ p of 27 groups of test data ( Table 2) was selected as the reference sequence, and the two parameters m and u 0 were analyzed as the comparative sequences. The original data matrix [X] can be derived as: x 10 x 20 . . .
In the above matrix, the first column of the original data matrix is the peak shear stress corresponding to each test, the second and third columns are the corresponding parameters m and u 0 , where x ij (i = 1,2, . . . 27; j = 0,1) represents the corresponding data of the i-th test on the j-th factor.
(2) Construction of the initial value matrix In order to eliminate the influence of dimension and transform the original data into relative values around 1, this paper adopts the method of initial value processing according to the following equation: The relevant data are substituted into Equation 22 to obtain the initial value matrix as shown below: (4) Construction of the gray relational coefficient matrix In the absolute difference matrix [∆], the maximum ∆ max and minimum difference ∆ min can be defined as: According to Equation 27, the gray relational coefficient matrix [ϑ] is expressed as: where ϑ ij = min +ω max ξ ij +ω max , and ω is the distinguishing coefficient within the range of [0,1], normally being 0.5.
(5) Calculation of the gray relational coefficient The gray relational coefficient γ reflects the effect degree of the influencing factors (comparative sequences) on the system characteristic factors (reference sequence). The gray relational coefficient γ j can be expressed as Equation 29: It can be calculated that the gray relational coefficient of m (γ 1 = 0.3560) is greater than the gray relational coefficient of u 0 (γ 2 = 0.3471), that is, the relation between m and the peak shear stress is closer than that between u 0 and peak shear stress.

CONCLUSION
The prediction of shear deformation of rock joints in the percolation environment is a vital component in geotechnical engineering. In this study, a new method for estimating the shear constitutive model of rock joints subject to high seepage pressure is presented in detail. This method is originated from the new joints constitutive model proposed by Lin et al. (2019b). The modeling results are validated and compared with the experimental results. And the effects of Weibull distribution parameters m and u 0 on the constitutive relationship and their response to normal stress and seepage pressure were also detailedly lucubrated as well. The following key conclusions can be drawn: (1) The joint shear stress-shear displacement curve in the percolation environment can be divided into four stages: elastic deformation stage, yield stage, post-peak deformation stage, residual deformation stage. Based on the damage statistics theory and the relationship between model parameters and seepage pressure, a statistical damage constitutive model considering seepage pressure was proposed.
(2) The model proposed in this study can accurately predict the whole shear process of rock joints under different seepage pressure. In addition, the physical meaning of the model parameters is clear and sufficient. All the model parameters can be directly determined from the laboratory direct shear tests, indicating the convenience of the model in the modeling process. (3) The mechanism of influence of normal stress and seepage pressure on Weibull parameters m and u 0 is different. As the normal stress increases, both Weibull parameters m and u 0 decrease, and as the seepage pressure increases, both Weibull distribution parameters m and u 0 increase. The Weibull distribution parameters m and u 0 decrease has different influence mechanisms on the peak shear stress. Gray relational analysis reveals that the effect of m on the peak shear stress is greater than the effect of u 0 .

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
HL, YC, and JL contributed conception and design of the study. SX organized the database and wrote the first draft of the manuscript. YW performed the statistical analysis. RC wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.