## 2 Characterization of Spatial Variability and Slope Stability Assessment Methods  
### Ordinary Kriging Interpolation  
In this study, the primary dataset comprises RQD measurements. The exploratory data analysis did not reveal a statistically significant large-scale drift in the RQD; therefore, a trend-based formulation, such as Universal Kriging, was not warranted. No independent, spatially co-extensive secondary variable was available, making co-Kriging inappropriate (the GSI employed later was derived from the RQD and is therefore not an independent covariate). In addition, the modeling domain was large and the numerical analysis used a high-resolution grid on the order of two million grid zones, so fully probabilistic alternatives (e.g., conditional simulation or random-field Monte Carlo coupling) would have required hundreds to thousands of realizations and were computationally impractical under our budget. Under these conditions, Ordinary Kriging provided an optimal linear unbiased estimator for the given sample set and yielded a reliable deterministic spatial model of RQD to support the subsequent mechanical parameterization and slope stability assessment [25]. The core methodology is based on the formulation and application of a semivariogram, mathematically expressed as follows:  
\[\gamma(h)=\frac{1}{2n(h)}\sum_{i=1}^{n(h)}\left(X(z_{i})-X(z_{i}+h)\right)^{2} \tag{1}\]  
where \(X(z_{i})\) and \(X(z_{i}+h)\) denote the observed values at positions \(z_{i}\) and \(z_{i}+h\), respectively; \(h\) represents the distance between these two points; and \(n(h)\) is the number of sampling point pairs separated by distance \(h\).  
In practical applications, an experimental semivariogram's data points are typically fitted to derive a mathematical model. Among the available models, the spherical model incorporating a nugget effect is the most widely utilized in practice. Its mathematical expression is as follows:  
\[\gamma(h)=\left\{\begin{matrix}N+C\left(\frac{3h}{2\beta}-\frac{h^{3}}{2\beta ^{3}}\right)&h\leq\beta\\ N+C&h>\beta\end{matrix}\right. \tag{2}\]  
in which \(N\) represents the nugget value, \(C\) the sill value, and \(\beta\) the range. When the distance \(h\leq\beta\), the spatial correlation decays nonlinearly with an increasing \(h\). For \(h>\beta\), the spatial autocorrelation becomes negligible, and \(\gamma(h)\) approaches \((N+C)\), corresponding to the overall sample variance.  
According to the theory of optimal estimation, the Kriging weights \(w\) must satisfy both the unbiasedness condition and the minimum variance criterion. Their matrix expression is as follows:  
\[\left[\begin{matrix}w_{1}\\ w_{2}\\ \vdots\\ w_{n}\\ \lambda\end{matrix}\right]=\left[\begin{matrix}\gamma_{11}&\gamma_{12}&\dots &\gamma_{1n}&1\\ \gamma_{21}&\gamma_{22}&\dots&\gamma_{2n}&1\\ \vdots&\vdots&\dots&\vdots&\vdots\\ \gamma_{n1}&\gamma_{n2}&\dots&\gamma_{nn}&1\\ 1&1&\dots&1&0\end{matrix}\right]^{-1}\left[\begin{matrix}\gamma_{01}\\ \gamma_{02}\\ \vdots\\ \gamma_{0n}\\ 1\end{matrix}\right] \tag{3}\]  
where \(n\) denotes the number of effective sampling points, \(\gamma_{0i}\) represents the semivariogram value between the point to be estimated and the \(i\)-th known point, and \(\lambda\) is the Lagrange multiplier. Solving this system of linear equations yields the optimal weight coefficients, enabling the predicted value of the point to be estimated to be expressed as follows:  
\[\hat{x}_{0}=\sum_{i=1}^{n}w_{i}x_{i} \tag{4}\]