Conditional Random Field Approach Combining FFT Filtering and Co-Kriging for Reliability Assessment of Slopes

1. Introduction

Rock and soil masses are natural materials with complex genesis. Owing to the multitude of structural planes and the highly uneven distributions of joints and fissures, critical mechanical parameters—such as cohesion and friction angle—exhibit marked spatial heterogeneity and multiscale variability [1,2,3,4]. Because limited field tests and sampling cannot fully capture this variability, substantial uncertainty surrounds assessments of overall slope stability [5,6]. For decades, engineers have largely relied on deterministic analyses, representing an entire study area with a single parameter set (or at most a few) to represent the entire area, thereby neglecting the spatial variability of material properties [7,8]. Such deterministic approaches fail to capture potential weak zones and the composite effects of alternating soft and hard strata that occur in real slopes, often under- or over-estimating safety risks and thereby compromising both the economy and reliability of design.

To overcome these deficiencies, numerical techniques such as the Random Finite Element Method (RFEM) [9,10] and the Random Finite Difference Method (RFDM) [11] have been widely adopted for complex slope-stability analyses. Anchored in random field theory, these methods characterize the spatial correlation, variability, and overall statistical distribution of geotechnical parameters within a probabilistic framework, thereby surmounting the limitations of deterministic approaches. Large-scale Monte Carlo simulations then produce reliability indicators—including failure probability, the distribution of safety factors, and the variability of potential slip surfaces—that furnish a more rigorous basis for risk evaluation and decision-making [12,13]. Most existing studies, however, still employ unconditional random fields (URFs) [10,11,13]; that is, they generate domain-wide random realizations using only prior statistics, without directly incorporating site investigation and testing data. This practice can result in “over-dispersion” or excessive conservatism, obscuring local strength characteristics and reducing the specificity and credibility of the reliability results.

Recent advances in geotechnical theory and information technology have underscored the pivotal role of field observations in producing scientifically sound reliability assessments [9,14]. Modern investigation techniques—such as borehole sampling, in-situ shear tests, and acoustic surveys—now yield reliable localized parameter data [15,16]. Accordingly, embedding these limited yet invaluable observations into random field models through conditional simulation has become a frontier research topic [17,18]. Appropriate conditioning calibrates the global random field, thereby improving both the accuracy and engineering relevance of slope-reliability analyses.

In response, this study proposes an integrated procedure for constructing conditional random fields (CRFs) by combining Fast Fourier Transform (FFT) filtering with co-Kriging interpolation. Using a representative three-bench slope case study, we conduct Monte Carlo simulations to compare slope failure patterns and safety-factor distributions for three models: the unconditional random field (URF), the independent random field (IRF), and the proposed CRF. These comparisons are carried out under varying bedding-dip angles and geotechnical parameter cross-correlation conditions. The specific objectives are to: (1) generate URFs via FFT filtering and examine how bedding-dip angle and parameter cross-correlation influence failure modes; (2) evaluate changes in stability results after removing cross-correlations (IRF); (3) construct CRFs with co-Kriging and quantify how field observations modify reliability outcomes; and (4) analyze the statistical distribution of safety factors, highlight differences among the random field models, and elucidate the technical advantages and practical value of CRFs in engineering reliability evaluations.

2. Method for Constructing CRFs

2.1. Generation of URFs via FFT Filtering

In geotechnical and rock-mechanics reliability analysis, both local spatial correlation and large-scale structural features of geomechanical parameters must be considered simultaneously [5,18,19]. An URF that satisfies a prescribed covariance function yields spatially correlated values at every point within the domain [1,20]. The FFT filtering technique provides an efficient means of generating such fields: by transforming the covariance function from the spatial domain into the frequency domain and filtering a white-noise field accordingly, one can rapidly construct a random field with the required correlation structure [21].

