## 3 Methods  
### Dynamic Fault Rupture Modeling in PFC2D  
In order to simulate an earthquake event in PFC2D (Particle Flow Code 2D), we developed a numerical workflow for fault dynamic rupture, as shown in Figure 3. This workflow outlines the process for simulating dynamic fault rupture using PFC2D. (1) The workflow begins with the creation of a bonded particle assembly, where particles are densely packed within the study domain. The study area domain is 120 km \(\times\) 180 km in size, and densely packed with particles of radii ranging from 200 to 300 m. In total, 93,318 particles are used. The particles are then bonded at their contacts, which mimics the cohesive properties of rockmass and faults. (2) After the creation of a bonded particle assembly, a pore-pipe network model is generated where a hydro-mechanical coupled scheme (Section 3.2) is implemented. (3) Traces of the major faults (Figure 1) are then inserted into the model using the smooth joint contact model [18], the properties of which are listed in Table 2. (4) We use the \"lock and release\" approach [19] for fault dynamic rupture. In this approach, a finite length of the fault at the hypocenter location is locked. (5) After this, the entire model domain is initialized uniformly with the in situ stress model, which is representative of the southeast region of the Korean Peninsula [20], produced from the field experiments by hydraulic fracturing and overcoming, as follows:  
\[\text{SH}(\text{MPa})=0.0305\;z(\text{m})+2.349 \tag{1}\]  
\[\text{Sh}(\text{MPa})=0.0178\;z(\text{m})+1.924 \tag{2}\]  
\[\text{Sv}(\text{MPa})=0.0266\;z(\text{m}) \tag{3}\]  
\begin{table}
\begin{tabular}{c c c} \hline
**Well No.** & **Name** & **Co-Seismic GWL (cm)** \\ \hline
4 & [PERSON], Rock & 38 \\
5 & [PERSON], Rock & 6 \\
6 & [PERSON], Rock & 1 \\
7 & [PERSON], Rock & 1 \\
8 & [PERSON], Rock & 12 \\
29 & [PERSON], Rock & 1 \\
30 & [PERSON], Rock & 4 \\
31 & [PERSON], Rock & 74 \\ \hline \end{tabular}
\end{table}
Table 1: Co-seismic GWL changes. The values are calculated from the record on 11 September and 12 September 2016.  
Figure 3: Workflow for fault dynamic rupture simulation using PFC2D v7.00.161.  
The maximum horizontal stress (SH) orientation known to be prevalent in the southeast region of the Korean Peninsula is N70\({}^{\circ}\) E [21]. The stress initialization is conducted by the PFC2D command \"BALL TRACTION\", so that the maximum stress is in the orientation of N70E throughout the entire model domain, as shown by the stress tensors in Figure 4. A pore pressure field is initialized uniformly throughout the entire model domain. (6) A quiet boundary condition is set to the boundaries to prevent the seismic wave from reflecting back to the model. This is to mimic a semi-infinite space. (7) A fault dynamic rupture is then simulated by instantaneously unlocking the previously locked part of the fault. A locking radius of 1400 m, corresponding to a fault length of 2800 m, is applied. The rupturing fault trace size is calibrated based on the fault scaling relations introduced in [22], to result in an event moment magnitude consistent with that of the 2016 GJEQ mainshock. (8) Upon fault unlocking, the strain energy previously stored along the locked fault trace is released and results in seismic wave initiation and propagation. Due to the damping of the particles, the seismic wave attenuates as it propagates. During this process, pore volume changes due to mechanical stress effects and pore pressure changes occur due to the poro-elasticity effect. Pore pressure evolutions at several selected locations are recorded to analyze their evolutions over time. The model is then run for 100 s to simulate the dynamic rupture process, capturing the fault's mechanical and hydraulic responses.