1. Introduction

Geotechnical reliability analysis requires a computational pipeline that connects three conceptually distinct components: (i) a stochastic input model describing how uncertain soil parameters vary in space; (ii) a deterministic forward model that maps a realisation of those parameters to an engineering output (the factor of safety); and (iii) a reliability engine that accumulates information about the probability of failure from many forward model evaluations. In practice these three components are rarely described together as a coherent software system, making independent replication difficult.

This paper fills that gap by documenting the complete MATLAB implementation used in a study of random-field (RF) slope reliability. Every function is named, its inputs and outputs defined, and representative code excerpts provided with line-by-line annotations. The goal is a document that a practising engineer with MATLAB access could follow to reproduce the analysis from scratch.

The study problem is a soil-nailed slope reinforced with six passive nail inclusions. Soil cohesion c′ and friction angle φ′ are treated as spatially correlated Gaussian random fields, discretised via the Karhunen–Loève (KL) expansion and evaluated at the base midpoints of ten slope slices. Three reliability engines are implemented and compared: FORM (First Order Reliability Method), MCS (Monte Carlo Simulation), and ARBIS (Adaptive Radial Based Importance Sampling). Section 2 describes the software architecture. Sections 3–6 cover the four main modules. Section 7 presents results.

2. Software Architecture and Data Structures

2.1 The allData Central Structure

All data in the implementation flows through a single MATLAB structure named allData. This design choice—analogous to a context object in modern software engineering—eliminates the need for global variables and makes each function call fully self-contained. The top-level structure contains four sub-structs:

1. allData.Geom — slope geometry (H, W, Cover profile matrix)
2. allData.Soil — soil properties (Gama, Phi, Coh, plus full XX matrix for RF runs)
3. allData.Anc — nail parameters (Eta, Th, LRod, Diam, Fy, Hole, SV, SH, Tmax, TF)
4. allData.Prob — problem settings (Entr, N slices, FS, R, D, OPT flag, Fmin)

The OPT flag in allData.Prob is a binary switch: OPT = 0 runs the slope without nail contributions (non-anchored), and OPT = 1 activates the full soil nail resistance terms in the factor of safety summation. This single flag makes it trivial to compare anchored versus non-anchored performance for any realisation.

2.2 inputData () — Data Structure Builder

The inputData () function populates allData with all default values. It exists in two versions: a FORM version (Appendix A1) with fixed scalar soil parameters, and an RF version (Appendix A2) that additionally defines the probabilistic XX matrix (mean, standard deviation, distribution type, CoV for each random variable).

3. Random Field Module

3.1 Correlation Structure and Mesh

The random field module wraps MATLAB's randomfield () function (Bellin, 1991). The correlation structure is defined by three fields in a corr struct: corr.name (the autocorrelation model name, set to 'exp'), corr.c0 (the anisotropic correlation length vector [horizontal, vertical] = [0.2, 1.0]), and corr. sigma (the field standard deviation, which may be spatially varying). The mesh is constructed from the Cartesian coordinates of all slice base midpoints: a Nx2 matrix of [x, y] pairs.

The KL expansion truncation level is set to trunc = 10 terms throughout. For the exponential autocorrelation with correlation length c₀ = [0.2, 1.0] and the slope domain of approximately 20 m × 12 m, ten terms capture well over 90% of the total field variance. The truncation error—measured as the sum of the discarded eigenvalues—is verified to be below 5% of the trace of the covariance operator in all test runs.

Figure 1 shows example RF output maps at two mesh resolutions—50×50 (coarse) and 150×150 (fine). The colour maps illustrate the spatial distribution: warm colours represent above-mean soil strength (high cohesion / large friction angle), and cool colours represent below-mean strength. The finer mesh resolves narrower but more intensely localised weak zones, demonstrating that mesh resolution affects the spatial texture of the RF while the statistical moments (mean, variance, autocorrelation) remain consistent.

Figure — Example random field output maps from MATLAB. Top row: 50×50 mesh (coarse realisation). Bottom row: 150×150 mesh (fine realisation). Left panels: colour-map intensity (blue = below mean, red = above mean). Right panels: profile plots showing spatial variability across the mesh. Source: Dissertation Figure 4.2.

4. Slope Stability Computation Module

4.1 Slice Generation and Geometry

