## 1 Introduction  
Slope rock masses are typically characterized by multiple joint sets and diverse lithologies, resulting in pronounced spatial heterogeneity [1, 2, 3]. First, the spatial distribution of rock mass properties--such as compressive strength, elastic modulus, and friction angle--varies due to mineral composition gradients, weathering processes, and tectonic deformation [4, 5]. Second, the development of structural planes, including joints, fissures, and faults, introduces additional anisotropy and weak-plane effects [6, 7]. These factors collectively lead to significant deviations in stress distribution and failure mechanismscompared with homogeneous assumptions. Neglecting or oversimplifying such spatial heterogeneity in slope stability analyses often results in underestimated hazard levels and inaccurately predicted failure extents, thereby compromising engineering safety and decision-making reliability [8,9].  
To address this, researchers have extensively explored methods for characterizing the spatial variability in rock mass mechanical parameters and integrating them into slope stability assessments. Currently, two primary approaches are prevalent: geostatistics [10] and random field theory [11]. Originating from mineral resource estimation, geostatistics quantifies the spatial correlation of rock mass parameters through variogram analysis [11,12], enabling the construction of spatial distribution models via interpolation techniques, such as Ordinary Kriging [13]. Its notable applications include the simulation of Rock Mass Rating (RMR) spatial variability by [PERSON] et al. [14], and the estimation of geological strength index (GSI) distributions by [PERSON] et al. [10].  
Alternatively, random field theory treats rock mass parameters as spatially correlated random variables, typically assuming normal or lognormal distributions [15]. Random parameter fields are generated using methods such as Cholesky decomposition [16], Karhunen-Loeve (K-L) expansion [17], and local average subdivision [18], and incorporated into finite element or finite difference models via Monte Carlo simulations [19]. Representative studies include [PERSON] et al.'s [4] investigation of weak-layer variability and seismic randomness effects and [PERSON] et al.'s [20] development of cross-correlated non-Gaussian random fields for stratified rock masses.  
While both approaches have established a solid foundation for the stability analysis of jointed rock masses [21,22], challenges remain, particularly for large-scale, highly heterogeneous open-pit slopes. Geostatistical methods typically treat the interpolated mechanical parameters as deterministic inputs [10], whereas random field models, though probabilistically robust, often entail significant computational costs, restricting their application to small- or medium-scale problems [23]. For example, [PERSON] et al. [24] performed a full Monte Carlo simulation (MCS) with 120,000 three-dimensional finite-difference realizations, each containing an 800-cell mesh; the complete run consumed \(\approx\) 306 CPU days (i.e., several hundred days of processor time), underscoring that reliability analyses of large study areas can demand enormous computational resources.  
This study addresses these limitations by investigating the southern slope of the Bayan Obo open-pit mine using a high-resolution geostatistical approach. The proposed methodology involves four sequential phases: (1) determining the representative elementary volume (REV) through three-dimensional discrete fracture network (DFN) simulations; (2) constructing a spatial distribution model of rock quality designation (RQD) based on borehole exploration data via Ordinary Kriging interpolation; (3) integrating the generalized Hoek-Brown failure criterion to convert the rock mass quality indices into the equivalent mechanical parameters, thereby generating a spatially variable mechanical parameter field; and (4) evaluating the slope stability using the point safety factor method within a numerical FLAC3D framework.