Understanding ice rubble strength and associated failure mechanics is important for a variety of engineering applications in marine ice environments, including the design and operation of coastal, offshore, subsea and floating structures. 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. These new results serve as an important guide for the development of improved numerical models. The discrete element method (DEM) is a direct modeling approach which has the potential to both describe and enhance understanding of the behavior of brittle granular materials, especially with regard to the evolution of damage towards failure. In this study we present results obtained from a newly developed model for the 3D 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 for ice rubble beams with nominal dimensions of 0.50m × 0.94m × 3.05m. 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 chosen based on representative freeze bond experiments. The ice rubble beam was then deformed by pushing a platen vertically downward through the center of the beam until failure occurred. For the numerical simulations presented here, two types of block size and shapes have been considered: 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. Results obtained for these two scenarios are compared with corresponding experimental test data. These results highlight that the DEM model is useful for estimating the flexural strength of the rubble, simulating the failure mechanism and for examining the extent to which the ice rubble beam failure is controlled by the strength of the freeze bonds. These results also provide valuable new insights regarding the importance of shape and size distribution of ice blocks on simulated ice rubble strength and failure behavior. Recommendations for future work are provided.