As part of the Mechanics of Ice Rubble project, recent experiments have been carried out to study the strength and failure behavior of ice rubble beams and the freeze bonds that form between individual ice blocks. In this study, we present results obtained from a newly developed model for the three-dimensional (3D) discrete element modeling (DEM) open-source code LIGGGHTS. The ice model contains normal and shear springs that operate between neighboring particles which are bonded or that overlap due to compressional stresses. Energy dissipation is accounted for by using a viscous damping model. Using this DEM model, medium-scale freshwater ice rubble punch tests have been simulated. Rubble specimens were generated by “raining” individual DEM ice pieces into a rectangular form and compacting the rubble mass to achieve the target porosity. Before the compacting pressure was removed, bonds between contacting blocks were introduced with parameter values based on representative freeze bond experiments. The rubble beam was then deformed by pushing a platen vertically downward through the center of the beam until failure occurred. Two types of block size and shapes have been simulated: cuboid blocks generated based on the size distribution of the actual rubble, and rubble blocks generated by image processing of actual blocks of broken ice used in the comparison experiments. The mechanism of flexural rubble failure for both cuboid block [s4.2] simulations and empirical block [s4.3] simulations is in line with experimental results; however, the empirical block simulations provide a significantly better estimation of the failure force.