## 2 Numerical Model and Methodology  
The motion of particles is governed by the well-established [PERSON]'s second law. For one particle, the translational and rotational equations are as follows:  
\[\mathbf{F}_{i}=m_{i}\frac{\mathrm{d}\mathbf{\omega}_{i}}{\mathrm{d}t} \tag{1}\]  
and  
\[\mathbf{T}_{i}=\mathbf{I}_{i}\frac{\mathrm{d}\mathbf{\omega}_{i}}{\mathrm{d}t}-(\mathbf{I}_{i} \cdot\mathbf{\omega}_{i})\times\mathbf{\omega}_{i} \tag{2}\]  
in which \(m_{i}\), \(\mathbf{I}_{i}\), \(\mathbf{\
u}_{i}\), and \(\mathbf{\omega}_{i}\) are the mass, moment of inertia, and translational and rotational velocities of the particle \(i\), respectively; \(\mathbf{F}_{i}\) and \(\mathbf{T}_{i}\) are external forces and moments exerting on it. The discrete element method (DEM), proposed by [PERSON] and [PERSON] (1979), is adopted to simulate the interaction between two objects. The DEM used in this study is a soft-sphere model based on spherical contact theories, and there are four contact scenarios:particle-particle, particle-hemispherical nose, particle-container inner surface, and particle-intruder outer surface. The contact models used in this study include Hertz for normal force and Mindlin for tangential force. The normal component of the contact force \(F_{n}^{c}\) of the first three contact scenarios can be expressed as  
\[F_{n}^{c}=\frac{4}{3}E^{*}R^{*\frac{1}{2}}\delta_{n}{}^{\frac{3}{2}} \tag{3}\]  
in which \(\delta_{n}\) is the overlap in normal direction. \(E^{*}\) is equivalent [PERSON]'s modulus, and \(R^{*}\) is equivalent radius of two objects in contact, which are defined as \(\frac{1}{E^{*}}=\frac{1-v_{1}^{2}}{E_{1}}+\frac{1-v_{2}^{2}}{E_{2}}\) and \(\frac{1}{R^{*}}=\frac{1}{R_{1}}+\frac{1}{R_{2}}\), respectively, where \(E_{1}\) and \(E_{2}\), \(R_{1}\) and \(R_{2}\), \(v_{1}\) and \(v_{2}\) are the Young's moduli, the radii, and the Poisson's ratios of two contacting particles. For the particle-intruder outer surface, the normal component of the contact force \(F_{n}^{c}\) can be expressed as  
\[F_{n}^{c}=8\sqrt{\frac{R^{*}}{27}}\alpha^{-\frac{3}{2}}E^{*}\delta_{n}^{\frac {3}{2}} \tag{4}\]  
in which \(\alpha\) is determined by the shape of the contact area and set to be 0.974 [37]. The Mindlin model is utilized for the tangential contact force [38]. When two contacting particles are subjected to an increasing tangential displacement, \(\delta_{t}\), the incremental tangential force \(\Delta T\) due to the incremental tangential displacement \(\Delta\delta_{t}\) depends not only on the loading history but also on the variation of the normal force. Therefore, \(\Delta T\) can be expressed as  
\[\Delta T=8G^{*}a\theta_{k}\Delta\delta_{t}+\left(-1\right)^{k}\mu\Delta N(1- \theta_{k}) \tag{5}\]  
where \(\Delta N\) is the incremental normal force, \(\mu\) is the coefficient of sliding friction, \(a=\sqrt{\delta_{n}R^{*}}\) is the radius of the contact area, and the parameter \(\theta_{k}\) depends on the loading status, e.g., when \(|\Delta T|<\mu\Delta N\), then \(\theta_{k}=1\), so \(\Delta T=8G^{*}a\Delta\delta\), which corresponds to a case in which there is no microslip, otherwise microslip occurs and  
\[\theta_{k}=\sqrt{1-\frac{T+\mu\Delta N}{\mu N}}\text{ if }k=0 \tag{6a}\] \[\theta_{k}=\sqrt[3]{1-\frac{\left(-1\right)^{k}\left(T-T_{k}\right)+2\mu\Delta N }{2\mu N}}\text{ if }k=1,2 \tag{6b}\]  
here, \(k\) = 0, 1, 2 represent loading, unloading, and reloading paths, respectively, and \(T_{k}\) is the historical tangential force at which loading or reloading commenced. The parameter \(G^{*}\) is defined as  
\[G^{*}=\frac{2-v_{1}}{G_{1}}+\frac{2-v_{2}}{G_{2}} \tag{7}\]  
in which \(G_{1}\) and \(G_{2}\) are the shear moduli of each particle. Also, the contact damping forces are added to dissipate the energy upon the contacts.