Before undertaking a reliability analysis, data from geological investigations, engineering experience, and laboratory tests must be integrated to derive the statistical characteristics of the rock–soil mass—namely, the mean, coefficient of variation, distribution type, cross-correlation coefficient, and correlation length [22,23]. These prior statistics underpin an accurate representation of the randomness and uncertainty inherent in geotechnical parameters. Because these parameters are non-negative and typically right-skewed, cohesion is commonly modeled with a log-normal distribution, whereas the friction angle is treated as normally distributed [17,24]. The relationship between their coefficients of variation $\delta$ and the logarithmic-space parameters is expressed as follows:

$$\delta ={\displaystyle \frac{\sigma }{\mu }}$$

(1)

$$\left\{\begin{array}{c}{\xi }^2=\ln (1+{\delta }^2)\\ \lambda =\ln (\mu )-{\displaystyle \frac{1}{2}}{\xi }^2\end{array}\right.$$

(2)

where $\sigma$ and $\mu$ are, respectively, the standard deviation and the mean of the parameter in real space, while ${\xi }^2$ and $\lambda$ represent the variance and the mean of the parameter in logarithmic space.

In view of the pronounced short-range correlation, long-range attenuation, and marked anisotropy of geotechnical parameters [25,26,27], an exponentially separable covariance function $C$ with explicit spatial structure is adopted. Its general form can be written as

$$C(x,y,z)={\sigma }^2\exp [-2({\displaystyle \frac{\left|\left.x^{\prime }\right|\right.}{dx}}+{\displaystyle \frac{\left|\left.y^{\prime }\right|\right.}{dy}}+{\displaystyle \frac{\left|\left.z^{\prime }\right|\right.}{dz}})]$$

(3)

where $dx$, $dy$, and $dz$ are the correlation lengths of the parameter along the three principal directions, and $x^{\prime }$, $y^{\prime }$, and $z^{\prime }$ are the spatial lags in a rotated coordinate system that aligns with the dominant geological fabric. These rotated distances are given by:

$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right)=R_x(\theta )R_y(\beta )R_z(\alpha )\left(\begin{array}{c}x\\ y\\ z\end{array}\right)$$

(4)

where $R_x(\theta )$, $R_y(\beta )$, and $R_z(\alpha )$ are the rotation matrices that correspond to counter-clockwise rotations through the angles $\theta$, $\beta$, and $\alpha$about the $x-$, $y-$, and $z-$axes, respectively.

To impose the desired covariance structure in the frequency domain and thus create a random field with the prescribed statistics, the frequency-domain power spectrum $G(u,v,w)$ is constructed as

$$G(u,v,w)=R(u,v,w)\circ \sqrt{T(u,v,w)}$$

(5)

where $u$, $v$, and $w$ are the discrete frequency indices that label grid points in the transformed domain; $\circ$ denotes the Hadamard (element-wise) product; $R(u,v,w)$ is the white-noise spectrum, whose grid-point values are independent and follow $N(0,1)$; and $T(u,v,w)$ is the power spectral density derived from the Wiener–Khinchin theorem. $R(u,v,w)$ and $T(u,v,w)$ are the three-dimensional FFTs of the white-noise field $R(x,y,z)$ and the covariance function $C$, respectively. Their explicit forms are given by:

$$R(u,v,w)={\displaystyle \sum _{m=0}^{N_x-1}{\displaystyle \sum _{n=0}^{N_y-1}{\displaystyle \sum _{p=0}^{N_z-1}R(x,y,z)e^{-2\pi i(\frac{um}{N_x}+\frac{vn}{N_y}+\frac{wp}{N_z})}}}}$$

(6)

$$T(u,v,w)={\displaystyle \sum _{m=0}^{N_x-1}{\displaystyle \sum _{n=0}^{N_y-1}{\displaystyle \sum _{p=0}^{N_z-1}C(x,y,z)e^{-2\pi i(\frac{um}{N_x}+\frac{vn}{N_y}+\frac{wp}{N_z})}}}}$$

(7)

where $N_x$, $N_y$, and $N_z$ denote the numbers of grid points along the three principal axes; $m$, $n$, and $p$ are the spatial-domain indices; and $i$ is the imaginary unit.

