In this thesis, an innovative probabilistic approach was developed and implemented to explicitly incorporate the spatial variability of the bedrock and the geotechnical parameters of the soil through random fields. The methodology is designed for the numerical modeling of landslides in natural slopes triggered by rainfall events and encompasses the entire process, from generating the finite element mesh - based on elevation data and bedrockdepth - to creating 2D and 3D random fields to represent the spatial variability of geotechnical parameters. The transient unsaturated flow was simulated by solving Richards equation numerically, while slope stability was assessed using Numerical Limit Analysis. The entire process was fully automated to enable the efficient execution of multiple simulations within a probabilistic framework based on the Monte Carlo method. Additionally, advanced post-processing tools were developed, including an algorithm for quantifying the volumes of displaced material and delineating the most landslide-prone zones, allowing for a detailed analysis of the temporal evolution of instability areas. The results highlight the importance of the combined effect of uncertainties associatedwith bedrock depth and morphology, as well as the spatial distribution of flow and shear strength parameters, whose influence was quantified throughout the simulations. By integrating safety factor estimates, pore pressure, failure probability, and associated risk, this methodology enables a comprehensive assessment of unsaturated slope stability and was successfully applied to a real case study.