A numerical model for direct phase-resolved simulation of nonlinear ocean waves propagating through fragmented sea ice is proposed. In view are applications to wave propagation and attenuation across the marginal ice zone. This model solves the full equations for nonlinear potential flow coupled with a nonlinear thin-plate formulation for the ice cover. A key new contribution is to modeling fragmented sea ice, which is accomplished by allowing the coefficient of flexural rigidity to vary spatially so that distributions of ice floes can be directly specified in the physical domain. Two-dimensional simulations are performed to examine the attenuation of solitary waves by scattering through an irregular array of ice floes. Two different measures based on the wave profile are used to quantify its attenuation over time for various floe configurations. Slow (near linear) or fast (exponential-like) decay is observed depending on such parameters as incident wave height, ice concentration and ice fragmentation.