Abstract: A probabilistic numerical approach for the stability analysis of soil and rock slopes is presented based on the polynomial chaos expansion. The slope is divided into a family of inclined slices. By using the Mohr-Coulomb yield criterion and according to the upper bound theorem in plasticity theory, a work-energy balance equation is proposed for each soil/rock slice with respect to the kinematically-admissible global velocity field of the slope. Geotechnical parameters are modelled as random variables, and the safety factor is treated as a functional of these random parameters. By using the polynomial chaos expansion and a stochastic Galerkin optimization process, the random safety factor and the corresponding critical failure surface are solved. Several numerical examples demonstrate that the proposed method has the advantage of accurately estimating not only the safety factor and the critical failure surface but also the failure probability of the slope.