## 2 Materials and Methods
### Flat Joint Model  
DEMs are controlled at the particle contact level. Particle contact constitutive models are implemented to represent a material macro-behavior under different scenarios. PFC 2D/3D [35] has a variety of built-in contact models, such as the Parallel Bond, Hertz, and Flat Joint (FJ) models. The FJ model was introduced by [PERSON] [15] and simulates the mechanical behavior of circular or spherical particles that might be bonded together at their contact points. In this paper, the FJ contact model is used because it is capable of reproducing the mechanical behavior of cemented materials [20].  
Figure 1a presents the FJ contact model for particle-particle contact behavior and the representation of the contact interface. The contact interface is the plane formed by 2 particles in contact with a length equal to 2-times the radius (\(R\)) of the smallest particle. The contact interface is split into \(N\) equal-length elements (\(e\)). Each element can sustain tensile (\(\sigma_{c}\)) and shear (\(\tau\)) stresses, calculated using Equation (3). The tensile (\(\overline{\sigma_{c}}\)) and shear (\(\overline{\tau}\)) strengths limits are set at the beginning of the simulation. The \(\sigma_{c}\) is calculated through the summation of normal forces (\(F_{F}^{n}\)) in the contact element and the area of the contact element (\(A_{e}\)), which is a function of R (Figure 1a). The shear stress (\(\tau\)) is calculated through the summation of tangential forces (\(F_{F}^{s}\)). The bonded state transitions to an unbonded state when either the \(\overline{\sigma_{c}}\) or \(\overline{\tau}\) strength limits are exceeded. The contact interface evolves from a fully bonded state to a fully unbounded (frictional) state (Figure 1b).  
\[\sigma_{c}=\frac{F_{e}^{n}}{A_{e}},\ \tau=\frac{\|F_{e}^{n}\|}{A_{e}},A_{e}= \frac{2R}{N} \tag{3}\]  
The bonded element state can sustain tensile and shear stress. The element breaks in tension if \(\sigma_{c}\) is greater than the tensile strength (\(\overline{\sigma}_{c}\)), Equation (4). The element breaks in shear if \(\tau\) is greater than to shear strength (\(\overline{\tau}\)) that is calculated with a Coulomb criterion,  
Figure 1: Flat Joint Model Description: (**a**) FJ particle-to-particle contact interface; (**b**) Rheological model components.  
Equation (5), and is a function of cohesive strength (\(\overline{c}\)), tensile strength (\(\overline{\sigma}_{c}\)) and friction angle (\(\overline{\phi}\)). These strength parameters are set at the beginning of the simulation.  
\[\sigma_{c}>\overline{\sigma}_{c} \tag{4}\]  
\[\tau>\overline{\tau}=\overline{c}+\overline{\sigma}_{c}\tan\overline{\phi} \tag{5}\]  
The unbonded element state is linearly elastic with a frictional slip. The tensile strength (\(\sigma\)) is zero, and the compressive strength (\(\widetilde{\sigma}\)) is a function of the element distance gap \(g_{s}\) and the normal stiffness of the unbonded element (\(k_{n}\)) (Equation (6)). Therefore, the shear strength (\(\overline{\tau}\)) follows the Coulomb criterion under compressive stress and is 0 otherwise (Equation (7)). In this equation \(\mu\) is the friction coefficient.  
\[\widetilde{\sigma}\Big{\{}\begin{matrix}0,&g_{s}\geq 0\\ -k_{n}g_{s},&g_{s}<0\end{matrix} \tag{6}\]  
\[\widetilde{\tau}\Big{\{}\begin{matrix}-\sigma\mu,&\sigma<0\\ 0,&\sigma=0\end{matrix} \tag{7}\]