## 2 Numerical Model  
As one of the explicit discrete element methods (DEMs), PFC2D (Particle Flow Code in two dimensions) of the company Itasca Consulting Group can observe the fracturing process indirectly, and has been widely used to simulate the failure behavior of rock materials [25]. In the PFC2D, the contact bond, parallel bond, and flat-joint model are the main bonding approaches to simulate the cementation between particles. Comparatively,the flat-joint model can mimic the grain interlocking of polygonal grain structure and replicate the mechanical parameters and crack propagation characteristics of real rocks effectively [26; 27; 28]. Therefore, the flat-joint model was applied to establish the numerical models in this study. In this study, the mechanical parameters of the surrounding rock were calibrated by the red sandstone from Linyi, Shandong Province, China, and the corresponding macroscopic mechanical parameters of experimental results and numerical results are shown in Table 1. The calibrated microscopic parameters of the numerical model are illustrated in Table 2, which have been validated in our previous work [29]. A numerical model with the geometry of \(20\times 15\) m was established by PFC2D, and two U-shaped caverns were set up in the middle of the model, as shown in Figure 1. In this section, the geometric dimensions of cavern 1 (left) and cavern 2 (right) are \(w_{1}=w_{2}=2.5\) m and \(h_{1}=h_{2}=3.75\) m, and the buried depth of the surrounding rock is 1000 m. According to the stress field of China, the horizontal and vertical principal stresses can be obtained by [30]:  
\[\left\{ \begin{array}{ll}\sigma_{v}=0.02532H+0.8177\\ \sigma_{h}=0.01766H+1.0583\end{array} \right. \tag{1}\]  
where \(H\) represents the buried depth of the cavern. According to Equation (1), the vertical stress \(\sigma_{v}=26.1\) MPa and the horizontal stress \(\sigma_{h}=18.7\) MPa are applied to the model through the particle boundary. To eliminate the reflected waves on the model boundary, the viscous boundary is set up [31]. When the disturbance wave source is relatively close to the cavern, the dynamic disturbance can be regarded as a cylindrical stress wave. Conversely, when the disturbance wave source is far from the cavern, the curvature of the incident wave can be neglected; thus, the dynamic disturbance can be simplified as a plane stress wave [32]. In this study, we take the blasting load as an example, and the incident stress wave is simplified to a triangular plane wave [33; 34]. When \(t<t_{r}\), it belongs to the loading stage of the incident wave, and when \(t_{r}<t<t_{s}\), it belongs to the unloading stage of the incident wave. The mathematical expression of the incident stress wave can be expressed as the following:  
\[P(t)=\left\{ \begin{array}{ll}0,&(t<0)\\ \frac{t_{r}}{t_{r}}P_{d},&(0\leq t<t_{r})\\ \frac{t_{r}-t_{r}}P_{d},&(t_{r}\leq t<t_{s})\\ 0,&(t\geq t_{s})\end{array} \right. \tag{2}\] \begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Mexological** & **Density** & **Elastic** & **Poisson’s** & **Uniaxial** & **Brazilian** \\
**Parameters** & **(kg/m\({}^{3}\))** & **Modulus** & **Poisson’s** & **Conpressive** & **Tensile** \\
**Ratio** & **Strength** & **Strength** & **(MPa)** & **Strength** \\ \hline Experimental & 2440 & 18.60 & 0.21 & 97.60 & 4.87 \\ results & 2440 & 18.74 & 0.21 & 97.38 & 4.91 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Macro-mechanical parameters of surrounding rock.