A growing number of tasks require subterranean locomotion and sensing, including infrastructure inspection and soil characterization. However, the design and control challenges of burrowing into flowable, complex substrates have limited the development of bioinspired burrowing devices. To begin to develop design and control principles for multifunctional burrowing robots, we explore the burrowing biomechanics of the Pacific mole crab, Emerita analoga, which uses its appendages to swim, crawl, and rapidly burrow into intertidal substrate. Experiments on physical robot analogs of E. analoga motivated the development of a large N-body granular simulation model built on the Project Chrono physics engine to quickly model and test morphological and control parameters. We find that design parameters such as relative limb lengths have a strong effect on depth reached and the angle of intrusion. Analysis of the connection between internal shape velocities and vertical intrusion velocity reveals a periodic relationship between the phase of limb rotation and the stages in the burrowing process. The active phases expand the burrow, raising the robot, and then excavate the substrate behind the robot, moving it deeper during excavation and continuing into the recovery phase. Upward slips occur when limbs are expanding the substrate during active power strokes, and effects are amplified when both limb groups are in the same phase. The simulation environment is parameterized and reconfigurable to facilitate further studies of diverse body morphologies and controllers, and results can guide the design of novel legged burrowing robots.