The Random Finite Element Method (RFEM) is a robust stochastic method for slope reliability analysis that incorporates the spatial variability of soil properties. However, the extensive computational time associated with the direct Monte Carlo simulation limits its practical application. To overcome this problem, this study investigates the use of machine learning (ML) models as surrogate models for the RFEM in both 2D and 3D contexts. It investigates the performance of three ML models in predicting slope stability by means of the factor of safety (FoS) based on a generated random field of the undrained shear strength. Additionally, a data augmentation technique is employed to improve performance. The models' performance is assessed for various slope cases, characterised by varying spatial variability. Two surrogate modeling approaches are employed: semi-surrogate modeling and full-surrogate modeling. In the semi-surrogate modeling approach, a small number of RFEM simulations are conducted for a specified case. The machine learning models are trained using the generated random fields as input data and the calculated factors of safety as output data. The mathematical models are then used to predict outcomes of FoS for a large number of random fields for the same specific slope case. In the full-surrogate modeling approach, many RFEM simulations are conducted for the training set, covering a range of spatial correlation lengths. Once trained, the full-surrogate models are ready for application to another different slope case without the need for any additional numerical simulation. The results indicate that the prediction accuracy of the ML models typically decreases for slope cases with smaller scales of fluctuation. Nonetheless, the FoS predictions by the best-performing semi-surrogate model are highly consistent with the results from RFEM simulations for the whole range of considered slope cases. In terms of predicting the probability of failure for 2D-modeled slopes, the accuracy is high, with relative errors within 10% across the cases considered. This level of accuracy is achieved using no more than 13% of the total number of realisations needed for RFEM analysis. Consequently, the computational time for reliability analysis involving 4000 realisations reduces from 67 hours using the RFEM to between 4 and 8 hours using a semi-surrogate model, with the time increasing as the spatial correlation length decreases. Predicting the p_f for 3D slopes using a semi-surrogate model showed larger errors, indicating a need for improvement. The full-surrogate models prove to be accurate for testing cases characterised by spatial correlation lengths within the training set's range. Notably, the best-performing full-surrogate model in 3D predicted the p_f within a relative error of 10% for two slope cases. This model performs a stochastic analysis of 4000 simulations within seconds, compared to 83 days of computational time required for RFEM reliability analysis.