## 2 DEM Simulation of Direct Shear Test
### Development of PFC3D Model  
The materials, dimensions, and boundary conditions in the numerical model were based on direct shear tests conducted by [PERSON] et al. [32]. Figure 2a shows the DEM model of this direct shear box for FRS with random fibers. The size of the direct shear box was chosen as 25 mm in height, 60 mm in length, and 60 mm in width (Figure 2c). The virtual shear box was enclosed by 10 rigid wall elements and divided horizontally into 2 equal halves (upper and lower). The lid of the shear box placed on top of the particles was a servo-controlled wall, and normal force was exerted on the particle assembly within the box. The rigid load plate moved up and down during compression and shearing. In addition to the 10 walls, 2 more walls (1 on the left and 1 on the right) were used for closing the horizontal surfaces on either side of the shear plane. This is a necessary step to avoid particles escaping from the box during shearing.  
Sand particles were modeled as rigid balls with conventional contact elements (Figure 2b). The grain size was slightly modified and scaled up to simulate the problem effectively. In the sand matrix, a normal distribution was used for a particle radius range from 1.0 to 1.4 mm. A fiber was modeled as a series of linked particles with a bond contact algorithm, where the strength of the bond is fairly high. Each ball element in a fiber had the same diameter. The ball diameter and number of balls could be  
Figure 1: Computational cycle of the discrete element method.  
varied to match the diameter and length of the fiber being simulated. Any arbitrary fiber shape also could be constructed by changing the arrangement of the balls. In this study, we used a straight virtual fiber (Figure 2d). The fibers could be located at horizontal, vertical, and oblique directions with respect to the shear plane. Considering the test situation, fibers were randomly oriented in a direct shear box during the calibration process (Figure 2c).  
The material friction coefficients could be adjusted to large and small values to achieve loose and dense setups, respectively. A very low friction coefficient was adopted to generate a dense assemblage. The assembly porosity of the test sample at the initial condition was controlled at 0.4, similar to the samples tested in the laboratory [32]. The sand-fiber system was then brought to a force equilibrium state by iteration. The friction coefficients were then reset to a specific value. A vertical displacement was then applied slowly to the top wall until the desired vertical normal loading was reached. This process formed a test specimen containing a homogeneous and well-connected granular assembly with a predetermined normal load.  
When the expected normal pressure (100 or 200 kPa) was achieved and the system reached equilibrium again, the shearing stage was started. Shearing was induced by moving the upper part of the shear box at a constant velocity in a positive x-direction, while the lower section was fixed. Note that the loading speed (0.01 m/s) needed to be extremely slow to guarantee that the test was  
Figure 2: Specimens for three-dimensional DEM simulation: (**a**) DEM model of FRS, (**b**) sand matrix and direct shear box, (**c**) random fibers and direct shear box, and (**d**) typical modelfor fibers in DEM.  
conducted under quasi-static loading conditions. During the whole shearing stage, the velocity of the top wall was adjusted using Iitasca's built-in servo-controlled functions to maintain a constant normal pressure. A horizontal velocity was applied to the upper box until the sample failed by shearing or the horizontal displacement reached 10 mm. The physical quantities, such as shear stress and horizontal displacements, were measured and recorded at the corresponding boundaries.