Significant commitment of cost, man power, and time are often required for characterizing particle morphology of granular materials based on X-ray micro-computed tomography (mu CT) scanning and image processing. As a result, only a small portion of particles can be scanned, leading to uncertainty in characterizing realistic particle shape. This paper develops a probabilistic method to characterize and simulate sand particle shape based on Sphere Harmonic (SH) analysis, where Nataf transformation is used to construct the joint probability distribution of SH shape descriptors and to generate their random samples for reconstructing virtual particles. The proposed approach was applied to characterizing and simulating particle shapes of Leighton Buzzard sand (LBS), and these simulated particles were used as clump templates for discrete element method to account for effect of realistic particle shape. It was shown that the proposed approach properly characterizes and simulates the particle shape of LBS in a probabilistic manner. (C) 2019 Elsevier B.V. All rights reserved.