*3.1. Discrete Element Method (DEM)*

The simulation part of this research is based on DEM, which was firstly proposed by Cundall and Strack [26]. This method considers two types of motion of a particle, translation and rotation, which are governed by Newton's second law of motion. The elastic contact force expression used in this work is the non-linear Hertz-Mindlin no-slip model [27], which is illustrated in Figure 2. The basic expressions are given in Equations (2) and (3). The former is the translational equation, which is composed of gravitational force, *mig*, contact force and viscous contact damping force, where *Kn*, *Kt*, γ*<sup>n</sup>* and γ*<sup>t</sup>* express the normal elastic constant, tangential elastic constant, normal damping constant and tangential damping constant, respectively. A particle with the mass of *mi* contacts with *K* particles, and the contact force between them depends on the deformation between particles, δ*n.* In the equation, *ui*, *vn* and *vt* represent the translational velocity, and the component of relatively velocity for the normal and tangential directions.

$$m\_i \frac{du\_i}{dt} = \sum\_{j=1}^{K} \left( \mathbf{K}\_n \delta n\_{ij} - \gamma\_n v\_{nij} \right) + \left( \mathbf{K}\_t \delta n\_{ij} - \gamma\_t v\_{tij} \right) + m\_i \mathbf{g} \tag{2}$$

Equation (3) represents the rotational movement of particles, where *M<sup>k</sup> <sup>r</sup>* and *M<sup>d</sup> <sup>r</sup>* are two torques, which are caused by a tangential force and rolling friction. A Coulomb-type friction law is used to express the friction between the two particles. *Ii* and ω*<sup>i</sup>* denote the moment of inertia and rotational velocity, respectively.

$$I\_i \frac{d\omega\_i}{dt} = \sum\_{j=1}^{K} (\mathcal{M}\_r^k + \mathcal{M}\_r^d) \tag{3}$$

Table 1 presents the formulas used to calculate the forces and torques between the particles. In this work, we use the open-source software LIGGGHTS 3.5.0 [28] implementation of DEM.

**Figure 2.** Depiction of forces acting particle *i* in contact with particle *j*.


**Table 1.** Detailed description of the parameter expressions in discrete element method (DEM).

In the expressions, *G*, *Y*, *e* and υ represent the Shear modulus, Young's modulus, coefficient of restitution and Poisson's ratio, respectively, while μ*<sup>r</sup>* and μ*<sup>s</sup>* represent the rolling and static friction coefficient, respectively.
