An Enhanced Ubiquitous-Joint Model for a Rock Mass With Conjugate Joints and Its Application on Excavation Simulation of Large Underground Caverns

A conjugate jointed rock mass (CJRM) is a rock mass with two sets of intersecting joints formed from intact rock under shear. Its mechanical properties and excavation-induced hazards of large underground caverns are different from those of common rock masses because of the unique geological origin thereof. To demonstrate numerically the excavation responses of CJRM, the ubiquitous-joint model is enhanced by consideration of the specific mechanical behaviors of the rock mass. In the enhanced model, CJRM is considered as the composite of columns of rock and two sets of weak planes of joints. The local coordinates, failure modes, and failure sequences of the rock columns and joints are redefined based on the composite characteristics of CJRM, and the failure criteria and plastic potential functions are accordingly modified. The enhanced model is verified numerically by triaxial compression tests and then employed to simulate the excavation of large underground caverns of a pumped storage power station in China. Results show that the modification of the local coordinate system, failure modes, and failure sequences made in the enhanced model is suited to the simulation of the mechanical behaviors of CJRM. Compared with the original ubiquitous-joint model, the enhanced model allows better predictions of the distribution of plastic zones and magnitudes of deformations in simulating underground excavations in CJRM and helps to assess the excavation-triggered hazards more accurately.