A three-dimensional inverse FFT is then applied to the filtered spectrum $G(u,v,w)$, transforming the signal back into physical space and producing a realization of the random field that conforms to the prescribed covariance structure $G(x,y,z)$:

$$G(x,y,z)={\displaystyle \frac{1}{N_x{N}_y{N}_z}}{\displaystyle \sum _{u=0}^{N_x-1}{\displaystyle \sum _{v=0}^{N_y-1}{\displaystyle \sum _{w=0}^{N_z-1}G(u,v,w)e^{-2\pi i(\frac{um}{N_x}+\frac{vn}{N_y}+\frac{wp}{N_z})}}}}$$

(8)

The filtered spatial field $G(x,y,z)$ is finally standardized to produce a standard normal field $S$, defined as

$$S(x,y,z)={\displaystyle \frac{G(x,y,z)-{\mu }_G}{{\sigma }_G}}$$

(9)

where ${\mu }_G$ and ${\sigma }_G$ are the mean and standard deviation of $G(x,y,z)$, respectively.

To incorporate the spatial cross-correlation between the cohesion field and the friction-angle field, the correlated standard-normal fields are constructed as follows:

$$S_{\phi }^{\prime }=\rho S_c+\sqrt{1-{\rho }^2}{S}_{\phi }$$

(10)

where $S_c$ and $S_{\phi }$ are the independent standard-normal fields produced in Eq. (9); $\rho$ is the cross-correlation coefficient between cohesion and friction angle; and $S_{\phi }^{\prime }$ denotes the standard-normal field for friction angle that is correlated with $S_c$.

Given the marginal distribution assumptions for cohesion and friction angle, the correlated standard-normal fields obtained from Eq. (10) are mapped to random fields of mechanical parameters that respect the specified priors, namely,

$$Z_{uc}=\exp ({\lambda }_c+{\xi }_c{S}_c)$$

(11)

$$Z_{u\phi }=\max ({\mu }_{\phi }+{\sigma }_{\phi }{S}_{\phi }^{\prime },{\phi }_{\min })$$

(12)

where $Z_{uc}$ represents the cohesion random field that follows a log-normal distribution, whereas $Z_{u\phi }$ denotes the friction-angle field modeled as a spatially correlated truncated-normal distribution. The parameter ${\phi }_{\min }$ is the empirically determined lower bound that guarantees $\phi$ remains non-negative across the domain. An illustrative realization of the fully spatially structured URFs generated by the FFT-filtering procedure is provided in Figure 1.

2.2. Construction of CRFs with the Co-Kriging Method

Co-Kriging is a multivariate geostatistical interpolation technique that exploits spatial cross-correlations between variables to obtain the statistically optimal estimate of multi-source data sets [28,29]. In this study, we refine the URFs with site-specific measurements through co-Kriging, thereby generating CRFs.

First, simulated values at the sampling locations are extracted from the URF and compared with their corresponding observations. The residuals of cohesion or friction angle, denoted $Z_{\Delta }$, are calculated and adopted as conditioning data:

$$Z_{\Delta }(x_i)=Z_o(x_i)-Z_u(x_i)$$

(13)

where $Z_o(x_i)$ is the measured cohesion (or friction angle) at sampling point $x_i$, and $Z_u(x_i)$ is the corresponding value extracted from the URF.

Using these residual data $Z_{\Delta }$, a residual field of the mechanical parameter is produced with the best linear unbiased estimator framework. Taking cohesion as an example, the estimated residual field ${\widehat{Z}}_{\Delta c}$ at an unsampled location $x_0$ is expressed as

$${\widehat{Z}}_{\Delta c}(x_0)={\displaystyle \sum _{i=1}^{n_c}{\epsilon }_i}{Z}_{\Delta c}(x_i)+{\displaystyle \sum _{j=1}^{n_{\phi }}{\delta }_j}{Z}_{\Delta \phi }(x_j)$$

(14)

