Abstract: A novel Gaussian process response surface method (GPRSM) is proposed, in which the Gaussian process regression algorithm is used to construct the relationship between the random variables and the response value of the performance function. Compared with the polynomial-based response surface method, the proposed method has high accuracy and efficiency for the reliability analysis with high-dimensional and highly nonlinear performance functions. This method can update the response surface dynamically by adding new training points. Meanwhile, to consider the spatial variability of soil properties, the random field is constructed by KL expansion and combined with the limit equilibrium method to evaluate the stability of slopes. The proposed GPRSM is used to build the surrogate model and used in Monte Carlo simulation for the failure probability of slopes, which reduces the calls to the slope stability analysis program while ensuring the calculation accuracy. Finally, the proposed method is applied to two case studies with explicit and implicit performance functions, respectively. Compared with those of other methods in the published papers, the validity and capability of the proposed method are proved.