## 1 Introduction  
In geotechnical engineering, embankment slope failure has always been an important issue within one of the most common failure modes. As one of the water-retaining structures, embankments play an important role in protecting the lives and properties further downstream during flood periods. The material of an embankment slope is always simplified as homogeneous, and the FOS is used to describe the failure mechanism in the conventional method. The limit equilibrium theory is used to evaluate slope stability in geochemical engineering [1, 2]. The shear strength parameters, e.g., cohesion and internal friction angle, are generally considered as deterministic parameters in the conventional method. [PERSON] et al. [3] investigated the effects of rainfall intensity on the unsaturated slope. [PERSON] et al. [4] studied the surficial stability of unsaturated expansive soil slope considering the hydro-mechanical effect. However, geotechnical material has obvious uncertainty and randomness as well as spatial correlation structure after long-term consolidation [5]. [PERSON] and [PERSON] [6] pointed out that there are three primary sources of geotechnical uncertainties: inherent soil, measurement errors and transformation uncertainties, and the vertical correlation length is less than the horizontal correlation length (\(l_{x}\)). Many previous studies have shown that the strength parameters, e.g., cohesion and internal friction angle, of geomaterials follow Gaussian distributions or log-normal distributions [7, 8]. So far, probabilistic analysis has been widely adopted in geotechnical engineering to describe the uncertainty of geomaterials [9, 10, 11, 12, 13, 14]. The reliability theory was used to investigate the structure failure in which the uncertainty of materials can be considered. The point estimation method was used to calculate the reliability index with the advantage of time saving and simplicity [15; 16; 17]. The first-order second-moment method was also used to study the reliability of geotechnical engineering [18; 19]. [PERSON] [20] studied the reliability index of two-layered slope stability considering the spatial variability of soil. The limit state function was used to calculate the failure probability and reliability index of engineering structures; however, the limit state function of the complex structures sometimes cannot be obtained easily. One of the obvious shortcomings of reliability theory is that the spatial correlation structure of the soil cannot be considered. In recent decades, more advanced probabilistic analysis methods have been developed to study the failure of soil structures caused by the uncertainty of geomaterials. [PERSON] and [PERSON] [21] analyzed the two neighboring strip footings by considering the spatial randomness of the soil using subset simulation coupled with the random field model. [PERSON] and [PERSON] [22] first studied steady-state seepage by combining random field theory and the finite element method (RFEM), and since then, a lot of work using RFEM has been conducted to study geotechnical engineering, including slope stability [23; 24; 25; 26]. This excellent work has attracted increasing attention and has been used to study various uncertainties in geotechnics. Based on different study cases, seepage analyses were also conducted using RFEM to describe the spatial variation [27; 28; 29; 30]. [PERSON] et al. [31] studied the multi-stage slopes using the bivariate RFEM and concluded that although the FOS obtained from deterministic results was close to the mean value based on random field results, there was still a high risk of instability. [PERSON] et al. [32] investigated reservoir slope failure considering the spatial variability of the soil properties and concluded that the spatial variability has a strong effect on the stability. In this study, the influence of spatial variability of the soil strength parameters, cohesion and internal friction angle, on the slope stability of embankments is studied by RFEM coupling with seepage.