where $Z_{\Delta c}(x_i)(i=1,2,\dots ,n_c)$and $Z_{\Delta \phi }(x_j)(j=1,2,\dots ,n_{\phi })$ are the cohesion and friction-angle residuals located at $x_i$ and $x_j$, respectively; ${\epsilon }_i$ and ${\delta }_j$ are the corresponding weight vectors. By imposing the unbiasedness condition and minimizing the estimation variance, these kriging weights are determined from the following system of linear equations:

$$\left[\begin{array}{cccc}{C}_{\Delta c}& C_{\Delta c\phi }& 1_{n_c}& 0_{n_c}\\ C_{\Delta \phi c}& C_{\Delta \phi }& 0_{n_{\phi }}& 1_{n_{\phi }}\\ 1_{n_c}^T& 0_{n_{\phi }}^T& 0& 0\\ 0_{n_c}^T& 1_{n_{\phi }}^T& 0& 0\end{array}\right]\left[\begin{array}{c}\epsilon \\ \delta \\ {\kappa }_1\\ {\kappa }_2\end{array}\right]=\left[\begin{array}{c}{C}_{x_0\Delta c}\\ C_{x_0\Delta \phi }\\ 1\\ 0\end{array}\right]$$

(15)

where $C_{\Delta c}$and $C_{\Delta \phi }$ are the auto-covariance matrices of the cohesion and friction-angle residuals (sizes $n_c\times n_c$ and $n_{\phi }\times n_{\phi }$); $C_{\Delta c\phi }$ and $C_{\Delta \phi c}$ are the cross-covariance matrices (sizes $n_c\times n_{\phi }$ and $n_{\phi }\times n_c$); $1_{n_c}$, $1_{n_{\phi }}$, $0_{n_c}$, $0_{n_{\phi }}$ are column vectors of ones or zeros of length $n_c$ or $n_{\phi }$, introduced to enforce the unbiasedness constraints; ${\kappa }_1$ and ${\kappa }_2$ are the Lagrange multipliers; $C_{x_0\Delta c}$ and $C_{x_0\Delta \phi }$ are the covariance vectors between the residuals at the sampled points and the target location $x_0$ for cohesion and friction angle, respectively (dimensions $n_c\times 1$ and $n_{\phi }\times 1$).

Assuming second-order stationarity, the auto--covariance function $C_{\Delta }(h)$ and cross-covariance function $C_{\Delta c\phi }(h)$ are written as

$$C_{\Delta }(h)={\displaystyle \frac{1}{N(h)}}{\displaystyle \sum _{i=1}^{N(h)}\left[Z_{\Delta }(x_i)-{\mu }_{\Delta }\right]}\left[Z_{\Delta }(x_i+h)-{\mu }_{\Delta }\right]$$

(16)

$$C_{\Delta c\phi }(h)={\displaystyle \frac{1}{N(h)}}{\displaystyle \sum _{i=1}^{N(h)}\left[Z_{\Delta c}(x_i)-{\mu }_{\Delta c}\right]}\left[Z_{\Delta \phi }(x_i+h)-{\mu }_{\Delta \phi }\right]$$

(17)

where $h$ is the lag distance—adjusted to the rotated coordinate system defined in Eq. (4); ${\mu }_{\Delta }$ is the mean residual of cohesion or friction angle; $N(h)$ is the number of sample pairs separated by lag $h$.

Finally, the CRF $Z_{CRF}$ is obtained by superposing the URF $Z_u$ with the estimated residual field ${\widehat{Z}}_{\Delta }$:

$$Z_{CRF}=Z_u+{\widehat{Z}}_{\Delta }$$

(18)

The CRF constructed herein simultaneously encapsulates the spatial variability of regionalized variables, their cross-correlations, and the constraints imposed by in-situ measurements, thereby furnishing a model that more faithfully represents in-service engineering conditions. Within a non-intrusive RFDM framework, the CRF parameters are mapped element-by-element onto the corresponding zones of a FLAC3D mesh. An iterative strength-reduction procedure is then executed to compute key response variables—most notably the safety factor of the slope and the maximum shear strain. These outputs form the basis for reliability assessment. The reliability index $\gamma$ can be expressed mathematically as

