## 2 Numerical Simulation Method
### Numerical Model  
The layered structure of shale imparts significant anisotropic characteristics, and the development position, thickness, and strength of the bedding are highly variable. Therefore, accurately reproducing the bedding structure of shale in numerical simulations is crucial.  
A layered rock block was generated using a MATLAB (R2023a) program to simulate the random distribution of shale bedding. Shale was sampled at different angles within the rock block, simulating the coring process of actual rock samples. This resulted in a shale model with varying sampling angles, as shown in the figure. This process reproduces the bedding distribution within the rock and allows for the simulation of rock sample characteristics under different experimental conditions by adjusting the sampling angle.  
Based on the constructed rock block, the model identified different regions using threshold segmentation, and the matrix and bedding characteristics were assigned accordingly. Using numerical simulation software, the partitioning of the matrix and bedding in the shale model was implemented, enabling a more accurate representation of the layered structure and anisotropic behavior of shale. Shale, a sedimentary rock with high heterogeneity, consists of clay minerals, quartz, organic matter, and other components exhibiting distinct mechanical responses. Figure 2 presents the layered structure model after meshing with the RFPA3D. The color variation in the RVE within the model reflects the distribution of the intensity heterogeneity. As mentioned in the previous section, the intensity parameters in RFPA3D follow a Weibull distribution function, effectively simulating the heterogeneous nature of rocks.  
Shale is subjected to various stress states during underground storage and related engineering processes. To investigate the anisotropic fracture characteristics of shale under different stress conditions, this study constructs uniaxial compression, uniaxial tension, and triaxial compression models of randomly layered shale:  
(1) Uniaxial compression:  
According to the method in Figure 2, we extracted four cylindrical models with typical bedding angles (0\({}^{\circ}\), 30\({}^{\circ}\), 60\({}^{\circ}\), 90\({}^{\circ}\)) from the complete model as one group. We extracted a total of three groups with different bedding distributions to ensure the universality of the results. The models are heterogeneous cylinders with a bottom radius of 25 mm and a height of 100 mm, with circular homogeneous plates of 5 mm thickness at both ends, consisting of approximately 370,000 elements. The loading is applied as a downward compressive axial load, with a stepwise loading rate of 0.02 mm per step. The specific model is shown in Figure 3.  
Figure 2: Three-dimensional shale numerical models.  
(2) Direct tension:  
The same model as the uniaxial compression test was used, but the upper and lower plates were removed. The model consists of approximately 340,000 elements. The loading is applied as an upward tensile axial load.  
(3) Triaxial compression:  
Considering the extended simulation time under confining pressure, a set of rectangular prism models was selected for triaxial compression simulations at four bedding angles (0\({}^{\circ}\), 30\({}^{\circ}\), 60\({}^{\circ}\), and 90\({}^{\circ}\)). Each model, with a base edge length measuring 50 mm and a height measuring 100 mm, with square caps at both ends, consists of approximately 480,000 elements. The confining pressures were set at 7.5 MPa, 15 MPa, and 22.5 MPa, while the axial loading conditions remained consistent with those in the uniaxial compression models. Notably, the bedding distribution of the rectangular models used for triaxial compression was identical to that of the first group in the uniaxial compression test.