Data-Driven Prediction of Stress–Strain Fields Around Interacting Mining Excavations in Jointed Rock: A Comparative Study of Surrogate Models

1. Introduction

The assessment of the stress–strain state (SSS) of fractured rock masses in the interaction zone of mining excavations is a central problem in geomechanics [1,2,3], directly affecting the safety and economic efficiency of engineering decisions [4,5]. The finite element method (FEM) remains the standard tool for solving this problem across a wide range of geotechnical applications [6,7,8]. However, in tasks requiring multiple recalculations—such as probabilistic analyses, parameter optimization, and the comparison of alternative design solutions within intelligent mining frameworks—direct use of FEM becomes computationally expensive. Every new configuration requires mesh regeneration and a full re-solution of the problem, and both meshing and computation remain the bottleneck even with automated workflows.

In this context, surrogate modeling based on machine learning (ML) has emerged as a promising alternative capable of rapidly approximating dependencies identified through FEM. Existing studies show successful applications of surrogate models primarily for predicting integral scalar indicators, such as tunnel stability coefficients in different geological conditions, bearing capacity of foundations, slope stability, and retaining structure performance [9,10,11,12,13]. In such works, the target variables typically represent single aggregated characteristics of system behavior.

In contrast, predicting full spatial fields of the SSS is a considerably more complex task, generally addressed using architectures designed for spatial data processing. One common approach represents geometry and output fields as images and applies convolutional neural networks (CNNs), including conditional generative adversarial networks (cGANs), which have demonstrated the ability to reproduce stress fields for complex geometries and boundary conditions [14,15,16]. Another approach represents the finite element mesh as a graph and uses graph neural networks (GNNs), such as MeshGraphNets, which can significantly accelerate FEM simulations in solid mechanics by operating directly on unstructured meshes [17,18,19]. Although highly accurate, these models require specialized data representations (images or graphs) and are often computationally demanding to implement in routine engineering practice.

To facilitate AI-assisted design in mining, the present study considers an alternative paradigm: encoding the physics of the problem directly into the feature space by designing a set of engineering (contextual) features that describe the mutual arrangement of nodes and excavations as well as the geometry of the system. This strategy enables the application of classical tabular ML algorithms—primarily tree-based ensembles—that are well suited for structured data and provide high computational performance [20]. However, the applicability and comparative efficiency of such algorithms specifically for reproducing full SSS fields around interacting mining excavations remain insufficiently investigated.

The aim of this study is to compare the effectiveness of surrogate models (Random Forest, LightGBM, CatBoost, MLP) for the rapid prediction of SSS fields around two interacting underground excavations based on FEM simulations.

The objectives of the study are as follows:

• To generate a parametric dataset of 1000 FEM simulations and develop an engineering feature space that reflects the geometry of the system and the spatial position of nodes relative to the excavations.
• To train and compare surrogate models in terms of accuracy ($R^2$, MAE) and computational performance (training time, inference speed, model size, and learning curves).
• To evaluate model interpretability using permutation feature importance and to conduct qualitative validation of SSS fields (including plastic zones) against reference FEM results.

2. Materials and Methods

The research methodology (Figure 1) included: parametric FEM modeling to generate a synthetic dataset of simulations; extraction of SSS fields and construction of engineering features; splitting the dataset into training/validation/test subsets with cumulative subsamples for learning curves; hyperparameter optimization using Optuna and model training (Random Forest, LightGBM, CatBoost, MLP); followed by a comprehensive evaluation of model accuracy and computational efficiency.

2.1. Numerical Modeling

The dataset used for model training was generated by conducting a series of 1000 parametric simulations in a geotechnical finite element software package. To automate geometry generation, modification of input parameters, and extraction of simulation outputs, an application programming interface (API) controlled by a Python script was employed.

The computational model (Figure 2) represents a 50 × 50 m rock mass domain containing two parallel horseshoe-shaped underground excavations. The model dimensions were selected to eliminate the influence of boundary effects. The mechanical behaviour of the rock mass was simulated using an elastoplastic constitutive model with a Hoek–Brown yield criterion, developed for fractured rock masses [21,22]. Equivalent Mohr–Coulomb parameters were derived using the procedure proposed by Hoek and Brown [23].