INTRODUCTION
Joints are often found in rock masses around underground caverns and usually occur in one or more sets in different directions, cutting the rock mass into blocky structures (Bandis et al., 1983;Jaeger et al., 2007;Wu et al., 2018). Jointed rock masses can be divided into several categories based on the number of joint sets, such as layered (or bedded) jointed rock masses with only one dominant set of joints, conjugate jointed rock masses (CJRM) with two dominant sets of joints, and other jointed rock masses with three or more sets of joints (Jaeger et al., 2007). Each set of joints affects the mechanical properties of the rock mass (Hoek and Brown, 1997); therefore, suitable descriptions of the constitutive relationships thereof are necessary for numerical stability analysis and hazard assessment of underground excavations (Agharazi et al., 2012;Ding et al., 2019).
Many studies have been conducted to assess the strength and deformation characteristics as well as the excavation responses of jointed rock masses. Theoretical investigations mainly focused on the prediction of the mechanical parameters of jointed rock masses based on superposition theory (Goodman et al., 1968;Jaeger et al., 2007), elasto-plastic damage theory (Chen et al., 2010;Yang et al., 2019), or displacement discontinuity method (Shen et al., 2016;Do and Wu, 2020). In experimental investigations, researchers obtained stressstrain relationships by conducting triaxial compression tests and validated the theoretical predictions through comparative analysis with the experimental results (Nova, 1980;Tien and Kuo, 2001;Singh et al., 2002;Tien et al., 2006;Chen et al., 2012;Chang et al., 2019). Numerical-analysis-based researches investigated the effects of the distribution properties (such as the orientation and spacing) and mechanical parameters of layered joints on the failure mode, plastic zones, and deformation of the surrounding rock masses (Adhikary and Dyskin, 1997;Park and Adachi, 2002;Wang and Huang, 2014;Zhou et al., 2016;Sainsbury and Sainsbury, 2017;Zhou et al., 2019;Yang et al., 2019;Zhou et al., 2021) as well as the internal forces and failure modes of the reinforcements (Hatzor et al., 2015;Gao et al., 2019) of underground caverns. Furthermore, many researchers proposed or modified constitutive models for layered jointed rock masses. For example, Jaeger (1960) proposed a theory based on the action of a single weak plane to predict the strength properties; Nova (1980) (Nova and Zaninetti, 1990) used the traditional Mohr-Coulomb (M-C) strength criterion to analyze the directional failures and this method has been used in the development of the ubiquitousjoint (U-J) model in Flac 3D (Itasca Consulting Group, 2011). Other researchers (Adhikary and Dyskin, 1998;Sitharam et al., 2001;Wang and Huang, 2009;Sainsbury and Sainsbury, 2017;Zhou et al., 2017;Das et al., 2019) further improved the constitutive equations of U-J model based on equivalent continuum methods.
These studies indicate that layered jointed rock masses have significant anisotropy of strength and deformation, are prone to shear and tensile failures on the weak planes of joints, and are prone to bending, shear, and tensile failures of the intact rock layers (Chen et al., 2012;Shen et al., 2016;Chang et al., 2019;Do and Wu, 2020). In addition, other investigations have proved that many engineering failure events triggered by underground excavation occur in such jointed rock masses (Chen et al., 2013;Ding et al., 2019;Do and Wu, 2020), and with increasing complexity of the joint sets, the excavation-induced disaster is usually greater (Jiang et al., 2014;Hatzor et al., 2015;Zhao et al., 2020). However, although the mechanical behaviors of layered jointed rock masses have been widely studied, little attention has been paid to CJRM which has more complexity of joint sets and few corresponding engineering case studies have been reported. CJRM is usually shaped into rock columns (intact columnar rocks) by conjugate joints. Geological research indicates that the conjugate joints are formed when two sets of in situ shear stress act, and since the direction of shear stress does not change significantly within a local area, the joints in the same set are roughly parallel (Deng et al., 2009;Ning et al., 2009;Wang et al., 2020). Besides, there are no weak interlayers in the joints due to the shear stress-induced formation process thereof, so conjugated joints can be regarded as two sets of parallel weak planes instead of weak interlayers. Thus, CJRM can be considered as the composite of columns of rock and two sets of weak planes of joints (Nova and Zaninetti, 1990;Wang and Huang, 2014).
The rock columns and weak planes of joints have different failure modes in the surrounding rock masses during underground excavations (Jia and Tang, 2008;. Chinese researchers (Jiang et al., 2014;Fan et al., 2017;Li et al., 2018;Jiang et al., 2019;Zhao et al., 2020) investigated the excavation-induced failure of the columnar jointed basalt of Baihetan hydropower station and ascertained the effects of joints on the spatially inhomogeneous distribution of damage zones, but the failure modes of the rock columns and joints are scarcely understood. Wang and Huang (2009; extended U-J model to include the failure mode and deformation performance for blocky rock masses with multiple sets of joints; however, their model did not consider the specific mechanical behaviors of the two sets of joints and the rock columns of CJRM, such as the bending failure of the rock columns and the effect of the angle between the two sets of joints. This study improves U-J model for such specific mechanical behaviors of CJRM and uses the enhanced model to analyze the stability of underground caverns during excavation. In the enhanced model, the local coordinate system is redefined and the related failure modes, failure sequences, failure criteria, and plastic potential functions are accordingly modified on the basis of U-J model in Flac 3D . The enhanced model is verified by numerical triaxial compression tests and used to investigate the excavation-induced deformation and plastic zones of the rock mass surrounding the underground caverns of a pumped storage power station in China. By comparing the simulation results arising from use of the enhanced model with those by the similar constitutive models, including the Mohr-Coulomb (M-C) model and U-J model, results show that the enhanced model provides a better prediction of the distribution of plastic zones and magnitudes of deformation when used to simulate underground excavations in CJRM.

PREPARATION FOR THE ENHANCED MODEL
Assumption of the Failure Process CJRM ( Figure 1A) contains two sets of joints (joint sets 1 and 2), which cut the rock mass into columns. The mechanical behaviors of such rock mass are determined by the joints and rock columns (Wang and Huang, 2014;Chang et al., 2019;Do and Wu, 2020). In the enhanced model, the assumed failures of CJRM include three categories: 1) the shear and tensile failures on the weak Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 744900 planes of joints, 2) shear and tensile failures of the rock columns, and 3) bending failure along the rock column axes. In addition, the assumed failure sequences between the weak planes and rock columns are such that 4) the weak planes fail before the rock columns, and 5) the weak plane at the greater distance from the stress state point to the corresponding strength envelope fails first if both weak planes fail in the same iteration process. The reasons for making such assumptions are as follows: firstly, in the cross-section plane of the rock columns ( Figure 1B), the rock mass is blocky so that the strength and deformability are mainly controlled by the joints (Singh et al., 2002;Singh and Rao, 2005). Assumption 1) ensures that material failures in such plane occur only on the weak planes of joints and that plastic deformation in the cross-section plane is mainly controlled by sliding deformation on the weak planes and tensile deformation perpendicular thereto. Secondly, in the axial direction of rock columns, the failures of the rock mass mainly occur in the rock columns. The failure modes include shear and tensile failures of the rock columns and bending failure (Adhikary and Guo, 2002;Neff et al., 2008) along the rock columns, as per assumptions 2) and 3), respectively. Thirdly, studies (Jiang et al., 2006;Hatzor et al., 2015) have shown that jointed rock masses usually fail on the joints first, followed by failure of the intact rocks, so assumption 4) guarantees a reasonable sequence of failure events between joints and rock columns.

