Uncertainties in parameters of landslide susceptibility models often hinder them from providing accurate spatial and temporal predictions of landslide occurrences. Substantial contribution to the uncertainties in landslide assessment originates from spatially variable geotechnical and hydrological parameters. These input parameters may often vary significantly through space, even within the same geological deposit, and there is a need to quantify the effects of the uncertainties in these parameters. This study addresses this issue with a new three-dimensional probabilistic landslide susceptibility model. The spatial variability of the model parameters is modeled with the random field approach and coupled with the Monte Carlo method to propagate uncertainties from the model parameters to landslide predictions (i.e., factor of safety). The resulting uncertainties in landslide predictions allow the effects of spatial variability in the input parameters to be quantified. The performance of the proposed model in capturing the effect of spatial variability and predicting landslide occurrence has been compared with a conventional physical-based landslide susceptibility model that does not account for three-dimensional effects on slope stability. The results indicate that the proposed model has better performance in landslide prediction with higher accuracy and precision than the conventional model. The novelty of this study is illustrating the effects of the soil heterogeneity on the susceptibility of shallow landslides, which was made possible by the development of a three-dimensional slope stability model that was coupled with random field model and the Monte Carlo method.