$$\gamma ={\displaystyle \frac{{\mu }_{\mathrm{FoS}}-F_{\mathrm{cr}}}{{\sigma }_{\mathrm{FoS}}}}$$

(19)

where ${\mu }_{\mathrm{FoS}}$ and ${\sigma }_{\mathrm{FoS}}$ are, respectively, the mean and standard deviation of the safety factor obtained from the RFDM simulations, and $F_{\mathrm{cr}}$ is the critical (target) safety factor. A comprehensive workflow of the entire procedure is presented in Figure 2.

3. Case Study

3.1. Model Configuration

To evaluate how the spatial variability of cohesion and friction angle—together with their cross-correlation—affects slope stability, we applied the proposed random-field framework to the slope model illustrated in Figure 3. The model measures 150 m (length) × 1 m (width) × 100 m (height). The engineered slope has a total height of 60 m and is divided into three benches, each 20 m high; the slope face dips at 60°, and every bench platform is 20 m wide.

Boundary conditions are specified as follows: normal displacements are constrained on the lateral faces, all degrees of freedom are fixed at the base, and the upper surface is free, being loaded only by self-weight in the vertical (-z) direction. The finite-difference mesh primarily comprises 1 m hexahedral zones, with a small number of tetrahedral zones of comparable size near geometric irregularities, giving a total of 10,751 zones. The prior statistical parameters adopted for the rock–soil mass are summarized in Table 1.

In the model, six investigation points are predefined and arranged in two boreholes—one on the uppermost bench and the other on the lowest bench. Each borehole contains three sampling points spaced 6 m apart. To quantify how these pointwise constraints affect slope stability, we first perform Monte Carlo simulations with URFs and IRFs generated via FFT filtering; we then incorporate the co-Kriging scheme to create CRFs. Finally, the failure patterns and safety factors obtained under the three random-field scenarios are compared.

3.2. Slope-Stability Analysis with URFs

Building on the model introduced in Section 3.1, we generate spatially complete random fields of the mechanical parameters in the frequency domain. A series of Monte Carlo experiments—with diverse combinations of geological orientation (bedding dip) and parameter cross-correlation—are carried out to systematically examine how spatial anisotropy and inter-parameter dependence influence failure mechanisms and safety factors.

Correlation lengths are set to 20 m in the x direction and 3 m in the z direction. To emulate the typical mismatch whereby local zones of low cohesion coincide with relatively high friction angles, the cross-correlation coefficient between the two parameters is fixed at −0.2. Five bedding-dip angles—15°, 30°, 45°, 60°, and 75°—are considered to capture anisotropy controlled by stratification; the resulting anisotropic parameter fields are illustrated in Figure 4.

For each bedding-dip angle, 200 independent $c-\phi$ URFs were generated, yielding a total of 1,000 realizations. After mapping these parameter fields to the numerical-model zones, the slope was driven to failure with the strength-reduction method (SRM), and the plastic-shear-strain distribution was used to delineate potential failure zones. As illustrated in Figure 5, the large result set was systematically classified according to the spatial relationship between the through-going plastic zones and the bench structure, enabling cross-comparison among dip angles. Nine representative failure modes (M1–M9) were ultimately distilled: Specifically, M1 is an overall failure penetrating all three benches; M2–M4 are local failures of the upper, middle, and lower single benches, respectively; M5–M7 are two-bench combined failures; M8 is a two-bench failure penetrating the upper and middle benches; and M9 corresponds to three benches failing independently.

To quantify how the bedding-dip angle governs both the distribution of failure modes and the slope-safety margin, the Monte Carlo results for each dip setting were post-processed to obtain the percentage share of every failure mode and the mean safety factor. The aggregate statistics for all 1,000 realizations serve as a benchmark. The compiled results are displayed in Figure 6.

Figure 6a shows that slope failures are dominated by either global, through-going collapse (M1) or single-bench, localized failures (M2–M4), which together account for more than 80 % of all realizations; complex composite modes are comparatively rare. The share of global failure (M1) follows a U-shaped trend with respect to bedding-dip angle, being high at shallow (15°) and steep (75°) dips but reaching its minimum (42 %) at 45°. In contrast, the lower-bench failure mode (M4) peaks at a dip of 45 % (15.5 %).