Definition of the Local Coordinate System
Similar to the principle of U-J model, the enhanced model is a completely equivalent continuum model (Mühlhaus, 1993;Adhikary and Dyskin, 1998;Sitharam et al., 2001) that the influence of the joints is smeared into the continuum description of the rock mass, so the distributions (such as the spacing and length) of the joints, except the orientation, are not necessarily defined explicitly (Agharazi et al., 2012;Wang and Huang, 2014;Zhou et al., 2021). In the enhanced model, the local coordinate system is redefined based on the orientations of the two sets of joints the better to describe the mechanical properties of CJRM.
One of the rock columns is taken as an example to illustrate the method for establishing the local coordinate system ( Figures  1B,C). First, the axial direction of the rock column is defined as the z′-axis; second, the direction perpendicular to z′-axis and in the plane of joint set 1 is defined as the x′-axis, and the direction perpendicular to z′-axis and in the plane of joint set 2 is defined as the y′-axis. Since the orientations of the two joint sets are not necessarily perpendicular, the angle between the x′-axis and y′-axis is between 0 and 180°on the cross-section of the rock column ( Figure 1D). The directions of the local coordinates x′-, y′-, and z′-axes are determined as described below.
Vectors J 1 and J 2 are defined as direction vectors normal to the weak planes of joint sets 1 and 2, respectively, and vectors x′, y′, and z′ are defined as the positive direction vectors of the local coordinate x′-, y′-, and z′-axes, respectively. Based on the perpendicularity relationships between the five vectors, the directions of the local coordinate axes are given by z′: J 1 · z′ 0 J 2 · z′ 0 z′ 1 , x′: J 1 · x′ 0 z′ · x′ 0 x′ 1 , y′: J 2 · y′ 0 z′ · y′ 0 y′ 1 , where J 1 and J 2 are calculated from the dip angles (Dip 1 and Dip 2 ) and dip directions (DD 1 and DD 2 ) of joint sets 1 and 2, respectively.
(2) FIGURE 1 | Method for establishing the local coordinate system of CJRM.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 744900 Matrix C is defined as the rotation tensor between the local and the global coordinate systems, and then the transformation of stress components between the two coordinate systems (Itasca Consulting Group, 2011) is given by where C is x cos x′, y cos x′, z cos y′, x cos y′, y cos y′, z cos z′, x cos z′, y cos z′, z

Generalized Stress and Strain Components
Weak planes of joints in the enhanced model are ubiquitous in the rock mass based on the theory of the equivalent continuum model (Mühlhaus, 1993;Adhikary and Dyskin, 1998;Sitharam et al., 2001). The stress and strain components at any point in the rock mass can be expressed as shown in Figure 2 without consideration of the spacing and length of the joints in the local coordinate system. On the weak planes of joint set 1, shear stress τ 1′ and shear strain c 1′ have the following forms: The equivalent forms are readily obtained for joint set 2: On the cross-sections of the rock columns, the similar forms are In what follows, the failures of rock mass are detected by the six components of the generalized stress vector: σ 1′1′ , σ 2′2′ , σ 3′3′ , τ 1′ , τ 2′ , and τ 3′ whose components of the corresponding generalized strain vector are ϵ 1′1′ , ϵ 2′2′ , ϵ 3′3′ , c 1′ , c 2′ , and c 3′ .

ENHANCEMENT OF U-J MODEL FOR CONJUGATE JOINTED ROCK MASS
On the basis of U-J model, the enhanced model for CJRM is proposed using the redefined local coordinate system and the generalized stress and strain components. The mechanical behaviors of CJRM, such as the number of joint sets as well as the failure modes and failure sequences of the joints and rock columns, are considered in the enhanced model and then implemented as a plug-in dynamic link library (.DLL) file into the finite difference code FLAC 3D using the user-defined constitutive models function.

