This paper develops an effective approach to evaluate the system reliability of slope with stochastic response surface method. Firstly, the proposed approach selects the representative slip surfaces from a large number of potential slip surfaces. For each representative slip surface, a stochastic response surface (SRS) is constructed to estimate its factor of safety with the Hermite polynomial chaos expansion. Then, the direct Monte Carlo simulation (MCS) is employed to calculate the system probability of slope failure. The minimum factor of safety for each random sample during the MCS is calculated with SRSs of representative slip surfaces. For illustration, the proposed approach is applied to estimate the system reliability of two slopes with multiple soil layers. The results show that the proposed approach can effectively identify the representative slip surfaces of the slope and yield the system probability of slope failure with reasonable accuracy. It can also evaluate the system reliability of slope at small probability levels, and provide an effective tool for system reliability analysis of slope considering correlated nonnormal soil parameters. In addition, the tedious procedure of calculating correlation coefficients between potential slip surfaces to determine the representative slip surfaces is avoided in the proposed approach. There may be multiple failure modes for a slope with multiple soil layers. If only one critical slip surface (e.g. critical deterministic slip surface) or insufficient representative slip surfaces are considered in system reliability analysis, the probability of slope failure would be underestimated.