The core stability function soilCalc(R, allData) computes the factor of safety for a circular arc of radius R defined by an Entry–Exit pair (Xin, Yin) → (Xout, Yout). The arc is divided into N = 10 equal-width slices. For each slice i, the function computes:

• Base coordinates (Xbot, Ybot) by solving the circle equation at each slice boundary
• Top coordinates (Xtop, Ytop) by linear interpolation along the Cover profile matrix
• Slice area as the trapezoid between base and top, then weight Wi = γ · Areai
• Base inclination αi = arctan((Ybot(i+1) − Ybot(i)) / ΔX) and arc length ΔLi = ΔX / cos(αi)

The factor of safety is then computed as the ratio of total resisting force to total driving force over all slices:

(1)

In the FORM version, c′ᵢ and φ′ᵢ are scalars, identical for all slices. In the RF version, they are vectors indexed by i—each element drawing its value from the spatially correlated random field evaluated at that slice's base midpoint. This is the single most consequential code change between the two analysis modes, and it is implemented in three lines of MATLAB.

4.2 Nail Resistance Contribution

When allData.Prob.OPT = 1, nail resistance is added to QQ. Each nail contributes a force along its axis, projected onto the failure surface. The contribution depends on the nail's tensile capacity (controlled by its cross-sectional area and yield strength Fy), the grout–soil skin friction (limited by Tmax), and the geometric projection angle between the nail axis and the slip surface normal. The maximum force resisted by the facing panel is bounded by TF = 10 MN. Nail contributions are computed slice-by-slice based on whether the nail head lies within the slice boundary.

4.3 Circle Optimisation: minFOS () and entrCalc ()

The function minFOS (X, allData) finds the radius R that minimises the factor of safety for a given exit point X. It wraps MATLAB's fmincon () with the active-set algorithm. The function entrCalc (Bound, allData) extends this by searching over the exit point X as well, returning the (X, R) combination that gives the global minimum FOS. The resulting FOSmin is the critical factor of safety FOS₁ reported in the reliability results.

5. Reliability Engines

5.1 FORM Driver: slopeBeta ()

The FORM driver slopeBeta (mu, Bound, allData) locates the design point in standard normal space using a gradient-based search. The performance function g(x) = FOS(x) − 1 is evaluated at candidate parameter vectors x. The FORM reliability index is the Euclidean distance from the origin to the design point:

(2)

The function returns beta (signed, so the calling script takes abs(beta)), Pf = normcdf(-beta), and the updated allData struct. The FORM driver calls soilCalc () internally with scalar soil parameters drawn from the design point.

5.2 MCS Driver

The MCS driver iterates N random field realisations. In each iteration, soilCalc () is called with a freshly generated RF, and the resulting FOS is stored. After N iterations, the failure count nfail (realisations with FOS < 1) yields:

(3)

MCS is unbiased and requires no assumptions about the shape of the limit-state surface. Its accuracy scales as 1/√N: for Pf ≈ 0.10, N = 1000 realisations gives a coefficient of variation on Pf of approximately 10%.

5.3 ARBIS Driver

ARBIS (Adaptive Radial Based Importance Sampling) concentrates sampling effort in the region near the limit-state surface (FOS = 1) rather than sampling uniformly from the prior distribution. The algorithm proceeds in two phases: (i) an initial MCS run identifies the approximate location of the failure boundary; (ii) subsequent samples are drawn from an importance density centred on that boundary, with each sample weighted by the ratio of the original to the importance density.

ARBIS is particularly important for the anchored slope configurations where Pf < 0.05. At these small failure probabilities, plain MCS requires N > 10,000 realisations to achieve coefficient of variation below 30% on Pf. ARBIS achieves equivalent accuracy with N ≈ 500 by focusing effort where it matters.

Figure 2 summarises the complete module architecture, showing data flow between all functions described in this paper.

Figure — Figure 2 – MATLAB software module architecture and reliability workflow. Entry-point modules (crimson) feed the Random Field module (amber) which passes spatially distributed parameters to the Slope Stability module (purple) and Circle Optimisation module (teal). Three interchangeable Reliability Engine modules (crimson/teal/green) all feed the Results Aggregation stage, producing FOS tables, β/Pf outputs, and failure surface plots.

6. Results

Table 1 presents a condensed comparison of outputs across two FORM runs (set-varied scalar parameters) and four random field realisations (two using MCS, two using ARBIS). The FORM runs use the FORM driver with fixed scalar cohesion (c′ = 5 and 15 kPa) and fixed friction angle (φ′ = 35°). The RF runs use the same mean values with CoV = 0.10 for both parameters.