Description of the Joints
The failure criterion for the weak planes of joints used in the enhanced model is a composite M-C criterion with a tension cutoff expressed in terms of σ i′i′ and τ i′ (where i 1, 2, 3), as illustrated in Figure 3.
When shear failures occur on the weak planes of joints, the stress states calculated using the initial elastic estimate are located in Domain 2 (Itasca Consulting Group, 2011). The corresponding failure criteria meet the following conditions: For joint set 1, For joint set 2, and where c J1 , c J2 , ϕ J1 , ϕ J2 , σ t J1 , and σ t J2 are the cohesion, angle of internal friction, and tensile strength of the weak planes of joints, respectively. For weak planes with non-zero internal friction angles, the maximum values of the tensile strengths are given by The potential functions g s J1 and g s J2 corresponding to a nonassociated flow rule have the following forms: For joint set 1, For joint set 2, where ψ J1 and ψ J2 are the angles of dilation of the weak planes. When tensile failures occur on the weak planes of joints, the stress states calculated from the elastic estimate are located in Domain 3. The corresponding failure criteria meet the following conditions: For joint set 1, For joint set 2, The potential functions g t J1 and g t J2 corresponding to an associated flow rule have the following forms: For joint set 1, For joint set 2, When no failures occur on the weak planes of joints, the stress states calculated using the elastic estimate are located in Domain 1 and are considered as the final incremental stress without any correction by way of potential functions. The corresponding failure criteria meet the following conditions: For joint set 1, For joint set 2, In addition, it should be noted that both of the weak planes may fail in the same iteration process under a certain stress state in the rock mass. The choice of potential functions for stress correction is determined by the failure sequence and affects the corrected stress state. In the enhanced model, the failure sequence is judged by the distances, d s Ji and d t Ji , from stress state points A and B to the corresponding envelopes, f s Ji and f t Ji , as illustrated in Figure 3. The greater the distance (either d s Ji or d t Ji ), the greater the amount by which the stress calculated by the elastic estimate exceeds the corresponding failure envelope; so, the weak plane with the greater value of d s Ji or d t Ji is more prone to failure than the other and it fails first. The distance functions have the forms: (20)

Description of the Rock Columns
The failure mode of the rock columns involve shear and tensile failures of the columns and bending failure along the columns as assumed (Adhikary and Guo, 2002;Neff et al., 2008). When shear and tensile failures occur, a similar procedure to that used in M-C constitutive model (Itasca Consulting Group, 2011) can be used to calculate the incremental stress and strain. This section mainly discusses the bending failure of the rock columns.  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 744900 Figure 4 illustrates the bending failure of a rock column in the local coordinate system. When a bending moment acts on the rock column, as shown in Figure 4, the rock column is subjected to compressive stress on the upper side and tensile stress on its lower side. With increasing bending moment, the tensile stress on the lower side increases until it exceeds the tensile strength, resulting in a bending failure on the lower side (Agharazi et al., 2012;Sainsbury and Sainsbury, 2017). However, the tensile stress is difficult to calculate accurately, as joints are smeared into the continuum description of the rock mass by an equivalent continuum method and have no boundaries and thickness (Agharazi et al., 2012;Wang and Huang, 2014;Zhou et al., 2021). Thus, the normal stress on the column cross-section σ′ 3′3′ is used as an approximate replacement for the tensile stress. The bending failure criterion has the form: where σ b 0 is the bending strength of rock columns. Generally, the bending strength is less than the tensile strength although it has similar forms of failure criterion in the z′-direction. The potential function g b corresponding to an associated flow rule has the form:

VERIFICATION OF THE ENHANCED MODEL Verification Settings
Numerical triaxial compression tests were conducted to verify the enhanced model. The samples used in the tests were cylinders with a diameter of 1 m and a height of 2 m. The compressive strength, plastic zones, and deformation of the samples were analyzed: contrasting simulations were conducted using the enhanced model and U-J model. The material parameters used in these simulations are listed in Table 1. The confining pressure used in all such simulated tests is 1 MPa. The difference of the triaxial compression test samples between use of U-J model and the enhanced model mainly lies in the consideration of the number of joint sets and local coordinate system. For the sample simulated using U-J model, only one set of joints is considered at a time ( Figure 5A); the three local coordinate axes are mutually orthogonal and meet the righthanded spiral criterion. By contrast, two sets of joints are considered at a time for the sample simulated using the enhanced model ( Figure 5B); the included angle of the local coordinate axes x′ and y′ is determined by the occurrence of joints, which ranges from 90°to 270°, and the coordinate axis z′ is perpendicular to the plane formed by axes x′ and y′.
During the numerical simulations, the compression stress and confining pressure acting on the samples are initially 0.1 MPa and gradually increased with an increment of 0.1 MPa. The confining pressure stops increasing when it reaches the target value of 1 MPa and then stays the same. The compression stress keeps increasing until the calculation no longer converges. The total compression stress of the final converged solution is taken as the compressive strength of the samples (Neff et al., 2008;Chang et al., 2019).

Results and Discussion
Strength Figure 6A shows the compressive strength of the samples simulated using the enhanced model and U-J model, considering different joint sets with varying dip angles. As shown in the figure, there is a significant difference between the compressive strength calculated using the U-J model and that calculated using the enhanced model. When U-J model is used to consider one of the joint sets, the relationship between the triaxial compressive strength and the dip angles of joints is a typically U-shaped curve, where the dip angles are the acute angles between the normal direction of the joints and the axis direction of the samples. The compressive strength is at a minimum when the dip angle of the joints is close to (π/4 + ϕ j /2) where ϕ j is the friction angle of the weak planes of joints (Jaeger et al., 2007). However, when the enhanced model is used to consider the two joint sets, the compressive strength is the lower value of those calculated by U-J model, considering one of the two joint sets, and the relationship may no longer be U-shaped curve.
It also can be seen from Figure 6A that the strength of the samples of CJRM in Table 1 is mainly controlled by one of the two sets of joints with dip angles ranging from 30°to 90°. The strength is controlled by joint set S2 when the dip angles range from 30°to 48.5°, by joint set S1 when the dip angles range from 48.5°to 90°, and by both of S1 and S2 when the dip angles are about 48.5°. The strength of CJRM behaves consistently with the theoretical prediction (superposition theory) of multiple sets of joints proposed by Jaeger (Jaeger et al., 2007), which indicates that the enhanced model is implemented correctly and suitable for strength simulation of CJRM.

Plastic Zones
Figures 6B-D show the performance of the enhanced model in terms of the simulation of plastic zones. The samples with dip angles of 35°, 48.5°, and 75°, respectively, are taken to study the failure modes of CJRM based on the behavior in compressive strength. As shown in the figures, the plastic zones mainly occur in the strike direction of joint set S2 ( Figure 6B) when the dip angles of the two sets of joints are 35°, in the strike direction of joint set S1 ( Figure 6D) when the dip angles are 75°, and in the strike directions of both the joint sets S1 and S2 ( Figure 6C) when the dip angles are about 48.5°. This phenomenon shows that with the change of the dip angle of the two sets of joints, the failures of the samples occur on different sets of joints. In general, the failure modes of the samples coincide with the compressive strength, indicating that the enhanced model can reasonably reflect the failure modes of CJRM.

Deformation
The deformations of the samples of CJRM with dip angles of 35°, 48.5°, and 75°, respectively, were calculated using the enhanced model. The horizontal components of the total deformations are presented to illustrate the contribution of joints to the total deformations ( Figures 6E-G). As shown in Figure 6E, the horizontal deformation of the sample is at a maximum in the strike direction of joint set S2 when the dip angles of joints are 35°. Similar results can be obtained when the dip angles of joints are 48.5°( Figure 6F) and 75°( Figure 6G), respectively, showing that the deformations of the samples increase in the strike directions of the joints controlling the strength and plastic zone distribution of the samples. Comparing the distribution of plastic zones of the samples, it can be inferred that the increased deformation in the strike direction of the joints is derived from the plastic flows of the plastic zones occurring on the weak planes of joints. In conclusion, the joints in CJRM exert an important influence on the deformation, which is consistent with our previous assumptions.

