The discrete element method (DEM) can provide detailed descriptions of sea ice dynamics that explicitly model floes and discontinuities in the ice, which can be challenging to represent accurately with current models. However, floe-scale stresses that inform lead formation in sea ice are difficult to calculate in current DEM implementations. In this paper, we use the ParticLS software library to develop a DEM that models the sea ice as a collection of discrete rigid particles that are initially bonded together using a cohesive beam model that approximates the response of an Euler-Bernoulli beam located between particle centroids. Ice fracture and lead formation are determined based on the value of a non-local Cauchy stress state around each particle and a Mohr-Coulomb fracture model. Therefore, large ice floes are modeled as continuous objects made up of many bonded particles that can interact with each other, deform, and fracture. We generate particle configurations by discretizing the ice in MODIS satellite imagery into polygonal floes that fill the observed ice shape and extent. The model is tested on ice advecting through an idealized channel and through Nares Strait. The results indicate that the bonded DEM model is capable of qualitatively capturing the dynamic sea ice patterns through constrictions such as ice bridges, arch kinematic features, and lead formation. In addition, we apply spatial and temporal scaling analyses to illustrate the model’s ability to capture heterogeneity and intermittency in the simulated ice deformation.