Figure 6b indicates that the safety factor (FoS) remains within a narrow band of 1.113 – 1.219 across all dip angles. The FoS for global failure (M1) is consistently lower than that for upper- (M2) and middle-bench (M3) failures, yet similar to that for lower-bench failure (M4). This pattern suggests that localized shear slips generally require a larger reduction factor to coalesce, whereas the lower bench—owing to stress concentration—tends to fail first and can trigger a transition from localized to global collapse (Figure 7). Moreover, the FoS associated with the global mode (M1) increases with bedding dip, implying that gently dipping strata oriented subparallel to the 60° bench face are more susceptible to slope instability.

3.3. Slope-Stability Analysis with IRFs

To isolate and quantify the influence of cohesion–friction-angle cross-correlation on slope behavior, we first examined the –0.2 case and then generated a second suite of random fields in which the cross-correlation coefficient was deliberately set to zero. In this IRF, the spatial fluctuations of cohesion and internal friction angle are entirely decoupled, eliminating the “weak-versus-strong” compensatory effect and enabling a clearer assessment of how correlation shapes failure patterns and safety factors. Figure 8 presents representative realizations of the mechanical-parameter fields produced under both correlation settings.

Figure 8 highlights how altering the cross-correlation coefficient markedly changes the heterogeneity of the mechanical-parameter fields. When the coefficient is –0.2, zones of low cohesion are typically offset by locally high friction angles, illustrating a “weak-strong” compensatory effect.

For the IRF analysis, 200 statistically independent realizations were generated for each of the five bedding-dip angles, yielding 1,000 simulations in total. All other model attributes—marginal statistics, correlation lengths, mesh, and boundary conditions—were identical to those described in Section 3.2. The results were post-processed to determine, for each dip angle, the proportion of the various failure modes and the mean safety factor (Figure 9a,b). Comparison matrices that summarize the differences relative to the URF case are presented in Figure 9c,d.

When the cross-correlation between cohesion and friction angle is removed, weak zones are characterized by simultaneously low values of both parameters, which reduces the local shear strength and renders the slope more susceptible to localized instability. As illustrated in Figure 9(c), the failure-mode distribution in the IRF undergoes a pronounced redistribution. In the absence of the compensatory effect, the proportions of the upper-bench (M2) and lower-bench (M4) localized-failure modes increase markedly—to 1.20 and 1.10 times those observed in the URF case, respectively. Because both M8 and M3 involve slip-out of the middle bench, the risk of failure in that bench also rises. By contrast, the frequency of the global through-going mode (M1) at steep bedding dips (60°–75°) decreases to 0.77 times that of the URF. Figure 9(d) further shows that the mean safety factor of the IRF declines overall, yet the reduction is modest (less than 2 %). These results indicate that, relative to the slight change in safety factor, the loss of parameter correlation exerts a far more significant influence on the distribution of failure modes.

3.4. Slope-stability analysis with CRFs

The preceding URF and IRF studies underscored how the spatial correlation structure of cohesion and friction angle shapes failure mechanisms. In practice, however, the sparse site-investigation data available for most projects are seldom exploited in traditional random-field frameworks. To make reliability assessments more targeted and credible, CRFs approach was adopted: the URFs generated in Section 3.2 were jointly conditioned on the six in-situ measurements by means of co-Kriging, thereby quantifying how field observations modify the reliability outcome. Figure 10 presents the resulting mechanical-parameter CRFs with a cohesion–friction-angle cross-correlation coefficient of −0.2.

For each of the five bedding-dip angles considered—15°, 30°, 45°, 60°, and 75°—200 realizations of the $c-\phi$ CRF with a cross-correlation coefficient of −0.2 were generated, yielding a total of 1,000 realizations. Under otherwise identical modeling conditions, the simulation outcomes were post-processed to obtain the proportion of each failure mode and the mean safety factor for every dip angle (Figure 11a,b). To highlight the influence of conditioning, the corresponding statistics were further organized into comparison matrices (Figure 11c,d).