Background of the Engineering Case
The case study involves the underground excavation simulation of a pumped storage power station in China. The underground caverns are located in a medium-coarse-grained granite ( Figure 7A). Two conjugate sets of joints, J1 and J2, occur therein ( Figure 7B). The axial directions of the powerhouse and transformer chamber are parallel and lie along the north-south direction ( Figure 7C). The in situ stresses in the rock mass are such that 1) the horizontal major principal stress σ 1 ranges from 12 to 18 MPa with the east-west direction; 2) the horizontal minor principal stress σ 3 ranges from 7 to 11 MPa with the northsouth direction; and 3) the vertical principal stress ranges from 6.6 to 9.3 MPa. The dimensions and excavation sequence of the underground caverns are shown in Figure 7D. The main support measures include 1) mortar anchors with lengths of 6 or 9 m at 1.5 m × 1.5 m grid spacings, 2) prestressed anchor cables with a length of 25 m at 4.5 m × 4.5 m spacings, and 3) steel fibreconcrete liners with a thickness of 150 mm sprayed onto the surface of the cavern walls. The mechanical properties of the rock mass are listed in Table 2 and those of the support structures are listed in Table 3.
The main engineering problems encountered in the underground caverns during the excavation are as follows: 1) with the downward excavation of the powerhouse and transformer chamber, the stress field of the surrounding rock mass evolves, resulting in joint slippage and cracking of the liners ( Figure 8A); 2) after the excavation of IV step is finished, the total deformation of the surrounding rock mass with a maximum value of 100 mm is much greater than the designed value of 50-60 mm. The monitored total deformations on some main cross-sections of the powerhouse are shown in Figure 8B. Compared with other underground engineering with similar lithology, in situ stress, size of excavation, and other conditions (Zhu et al., 2008;Zhu et al., 2010), the most prominent geological feature of the present engineering is that the conjugate joints in the rock mass are relatively well-developed. Therefore, it can be inferred that one of the possible factors inducing the aforementioned engineering problems is the strength reduction of the rock mass due to the conjugate joints, making it necessary to use the enhanced model Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 744900 to study the plastic zones and deformation characteristics of the surrounding rock mass.

Preparation for the Simulation
In the excavation simulation of the underground caverns, the excavation sequence is I through VI, as shown in Figure 7D. The anchors and cables are simulated by cable elements and the liners are simulated by shell elements (Itasca Consulting Group, 2011). The parameter settings of the anchors, cables, and liners are consistent with those listed in Table 3. The enhanced model, as well as M-C model and U-J model used as contrasting models, are each used to calculate the distribution of plastic zones and deformation of the surrounding rock mass. The rock mass properties used in the simulation are those listed in Table 2.  M-C model and U-J model are selected as contrasting models for the enhanced model due to both the similarities and differences between them. The similarities lie in the failure criteria and the potential functions for joints and intact rocks, i.e., they are all derived from M-C yield criterion. The differences lie in the consideration of the numbers of joint sets and the enhancements in the local coordinate system, failure modes, and failure sequences of the enhanced model. In terms of the numbers of joint sets, neither joint set J1 nor J2 can be considered in M-C model, either J1 or J2 can be considered in U-J model, but both J1 and J2 can be considered in the enhanced model. Thus, by comparing the results calculated using the three models, the effects of joints on strength and deformation behaviors of the surrounding rock mass and the responses of the enhanced model in underground excavation simulations can be determined.

Distribution of Plastic Zones
The distributions of plastic zones calculated by the three models are presented in Figure 9: when U-J model is used to simulate either joint set J1 or J2, the distributions of plastic zones are mainly concentrated on the upper left and lower right corners of the powerhouse (when considering J1) or on the lower left and upper right corners (when considering J2). By contrast, when the enhanced model is used to simulate both joint sets J1 and J2, the plastic zones are distributed all around the powerhouse. The notable difference in the distributions of plastic zones indicates that the plastic zones are significantly affected by the joints and it is necessary to fully consider the effects of the two sets of joints in the simulation of excavation-induced failure in CJRM.
The volumes of plastic zones calculated using the enhanced model are significantly larger than those calculated using the U-J model (Figure 9): the difference in the volumes is quantified by comparing the normalized volumes calculated using the three models ( Figures 10A,B). Compared with the plastic zone volumes calculated using M-C model considering no joints, those calculated using U-J model considering either joint set J1 or J2 are almost 1.8 times greater, and those calculated using the enhanced model considering both joint sets J1 and J2 are 2.5 times greater ( Figure 10A). As the excavation proceeds, the growth rates of plastic zone volumes calculated using each of the three models become significantly different ( Figure 10B). From the first excavation step to the final excavation step, the plastic zone volumes increase 1-2.6 times when using M-C model, 1.8-5 times when using U-J model, and as many as 2.2-7 times when using the enhanced model.
The proportional occurrence of the failure types of the rock mass is different when using different models. As shown in Figures 10C,D, when M-C model is used, the failure of the rock mass is mainly tensile, and the proportions of shear failure and tensile failure are 20 and 80%, respectively. On the contrary, when the U-J model or the enhanced model is used, the rock mass is mainly subjected to shear failure. The proportions of shear failure and tensile failure calculated using U-J model are 60 and 40%, respectively, while the proportions calculated using the enhanced model are 90 and 10%, respectively.

