## 2 Numerical Modeling
### Bond Model  
The bonding method between particles was based on the parallel bond method [22]. Figure 1 shows two disk-shaped particles bonded by a cuboid element in the two-dimensional model [13]. Each bond was characterized by the thickness (\(h_{b}\)), width (\(2R_{b}\)), length (\(l_{b}\)), radius factor (\(\lambda_{R}\)), bond Young's modulus (\(E_{b}\)), and the normal stiffness (\(k_{b,n}\)) and shear stiffness (\(k_{b,t}\)) per unit area. Here, the bond thickness (\(h_{b}\)) and the disk-shaped particle thickness (\(h_{d}\)) exist in the z-direction, and both values are identical. The width and length of the bond could be expressed as  
\[R_{b}=\lambda_{R}\text{min}\ \big{\{}r_{i},r_{j}\big{\}} \tag{5}\]  
\[l_{b}=r_{i}+r_{j} \tag{6}\]  
where the radius factor (\(\lambda_{R}\)) ranged from 0 to 1. \(k_{b,n}\) and \(k_{b,t}\) depend on the bond Young's modulus (\(E_{b}\)), bond length (\(l_{b}\)), and the ratio of the normal to shear stiffness (\(\lambda_{ns}\)), \(k_{b,n}\) and \(k_{b,t}\) could be expressed as  
\[k_{b,n}=\frac{E_{b}}{l_{b}} \tag{7}\]  
\[k_{b,t}=\frac{k_{b,n}}{\lambda_{ns}} \tag{8}\]  
The forces (\(F_{b}\)) and torques (\(M_{b}\)) acting on the particles connected with a bond caused a relative displacement and rotation of those particles. \[F_{b}=F_{b,n}+F_{b,t} \tag{9}\]  
\[M_{b}=M_{b,n}+M_{b,t} \tag{10}\]  
Before the start of the simulation, the forces (\(F_{b}\)) and moments (\(M_{b}\)) in the bond were initialized to zero. The relative displacement and rotation of the particles were governed by the linear elastic material and calculated as follows:  
\[\Delta F_{b,n}(t)=k_{b,n}S_{b}\Delta t\Delta u_{n} \tag{11}\]  
\[\Delta F_{b,t}(t)=-k_{b,t}S_{b}\Delta t\Delta u_{t} \tag{12}\]  
\[\Delta M_{b,tu}(t)=-k_{b,t}l_{b}\Delta t[\Delta{\mathbf{x}}_{n},\ 0] \tag{13}\]  
\[\Delta M_{b,bn}(t)=-k_{b,n}\Delta t\Big{[}l_{b,n}\Delta{\mathbf{v}}_{t},\ l_{b,z} \Delta{\mathbf{w}}_{z}\Big{]} \tag{14}\]  
\[l_{b}=r_{i}+r_{j} \tag{6}\]0 Two particles connected with a bond. where \(S_{b}=2R_{b}h_{b}\) is the cross-sectional area of the bond, and \(\Delta u_{n}\) and \(\Delta u_{t}\) represent the parallel and perpendicular components to the bond axis, respectively. The subscripts \(tw\) and \(bn\) indicate a twisting and a bending, respectively. \(\Delta w_{n}\) and \(\Delta w_{t}\) represent the normal and tangential components, respectively. \(I_{b,x}=I_{b,y}=\frac{1}{b}h_{b}^{3}R_{b}\), \(I_{b,z}=\frac{2}{3}h_{b}R_{b}^{3}\), and \(J_{b}=I_{b,y}+I_{b,z}\) represent the moments of inertia. \(J\) denotes the polar moment of inertia. The torque could be expressed as follows:  
\[l_{b}=r_{i}+r_{j} \tag{6}\]1  
In the 2-D model, \(\mathbf{F}_{b,t}\) contributed to the particle's rotation around the \(z\)-axis. The twisting moment \((\mathbf{M}_{b,tw})\) was zero, and only the z-component of the bending moment \((\mathbf{M}_{b,bn})\) was not zero. To simulate the ice fracture, we used the elastic beam theory. The normal bond strength (\(\sigma_{b}\)) and the shear bond strength (\(\tau_{b}\)) could be expressed as a maximum value as follows [22]:  
\[l_{b}=r_{i}+r_{j} \tag{6}\]2  
\[l_{b}=r_{i}+r_{j} \tag{6}\]3  
where \(n\) is a unit normal vector, \(\mathbf{H}\) is a \(3\times 3\) projection matrix from the 3-D space to the \(x\)-\(y\) plane, \(\mathbf{H}_{11}=\mathbf{H}_{22}=1\), and all remaining elements are zero. The negative sign in Equation (16) indicated that only tensile normal stress could cause bond failure. The normal stress acting on the bond element includes the bonding moment due to the relative rotation of the particle and the normal force. The bond began to break if at least one of the maximum stress components acting on that bond exceeded the bond strength as follows:  
\[l_{b}=r_{i}+r_{j} \tag{6}\]4