## 2 Numerical Simulation Model  
### GPU-Based Discrete Element Method  
The discrete element method (DEM) was first introduced by [PERSON] and [PERSON] [46] to study the behavior of granular assemblies. By introducing bonds, DEM was successfully used to simulate the behaviors of cohesive materials [47], e.g., soil and rock. For bonded particle models, materials are assembled from many mass particles linked to other particles via springs. Contact forces are generated when particles penetrate, slide, and rotate with other contact particles. The sum of contact force, shear force, external force, and gravity on a particle constitutes the combined force of the particle. The dynamic evolution of this particle system can be simulated by integrating equilibrium equations using [PERSON]'s physics laws and explicit time algorithms [46; 47]. In the explicit time algorithm, the step time is set very small to keep numerical stability. In addition, the velocity and acceleration of a particle are very low so that they can be assumed to be constant within a time step, allowing the calculation of the displacement of each particle. For quasi-static problems, the adaptive damping scheme can be applied to achieve equilibrium and can be written as follows:  
\[f_{i}=f_{i}-c\cdot sign(\dot{u})|f_{i}| \tag{1}\]  
where \(c\) is the local damping coefficient, \(f\) and \(\dot{u}\) are the force vector and velocity vector for the particle \(i\), respectively. Two contact particles are linked by a breakable bond. In this work, the widely used enhanced parallel bonded model was used for modelling rocks. In this contact model, the bond is represented as a finite area beam and includes seven parameters: normal and shear stiffness per unit area \(\overline{k}_{n}\), and \(\overline{k}_{s}\), tensile strength \(\sigma_{t}\)cohesion strength \(c\), friction angle \(\varphi\), a moment contribution coefficient \(\beta\), and bond-radius multiplier \(\lambda\), which defines the bond radius \(\overline{R}=\lambda min(R_{A},R_{B})\); with \(R_{A}\) and \(R_{B}\) being the radii of the two bonded particles. The relative motions between two contact particles cause increments of contact forces and moments due to the contact stiffness. The particle movements and the resultant forces and moments follow [PERSON]'s law of motion. The change of contact forces and moments due to the relative particle movements are determined, respectively, by Equations (2) and (3):  
\[\Delta\overline{F_{n}}=-\overline{k}_{n}\overline{A}\Delta U_{n},\ \Delta \overline{F_{s}}=-\overline{k}_{s}\overline{A}\Delta U_{s} \tag{2}\]  
\[\Delta\overline{M}_{n}=-\overline{k}_{n}\overline{J}\Delta\theta_{n},\Delta \overline{M}_{s}=-\overline{k}_{s}\overline{I}\Delta\theta_{s} \tag{3}\]  
where \(\overline{F_{n}}\), \(\overline{F_{s}}\), and \(\overline{M}_{n}\), \(\overline{M}_{s}\) are the contact forces and moments at the contact zone center, respectively, in the normal (n) and shear (s) directions; \(U_{n}\), \(U_{s}\) and \(\overline{M}_{n}\), \(\overline{M}_{s}\) are the relative displacements and rotations between the two bonded particles, respectively, in the normal (n) and shear (s) directions; and \(\overline{A}\), \(\overline{I}\) and \(\overline{J}\) are, respectively, the area, moment of inertia, and polar moment of inertia of the bond cross-section, determined by Equation (4):  
\[\overline{A}=\pi\overline{R}^{2},\ \overline{I}=\frac{1}{4}\pi\overline{R}^{ 4},\ \overline{J}=\frac{1}{2}\pi\overline{R}^{4} \tag{4}\]  
Elements are originally bonded to each other.