Figure 11(c) shows that the high-strength constraint at the mid-depth sampling point of Borehole 2 suppresses the occurrence frequency of the M1 failure mode to 0.88 times its baseline value, whereas the low strength measured at the upper sampling point of Borehole 1 promotes a 1.21-fold increase in the proportion of M2. Notably, the “low-over-high” vertical strength profile observed in Borehole 2 imparts a marked dip-angle dependence to M4. This constraint narrows the extension angle available to potential slip surfaces and favors their penetration through the steeply dipping strata; consequently, the share of M4 rises sharply at steep dips (60° and 75°) but declines at gentle dips (15° and 30°). The corresponding increases in the mean factor of safety for M4 under gentle dips—1.038 and 1.023 times, respectively, as plotted in Figure 11(d)—corroborate this trend.

Moreover, the proportion of M1 still follows a “U-shaped” distribution. From a mechanical standpoint, both gently and steeply dipping strata easily develop local step-like rock bridges that delay instability. During the strength-reduction analysis, an incipient slip surface that cannot breach these local steps eventually propagates along zones of concentrated shear stress, leading to overall failure.

To further quantify the effects of dip angle and the spatial random-field model on slope stability, the distribution maps of the factor of safety are provided in Figure 12.

Figure 12 systematically displays the kernel-density curves of the FoS obtained from three spatial random-field models— URF, IRF, and CRF—evaluated at five bedding-dip angles, together with their key statistics ($\mu$ and $\sigma$). For direct comparison, Figure 12(d) superimposes the aggregate histograms of the three models and their corresponding Gaussian fits.

The results reveal that both URF and IRF distributions are mildly right-skewed, with the IRF exhibiting a noticeably heavier right tail, whereas the CRF curve is nearly symmetric and shows well-converged tails. Relative to the URF baseline, the CRF distribution shifts rightward by 0.010 and becomes narrower, while the IRF shifts leftward by 0.006 and widens; the URF lies between these two extremes. These shifts indicate that, once the negative cross-correlation between cohesion and friction angle is removed in the IRF, the “weak–weak” co-location mechanism of mechanical parameters not only lowers the mean FoS but also enlarges its dispersion ($\sigma$ rises by 10.26 %). Conversely, incorporating field observations in the CRF effectively contracts the predictive uncertainty ($\sigma$ declines by 15.38 %), an effect that is most pronounced at low dip angles where local weakness exerts stronger control.

Assuming a short service life for the slope and negligible consequences upon failure, a critical FoS of 1.1 is adopted. Based on the fitted distributions, the calculated failure probabilities for the URF, IRF, and CRF are 7.49 %, 12.30 %, and 2.30 %, respectively, yielding reliability indices of 1.452, 1.152, and 2.011. These figures substantiate the substantial improvement in slope reliability achieved by incorporating observation-based constraints.

4. Discussion

This study systematically compares slope stability outcomes using three different random field descriptions: URF, IRF, and CRF. The comparison highlights the strong influence of spatial variability, cross-correlation, and fields measurements on both failure mechanisms and factor-of-safety (FoS) statistics.

(1)

Failure mechanisms. Incorporating field observations substantially changes the evolution of the incipient slip surface. In the CRF model the probability of through-going failure (Mode M1) is markedly lower than in the URF, whereas modes characterised by local sliding or multi-step combinations (e.g. M2, M8) become much more frequent. Borehole data locally update the material strength, which in turn inhibits the formation of continuous low-strength bands and mitigates the risk of a large, penetrating slip surface. At the same time, the revealed weak zones trigger shear instability in the upper benches. This behaviour is consistent with the shallow, small-scale failures commonly observed in practice.

(2)

Statistical characteristics of the FoS. The CRF increases the mean FoS slightly (≈ 0.01) while reducing its variance by roughly 15 %, yielding a markedly more concentrated distribution with thinner tails. This convergence is most pronounced at intermediate bedding dips of 15°–45°, underscoring the decisive role of observation constraints in controlling local weakening zones and illustrating a robustness effect arising from the inclusion of field information together with parameter cross-correlation..

