Numerical Investigation of Influence Factors on Underground Powerhouse Using an Anisotropic Ubiquitous Joint Model

1  Introduction

Pumped-storage hydroelectric facilities are essential components of modern energy systems, providing balance to grid demand through peak shaving, load leveling, and the integration of intermittent renewable energy sources such as wind and solar. These facilities function by transferring water between upper and lower reservoirs, converting gravitational potential energy into electricity during periods of peak demand and reversing the process to store surplus energy. At core of these systems is the underground powerhouse—a complex engineering structure typically excavated in anisotropic rock masses. The stability of the surrounding rock during and after excavation is paramount, as it directly impacts construction safety, operational longevity, and economic viability. However, optimizing the orientation of the powerhouse axis orientation presents a significant challenge due to the interplay of geological discontinuities (e.g., joints, faults) and in-situ stress anisotropy. Joint geometry—specifically dip angle (inclination from horizontal) and dip direction (compass orientation of steepest descent)—along with stress redistribution induced by excavation, govern deformation patterns and failure mechanisms such as shear slippage or tensile fracturing [1]. This study aims to address this complex issue by systematically analyzing how these factors influence rock mass stability, with the ultimate goal of establishing actionable criteria for axis alignment in geologically intricate environments.

To assess the stability of underground powerhouses, researchers have traditionally relied on two primary approaches: physical model testing and numerical analysis. Physical models, such as large-scale 3D geomechanical setups or acoustic emission-based monitoring systems [2], provide tangible insights into rock behavior by replicating excavation processes under controlled conditions. For instance, Zhu et al. [3] utilized physical models to demonstrate that high in-situ stress concentrations trigger spalling failures in deep tunnels, while Liang et al. [4] employed microseismic seismic monitoring to correlate crack evolution with stress redistribution. Despite their observational strengths, physical methods face inherent limitations: high costs, lengthy preparation times, and difficulties in scaling to field dimensions. Tiwari and Latha [5] further highlighted that simplified laboratory conditions often fail to capture the heterogeneity of natural joint networks, leading to overly optimistic stability predictions. In contrast, numerical methods—such as finite element analysis (FEM), discrete element modeling (DEM), and finite difference codes—offer scalable, cost-effective alternatives for simulating complex geological conditions. Recent advancements include Chen et al.’s [6] DEM-based investigation of spalling failure mechanisms and Meng et al.’s [1] discrete fracture network (DFN) models for hydraulic fracturing analysis. Nevertheless, existing numerical studies often isolate individual factors (e.g., stress orientation or joint persistence), neglecting the synergistic effects of joint dip angle, dip direction, and stress deviation—a gap this study aims to bridge. A novel DFN-DEM approach has been proposed to simulate the long-term behavior of deep resources, emphasizing the importance of mechanical calculations alongside hydraulic assessments [7]. This integrated methodology allows for a more comprehensive understanding of the interactions between fractures and the surrounding rock mass, which is crucial for effective resource extraction and stability analysis. In the context of underground excavations, the Discrete Element Method (DEM) modeling framework has been utilized to analyze the stability of large underground structures. Specifically, the integration of Discrete Fracture Networks (DFN) with DEM has been instrumental in evaluating structurally controlled failure mechanisms in blocky rock masses [8]. This approach offers valuable insights into the failure patterns that may arise during excavation, thereby informing design and operational strategies. Further advancements in the field include the development of a unique DEM model that incorporates the weakening of fracture strength, particularly in the context of borehole stability in naturally fractured rocks [9]. This model addresses the complexities associated with drilling mud interactions and the mechanical integrity of boreholes, which are critical for successful drilling operations in deep geological formations.

This research adopts Flac3D, a robust finite difference code widely recognized for its ability to handle large-scale, nonlinear geomechanically problems involving anisotropic materials and iterative stress corrections. In contrast to conventional Finite Element Method (FEM), which often struggles with complex joint configurations, or DEM, which requires excessive computational resources for fine-scale fracture modeling, Flac3D strikes a balance between accuracy and efficiency. Its explicit solution scheme enables dynamic simulation of excavation sequences, while the embedded ubiquitous joint model allows for the explicit representation of weak planes with user-defined orientation and strength parameters. Leveraging field data from a 120 MW pumped-storage project in sericite phyllite—a foliated rock characterized by well-developed joint sets with dip angles ranging from 30° to 40° and high horizontal stresses between 13.09 and 17.77 MPa—this study constructs a 3D numerical model encompassing 687,917 elements. Key parameters, including joint cohesion and friction angle, are derived from in-situ triaxial tests and geological surveys. The analysis systematically varies three factors: joint dip angle (30°–70°), dip direction (0°–90° relative to the powerhouse axis), and in-situ stress deviation angle (0°–60°), to quantify their impacts on displacement and plastic zone development.

The study’s objectives are threefold: to identify optimal joint orientations that minimize deformation and rock damage, to compare the predictive capabilities of physical and numerical methods in capturing joint-controlled failure mechanisms, and to propose axis alignment criteria tailored to anisotropic rock masses. By integrating theoretical models with practical engineering data, this work not only advances the understanding of joint-stress interactions but also provides a decision-making framework for designers and engineers. For instance, preliminary results indicate that aligning the powerhouse axis within 40° of the joint dip direction reduces displacement by 35%–40%, while stress deviations beyond 30° significantly amplify plastic zones—findings that have direct implications for site selection and support design. Ultimately, this research underscores the importance of holistic stability analysis in mitigating excavation risks and enhancing the sustainability of pumped-storage hydropower projects.

A pumped-storage power station utilized in this numerical study. The station is located in the southern region of a mountainous area, precisely at the juncture of low mountains and highlands. The outlet of the power station is initially designed to be linked to a substation, offering convenient access for electricity supply and reception. The conditions for receiving and sending electricity are favorable. Once completed, the power station will primarily perform functions such as grid peak shaving, valley filling, energy storage, frequency regulation, phase modulation, and emergency accident standby [6,7,10,11]. The pumped-storage power station is a daily regulated pure pumped-storage power station with an installed capacity of 120 MW. The hydraulic engineering project mainly comprises upper and lower reservoirs, water supply systems, underground powerhouses, and switching stations, among other components.

The upper reservoir is located upstream of a certain village and is formed by damming a valley. The water supply system mainly consists of the inlet and outlet structures of the upper reservoir, the water intake tunnel, the water pressure chamber, the water intake shaft, the lower horizontal water tunnel, the water steel branch pipe, the high-pressure steel support pipe, the tailwater branch pipe, the tailwater pressure chamber, the tailwater tunnel, and the inlet/outlet structures of the lower reservoir, etc. The dam site of the lower reservoir currently has a simple village road passing through, which can be connected to the surrounding road network via a county road and a provincial highway. The nearest wharf to the project is a port located at the confluence of the river and the Yangtze River.

