Abstract: A novel method is put forward for the random field discretization based on the Chebyshev-Galerkin-KL (Karhunen-Loève) expansion, followed by the derivation of equations for the proposed method. By means of Python language, an efficient program is exploited for automatically calculating the slope sliding volume and identifying the slope failure mode. The proposed method is validated through an unsaturated slope example subjected to water rising. The results indicate that the proposed method for the random field generation provides a new way to solve the Fredholm integral equation of the second kind, which can accurately characterize the spatial variability of geotechnical parameters. The Python-based program for risk estimation is decoupled from the random finite element calculations, which ensures the slope risk estimation with sufficient efficiency and promotes the reduction of the required time to predict the landslide risk. In addition, the obtained results from the unsaturated slope example show that a lower water rising velocity and a greater maximum water level will lead to the decrease of slope stability. The vertical spatial variability of geotechnical parameters has marginal effects on the safety factor of slopes. However, the sliding volume may be significantly affected. Attention should be paid to the negative cross-correlation between shear strength parameters when conducting reliability analysis of slope stability under water rising. Otherwise, the slope stability will be underestimated.