Soil spatial variability has essential influence on the reliability of geotechnical structures. Karhunen–Loève (KL) series expansion is an effective approach to characterize such features of soil properties. Acquiring solution for Fredholm integral equation of the second type is a necessary prerequisite; however, the corresponding analytical expressions are only available for limited circumstances. To overcome this challenge, a newly proposed method called Chebyshev–Galerkin–KL expansion was developed to discretize the random fields of soil parameters, from which the approximated eigenvalues and eigenfunctions can be obtained using the Chebyshev orthogonal polynomials of the second kind combined with Galerkin technique. Application of the proposed approach is illustrated through reliability analysis of an unsaturated slope example under different rainfall patterns, where the uncertainty in selection of a “best” soil-water characteristic curve (SWCC) model and statistical uncertainties in SWCC model parameters are taken into account. Results show that the developed approach is feasible to generate random fields with sufficient accuracy. Under a constant rainfall duration, the Advanced pattern may lead to shallow landslide with the highest probability, followed by Intermediate and Delayed. It should be noted that Bayesian inference and determination of optimal SWCC model should be carried out prior to reliability analysis. Otherwise, the landslide risk level would be exaggerated.