The near field area of the project is situated in the hinterland, characterized by highly undulating terrain. The landforms in this region range from basins to mountains. The landscape is predominantly composed of mountains and hills, with limited areas of plains (intermountain valleys) and a relative elevation difference of 100 to 200 m. Hills are mainly distributed on the exterior of low mountains and the interior of basins and valleys, and the slope is generally around 20 degrees. Hills with a relative height of less than 100 m are mainly scattered on the inner side of the mountain basin, valley and on the basin floor and valley floor, lacking a definite extension direction. Most of the hills are rounded, and the slope is generally less than 15 degrees, with a maximum not exceeding 20 degrees.

The middle and low mountains are divided into erosion-accumulation low mountains and denudation high mountains. The erosion-accumulation low mountains rise 500 to 1000 m and are the most extensively distributed mountain type, and the slope of the mountain is mostly 25 to 30 degrees. In contrast, denudation high mountains exceed an elevation of over 1000 m, the denudated mountains are mainly distributed among three parts of the mountain range in strip and sheet forms, and the distribution direction is consistent with the main structural lines, predominantly oriented in the near east-west, northeast, and northwest directions.

The lowest drainage surface in the project area is the river. The south side features a middle and low mountain landform with a slope of 35 to 50 degrees, and the peak elevation of the topographic watershed exceeds 1000 m. The project area is located south of the river, at the junction of the middle mountain and the low mountain. The terrain is undulating, and the gully incision is deep. The upper and lower reservoirs of the power station are located in the upper and lower reaches of the Huixi River gully, which is deep and generally trends N30°E.The river near a village turns nearly east-west to join another river close to the same village.

The project site is situated within a fault block area characterized by interspersed tectonic units and crisscrossing fault zones, which have experienced multi-stage tectonic movements. The faulted basin developed since the early Tertiary period, accompanied by the eruption of alkaline magma. Tectonic activity in this region was particularly vigorous during the early Himalayan orogeny. The fault structure is well-developed throughout the area, with variations in fault distribution patterns across different regions. The fault zone acts as a boundary. In the western region, the predominant fault orientations are NNW and EW, with NE trending faults also present. Conversely, in the eastern area, the primary fault orientation is NE, followed by NNE and ENE trending faults.

There are no active faults within 25 km of the near field area, and recorded earthquakes have magnitudes of less than 4.7. The horizontal design ground motion acceleration of the project site exceeds 43 to 45 gal with a probability of 10% in 50 years. The basic intensity of the earthquake is classified as VI, and the structural stability of the site area is favorable, which is demonstrated in Fig. 1.

Figure 1: Picture of regional faults and geological structures

The in-situ stress data were acquired through hydraulic fracturing tests conducted at strategically selected boreholes within critical zones of the underground powerhouse [12–14]. Vertical boreholes, drilled to depths exceeding 140 m, and complementary horizontal boreholes extending approximately 35 m, ensured comprehensive coverage of stress conditions across varying elevations and geological strata. Measurements of fracture pressures, crack orientations, and closure stresses were rigorously analyzed to determine the three-dimensional stress magnitudes and directional trends, including both horizontal and vertical components. These results were cross-validated against theoretical overburden stress calculations to confirm consistency with expected geological behavior. As a globally recognized in-situ testing methodology, hydraulic fracturing provided direct, high-resolution stress profiles, ensuring reliability for numerical modeling and design optimization in complex underground environments.

According to the geological survey, the names and meanings of each parameter are presented in Table 1, whilethe parameters of the anisotropy joint model are specified as Table 2.

2  Methodology

2.1 Initilization of In-Situ Stress

The Mohr-Coulomb criterion, initially proposed by the German engineer Christian Otto Mohr in 1900, states that the shear failure of materials primarily occurs when the shear stress on a specific plane reaches its ultimate strength. This failure issignificantly influenced by the normal stress on that plane.

In the process of model construction and analysis, for thoroughly investigate the impact of dip direction and dip angle on the stability of surrounding rock in underground powerhouses, we maintain a constant position for the main axis of the powerhouse while systematically altering the dip direction and dip angle of the joint planes [15]. In most cases, the orientation of the maximum principal stress in rock is influenced by the stress in all three axes [4]. However, during excavation of an underground powerhouse, the ground can be treated as a semi-infinite plane, and the direction of the maximum principal stress becomes vertical to the ground.

As is demonstrated in Fig. 2, in a triaxial stress state, the intermediate principal stress σ2 is often neglected to enable the determination of the maximum or critical shear stress at a point. This is achieved by constructing a stress circle using the maximum σ1 and minimum σ3 principal stresses. The ultimate stress circle represents the failure condition of the material. In the plane stress state, the stress tensor of soil can be expressed as follows:

[σij]=[σxτxy0τyxσy0000](1)

Figure 2: Angle of in-situ stress in 3D dimension and 2D dimension: (a) 3D in-situ stress; (b) 2D in-situ stress

Here, σx, σγ, and τxγ are the stress components that define the stress state at a point in the soil [16]. The stress component on any inclined plane within a soil microelement corresponds to coordinates of a point on the Mohr stress circle. As the outer normal vector of the inclined plane rotates counterclockwise by an angle α counterclockwise and the corresponding point on the Mohr circle rotates clockwise by 2α clockwise, which is demonstrated in Fig. 3.

Figure 3: Mohr-Coulomb circle of principle and minimum stress

Given the known maximum σ1 and minimum σ3 principal stresses, the normal stress σn and shear stress τ on an inclined plane can be formulated as:

σn=σ1+σ32+σ1−σ32cos⁡2α(2)

τ=σ1−σ32sin⁡2α(3)

Since the normal stress σn and shear stress τ on the inclined plane can be linearly expressed in terms of the maximum σ1 and minimum σ3 principal stresses, the gradients of normal stress ∇σn and shear stress ∇τ on the inclined plane can also be linearly expressed through the gradients of the maximum ∇σ1 and minimum ∇σ3 principal stresses [5,17–19]:

∇σn=∇σ1+∇σ32+∇σ1−∇σ32cos⁡2α(4)

∇τ=∇σ1−∇σ32sin⁡2α(5)

This framework offers a solid foundation for understanding and predicting the shear failure of materials under complex stress conditions, making it an essential tool in geotechnical engineering for analyzing the stability of surrounding rock. This application of the Mohr-Coulomb criterion, along with other relevant literature, underscores the importance of considering both normal and shear stresses in the assessment and design of geotechnical projects. The method to realize this function in FLAC3D is outlined below.