Displacement Characteristics
The displacement contours ( Figure 11) calculated using each of the three models are similar in terms of the deformation mode:  FIGURE 10 | Comparison of (A) the volume ratios of total plastic zones (V tot /V tot mo ), (B) the growth ratios of plastic zone volumes (V tot /V 1st mo ), (C) the volume ratios of shear failures (V sh /V tot ), and (D) the volume ratios of tensile failures (V ten /V tot ) at each excavation step. V tot is the total failure volume (including volumes of shear and tensile failures); V tot mo is the total failure volume calculated using M-C model; V 1st mo is the failure volume calculated using M-C model at the first excavation step; V sh is the shear failure volume; V ten denotes the tensile failure volume.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 744900 11 the deformation of the surrounding rock mass develops towards the opening and the largest deformation occurs on the sidewalls of the powerhouse upon completion of the excavation. The deformation characteristics are consistent with those reported elsewhere (Jia and Tang, 2008;Hatzor et al., 2015), indicating the correctness of the enhanced model in terms of the deformation response. However, there are some differences among the deformation magnitudes calculated using the three models. Compared with the displacement magnitude contour calculated using M-C model ( Figure 11A), those calculated using U-J model ( Figures 11B,C) are increased on the sidewalls (especially where the arrows indicate in the figures) and the magnitude calculated using the enhanced model ( Figure 11D) is increased on almost all of the sidewalls (especially where also found to have been increased when using U-J model, as the arrows indicate in the figure).
For further study, the magnitudes of the total displacements calculated using the three models at each excavation step are compared ( Figure 12). The total displacements calculated using the three models contain both of the elastic and plastic components. The elastic components are almost the same since they are calculated under the same stress condition, but the plastic components are different since they are derived from plastic flows of different plastic zones. As shown in Figure 12, the displacement magnitudes of the monitored points calculated using the enhanced model are generally similar to, but slightly larger than, those calculated using M-C model and U-J model. This indicates that the total displacements are mainly controlled by the elastic components and the plastic components which are derived from plastic flows of plastic zones entailing only minor corrections to the total displacements.

