Abstract: The existing geometries of the slopes in slope reliability analysis considering spatial variability of soil properties are relatively small. An efficient approach based on the multiple response-surface and subset simulation is proposed for solving slope reliability problems involving relatively large slope geometries. An example of reliability analysis of two-layered heterogeneous clay slope with the height of 24 m is presented to demonstrate the effectiveness of the proposed method. The effect of marginal probability distributions, namely Gaussian, lognormal, Extvalue I, Gamma and Beta on slope reliability is investigated. The results indicate that the proposed approach possesses the following advantages: (1) it can properly evaluate the slope reliability at low-probability levels (i.e., 10^-9 ~ 10^-4) in spatially variable soils; (2) it effectively solves slope reliability problems involving relatively large slope geometries; (3) it greatly improves the computational efficiency in parametric sensitivity analysis, and provides an effective way to investigate the effects of statistics (e.g., probability distribution, scale of fluctuation) on the slope reliability. Additionally, the marginal probability distributions of soil properties significantly affect the slope reliability. The commonly-used Gaussian and lognormal distributions may overestimate and underestimate the probability of slope failure, respectively.