## 3 Model Establishment  
### Determination of Parameters  
First, cylindrical specimens of sandstone-like rocks were fabricated. The dimensions were set to a height of 100 mm and a diameter of 50 mm, according to the recommendations of the International Society of Rock Mechanics, as shown in Figure 5a. The compression rate was set to 0.01 mm/min. The experiment stopped when the specimen was obviously broken and the loading curve exceeded the peak value. Other variables were kept constant during the test, and no confining pressure was set. Second, a three-dimensional uniaxial compression numerical model was developed with the same test dimensions. Under the same mesh size, compared with other shapes (wedge, hexahedron, etc. ), the tetrahedral mesh could divide the model more finely, resulting in finer simulation results. In addition, since the simulated crack propagation and mode depended on the initial finite mesh, the tetrahedral mesh could have produced oblique failure, so the crack generation effect was better. Consequently, we referred to previous similar studies ([PERSON] et al. [25], [PERSON] et al. [26], etc.) and finally selected the tetrahedral mesh. In the numerical simulation of the uniaxial compression of a cylinder, there were a total of 121,450 elements, including 76,230 cohesive elements (COH3D6) and 45,220 solid elements (C3D4). The use of the finite element plug-in to insert the cohesive element, and to set the loading plate and the stiffness of the loading plate was much greater than the stiffness of the specimen, to ensure that the load could be applied uniformly on the top of the specimen, as shown in Figure 5b. The lower loading plate was fixed, and the upper loading plate was moved downward with a speed of 0.1 mm/s, until the specimen was  
Figure 4: The schematic diagram of crack evolution: (**a**) schematic representation of tensile shear crack evolution; (**b**) schematic representation of compressive shear crack evolution; (**c**) tension crack evolution; and (**d**) failure of the removal of the cohesive element. damaged. If the loading speed of the numerical simulation is the same as the test, it is often necessary to set a smaller time step to ensure the convergence of the calculation, and a longer calculation time is required. Therefore, without affecting the simulation effect, the loading speed of the simulation was set to be ten times higher than the loading speed of the test. There are two methods for determining the parameters of the cohesive element. The first method is error reduction: the elastic modulus, Poisson's ratio, and other rock mass parameters are obtained through experiments. The creation of two uniaxial compression models that the cohesive elements have been inserted. Numerical simulations are carried out under the same boundary conditions. When the difference between the two stress fields tends to be stable, it is considered that the parameters are determined. The second method is the experimental comparison method: multiple sets of parameters are used for simulation, and the numerical simulation is compared with the experimental results. If the two are similar, the parameters used for the numerical simulation are considered feasible. In this study, we used the second method. In order to obtain accurate fracture damage parameters in the shear numerical simulation of joint planes, we referred to the research of [PERSON] et al. [25], and the uniaxial compression test was simulated by the method of the global embedding of the cohesion element, which was used to calibrate the parameters of cohesive elements. The macroscopic physical and mechanical parameters, such as [PERSON]'s modulus and Poisson's ratio, obtained in the uniaxial compression test, were brought into the numerical model as solid element parameters. The quadratic nominal stress criterion (QUADS) was used for the material behavior of the cohesive elements, to judge the onset of damage evolution, and the BK law was selected for damage propagation.