The reliability analysis of a spatially variable soil slope should be considered as a system reliability problem. The multiple response surfaces method (MRSM) provides a rational vehicle to evaluate it efficiently and accurately. This study firstly determines the reasonable response surface function for safety factor with considering spatial variability based on single slip surface. Then the approximate linear relationship between the accuracy of reliability analysis and that of random field discretization is explored. Using the MRSM on a large amount of potential slip surfaces, several representative slip surfaces (RSSs) are identified with the assessment of slope system reliability. The proposed method is illustrated and verified using two spatially variable soil slopes. The results indicate that as the spatial variability becomes stronger, the slope reliability becomes more systematic, and the failure probability obtained using single slip surface significantly underestimates the system failure probability. It is possible to ensure the accuracy of reliability analysis by predefining a proper value of the accuracy of random field discretization, aiming at reducing the number of unimportant random variables. In addition, choosing a proper response surface function helps improve both efficiency and accuracy of the proposed method. In addition, MRSM can simultaneously calculate the failure probability of all potential slip surfaces and the system failure probability, and identify RSSs, providing a guidance to make slope remedy measure.