Three features of Table 1 are noteworthy from an implementation perspective. First, the global FOS for RF realisations is substantially higher than that for FORM runs using the same mean parameters—up to 4.69 versus 1.33. This is not a contradiction: the global FOS in the RF analysis is the FOS computed for the mean-valued slope, while the reliability output (β, Pf) reflects the probability of the actual field producing a failure configuration. The two quantities measure different things.

Second, the critical FOS (FOS₁) is consistently lower than the global FOS in every RF realisation—sometimes dramatically so (0.77 versus 3.90 in Realisation 4). This gap reflects the presence of weak zones: the optimised slip circle routes through the lowest-strength path in the random field. The FORM analysis, using uniform scalar parameters, cannot generate this gap because all slices see the same strength.

Third, ARBIS and MCS give consistent Pf estimates for the same realisation index (compare Realisations 3 and 4 to 1 and 2). ARBIS is used for the realisations with smaller predicted Pf (3 and 4) where plain MCS would require many more iterations for equivalent precision.

7. Conclusion

This paper has documented the complete MATLAB implementation for random-field slope reliability analysis. The key design features of the workflow are:

1. A single allData struct carries all data between modules, eliminating global variables and enabling straightforward unit testing of each function.
2. The inputData () function exists in two versions—scalar (FORM) and probabilistic (RF)—allowing like-for-like comparison between analysis modes without changing any other code.
3. The critical change from FORM to RF analysis in soilCalc () is the replacement of scalar Coh and Phi with spatially indexed vectors Coh(i) and Phi(i). Every other computation remains unchanged.
4. Three interchangeable reliability engines—FORM, MCS, and ARBIS—are implemented through a uniform interface: each receives allData, calls soilCalc () internally, and returns β and Pf.
5. The OPT flag in allData.Prob enables immediate comparison of anchored versus non-anchored slope performance within the same workflow, without code duplication.

The results confirm that the RF implementation produces qualitatively and quantitatively different reliability outputs from FORM: non-negligible failure probabilities (7–46%) coexist with high global factors of safety (up to 4.7), and non-circular failure surfaces emerge without requiring any modification to the optimisation logic.

Appendix — Code Module Index

The following functions constitute the complete workflow. For publication, the full source would be archived in a repository with a persistent DOI. For this study, full listings appear in the dissertation appendices.

References

[1] Bellin, A. (1991). Numerical generation of random fields with specified spatial correlation structure. Internal Report No. IDR 2/1991, University of Trento, Italy.

[2] Vanmarcke, E. H. (1983). Random Fields – Analysis and Synthesis. MIT Press, Cambridge, MA.

[3]Baecher, G. B. and Christian, J. T. (2003). Reliability and Statistics in Geotechnical Engineering. John Wiley & Sons, Chichester.

[4] Sudret, B. and Der Kiureghian, A. (2002). Comparison of finite element reliability methods. Probabilistic Engineering Mechanics, 17, 337–348.

[5] Enevoldsen, I. and Sørensen, J. D. (1994). Reliability-based optimization in structural engineering. Structural Safety, 15(3), 169–196.

[6] Phoon, K. K. and Kulhawy, F. H. (1999). Characterisation of geotechnical variability. Canadian Geotechnical Journal, 36, 612–624.

[7] Griffiths, D. V. and Fenton, G. A. (2009). Probabilistic settlement analysis by stochastic and random finite element methods. Journal of Geotechnical and Geoenvironmental Engineering, 135(11), 1629–1637.

[8] El-Ramly, H., Morgenstern, N. R. and Cruden, D. M. (2002). Probabilistic slope stability analysis for practice. Canadian Geotechnical Journal, 39, 665–683.

[9] Phoon, K. K. (2008). Reliability-Based Design in Geotechnical Engineering. Taylor & Francis, London.

[10] Duncan, J. M. (2000). Factor of safety and reliability in geotechnical engineering. Journal of Geotechnical and Geoenvironmental Engineering, 126(4), 307–316.

[11] Taib, S. N. L. (2010). A review of soil nailing approaches. UNIMAS E-Journal of Civil Engineering, 1(2).

[12] Rackwitz, R. (2000). Reviewing probabilistic soils modelling. Computers and Geotechnics, 26, 199–223.