## 1 Introduction  
Rockfill embankments are widely constructed around the world due to their low cost and strong adaptability to complex topographical and geological conditions [1; 2; 3; 4; 5; 6; 7; 8; 9]. While playing an important role, accidents that may occur in rockfill embankments also pose potential threats to the environment and society. The San Fernando dam suffered a 6.6 magnitude earthquake in 1971, resulting in slope instability and extensive collapse [10]. The Teton dam experienced piping during its initial impounding in 1976, ultimately resulting in 11 fatalities and USD 2 billion in economic losses [11]. Slope stability analysis is a major aspect of the safety evaluation of water-retaining embankments [12]. Since rockfills are incohesive materials, only friction angle is considered as a shear strength parameter in the slope stability analysis of rockfill embankments [13; 14]. Due to inherent variability, test errors, and other unknown factors, the values of the friction angle of rockfills are uncertain and the uncertainty of the friction angle can be quantified by its variance (a measure of the dispersion extent of friction angle around its mean value). In addition,the shear strength of rockfill materials is closely related to confining pressure, specifically, the friction angle of rockfills decreases with an increase in confining pressure [15; 16; 17; 18; 19; 20; 21]. Recently, it was found that confining pressure can significantly affect the mean value and variance of the friction angle of rockfill materials which decrease with increasing confining pressure [16]. Since the confining pressure (i.e., the minor principal stress) spatially varies within a rockfill embankment, this finding indicates that the mean value and variance of the friction angle of rockfills vary from one location to another in the dam as well. For high rockfill embankments, the stress-dependent spatial variability of the friction angle is expected to be significant and may considerably affect the slope stability. Probabilistic slope stability analysis provides a reasonable way to quantitatively characterize stochastic uncertainties and their effects on slope stability [22; 23; 24; 25; 26; 27; 28; 29]. However, for the complex formulations and algorithms as well as the requirement for a large quantity of samples to accurately determine the probability distributions and statistical parameters, the practical application of probabilistic slope stability analysis is not convenient. In practice, the slope ratio of rockfill embankments is commonly designed with safety factors of slope stability evaluated using the Limit Equilibrium Method (LEM) or numerical method, and the computed safety factors must be higher than the allowable safety factors [30; 31; 32; 33; 34]. The LEM considers the analysis domain as a rigid body, based on the principle of static equilibrium to calculate the safety factor. It has the characteristics of simplicity and high efficiency, including specific methods such as the Bishop's simplified method, the Janbu method, and the Morgenstern-Price method [35]. Numerical methods such as the finite element method, finite difference method, and boundary element method can more accurately take into account the topographical and geological conditions as well as the deformation behavior of rockfills, but at a higher computational cost [36]. As discussed above, the effect of stress-dependent spatial variability in the friction angle of rockfills should be investigated for the slope stability evaluation of rockfill embankments by considering the adverse variation towards low values of friction angle based on the distribution of confining pressure in rockfill embankments.