This paper aims to propose an efficient approach for evaluating the system reliability of multi-layered soil slopes using representative slip surfaces and multiple stochastic response surfaces (SRSs). First, the representative slip surfaces are identified from a large number of potential slip surfaces. For each representative slip surface, a stochastic response surface using the Hermite polynomial chaos expansion is constructed to estimate its factor of safety (FS). Second, direct Monte-Carlo simulations are performed to compute the system failure probability of the slope, of which the minimum FS for each random sample is calculated using SRSs of representative slip surfaces. Finally, a three-layered clay slope is investigated to demonstrate the effectiveness of the proposed approach. The results indicate that the proposed approach can effectively identify the representative slip surfaces of multi-layered soil slopes and produce accurate system failure probability which is commonly at relatively low levels. In addition, the proposed approach does not need to calculate the correlations between different potential slip surfaces for identification of the representative slip surfaces. The system failure probability of a multi-layered soil slope could be significantly underestimated if only the critical slip surface or insufficient representative slip surfaces are used.