The physical and mechanical properties adopted in the model are presented in Table 1. Input parameters for such constitutive models—including deformability characteristics and shear strength parameters—are typically derived from laboratory tests on rock specimens [24,25]. The computational domain was discretized using 15-node triangular finite elements with local mesh refinement around the excavation contours.

To identify nodes that reached the limiting shear state, a dimensionless relative shear stress parameter ${\tau }_{rel}$ was used:

$${\tau }_{\mathrm{rel}}={\displaystyle \frac{{\tau }_{\mathrm{mob}}}{{\tau }_{\max }}}$$

(1)

where ${\tau }_{\mathrm{mob}}$ is the magnitude of mobilized shear stresses along the local shear plane, and ${\tau }_{\max }$ is the ultimate shear strength computed according to the Mohr–Coulomb criterion: ${\tau }_{max}=c+{\sigma }_n^{\text{'}}\mathit{tan}\phi$.

Within this study, the parameters listed in Table 1 were kept constant, while six global parameters governing the geometry and relative position of the excavations were varied: horizontal spacing and vertical offset between excavation centres, as well as the width and height of each excavation. Parameter values were sampled from predefined uniform distributions. The variation ranges are presented in Table 2.

For each simulation, the following data were extracted from all nodes of the finite element mesh:

Basic features: nodal coordinates (X, Y).

Target variables: six components of the stress–strain state, including effective normal stresses (${\sigma }_{xx}$, ${\sigma }_{yy}$), shear stresses (${\tau }_{xy}$​), strain components (${\epsilon }_{xx}$, ${\epsilon }_{yy}$, ${\gamma }_{xy}$​), as well as total displacement ($U_{tot}$​).

2.2. Feature Space and Its Construction

Since the original nodal coordinates (X, Y) do not provide ML models with explicit information about the geometry of the computational domain, and the global parameters alone are insufficient to capture dependencies explaining the target variables, a set of engineering (contextual) features was developed. These features supply the model with information about the spatial position of each node relative to the excavation contours and the overall system geometry. The full set of input features used for training is summarized in Table 3.

The average excavation width ($\mathrm{mean}\_\mathrm{width}$) and height ($\mathrm{mean}\_\mathrm{width}$) were computed as:

$$\mathrm{mean}\_\mathrm{width}={\displaystyle \frac{\mathrm{widt}{\mathrm{h}}_1+\mathrm{widt}{\mathrm{h}}_2}{2}},$$

(2)

$$\mathrm{mean}\_\mathrm{height}={\displaystyle \frac{\mathrm{heigh}{\mathrm{t}}_1+\mathrm{heigh}{\mathrm{t}}_2}{2}},$$

(3)

where${ \mathrm{width}}_1$, ${\mathrm{width}}_2$​ are the widths of the first and second excavation, and ​ ${\mathrm{height}}_1{, \mathrm{height}}_2$ are their heights.

The aspect ratios (${\mathrm{aspect}}_1$, ${\mathrm{aspect}}_2$) and area ratio ($\mathrm{area}\_\mathrm{ratio}$) were defined as:

$${\mathrm{aspect}}_1={\displaystyle \frac{{\mathrm{width}}_1}{{\mathrm{height}}_1}},$$

(4)

$${\mathrm{aspect}}_2={\displaystyle \frac{{\mathrm{width}}_2}{{\mathrm{height}}_2}}$$

(5)

$$\mathrm{area}\_\mathrm{ratio}={\displaystyle \frac{{\mathrm{width}}_1\times {\mathrm{height}}_1}{{\mathrm{width}}_2\times {\mathrm{height}}_2}}.$$

(6)

The normalized distance ($\mathrm{dist}\_\mathrm{norm}$) and shift ($\mathrm{shift}\_\mathrm{norm}$) were defined as:

$$\mathrm{dist}\_\mathrm{norm}={\displaystyle \frac{\mathrm{distance}}{\mathrm{mean}\_\mathrm{width}}},$$

(7)

$$\mathrm{shift}\_\mathrm{norm}={\displaystyle \frac{\mathrm{vertical}\_\mathrm{shif}}{\mathrm{mean}\_\mathrm{height}}}.$$

(8)

The normalized signed distance ($\mathrm{Signed}\_\mathrm{Dist}\_\mathrm{Norm}$) was given by:

$$\mathrm{Signed}\_\mathrm{Dist}\_\mathrm{Norm}\left(P\right)=\mathrm{sign}\left(P\right)\cdot {\displaystyle \frac{\mathrm{dist}\left(P,B\right)}{W_{\mathrm{mean}}}},$$

(9)

where $\mathrm{Signed}\_\mathrm{Dist}\_\mathrm{Norm}\left(P\right)$ is the feature value for node $P$; $B$ is the set of points forming the excavation contour; $W_{\mathrm{mean}}$ is the mean excavation width; and $\mathrm{sign}\left(P\right)$ is a function returning +1 for rock mass nodes and –1 for excavation nodes.

The overlap index ($\mathrm{Overlap}\_\mathrm{Index}$) was defined as:

$$\mathrm{Overlap}\_\mathrm{Index}\left(P\right)=\exp \left(-{\displaystyle \frac{d_1\left(P\right)}{W_{\mathrm{mean}}}}\right)\cdot \exp \left(-{\displaystyle \frac{d_2\left(P\right)}{W_{\mathrm{mean}}}}\right),$$

(10)

where $d_1\left(P\right)$, $d_2\left(P\right)$ are the minimum distances from node $P$ to the nearest points on the contours of the first and second excavation, respectively, and $W_{\mathrm{mean}}$​ is the mean excavation width.

The curvature feature ($\mathrm{Curvature}$) was defined as:

$$\mathrm{Curvature}\left(P\right)=\mathsf{\kappa }\left(b^{*}\left(P\right)\right)\cdot \exp \left(-{\displaystyle \frac{\mathrm{dist}\left(P,B\right)}{\sigma }}\right)$$

(11)

where $b^{*}\left(P\right)=\mathit{arg}\underset{b\in B\mathrm{dist}}{\mathit{min}}\left(P,b\right)$ is the point on the excavation contour $B$ closest to node $P$; $\kappa \left(b\right)$ is the local curvature of the contour at point $b$, estimated by approximating a local segment with a circular arc; and $\sigma$ is a parameter controlling the decay rate of curvature influence.

The vertical projection feature ($\mathrm{Vertical}\_\mathrm{Projection}$) was defined as:

$$\mathrm{Vertical}\_\mathrm{Projection}\left(P\right)=\mathrm{clip}\left(\mathrm{sgn}\left(\mathsf{\Delta }y\left(P\right)\right)\cdot f\left({\displaystyle \frac{\left|\mathsf{\Delta }y\left(P\right)\right|}{r\left(P\right)}}\right),-1,1\right),$$

(12)

where $\mathsf{\Delta }y\left(P\right)$ is the vertical offset from node $P$ to the nearest point on the excavation contour; $r\left(P\right)$ is the Euclidean distance to the same point; $f:\left[0,\mathrm{\infty }\right)\to \left[\mathrm{0,1}\right]$ is a monotonically increasing function; and $\mathrm{clip}\left(z,-\mathrm{1,1}\right)$ is a function that limits the final value to the range [–1, 1].

The density feature ($\mathrm{Density}\_\mathrm{Excavated}\_\mathrm{Distances}$) was computed as:

$$\mathrm{Density}\_\mathrm{Excavated}\_\mathrm{Distances}\left(P\right)={\left({\displaystyle \frac{1}{\left|E\right|}}\sum _{e\in E} \mathrm{dist}\left(P,e\right)\right)}^{-1},$$

(13)

where $\mathrm{Density}\_\mathrm{Excavated}\_\mathrm{Distances}\left(P\right)$is the feature value for node $P$; $E$ is the set of nodes belonging to the excavated domain; and $\left|E\right|$ is the cardinality of this set.

Visualizations of the engineering feature fields for a test example are shown in Figure 3.

After constructing the feature set, each FEM simulation was treated as a separate sample. The collection of files was randomly shuffled once and split at the file level into training (70%), validation (15%) and test (15%) subsets. Based on the training subset, cumulative subsamples were formed as follows: for k ∈{1–10, 16, 24, 32, 40, 56, 72, 80, 160, 240, 320, 400, 480, 560, 640, 700}, datasets were created containing the first k training files. Each subsequent subsample included all samples from the previous one, which enabled consistent construction of learning curves.

2.3. Machine Learning Algorithms Considered

A separate surrogate model was trained for each of the seven target parameters. The study compared four regression algorithms widely applied in engineering prediction tasks:

Random Forest (RF). An ensemble method based on bagging [26]. The model consists of T independent decision trees, each trained on a bootstrap sample, while only a random subset of features is considered when splitting a node. The final prediction is obtained by averaging the outputs of all trees in the ensemble, as shown in Equation (14):

$$\widehat{y}\left(x\right)={\displaystyle \frac{1}{T}}\sum _{t=1}^T{h}_t\left(x\right).$$

(14)

Gradient Boosting. An ensemble method that constructs a composition $F_T\left(x\right)$ by sequentially adding $T$ weak learners trained to approximate pseudo-residuals [27]. At each iteration, a new model $h_t\left(x\right)$ is fitted to a vector that is anti-parallel to the gradient of the loss with respect to the current ensemble prediction. For the mean squared error (MSE), pseudo-residuals coincide with standard residuals, and the ensemble is updated according to Equation (15):

$$F_{\left(t\right)}\left(x\right)=F^{t-1}\mathrm{ }\left(x\right)+\eta \cdot h_t\left(x\right),$$

(15)

where $\eta$is the learning rate. Two high-performance implementations of gradient boosting were used in this work:

LightGBM [28], which accelerates tree construction via a histogram-based split-search algorithm;

CatBoost [29], which employs ordered boosting to mitigate prediction bias and provides built-in tools for handling categorical variables.

Multilayer Perceptron (MLP). A feed-forward fully connected neural network. Each layer performs an affine transformation of the input vector followed by a nonlinear activation (e.g., ReLU). The network parameters are optimized via backpropagation [30] and gradient-based optimization methods. The loss function includes an L2-regularization term , as defined in Equation (16):

$$\mathcal{L}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left(y_i-\widehat{y}\left(x_i\right)\right)}^2+\mathsf{\lambda }\sum _l‖W_l‖\stackrel{2}{F}.$$

(16)

For tree-based models (RF, LightGBM, CatBoost), input features were not normalized because decision-tree-based algorithms are inherently invariant to monotonic transformations of the feature space. For MLP, all features were normalized using StandardScaler. The normalization parameters were computed exclusively on the corresponding training subset and subsequently applied to both validation and test subsets, preventing information leakage from the test data into the training process.

2.4. Training Methodology and Evaluation Criteria

For each algorithm and each target variable, hyperparameter optimization was carried out using the Optuna framework [31], which implements Bayesian optimization methods. The optimization objective was the coefficient of determination ($R^2$) on the validation set. A uniform computational budget was imposed in this study: 2 hours of hyperparameter search for each algorithm–target pair.

In addition, the following assumption was adopted: hyperparameter optimization was performed once on the full training set, and the resulting optimal configuration was then used to train models on all smaller subsamples. This strategy was chosen to standardise the comparison conditions and reduce the total computational time. However, it should be noted that each subsample may in principle have its own locally optimal set of hyperparameters.

The experiments were conducted on a workstation running Windows 11 Pro with an AMD Ryzen 5–class CPU (6 cores, 12 threads), 32 GB DDR4 RAM, an NVIDIA GeForce RTX 4070 GPU (12 GB VRAM) and a 512-GB NVMe SSD. The GPU was used for training and inference of the MLP (PyTorch) and LightGBM models, while Random Forest and CatBoost were executed on the CPU. The MLP models were trained in a full-batch regime, so the maximum size of the training set for this architecture was limited to 400 simulation cases by the available hardware resources.

All computations were performed in 32-bit precision (single precision), as increasing the numerical precision did not improve the results but led to higher demands on computational resources.

The performance of the surrogate models was assessed and compared using five criteria:

• Global predictive accuracy: The primary metric was the coefficient of determination ($R^2$) computed on the held-out test set. The dependence of $R^2$ on the size of the training data was analysed to construct learning curves and to identify the saturation point of accuracy [26].
• Accuracy in critical zones: To evaluate the ability of the models to predict material behaviour near the limit state, an additional analysis was carried out. For a series of simulations with fixed excavation geometry but varying distance between the excavations, the model predictions were compared with reference values of ${\sigma }_{xx}$, ${\sigma }_{yy}$ and ${\tau }_{xy}$ at nodes that had entered the plastic state (${\tau }_{\mathrm{rel}}=1$. The evaluation metric was the mean absolute error (MAE).
• Computational efficiency: Both training and inference (prediction) speed were assessed:
  • training time required to fit the model on training sets of different sizes;
  • inference time required to generate predictions of the stress–strain fields for a single simulation case from the test set. For the MLP and LightGBM models, these computations were performed on the GPU, whereas RF and CatBoost used the CPU.
• Model compactness: The physical size of the trained model file on disk (in megabytes) was measured as an indicator of storage and deployment requirements.
• Interpretability: The contribution of each input feature to the final prediction was quantified using the Permutation Feature Importance (PFI) method [26]. The method consists in measuring the drop in $R^2$ after random permutation of the values of a given feature on the test set.

3. Results

Figure 4 shows the dependence of $R^2$ on the training set size for each of the seven predicted stress–strain parameters, with separate curves for RF, LightGBM, CatBoost and MLP. The corresponding maximum $R^2$ values are summarised in Table 4.

The analysis demonstrates the high predictive capability of the developed surrogate models for the stress components (${\sigma }_{xx}$, ${\sigma }_{yy}$, ${\sigma }_{xy}$) and the total displacement $U_{tot}$. For these quantities, the tree-based models exhibit a consistent increase in $R^2$ as the size of the training set grows, with the accuracy approaching saturation at approximately 300–400 FEM simulations. When trained on the full dataset, high final $R^2$ values are achieved, in particular 0.967 for ${\sigma }_{yy}$ (LightGBM, CatBoost) and 0.998 for $U_{tot}$ (LightGBM, Random Forest, CatBoost).

A comparison of the algorithms shows that LightGBM, Random Forest and CatBoost provide similarly high performance for the stress components and $U_{tot}$. LightGBM yields the best results for ${\sigma }_{xx}$​ and ${\epsilon }_{xx}$, Random Forest performs best for ${\tau }_{xy}$ and ${\gamma }_{xy}$, while for ${\sigma }_{yy}$, ${\epsilon }_{yy}$ and $U_{tot}$ the differences between the three tree-based models are negligible. The MLP model is consistently less accurate than the tree-based algorithms in all cases.

The lowest $R^2$ values are obtained for the strain components, particularly for ${\epsilon }_{xx}$, where the maximum $R^2$ is 0.601.

Figure 5 presents the dependence of model training time on the number of training samples. The plots reveal fundamental differences in the scalability of the considered algorithms. The maximum training times are summarised in Table 5.

For the gradient-boosting models, the increase in computational cost is the most predictable and efficient. The training time of LightGBM on the full dataset is on the order of ≈130–170 s, while CatBoost requires ≈120–270 s depending on the target parameter.

Random Forest proves to be the most computationally demanding algorithm: for some target parameters, the maximum training time reaches ≈7.8–8.1 × 10³ s, i.e. roughly an order of magnitude higher than for the gradient-boosting models.

The MLP model trains faster than the other algorithms for small and medium training set sizes; however, after a certain threshold (around 160–240 files) the training time increases by an order of magnitude and reaches ≈1.2–2.6 × 10³ s at the maximum dataset size. This behaviour is likely related to the limitations of the available hardware resources and specifics of the training implementation.

In addition to accuracy and training time, an important practical aspect is the compactness of the trained models, i.e. the physical size they occupy on disk. Figure 6 shows the dependence of model size on the number of training samples, and the corresponding maximum values are reported in Table 6.

The MLP model is the most compact and is practically insensitive to the size of the training set: the file sizes remain in the range of ≈0.3–0.6 MB for all target parameters.

The gradient-boosting models differ substantially in compactness. The size of LightGBM models stabilises at ≈23–25 MB, whereas for CatBoost the model size after an initial growth reaches ≈70–160 MB, depending on the target parameter.

In contrast, the Random Forest models exhibit a pronounced and almost linear increase in size with the number of training samples, reaching several gigabytes on the full dataset.

A key indicator governing the applicability of surrogate models for rapid assessment tasks is the inference speed, i.e. the time required to generate predictions for new data. Figure 7 compares the prediction times of models trained on datasets of different sizes with the reference time of a single FEM simulation. The maximum inference times are summarised in Table 7.

The results indicate a substantial (one to two orders of magnitude) speed-up in prediction for most models compared with direct numerical simulation (4.86 s).

The highest and most stable prediction speed is achieved by the MLP: its inference time does not depend on the training set size and is ≈0.004–0.007 s, corresponding to a speed-up of about 700–1200 times relative to FEM.

The gradient-boosting models are also efficient. For LightGBM, the maximum inference time lies in the range ≈0.14–0.29 s, providing a speed-up of approximately 15–30 times. For CatBoost, the values are ≈0.07–0.15 s, which corresponds to a speed-up of about 30–70 times.

For Random Forest, a pronounced dependence of inference time on model size and complexity is observed. For some target parameters, the prediction time remains at 0.1–0.4 s (up to ≈40-fold speed-up), but for the shear components (${\tau }_{xy}$, ${\gamma }_{xy}$) it increases to 8–16 s, making the inference comparable to or even slower than a direct FEM simulation and severely limiting the practical applicability of this algorithm.

To assess the ability of the models to correctly predict material behaviour near the limit state, an additional analysis was performed. The prediction accuracy was evaluated at nodes that had entered the plastic regime according to the volumetric Hoek–Brown (HB) criterion or had reached the limiting shear state on the excavation boundary (${\tau }_{rel}=1$). The metric used was the mean absolute error (MAE), averaged over the three stress components.

Figure 8 shows the dependence of MAE on the distance between the excavations C. The training data covered the range $C\in \left[\mathrm{1,8}\right]$ m; therefore, the points at $C < 1$ m and $C > 8$ m characterise the behaviour of the models under extrapolation.

Within the training range ($C\in \left[\mathrm{1,8}\right]$ m), the gradient-boosting models provide the best accuracy: the average MAE is 197.3 kPa for LightGBM and 217.1 kPa for CatBoost. Random Forest yields a slightly higher error level (245.52 kPa), and the worst result is obtained with MLP (358.2 kPa).

The behaviour of the models outside the training range shows the expected degradation in accuracy. When extrapolating to small distances ($C < 1$ m), where stress concentrations and non-linear effects are most pronounced, a sharp increase in error is observed for all models; a similar trend is seen for $C > 8$ m.

To quantify the relative contribution of the input features to the model predictions, the Permutation Feature Importance (PFI) method was applied. The method consists in measuring the drop in $R^2$ after randomly permuting the values of a given feature on the test set; this drop is directly proportional to the importance of the feature. The results are shown in Figure 9, where the numerical values correspond to the normalised feature importance (in %).

The heatmaps reveal a high degree of consistency among the tree-based models (RF, LightGBM, CatBoost). For these algorithms, the main share of importance (on average about 70–75%) is attributed to the engineered features ($V\mathrm{ertical}\_\mathrm{Projection}$, $\mathrm{Signed}\_\mathrm{Dist}\_\mathrm{Norm}$, $\mathrm{Curvature}$, $\mathrm{Overlap}\_\mathrm{Index}$, etc.), whereas the raw coordinates X and Y account for the remaining ≈25–30%. For the normal components ${\sigma }_{xx}$, ${\sigma }_{yy}$, ${\epsilon }_{xx}$, ${\epsilon }_{yy}$, the contribution of the engineered features reaches 80–90%, while for the shear components and total displacement (${\tau }_{xy}$, ${\gamma }_{xy}$, $U_{tot}$) their share is about 55–60%. At the same time, global parameters such as $\mathrm{mean}\_\mathrm{width}$, $\mathrm{mean}\_\mathrm{height}$, $\mathrm{aspect}1/2$, $\mathrm{dist}\_\mathrm{norm}$, $\mathrm{shift}\_\mathrm{norm}$ and $\mathrm{area}\_\mathrm{ratio}$ show zero importance for all models.

In the updated setting, the MLP model preserves the general tendency of dominance of the engineered features, but redistributes importance differently. For most normal components (${\sigma }_{xx}$, ${\sigma }_{yy}$, ${\epsilon }_{xx}$, ${\epsilon }_{yy}$) and for the total displacement $U_{\mathrm{t}\mathrm{o}\mathrm{t}}$, the key predictor becomes $\mathrm{Overlap}\_\mathrm{Index}$, whose contribution reaches ≈45–55% of total importance, while the shares of $\mathrm{Signed}\_\mathrm{Dist}\_\mathrm{Norm}$ and $\mathrm{Vertical}\_\mathrm{Projection}$ are noticeably lower than in the tree-based models (typically 3–15%). For the shear components (${\tau }_{xy}$, ${\gamma }_{xy}$), the distribution is more balanced: X, $\mathrm{Vertical}\_\mathrm{Projection}$ and $\mathrm{Overlap}\_\mathrm{Index}$ all retain a substantial influence.

To qualitatively assess the ability of the surrogate models to approximate the spatial distribution of the stress–strain fields, a visual comparison was performed for a randomly selected test case. As an illustration, Figure 10 shows triptychs for ${\sigma }_{xx}$ for each algorithm. Each triptych includes the reference field from the FEM simulation, the field predicted by the model, and the map of absolute errors $\left|Y_{pred} - Y_{true}\right|$. To improve readability and suppress the influence of rare outliers, the colour scale on the error maps is clipped at the 99th percentile of the error distribution.

The obtained absolute error maps are generally consistent with the $R^2$ values (Table 4). The largest discrepancies for all models are localised in the immediate vicinity of the excavation boundaries, whereas in the bulk of the rock mass the errors remain relatively small. For LightGBM and Random Forest, the errors are mainly concentrated near the contours, with Random Forest providing a slightly more accurate reproduction of the fields; CatBoost is characterised by artefacts in the form of radially diverging “rays” emanating from the excavations. In the case of MLP, a fundamentally different pattern is observed: for ${\sigma }_{xx}$, the model fails to reproduce the detailed structure of the target field, and the predicted distributions are noticeably smoother than those of the tree-based models, which leads to an underestimation of local stress extrema around the excavations.

4. Discussion

The comparison shows that, in terms of the coefficient of determination $R^2$, the best performance is achieved by the gradient-boosting models (LightGBM, CatBoost) and Random Forest [32]. At the same time, the qualitative analysis of the fields revealed differences in the error patterns: for gradient boosting (CatBoost in particular), numerical artefacts were observed (radial “rays”, non-physical elongated zones of concentration), whereas Random Forest more frequently reproduced field patterns closer to the reference FEM solutions. The MLP exhibited smoothing, accompanied by a loss of local extrema.

From the viewpoint of computational characteristics, fundamental trade-offs were identified. LightGBM and CatBoost train rapidly on the full dataset (≈130–170 s and ≈120–170 s, respectively), with model sizes of ≈18–25 MB (LightGBM) and ≈70–160 MB (CatBoost), and stable inference times of ≈0.14–0.29 s and ≈0.07–0.15 s. Random Forest proved to be the most resource-demanding: training time reaches ≈7.8–8.1 × 10³ s, model size grows to tens of gigabytes, and inference time for some target variables increases to 8–16 s, making these configurations comparable to or slower than a direct FEM run (4.86 s), which negates the advantage of using surrogates for real-time mining monitoring tasks The MLP provides the lowest inference latency (≈0.004–0.007 s, a speed-up of about 700–1200 times relative to FEM) with very compact models (≈0.5 MB), but it is less accurate than the tree-based models.

The permutation feature importance analysis showed that the tree ensembles rely primarily on engineering features , validating the proposed physics-informed approach. For the vertical components (${\sigma }_{xx}$, ${\epsilon }_{yy}$), $\mathrm{Vertical}\_\mathrm{Projection}$ dominates (in total ≈50–65%); for ${\sigma }_{yy}$, ${\epsilon }_{yy}$, $\mathrm{Curvature}$ and $\mathrm{Signed}\_\mathrm{Dist}\_\mathrm{Norm}$ are particularly important; for ${\tau }_{xy}$, ${\gamma }_{xy}$, the contribution of X is large; and for $U_{\mathrm{t}\mathrm{o}\mathrm{t}}$, Y and $\mathrm{Vertical}\_\mathrm{Projection}$ play a leading role. The high importance of geometric and positional features is typical for stability assessment problems and is consistent with previous studies [11,13,33].

The dependence of accuracy on the data volume exhibits saturation for most targets at ≈300–400 simulations (30–40% of the dataset); further increases in $N$ yield only limited improvements in R². For the strains, especially ${\epsilon }_{xx}$, the achieved values are lower, which indicates the need to enrich the feature space with descriptors more sensitive to local deformation effects. The analysis of accuracy in plastic zones (according to the Hoek–Brown criterion) confirmed that prediction errors increase significantly when the distance between excavations C goes beyond the training range ($C\in \left[\mathrm{1,8}\right]$ m). This highlights the limited ability of the models to approximate behavior governed by the HB criterion outside the training domain—a critical consideration for safety assessments in mining design when solving extrapolation problems.

From a practical point of view, for rapid assessment and multi-scenario analysis, CatBoost and LightGBM are optimal in terms of the combination of accuracy, training time and model size, whereas MLP provides the lowest inference latency provided that adequate regularisation is applied. Although Random Forest can reproduce the stress–strain fields more faithfully than gradient boosting, its use is not justified for large datasets.

The practical advantages of the surrogate models are expected to be most significant for three-dimensional layouts, where a single FEM simulation may require hours or days. In such cases, surrogates can replace FEM in tasks involving a large number of forward evaluations, while FEM is retained for detailed verification. Typical examples include probabilistic analyses of pillar stability, optimization of mining excavation layouts, and routine screening of design alternatives within Digital Twin frameworks. For these applications, increasing the dimensionality of the problem primarily affects the cost of generating the training dataset, whereas the prediction time of the LightGBM, CatBoost, and MLP models remains essentially constant.

5. Conclusions

This study has demonstrated and systematically evaluated the applicability of classical tabular machine-learning algorithms for constructing surrogate models capable of predicting, with high accuracy and speed, the full stress–strain fields in a rock mass around a system of two interacting mining excavations.

It has been confirmed that tree ensembles (Random Forest, LightGBM, CatBoost) provide high predictive accuracy, reaching $R^2>0.967$ for the stress components and $R^2>0.998$ for the total displacements. The most important factor for model accuracy is the use of engineering features, whose importance is consistent with the physical pattern of stress distribution in jointed rock masses.

The comparative analysis revealed fundamental trade-offs between the algorithms suitable for intelligent design:

• The MLP offers the lowest inference latency (≈0.004–0.007 s, a speed-up of 700 to 1200 times relative to FEM) with compact models, but produces smoother fields with partial loss of local features.
• LightGBM and CatBoost provide the most balanced combination of accuracy and computational cost ($R^2$ close to Random Forest, speed-up of 15 to 70 times relative to FEM), making them ideal for iterative design optimization.
• Random Forest yields fields closest to reference FEM solutions but is too resource-demanding for large-scale deployment.

It has been established that model accuracy saturates at ≈300–400 parametric simulations, allowing for the optimization of computational costs during data preparation. The analysis of the mean absolute error (MAE) in plastic zones showed that within the training range ($C\in \left[\mathrm{1,8}\right]$ m), LightGBM attains the lowest error values. Outside this range, the error increases systematically, reflecting the limitations of data-driven models in extrapolation.

Overall, the proposed approach represents an effective tool for rapid stress–strain assessment at the design stage, substantially reducing computational costs. These properties make the proposed surrogates suitable as core components of probabilistic and optimization-based design workflows for underground mining. Large ensembles of scenarios can be evaluated using surrogate models, while detailed FEM analyses are reserved for critical cases. This division of roles allows engineers to account for parameter uncertainty and explore alternative layouts systematically, supporting the transition towards AI-assisted intelligent mining.