DISCUSSION
Previous work has assessed the use of various constitutive models for layered jointed rock masses (Adhikary and Dyskin, 1998;Sitharam et al., 2001;Sainsbury and Sainsbury, 2017;Zhou et al., 2017) and their strength and deformation performance in the simulation of underground excavations (Jiang et al., 2006;Jia and Tang, 2008;Das et al., 2019;Ding et al., 2019). These studies were mainly focused on rock masses with a set of joints, while little attention has been paid to the mechanical characteristics of CJRM and their excavation responses in underground engineering. In fact, CJRM or similar rock masses often appear in engineering practice. For example, conjugate sets of joints can be found among the multiple sets of joints in the rock mass at Baihetan hydropower station (Jiang et al., 2014;Fan et al., 2017;Li et al., 2018;Jiang et al., 2019;Zhao et al., 2020). In this study, U-J model was enhanced by considering the specific mechanical behaviors of CJRM based on the equivalent continuum method (Mühlhaus, 1993;Adhikary and Dyskin, 1998;Sitharam et al., 2001). The enhanced model was verified by numerical triaxial compression tests and used to simulate the excavation responses of the underground caverns of a pumped storage power station in China. The comparison of results calculated using the M-C model, U-J model, and enhanced model shows that the mechanical behaviors of jointed rock masses are greatly affected by the number of joint sets. In the triaxial compression tests, the strength of the samples of CJRM is the lower one of the two strengths of the corresponding layered jointed rock masses. These findings are consistent with the superposition theory of multiple sets of joints proposed by Jaeger et al. (2007). In addition, because of the modification of the local coordinate system, failure modes and failure sequences based on the characteristics of CJRM, the enhanced model shows a significant difference from M-C model and U-J model in the simulated excavation of underground caverns in CJRM. The mechanical responses, such as the distribution of plastic zones and displacement of the surrounding rock mass, are significantly underestimated by the M-C model or U-J model but are reasonably estimated by the enhanced model, suggesting that the enhanced model can be used to simulate the mechanical behavior of underground excavations in CJRM.
The prediction of the excavation-induced plastic zones (also considered as excavation-disturbed zones, EDZ) of the surrounding rock mass is of great significance to the design of support measures and the assessment of related hazards of underground caverns (Parise and Lollino, 2011;Chen et al., 2013;Zhao et al., 2020). For the present engineering case, the plastic zone volumes calculated using the enhanced model are 2.5 and 1.4 times greater than those, respectively, calculated using the M-C model and U-J model. Moreover, the failure modes of the surrounding rock mass calculated using each of the three models are quite different. The failure mode predicted using M-C model involved mainly tensile failure, which accounted for 80% of all failures; however, the failure modes predicted using the U-J model and the enhanced model involved mainly shear failure, accounting for 60 and 90% of all failures, respectively. Comparing the failure modes and failure scales predicted using the three models with the practical engineering problems, including joint slippage and cracking of the liners ( Figure 8A), it can be found that the enhanced model provides more accurate predictions than M-C model and U-J model in the simulation of underground excavations in CJRM.
The total deformations of surrounding rock masses are composed of elastic and plastic parts (Chen et al., 2013;Wang and Huang, 2014). For the three models used in the present engineering case, there is no difference in the calculation of elastic deformations, but a significant difference in those of the plastic deformations. Compared with the M-C model, the enhanced model and U-J model add deformation corrections of the joints to the total deformations, i.e., the U-J model considers one of the two sets of joints and the enhanced model considers both sets of joints. The results showed that the deformation patterns and magnitudes calculated using the three models differed (albeit not to any significant extent), indicating that the total deformations are mainly controlled by the elastic deformations and the plastic deformations made only minor corrections thereto. Comparing the calculated deformations with those measured in situ, the deformation calculated using the enhanced model is closer to those monitored, even if there remained certain differences between them. In future work, the prediction of the deformation of CJRM can be further modified by considering the stiffness of joints, thus improving such simulations (Wang and Huang, 2009;Chen et al., 2012).

CONCLUSION
CJRM or the similar rock masses often appear in engineering practice, but their specific mechanical performance is much less explored in numerical simulation (Wang and Huang, 2014;Yang et al., 2019;Zhao et al., 2020). This study enhanced U-J model by consideration of the particular mechanical behaviors of CJRM based on an equivalent continuum method and used the enhanced model to study the excavation responses of the large underground caverns in CJRM. It is obtained from the study that 1) The modification of the local coordinate system, failure modes, and failure sequences based on the structural and Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 744900 mechanical characteristics of CJRM is suitable for the simulation of the strength, deformation, and distribution of plastic zones of the rock mass. The presence of the conjugate joints greatly affects the performance of CJRM, and the effect of each set of joints on the mechanical behaviors is consistent with the superposition theory of general jointed rock masses.
2) The underground excavation in CJRM causes a much larger scale of plastic zones of the surrounding rock mass compared with those in unjointed or layered jointed rock masses. The failure modes of the underground caverns in CJRM involve mainly shear failure of the joints, which accounts for about 90% of all failures and is manifested as the joint slippage and cracking of the liners in the practical engineering.
3) The deformation of the underground excavations in CJRM is larger than those in unjointed or layered jointed rock masses. This is partly reflected in the plastic flow generated by the larger scale of the plastic zones in the calculation using the enhanced model; however, the deformation component derived from joint slippage is not considered, so the calculated deformation is smaller than that in practice. This weakness of the enhanced model can be improved in future work by considering the stiffness reduction of CJRM.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
XL designed the research. CW processed the corresponding data. XL and CW wrote the first draft of the manuscript. QS and HL helped to organize the manuscript. JC revised and edited the final version.