(3)

Consistency of input parameter fields (Figure 13). Goodness-of-fit tests show that the generated cohesion field (log-normal) and friction-angle field (truncated normal) reproduce the prescribed marginal distributions well (mean Kolmogorov–Smirnov statistic D < 0.035). The sample mean of the target cross-correlation coefficient (ρ = −0.2) deviates by only 0.007, indicating that the framework accurately captures the prior, although the sample standard deviation remains about 0.08. Such “accurate mean yet noticeable fluctuation” is expected when a finite domain, long correlation length and nonlinear distribution mapping coexist, and it exerts only a minor influence on Monte Carlo estimates of failure modes and FoS.

To further reduce dispersion in the cross-correlation coefficient, two technical options are available. (i) Refine the FFT spectral resolution (or use a closed-form transform) during random-field generation so that the target correlation is more closely achieved in standard-normal space. (ii) Apply a threshold filter (e.g., discard samples with correlation coefficient deviating by more than ±2$\sigma$) to reduce fluctuations and improve the stability of the reliability analysis.

Using a 13th-generation Intel® Core™ i7-13700F CPU, generating 1 000 CRFs required ≈ 90.4 s, whereas solving the FoS for a single field demanded ≈ 10 min—confirming that the computational bottleneck lies in the numerical-simulation step. Because site-investigation data are both time-consuming and expensive to obtain, the CRF approach makes full use of the available measurements, thereby avoiding the safety-margin misestimation often associated with URFs. Although conditioning slightly perturbs the global parameter distributions (mean K–S statistic D rises by < 2 %), the CRF attains markedly higher accuracy and relevance in reliability assessment at practically the same computational cost. The resulting FoS distribution is narrower with lighter tails, delivering more robust designs that support rapid decisions and on-the-fly adjustments — a clear practical advantage.

Now that the CRF framework’s effectiveness for slope reliability has been demonstrated, future work can focus on integrating this approach with optimized borehole-investigation planning. Failure-mode probabilities derived from CRFs can quantify the slope’s sensitivity to specific zones, thereby guiding strategic borehole placement. This reliability-sensitive drilling scheme pinpoints areas where geotechnical uncertainty is both high and influential, reducing investigation costs while improving information efficiency. Newly acquired borehole data can then update the priors and refine the CRF [30], closing the loop of “survey design → data acquisition → parameter updating → field-model correction”. Such a closed-loop workflow would progressively enhance both the accuracy and robustness of slope-stability assessments, shifting practice from passive safety evaluation to proactive risk management. An active investigation-optimisation strategy grounded in CRFs therefore represents a promising avenue for future research.

5. Conclusions

To quantify the influence of spatial variability in geomechanical parameters on slope reliability, this study proposes a CRF framework that couples FFT filtering with co-Kriging conditioning. A systematic comparison with URFs and IRFs over a range of bedding-dip angles yields three principal conclusions:

(1)

Cross-correlation governs failure mechanisms. A negative cohesion–friction-angle correlation (ρ = −0.2) produces a pronounced strength-compensation effect, whereas assuming independence (IRF) accentuates “weak-weak” spatial co-location. Consequently, the IRF increases the local-failure probability, reduces the mean FoS by 0.006, enlarges its standard deviation by 10.26 %, and raises the probability of low-FoS events (FoS < 1.1) from 7.49 % to 12.30 %.

(2)

Observation constraints optimise the failure-mode distribution. By suppressing the formation of extreme local weak zones, the CRF reduces the probability of through-going failure (Mode M1) by an average of 12 % and increases the incidence of local or multi-step failures (e.g. M2, M8), yielding patterns that more closely match field observations.

(3)

The CRF markedly enhances FoS robustness. Relative to the URF, the CRF produces a tighter FoS distribution with lighter tails: the mean FoS rises by 0.010, the standard deviation decreases by 15.38 %, and the probability of low-FoS events falls to 2.30 %. These improvements provide more reliable and actionable guidance for engineering design.