In Flac3D, the process begins by defining a function that calculates stress components and their gradients based on the input parameters, which include the maximum and minimum ground stresses, their gradients, and the deflection angle. The function is then initialize with key properties such as dip, dip direction, joint cohesion, and other relevant factors [20,21], while setting the initial conditions like zero displacement and velocity. The algorithm iteratively evaluates the stress distribution across a range of angles, solving the model for each angle to generate a stress distribution curve that demonstrates how stress varies with changes in the deflection angle.

The Mohr-Coulomb criterion is a fundamental principle in geotechnical engineering by linking theoretical shear strength concepts to practical stability assessments [3,22]. At its core, the criterion quantifies how materials resist failure through a balance of inherent cohesion and stress-dependent friction—a concept that transitions from theory to practice when integrated into computational tools like FLAC3D. Engineers begin by defining site-specific parameters such as in-situ stresses, joint orientations (dip and dip direction), and material properties derived from laboratory or field tests. These inputs anchor the numerical model in physical reality, enabling the software to simulate stress redistribution around structures like underground powerhouses. As FLAC3D iteratively calculates stress states across varying geological conditions, it applies the Mohr-Coulomb threshold to identify potential failure zones, effectively translating equations into visual risk maps. For instance, by analyzing how changes in joint angles influence displacement patterns, engineers can determine optimal orientations for excavations or reinforcements, ensuring that designs operate safely within the material’s strength limits. This seamless integration of theory and technology not only predicts failure mechanisms but also empowers data-driven decisions, transforming abstract concepts into actionable strategies for mitigating risks in slopes, tunnels [12–14], and other critical infrastructures [23–25]. By bridging mechanistic understanding with practical validation, the Mohr-Coulomb criterion remains indispensable for balancing safety, efficiency, and innovation in geotechnical engineering.

2.2 Transverse Isotropic Ubiquitous Joint Model

There are two methods for solving joint plane problems that have emerged in recent years: the discrete fracture model method and the equivalent continuous model method. In FLAC3D, the equivalent continuous model is employed to address these issues. This research takes into account a transverse isotropy joint model to effectively manage the analysis. In practical scenarios, failures can occur in both the matrix and the joint planes, depending on which section presents a greater risk. This is closely related to the inclination and orientation of the joint planes, as illustrated in Fig. 4.

Figure 4: Joint plane dip angle and direction: (a) joint set distribution; (b) dip angle and direction

The transverse isotropy joint model is grounded on the Mohr-Coulomb yield criterion. When subjected to external forces, the total elastic strain increment in a transverse isotropy jointed rock mass model is partitioned into contributions from the intact rock and three sets of joints. In the global coordinate system, this can be expressed as follows:

Δεij=ΔεijorR+ΔεijJ1+ΔεijJ2+ΔεijJ3(6)

where Δεij is total strain increment of the equivalent continua and ΔεijorR is the intact rock strain increment in transverse isotropy model, and ΔεijJ1,ΔεijJ2,ΔεijJ3 are the strain increments in the joint sets. Two parts are introduced to get these strain increments.

2.2.1 Transverse Isotropy Elastic Rock

Assuming the intact rock is transverse isotropy, the incremental strain in the rock can be related to the increment of stress in the global coordinate system as follows [13]:

{τxx=c11εxx+c12εyy+c13εzzτyy=c12εxx+c22εyy+c23εzzτzz=c13εxx+c23εyy+c33εzzτyz=2c44εyz,τxz=2c55εxz,τxy=2c66εxy(7)

where cij is the elastic compliance tensor of the intact rock.

As C is block-diagonal, its inverse is found to be:

s11=(c22c33−c232)/D,s22=(c11c33−c132)/D,s33=(c11c22−c122)/D(8)

s12=(c12c23−c12c33)/D,s13=(c12c23−c22c13)/D,s23=(c12c13−c11c23)/D(9)

where D is the determinant of the upper-left-corner block of C, that is:

D=det[c11c12c13c12c22c23c13c23c33](10)

The inverse form of Hooke’s law for a transverse isotropy material is therefore:

[εxxεyyεzz2εyz2εxz2εxy]=[s11s12s13000s12s22s23000s13s23s33000000s44000000s55000000s66][τxxτyyτzzτyzτxzτxy](11)

which can be written as:

εij=SijmnorRτij(12)

The increment of stress in the global coordinate system as follows:

Δεij=SijmnorRΔτij(13)

Finally, we can write Hooke’s law for a transverse isotropy material as:

{εxx=+1Exτxx−vyxEyτyy−vzxEzτzzεyy=−vyxEyτxx+1Eyτyy−vzyEzτzzεzz=−vzxEzτxx−vzyEzτyy+1Ezτzzεyz=12Gyzτyz,εxz=12Gxzτxz,εxy=12Gxyτxy(14)

2.2.2 Embeded of Weak Plane

For joint plane, we can use strain energy to get its displacement.

For a representative volume of V subjected to a stress increment of Δτij with corresponding strain increment of Δεij, the increment of strain energy due to deformation along joint set k, ΔEJk, can be mathematically written as:

ΔEJk=σijΔεijJkV(15)

where σij is average stress. Assuming that joints are persistent and planar in the model, the tractions on joint, τij(i≠j), can be calculated as follows:

τj=σijni(16)

where ni is normal unit vector of joint. The work done by traction on the kth joint plane can be calculated as:

ΔwJk=τjkΔδjkAk(17)

where Δδjk is elastic displacement along kth joint plane.

According to the principle of conservation of energy, the part of the elastic strain energy stored in the system due to elastic displacement along a given joint set (Eq. (18)), is equal to the work done by surface traction on the joint plane as defined in Eq. (20). By substituting Eq. (19) into (20), the strain increment caused by elastic displacement along kth joint plane is given by [15]:

ΔεijJk=niΔδjk1lk(18)

where lk is the dimension in the direction normal to the joint plane kth.

In the global coordinate system, the joint displacements can be related to stresses acting on the joint plane by using joint elastic compliance matrix D as follow:

δj=Djmτm(19)

Substituting Eq. (19) in (22) and an incremental form of Eq. (22) into (21), yields:

ΔεijJk=nikDjmknnk1lkΔσmn(20)

From Eq. (23) the elastic compliance of joint set kth can be expressed in the form of a fourth-order tensor, CijmnJk, as follows:

CijmnJk=nikDjmknnk1lk(21)

Following the same steps for all joint sets in the model, the equivalent compliance tensor of the model, CijmnEq can be theoretically determined by summation of the compliance tensors of model elements in the global coordinate system:

CijmnEq=CijmnorR+CijmnJ1+⋯+CijmnJk(22)

In summary, the total displacement Δεij can be scribed as:

Δεij=CijmnEqΔτmn(23)

2.2.3 Consideration of Powerhouse Arrangement

When discussing the impact of the horizontal deviation angle of in-situ stress on the stability of the surrounding rock in underground powerhouses, we examine the influence of the angle between in-situ stress direction and the axis of the underground powerhouse on rock stability. Although the axis of the underground powerhouse axis can be adjusted for optimal results, as illustrated by the yellow lines in Fig. 5, however, in our research, we maintain a fixed position for the underground powerhouse axis while systematically varying the horizontal deviation angle of the in-situ stress. We assign stress fields corresponding to different horizontal deviation angles to the model, allowing us to derive the stability variation patterns of the surrounding rock in response to changes in the horizontal deviation angle.

Figure 5: Relationship of the three factors and the axis of main powerhouse position

The strength parameters, dip direction, and dip angle of joint planes have a significant impact on the stability of surrounding rock in underground powerhouses. In practical engineering applications, the strength parameters of joint planes are generally fixed. However, the dip angle of joint planes may vary depending on the site selection for different underground powerhouses, directly influencing the deformation and stability of the surrounding rock. Additionally, the angle between the main axis of the underground powerhouse and the direction of the joint planes also plays a critical role in determining the stability of the surrounding rock. Variations in the dip direction can lead to substantial differences in stress distribution and failure patterns. The joint dip angle can change along the blue arc, while the in-situ direction and joint dip direction can vary along the red arc, as demonstrated in Fig. 5.

2.3 Flow Chart of the Nurmerical Analysis

This article divides four parts to make this research: Engineering geological analysis. Model Building, Numerical Analysis and Conclusion. The flowchart of overall technology is as Fig. 6.

Figure 6: Flow chart for getting the regularity between axis of underground powerhouse and joint dip angle, direction, in-situ stress direction in this study

The discrete fracture model and the equivalent continuous model have been utilized in recent years to address joint modeling. The discrete fracture model represents the distribution of fractures within a rock mass as a collection of finite-sized discrete fissures, depicted as line segments in two-dimensional space and flat disks in three-dimensional space. In rock mechanics, this model is applied to discrete elements that consist of fractures and intersections. An intersection occurs when one fracture intersects with another object, such as a straight line, a polygonal plane, a convex set of plane polygons, or another fracture. This model can effectively simulate the mechanical behavior of joint fractures within a rock mass.

The distribution of joints within a rock mass is highly complex, rendering it impractical to address engineering problems based solely on the mechanical and deformation characteristics of an individual joint. Unlike the discrete fracture model, the elastic matrix of the equivalent continuous model relies solely on the properties of the rock mass and joints, without taking element division into account. The elastic matrix is affected by parameters such as joint connectivity, spacing, width, and their random distribution characteristics, all of which must be considered. In practical applications, an incremental method can be utilized to effectively capture nonlinear characteristics. In Flac3D, the equivalent continuous modeling approach is adopted to simulate jointed rock masses using the Mohr-Coulomb and anisotropic ubiquitous joint models. This method treats the fractured medium as a homogenized continuum with equivalent mechanical properties, rather than explicitly modeling individual joints. The Mohr-Coulomb model effectively captures shear failure mechanisms in isotropic materials by balancing computational accuracy with practicality, particularly in large-scale underground excavations where explicit joint modeling would be prohibitively complex. In the case of anisotropic ubiquitous joint models, the equivalent continuum framework incorporates directional strength weakening through orientation-dependent parameters (e.g., dip, dip direction, joint cohesion), enabling a realistic simulation of stress-driven anisotropy caused by pervasive, closely spaced joints.

FLAC3D excels in such analyses due to its robust finite difference algorithm, which effectively manages nonlinear material behavior and large deformations characteristic of jointed rock. Its explicit solution scheme ensures numerical stability, even in highly heterogeneous or fractured media, while the built-in constitutive models (e.g., Mohr-Coulomb, ubiquitous joint) facilitate streamlined parameterization. Furthermore, FLAC3D’s capability to visualize stress redistribution and plastic zones in three dimensions provides critical insights for optimizing excavation layouts and reinforcement designs, making it a preferred tool for geomechanical stability assessments in complex geological environments. In this case, we utilize FLAC3D to analyze these models. In Flac3D, the analysis utilizing the Mohr-Coulomb and ubiquitous joint models is conducted as follows:

(1)   Initialize the model geometry and mesh using the model configure command.

(2)   Activate the desired constitutive model—model Mohr-coulomb for isotropic shear failure or model ubiquitous-joint for anisotropic behavior influenced by pervasive joints—and assign material properties (e.g., cohesion, friction angle, tensile strength, and joint orientation parameters) via the property command.

(3)   Define the initial in-situ stress field (initial stress) and establish boundary conditions to replicate geological constraints.

(4)   Simulate excavation stages by removing elements with model null, then execute the iterative solver using solve to compute stress redistribution and failure mechanisms.

(5)   Visualize results, such as displacement fields, plastic zone distributions, and joint shear failure patterns, using plotting or history commands to enable a comprehensive evaluation of stability.

3  Result

3.1 Development of the Geological Model

To ensure the practical applicability of this model, the paper incorporates detailed information from a specific underground powerhouse project for analysis. The surrounding rock in the underground powerhouse area consists of fresh sericite phyllite, which is characterized by well-developed foliation and joint sets. These joints are primarily aligned with the foliation and intersect the powerhouse axis at angles ranging from 35° to 50°, with dip angles between 30° and 40°. Together with two additional joint sets, these features create potential instability blocks in certain areas of the cavern. Overall, the rock mass surrounding the underground powerhouse is relatively intact, although some local areas exhibit reduced integrity. Within the depth range of the powerhouse, the maximum horizontal principal stress values range from 13.09 to 17.77 MPa, oriented between NW4° and NW13°, suggesting a possibility of stress-induced failure in the rock mass along the tunnel walls. To fully reflect the topographical, geological, and geomorphological characteristics, the computational range of the model is determined as follows based on the analyzed data:

(1)   The bottom elevation of the model is set at 148.5 m, while the top elevation is set at 203.7 m.

(2)   The model extends in the length direction from −40.00 to 142.5 m, and in the width direction from −12.5 to 132.5 m.

(3)   To account for the stability of the surrounding rock mass, the vertical boundaries on both sides are extended. The length extends from −100 to 200 m, while the width extends from −100 to 180 m.

(4)   In the vertical direction, it extends from 88.7 to 278.7 m.

This geological model measures a total length of 300 m, a width of 280 m, and a height of 190 m. It comprises 687,917 elements and 445,215 nodes. The underground powerhouse model is depicted in Fig. 7.

Figure 7: Overall modeling of underground powerhouse model for analysis

In order to enhance the analysis of the influence of three key factors on the stability of surrounding rock, this paper establishes observation points at various elevations within the plant model, as illustrated in Fig. 8 below, in which Ci is the apex of hole, Ai, Bi, Di, Ei is pore wall of the model.

Figure 8: Observation points

The observation points (Ci, Ai, Bi, Di, Ei) are strategically positioned to address critical structural and geological risks within the underground powerhouse model. Ci, located at the apex of the cavern (e.g., coordinates such as (0, 10, 203.65)), is chosen to monitor vertical stress concentrations and potential roof instability, particularly relevant due to the intersecting foliation and joint sets that form instability blocks. Points Ai, Bi, Di, Ei are placed along the cavern walls at varying elevations (e.g., 168 and 200 m) to capture stress redistribution and deformation near zones where joints intersect the powerhouse axis at angles of 35°–50°, which are prone to shear or tensile failure. The elevation range of these points spans from 88.7 to 278.7 m vertically across the model; this ensures comprehensive coverage of stress gradients influenced by a non-uniform horizontal principal stress field (13.09–17.77 MPa oriented NW). Redundant placements, such as multiple points, account for the heterogeneity of the sericite phyllite rock mass and validate spatial variations in stress scatter caused by foliation and joint systems. Furthermore, observation points near extended boundaries of the model (e.g., x = 0, y = 10) facilitate an assessment of boundary effects and stress equilibrium thereby ensuring computational accuracy. Collectively, these strategically placed points enable a thorough evaluation of stress-induced failures alongside joint interactions and overall rock mass integrity; this aligns with our objective for modeling real-world geological complexities as well as mechanisms underlying instability.

3.2 Initialization of the Ground Stress

Due to the inherent non-uniformity of the material, direct generation of stress may not accurately represent the actual conditions. In this study, there is a notable scatter in the range of initial stress, indicating that both the magnitude and direction of stress can vary significantly across different locations and depths. This non-uniformity may arise from the natural heterogeneity of the rock mass, complexities within geological structures, and influences from historical geological events. To address this issue, this study employs an approach utilizing the scatter in initial stress ranges as follows: (1) Calculation of stress in elastic materials: To achieve rapid convergence, an elastic parameter is initially applied for calculations. Subsequently, various elastic parameters are utilized after completing these calculations to determine the stress state of the rock body under elastic conditions; (2) Calculation of stress in elastic materials: Furthermore, the material model is changed to a plastic model, and the strength parameter is calculated to obtain the stress; (3) Final correction of initial stress: Based on this, the stress obtained from the test is applied to the rock body. Since the stress can only be applied in the XYZ coordinate direction, it is necessary to perform transformation calculations. According to the direction of the principal stress, the elastic stress tensor conversion formula is used to calculate and determine the XYZ direction stress and its gradient change, and then ‘zone initialize stress’ method is used to update the stress.

As illustrated in Fig. 9, it is evident that initial stresses at identical elevation levels remain consistent across all directions. Therefore, when adjusting for dip and strike variations on joint surfaces, it becomes unnecessary to consider any coupling effects between these two factors regarding induced stresses. The impact of inclination angles related to stress levels on underground powerhouse axis layout will be further discussed in Section 4.3 of Section 4.

Figure 9: Stress distribution of initialization of in-situ stress: (a) minimum principal stress of initialization of in-situ stress; (b) maximum principal stress of initialization of in-situ stress

3.3 Verification of Observation Points and This Research

In Fig. 8, we arranged 20 observation points to investigate these three factors. We now proceed to verify the rationality and representativeness of this arrangement. Based on Table 2, specific engineering parameters for the project are provided. Utilizing the anisotropic joint model in Flac3D, we conducted stepwise excavation to obtain the distribution of maximum principal stress, volumetric stress, and total displacement at the model’s crown (elevation 203.65 m) and sidewalls (elevations 200 and 168 m), as shown in Fig. 10.

Figure 10: Displacement and stress distribution: (a–c) are zone maximum principle stress of 203.65, 200, 168 elevation, (d–f) zone volumetric stress of 203.65, 200, 168 elevation, (g–i) zone displacement magnitude of 203.65, 200, 168 elevation

It can be observed that displacement and stress are primarily concentrated in the sidewalls and crown of the main powerhouse, confirming that the arrangement of observation points is both reasonable and representative.

Additionally, we have validated the feasibility and significance of this study. As shown in Fig. 11, different joint orientations, dip angles, and horizontal deflection angles of in-situ stress significantly affect the stability of the surrounding rock in the underground powerhouse. Although it is challenging to alter the orientation, dip angle, and direction of in-situ stress in practical engineering, the position and alignment of the underground powerhouse can be adjusted. By applying the conclusions of this study, the stability of the surrounding rock in the underground powerhouse can be enhanced. Through in-depth analysis of data on geological structure, hydrogeological conditions, and rock stability, we can more fully assess the suitability of potential sites and provide a scientific basis for future engineering construction. Therefore, this research is both necessary and valuable can provide a reference for the selection of axis and location of underground powerhouse.

Figure 11: Displacement change by different factors of the model of 203.65 m elevation: (a) displacement change by dip angle; (b) displacement change by dip direction; (c) displacement change by in-situ stress direction

4  Discussion

4.1 Joint Plane Dip Angle

The dip angle of the joint plane plays a crucial role in determining the stability of the surrounding rock. This angle influences the likelihood of sliding instability, particularly when it is steep and closely aligned with the direction of gravity, thereby increasing the risk of sliding. Additionally, the dip angle affects the shear strength of the rock mass; larger dip angles may lead to a reduction in shear strength, consequently heightening instability risks. Furthermore, variations in dip angles can modify internal stress distribution within the rock mass, with steeper angles potentially causing stress concentration that could result in rock mass failure.

The interplay between the dip angle and other structural planes—such as bedding planes and fault planes—also significantly impacts surrounding rock stability. Multiple intersecting structural planes with unfavorable dip angles can create unstable rock blocks. Therefore, it is imperative to thoroughly consider the dip angle of joint planes during engineering design and stability analysis for underground powerhouses. In subsequent discussions, we will delve into an in-depth examination of how the dip angle of joint planes correlates with displacement within surrounding rocks.

From Fig. 12a,b, it can be inferred that as the earth stress inclination increases, the deformation reduces suddenly before possibly increasing gradually, however the rise is less than the decrease. The increasing shift of failure mechanisms is the cause of this nonlinear response. At lower stress inclinations, shear slip primarily occurs along the joint plane because normal stress confinement is diminished, causing abrupt deformation. The alignment between the joint plane orientation and the major stress direction changes the stress partitioning as the inclination rises, favoring elastic strain buildup over plastic slip and leading to mild increments of deformation. The minima eventually emerge at the front position as displacement increases. When the joint plane inclination is 70°, the minima appear in Bi; when it is 55°, the minima appear in Ci; when it is 40°, the minima appear in Di; and when it is 5°, the minima appear in Ai and Ei. The model’s displacement is at its ideal value when the joint plane’s inclination falls between 65° and 70°. For low-dip joints (less than 40°), stress concentration zones mostly form at the junctions of joint planes and excavation boundaries, especially close to the arch shoulders (Ai and Ei regions). while shifting toward the sidewalls (Bi and Ci regions) as joint inclination increases. This spatial redistribution of stress concentrations correlates with the observed displacement minima positions, as high-dip joints redirect compressive stresses deeper into the rock mass rather than concentrating them near excavation surfaces.

Figure 12: Variation of displacement and plastic zone with joint dip angle: (a) original diagram of displacement variation with joint dip angle; (b) post-processing diagram of displacement variation with joint dip angle; (c) plastic zone with joint dip angle

Relationship between the dip angle of the joint plane and changes in the plastic zone aligns with the shoulder-shaped curve in Jaeger’s single weak plane criterion, as is demonstrated in Fig. 12c. The mechanical basis for this correspondence lies in the critical stress ratio required for joint activation—high-dip joints (55°–70°) approach the frictional lock-up angle under typical in-situ stress conditions, thereby restricting plastic zone expansion through interlocking effects. This further indicates that high dip angle structures result in less damage to the surrounding rock during the excavation of the underground powerhouse, whereas low dip angle structures lead to a larger volume of rock damage. The intensified damage at low angles arises from extended stress concentration belts along sub horizontal joints, facilitating tensile crack propagation parallel to the maximum principal stress trajectory near excavation peripheries.

4.2 Joint Plane Dip Direction

The stability of the surrounding rock depends on the joint plane’s dip direction since it establishes the main sliding direction. Sliding instability is more likely if the dip direction aligns with the major direction of in-situ stress. The redistribution of stress within the rock mass is also influenced by the dip direction. The stability of the nearby rock is impacted by the different stress concentration and release patterns caused by different dip directions. The geometry and stability of rock blocks are also influenced by the angle formed by the dip direction and other structural planes. The risk of collapse is increased by the formation of unstable wedge-shaped or columnar rock blocks caused by specific combinations of dip directions. Finally, understanding the dip direction helps optimize engineering design by selecting appropriate excavation directions and support measures to minimize negative impacts on rock stability.

Therefore, the dip direction of the joint plane must be thoroughly considered in engineering design and stability analysis. Now we take joint plane dip direction into consideration.

The relationship between the dip direction of the joint plane and the deformation of the surrounding rock, as demonstrated in Fig. 13a,b, indicates that stability is enhanced when the dip direction aligns parallel to the underground powerhouse axis or forms an angle of less than 40°. This advantageous behavior results from the kinematic compatibility between principal stress trajectories and joint orientation: low-angle or parallel configurations suppress block rotation and slip by keeping joint planes subparallel to the direction of maximum compressive stress, which minimizes shear stress mobilization along discontinuities. When the dip direction exceeds 40° relative to the powerhouse axis, the misalignment causes uneven stress distribution. The resulting torque impact promotes increasing joint dilatation and shear band formation, jeopardizing chamber stability. Fig. 13c further reveals that stress concentration effects intensify with distance from the excavation boundary, creating distinct damage progression patterns. Even in favorable dip orientations (parallel or less than 40°), localized deformation along the chamber periphery increases slightly, but the amplitude is still limited because of stress arching effects that shift loads into deeper rock masses. Crucial stress concentration zones, which correlate with observed asymmetric deformation patterns, notably appear as tensile-dominated regions at the downstream shoulder (for dip directions 30°–40°) and transform into compressive shear belts close to the upstream sidewall when dip directions surpass 50°. The relationship between dip direction and plastic zone evolution exhibits threshold-controlled characteristics: By allowing favorably oriented joints to interlock, configurations below 40° prevent significant damage, but angles above 40° cause the volume of the plastic zone to increase quickly and linearly. The critical angle for stress-induced joint reactivation is represented by this transition, where a greater misalignment between the joint planes and the highest primary stress lowers frictional resistance and permits progressive shear failure. The development of continuous shear surfaces connecting several joint sets is the cause of the increased damage at high angles, especially in areas where favorably aligned discontinuities and excavation-induced tensile pressures meet. Plastic zones exhibit a distinctive spatial progression, beginning at stress concentration points along the crown and inverting for parallel dip directions. As dip direction angles increase, the zones migrate toward sidewalls, and when angles exceed 60°, they form interconnected failure wedges. These mechanisms collectively explain the nonlinear stability degradation observed with increasing dip direction angles.

Figure 13: Variation of displacement and plastic zone with joint dip direction: (a) original diagram of displacement variation with joint dip direction; (b) post-processing diagram of displacement variation with joint dip direction; (c) plastic zone with joint dip

4.3 Direction of In-Situ Stress Levels

The horizontal deviation angle of in-situ stress is a critical factor in determining the stability of the surrounding rock, due to its multifaceted influence on the stress distribution and deformation behavior within the rock mass. This deviation angle directly affects the concentration of stresses, for larger angles lead to increased stress concentrations that heighten the likelihood of rock mass failure. When the horizontal deviation angle is significant, it can intensify the stress concentrations, particularly in localized areas, thereby elevating the risk of structural instability. Moreover, the horizontal deviation angle plays a pivotal role in altering the deformation patterns of the surrounding rock. Furthermore, the interplay between the horizontal deviation angle and the orientation of joints and fractures within the rock mass is a key factor that can significantly influence the stability of the surrounding rock. A larger horizontal deviation angle can align more closely with these pre-existing discontinuities, increasing the likelihood of instability and failure along these planes. This interaction underscores the importance of accurately assessing and considering the horizontal deviation angle in the context of rock mass behavior. From an engineering perspective, a comprehensive understanding of the horizontal deviation angle is indispensable for optimizing the design and construction processes.

Then we take deviation angle of in-situ stress levels into consideration.

The relationship between the horizontal deviation angle of in-situ stress and the deformation of surrounding rock, as demonstrated in Fig. 14a,b, exhibits a distinctive “increase-decrease-increase” trend. This triphasic behavior arises from the evolving interaction between principal stress trajectories and excavation-induced stress redistribution. The deformation minima observed at 0° and 60° deviation angles correspond to configurations where the powerhouse axis aligns either parallel to (0°) or approximately 30° oblique to (60°) the maximum principal stress direction. At these critical orientations, the kinematic constraints imposed by the stress field geometry suppress shear strain localization through optimized stress transfer mechanisms. Specifically, the 0° alignment maximizes coaxially between excavation geometry and principal stresses, minimizing differential strain accumulation, while the 60° configuration induces beneficial stress rotation that activates joint interlocking effects. Conversely, the deformation maximum between 30° and 40° reflects a stress-path instability threshold where progressive misalignment induces asymmetric stress partitioning. This transitional regime promotes strain localization through two concurrent mechanisms: (1) tensile stress amplification along the upstream arch shoulder due to principal stress rotation exceeding the rock mass dilatancy angle, and (2) shear band nucleation at the downstream sidewall resulting from critical mobilization of joint friction angles.

Figure 14: Variation of displacement and plastic zone with angle of direction of in-situ stress: (a) original diagram of displacement variation with angle of direction of in-situ stress; (b) post-processing diagram of displacement variation with angle of direction of in-situ stress; (c) plastic zone variation with angle of direction of in-situ stress

Fig. 14c further elucidates that the horizontal deviation angle governs nonlinear deformation progression through stress concentration zone migration. As the angle increases from 0° to 30°, stress concentrations shift from symmetric arch abutment zones to asymmetric crown-sidewall junctions, triggering rapid deformation escalation. Beyond 30°, the development of self-equilibrating compressive arches along favorably oriented discontinuities mitigates strain accumulation, explaining the subsequent deformation decline. The minimal deformation states at 0° and 60° ± 5° are mechanically distinct: the former achieves stability through optimal principal stress alignment, while the latter benefits from stress shadow effects created by oblique joint plane orientations that deflect compressive stresses into intact rock bridges.

Plastic zone evolution demonstrates threshold-dependent behavior, with negligible volume changes below 20° deviation transitioning to rapid expansion beyond this critical angle. This bifurcation corresponds to the transition from contained brittle fracture at low angles (dominated by tensile crack propagation along the vault) to pervasive shear failure at higher angles (characterized by through-going shear surfaces connecting the invert and sidewalls). The stress concentration zones exhibit progressive spatial reorganization—initiating as tensile-dominated regions near the downstream spandrel at 0°–20°, evolving into compressive shear belts along the upstream haunch at 30°–40°, and finally stabilizing as distributed microcrack networks at 60°. These spatial patterns correlate with observed deformation minima, as the 60° configuration promotes stress redistribution into deeper rock masses through inclined stress arching mechanisms, effectively decoupling near-excavation damage from global stability. The combined effects of stress rotation kinematics and damage zone localization explain the dual stability optima observed in the deviation angle-deformation relationship.

4.4 Limitations

Certain geomechanical assumptions, notably those pertaining to intact rock strength, joint surface characteristics, and in-situ stress regimes, are necessary to draw the conclusion that joint planes with dip angles between 65° and 70° produce ideal displacement circumstances. The assumption of constant friction angles along discontinuities is a crucial drawback. For example, the so-called “optimal” dip angle range may go downward if joint surfaces show lower friction angles (such as clay-filled or weathered interfaces). The interlocking effects seen under greater friction settings may be negated in such circumstances, as even high-dip joints may undergo shear mobilization as a result of decreased frictional resistance. The ideal dip angle range, on the other hand, may be extended beyond 70° by rougher joints with higher friction angles because asperity resistance may preserve stability at steeper orientations. This is consistent with the critical slip angle for discontinuities, which is directly impacted by changes in frictional strength according to Bayerle’s law. This generalization is further complicated by joint persistence. The study ignores situations where discontinuous fractures predominate because it implicitly assumes fully permanent joints. In contrast to field settings with interconnected fracture networks, short, non-persistent joints may not be able to produce continuous failure surfaces at high dip angles (65°–70°), artificially increasing stability. Even at ideal dip angles, stress concentrations at fracture ends may cause tensile crack propagation in rock masses with high joint persistence, paradoxically enhancing deformation. This phenomena is especially significant for brittle rock types, whose failure processes are governed by fracture toughness rather than only frictional slide. Furthermore, in contrast to the stability assumptions of the single-plane model, the interaction between numerous joint sets—which is typical in actual geological settings—may produce kinematic freedom for block rotation at high dip angles. The influence of confining pressure gradients introduces another layer of complexity. The conclusion assumes homogeneous stress redistribution post-excavation, whereas in stratified or anisotropic stress fields, high-dip joints may align with principal stress rotations, reactivating slip surfaces that would otherwise remain stable. For example, in deep excavations with significant horizontal stress anisotropy, joints dipping 65°–70° parallel to the intermediate principal stress direction could become prone to buckling failures—a mechanism not captured in conventional slip-based analyses. Furthermore, time-dependent behaviors such as creep in weak interlayers or stress corrosion cracking in intact rock bridges could progressively degrade the perceived optimality of high-dip structures, particularly in humid environments or seismically active regions. These factors collectively suggest that while the identified dip angle range provides useful guidance under idealized conditions, field applications require rigorous site-specific analysis incorporating 3D joint network geometry, scale effects, and lithological heterogeneity.

5  Conclusion

(1)   During underground powerhouse excavation, a high inclination plane will result in less damage to the surrounding rock than a low inclination plane, which causes more damage to the surrounding rock. The connection between the joint plane’s dip angle and the deformation of the surrounding rock in this subterranean powerhouse shows that dip angles between 60° and 70° increase the stability of the surrounding rock. Low dip angles cause the surrounding rock to become unstable, whereas high dip angles help the rock stay stable.

(2)   The surrounding rock will be more stable when the joint dip direction is parallel to or at a slight angle to the underground powerhouse’s axis. In this subterranean power plant, the relationship between joint dip direction and surrounding rock deformation is demonstrated that stability is optimized when the dip direction is parallel to the powerhouse axis or forms an angle of less than 40°.

(3)   The surrounding rock will be more stable when the underground powerhouse’s axis is parallel to or slightly angled with the direction of highest primary stress. The horizontal deviation angle of in-situ stress and the surrounding rock deformation in this subterranean powerhouse exhibit an increase-decrease-increase pattern, with the least amount of deformation taking place between 0° and 60°. Deformation reaches its maximum between 30° and 40°. When the powerhouse axis aligns at 60° ± 5°, minimizes deformation and plastic zone volume, or is parallel to or forms a minor angle with the greatest principal stress direction, stability is maximized.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: Junfei Dai: The main author of the article and the provider of ideas; Weijiang Chu: Data provider, experimental processing, data processing; Xiangming Ge: Experimental processing, data processing; Qingxiang Meng: Paper reviser and refiner. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: Not applicable.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

References

1. Meng Q-X, Wang H-L, Xu W-Y, Chen Y-L. Numerical homogenization study on the effects of columnar jointed structure on the mechanical properties of rock mass. Int J Rock Mech Min Sci. 2019;124(3):104127. doi:10.1016/j.ijrmms.2019.104127. [Google Scholar] [CrossRef]

2. Bizjak KF, Petkovšek B. Displacement analysis of tunnel support in soft rock around a shallow highway tunnel at Golovec. Eng Geol. 2004;75(1):89–106. doi:10.1016/j.enggeo.2004.05.003. [Google Scholar] [CrossRef]

3. Zhu Q, Li T, Ran J, Du Y, Zhang H, Jiang H. Model test on creep deformation and failure characteristics of soft rock roadways. Eng Fail Anal. 2022;141:106670. doi:10.1016/j.engfailanal.2022.106670. [Google Scholar] [CrossRef]

4. Liang Z, Xue R, Xu N, Dong L, Zhang Y. Analysis on microseismic characteristics and stability of the access tunnel in the main powerhouse, Shuangjiangkou hydropower station, under high in situ stress. Bull Eng Geol Environ. 2020;79(6):3231–44. doi:10.1007/s10064-020-01738-6. [Google Scholar] [CrossRef]

5. Tiwari G, Latha GM. Shear velocity-based uncertainty quantification for rock joint shear strength. Bull Eng Geol Environ. 2019;78(8):5937–49. doi:10.1007/s10064-019-01496-0. [Google Scholar] [CrossRef]

6. Chen L, Jin A, Wu S, Chu C, Li X. Numerical study on spalling failure of rock surrounding deep buried tunnel based on DEM. Comput Geotech. 2022;145(3):104653. doi:10.1016/j.compgeo.2022.104653. [Google Scholar] [CrossRef]

7. Tiedtke F, Konietzky H, Magri F. A novel DFN-DEM approach to simulate long-term behavior of crystalline rock under effects of glacial climate conditions. Deep Resour Eng. 2024;1(1):100002. doi:10.1016/j.deepre.2024.100002. [Google Scholar] [CrossRef]

8. Wang M, Cai M. Simulation of time-dependent response of jointed rock masses using the 3D DEM-DFN modeling approach. Int J Rock Mech Min Sci. 2025;188(1):106062. doi:10.1016/j.ijrmms.2025.106062. [Google Scholar] [CrossRef]

9. Wei Y, Feng Y, Tan Z, Yang T, Li X, Dai Z, et al. Borehole stability in naturally fractured rocks with drilling mud intrusion and associated fracture strength weakening: a coupled DFN-DEM approach. J Rock Mech Geotechnical Eng. 2024;16(5):1565–81. doi:10.1016/j.jrmge.2023.07.012. [Google Scholar] [CrossRef]

10. Ding X, Niu X, Pei Q, Huang S, Zhang Y, Zhang C. Stability of large underground caverns excavated in layered rock masses with steep dip angles: a case study. Bull Eng Geol Environ. 2019;78(7):5101–33. doi:10.1007/s10064-018-01440-8. [Google Scholar] [CrossRef]

11. Dong M, Zhang F, Lv J, Fei Y, Li Z. Study of stability influencing factors of excavated anti-dip rock slope. KSCE J Civil Eng. 2020;24(8):2293–303. doi:10.1007/s12205-020-1412-4. [Google Scholar] [CrossRef]

12. Meng Q, Xue H, Zhuang X, Zhang Q, Zhu C, He B, et al. An IFS-based fractal discrete fracture network for hydraulic fracture behavior of rock mass. Eng Geol. 2023;324:107247. doi:10.1016/j.enggeo.2023.107247. [Google Scholar] [CrossRef]

13. Jaeger JC, Cook NG, Zimmerman R. Fundamentals of rock mechanics. Hoboken, NJ, USA: John Wiley & Sons, Inc.; 2009. [Google Scholar]

14. Li X, Chen J, Ma C, Huang L, Li C, Zhang J, et al. A novel in-situ stress measurement method incorporating non-oriented core ground re-orientation and acoustic emission: a case study of a deep borehole. Int J Rock Mech Min Sci. 2022;152(6):105079. doi:10.1016/j.ijrmms.2022.105079. [Google Scholar] [CrossRef]

15. Li Z, Liu H, Dai R, Su X. Application of numerical analysis principles and key technology for high fidelity simulation to 3-D physical model tests for underground caverns. Tunnelling Undergr Space Technol. 2005;20(4):390–9. doi:10.1016/j.tust.2005.01.004. [Google Scholar] [CrossRef]

16. Ma H-P, Daud NNN, Yusof ZM, Yaacob WZ, He H-J. Stability analysis of surrounding rock of an underground cavern group and excavation scheme optimization: based on an optimized DDARF method. Appl Sci. 2023;13(4):2152. doi:10.3390/app13042152. [Google Scholar] [CrossRef]

17. Malmgren L, Saiang D, Töyrä J, Bodare A. The excavation disturbed zone (EDZ) at Kiirunavaara mine, Sweden—by seismic measurements. J Appl Geophy. 2007;61(1):1–15. doi:10.1016/j.jappgeo.2006.04.004. [Google Scholar] [CrossRef]

18. Qian D, Zhang N, Pan D, Xie Z, Shimada H, Wang Y, et al. Stability of deep underground openings through large fault zones in argillaceous rock. Sustainability. 2017;9(11):2153. doi:10.3390/su9112153. [Google Scholar] [CrossRef]

19. Schmalholz SM, Moulas E, Plümper O, Myasnikov AV, Podladchikov YY. 2D hydro-mechanical-chemical modeling of (de)hydration reactions in deforming heterogeneous rock: the periclase-brucite model reaction. Geochem Geophy Geosy. 2020;21(11):e2020GC009351. doi:10.1029/2020gc009351. [Google Scholar] [CrossRef]

20. Zhang H, Meng X, Yang G. A study on mechanical properties and damage model of rock subjected to freeze-thaw cycles and confining pressure. Cold Reg Sci Technol. 2020;174:103056. [Google Scholar]

21. Zhang Q, Liu C, Duan K, Zhang Z, Xiang W. True three-dimensional geomechanical model tests for stability analysis of surrounding rock during the excavation of a deep underground laboratory. Rock Mech Rock Eng. 2020;53:517–37. [Google Scholar]

22. Zhu G-Q, Feng X-T, Zhou Y-Y, Li Z-W, Fu L-J, Xiong Y-R. Physical model experimental study on spalling failure around a tunnel in synthetic marble. Rock Mech Rock Eng. 2020;53(2):909–26. doi:10.1007/s00603-019-01952-z. [Google Scholar] [CrossRef]

23. Eberhardt E. Numerical modelling of three-dimension stress rotation ahead of an advancing tunnel face. Int J Rock Mech Min Sci. 2001;38(4):499–518. doi:10.1016/s1365-1609(01)00017-x. [Google Scholar] [CrossRef]

24. Jia C, Li S, Fan C, Rong H, Yang L, Pu Z. Numerical simulation of roadway deformation and failure under different degrees of dynamic disturbance. Sci Rep. 2022;12(1):20017. doi:10.1038/s41598-022-24128-2. [Google Scholar] [PubMed] [CrossRef]

25. Jiang Y, Zhou H, Lu J, Gao Y, Zhang C, Chen J. Analysis of stress evolution characteristics during TBM excavation in deep buried tunnels. Bull Eng Geol Environ. 2019;78(7):5177–94. doi:10.1007/s10064-019-01466-6. [Google Scholar] [CrossRef]