## 2 Methodology
### Mohr-Slip Joint Model of Discrete Element Method  
DEM was proposed by [PERSON] [37], which regards the research object as a collection of multiple rigid or deformable units. In this study, the Mohr-slip joint model is adopted, which assumes a linear relationship between stress and displacement, as shown in Figure 2. Once the stress of the sub-contact element exceeds the strength threshold, the sub-contact  
Figure 1: Mesoscopic image of silty mudstone under a polarizing microscope.  
element will undergo tensile or shear failure. The specific principle is as follows. In the normal direction of the sub-contact element, it is assumed that the normal stress-displacement obeys a linear relationship, which can be expressed as:  
\[\Delta\sigma_{n}=-k_{n}\cdot\Delta u_{n}, \tag{2}\]  
where \(\Delta\sigma_{n}\) is the normal stress increment, \(\Delta u_{n}\) is the normal displacement increment, and \(k_{n}\) is the normal stiffness. In addition, the threshold of normal tensile stress is the tensile strength \(f_{\text{I}}\). When the normal tensile stress exceeds the tensile strength, \(\sigma_{n}=0\).  
Similarly, the shear stress-displacement relationship in the shear process is controlled by the shear stiffness \(k_{\text{s}}\), which can be expressed as:  
\[\Delta\tau_{\text{s}}=-k_{\text{s}}\cdot\Delta u_{\text{s}}^{e}, \tag{3}\]  
where \(\Delta\tau_{\text{s}}\) is the shear stress increment, and \(\Delta u_{\text{s}}^{e}\) is the elastic part of the shear displacement increment. Then, the maximum shear stress \(\tau_{\text{max}}\) on the sub-contact element can be expressed as:  
\[|\tau_{\text{s}}|\leq c+\sigma_{n}\cdot\tan\varphi=\tau_{\text{max}}, \tag{4}\]  
where \(c\) and \(\varphi\) are the cohesion and internal friction angle of the sub-contact element, respectively. When the shear stress \(\tau_{\text{s}}\geq\tau_{\text{max}}\), the residual strength of the sub-contact element can be expressed as:  
\[\tau_{\text{s}}=\sigma_{n}\cdot\tan\varphi_{\text{r}}, \tag{5}\]  
where \(\varphi_{\text{r}}\) is the residual internal friction angle of the sub-contact element.  
According to the XRD test, the silty mudstone studied in this paper is mainly composed of quartz, feldspar, and cement (e.g., non-clastic minerals), with contents of 35%, 33%, and 32%, respectively. Therefore, the grain-based heterogeneous model should contain three components: quartz, feldspar, and cement, as well as six types of contacts: quartz-quartz contact, quartz-feldspar contact, quartz-cement contact, feldspar-feldspar contact, and cement-cement contact. The FISH language built into Iasca  
Figure 2: Schematic diagram of the Mohr-slip joint model [38].  
software is utilized to construct an algorithm to generate the grain-based discrete element (GB-DEM) model. The distribution of material components within the model follows a uniform distribution, and the generation of Voronoi polygons uses the built-in algorithm in the UDEC software. The contents of each material in the numerical model are consistent with the experimental results, as shown in Figure 3. The size of the numerical model is consistent with the test conditions, with the standard cylindrical model having a height of 100 mm and a diameter of 50 mm; the Brazilian disk has a diameter of 50 mm and a thickness of 25 mm, according to the International Society for Rock Mechanics Suggested Methods. The range of loading rates used for uniaxial compression and Brazilian disk tests in the existing literature [35, 36, 38, 39] is 0.01 m/s to 0.1 m/s. Since the default time step size in the UDEC software is about \(10^{-8}\) to \(10^{-7}\) s/ step, a loading rate of 0.1 m/s is equivalent to loading \(10^{-6}\) to \(10^{-5}\) mm per time step. It means that loading 1 mm requires \(10^{5}\) to \(10^{6}\) time steps, indicating that the specimen is under quasi-static loading conditions. Therefore, considering computational accuracy and efficiency, the loading rate is taken as 0.01 m/s.