Development of neural network-based surrogate models in geotechnics: a state-of- the-art review

Geotechnical engineers must make informed decisions based on comprehensive analyses of complex structures. These decisions directly influence the final design implemented in construction or manufacturing. Advances in computational engineering and the increasing availability of computational power now allow for high-fidelity numerical simulations to assess structural behaviour under various loading conditions. However, the growing computational cost and time required for these simulations hinder the desired shift towards efficient sustainable design processes, as formulated in the United Nations Sustainable Development Goal 9 – Industry, Innovation, and Infrastructure (United Nations, 2015). Surrogate modelling is an approach that seeks to develop fast mathematical approximations of complex and computationally expensive problems. Rather than allocating significant computational resources to solve a single specific problem, a surrogate model is constructed to provide efficient predictions across a wide range of scenarios (Forrester et al., 2008).

Surrogate models provide a powerful solution for approximating computationally expensive simulations or mapping field measurements, enabling fast and accurate predictions in unexplored design spaces (Zheng et al., 2024). The complexity of modern numerical models, which incorporate multiple boundary conditions and constitutive laws with numerous parameters, demands substantial computational resources, limiting the possibility of obtaining real-time responses through back-analysis, a process aimed at understanding the behaviour of a structure based on measured features (Walton and Sinha, 2022; Zhao et al., 2021; Xing et al., 2025; Li et al., 2021). Complex constitutive models, such as the hypoplastic model developed by von Wolffersdorff (1996), often require preliminary parameter calibration based on element tests. Simplifying the model can reduce computational cost, but the reliability of numerical predictions must be carefully considered when deciding on the model setup (Wiles, 2006). Traditional modelling approaches still struggle to fully capture the combined influence of multiple factors on the non-linear behaviour of geomaterials (Cui et al., 2022; Zhao et al., 2023). In this context, neural network-based surrogate models emerge as a promising solution, offering a combination of accuracy and efficiency suitable for real-time back-analysis (Qi and Fourie, 2018; Fathi et al., 2024; Zhao et al., 2024; Wei et al., 2023; Chen et al., 2025) and structural health monitoring (SHM) (Mahar et al., 2025; Chen et al., 2025; Hlaing et al., 2023; de Nolasco Santos et al., 2024).

The reduced time required to obtain information about the behaviour of the primary system makes surrogate models highly advantageous across multiple engineering fields. They have been applied in preliminary design (Tang et al., 2024; Heidebrecht and MacManus, 2019), sensitivity analysis ( Yi et al., 2025 ; Aufiero and Janson, 2022; Mohamed et al., 2025), reduced-order modelling (Hu et al., 2024; Zhao, 2021; Zancanaro et al., 2021; Yang et al., 2026) and design optimisation (Jiang et al., 2020; Yondo et al., 2018; Jiao et al., 2025; Chang et al., 2026; Giacomini and Huerta, 2026). Moreover, once trained, these models are user-friendly, enabling non-specialist users to perform tasks such as short-notice maintenance checks by simply providing new input data, without requiring expert knowledge.

This paper provides a state-of-the-art review of neural network-based surrogate models in geotechnical engineering, covering the entire development cycle, from data set generation and sampling strategies to model evaluation and performance assessment. Various types of surrogate models are examined and compared, along with their potential applications. The study explores state-of-the-art approaches for data sampling and hyperparameter tuning, as well as key neural network architectures suited to geotechnical problems. Furthermore, performance metrics used for model assessment are discussed. Finally, the paper highlights the comparative advantages of different methods and outlines future research directions for neural network-based surrogates in geotechnical engineering.

In the following, Section 2 presents various approaches for neural network-based surrogate modelling. Section 3 elaborates on the design of experiment (DoE) for data generation and preprocessing. Sections 4 and 5 discuss the sampling and tuning methods, while Section 6 addresses activation functions. In Section 7, evaluation metrics are described. Sections 8 and 9 conclude the paper by summarising the discussed topics and giving recommendations for future research.

Surrogate models are commonly employed to augment the outcomes of computationally expensive simulations performed across a range of input parameters defined by a specific strategy. Two fundamental assumptions, often at odds with each other, underpin surrogate modelling: (1) the developed model significantly accelerates analysis and (2) it maintains sufficient accuracy, even when predicting beyond known data points. Surrogate models function as curve-fitting tools, learning the input–output mapping from a limited data set generated by a primary source (e.g. an expensive simulation) and subsequently generalising to new inputs (Forrester et al., 2008), see Figure 1. For less complex tasks, surrogate modelling can be seen as a non-linear inverse problem for which the goal is to identify a function f that provides input–output mapping based on a limited number of data points (Queipo et al., 2005). In general, the available training data cannot provide sufficient information to uniquely identify f. As a result, multiple surrogate models $fi^$ may be consistent with the available data and exhibit comparable predictive accuracy. This is illustrated in Figure 2, where both surrogate functions achieve similar accuracy with respect to the original data points.

When distinguishing surrogate models based on the origin of the training data, two groups can be identified: i) models developed using synthetic data generated deterministically from physical laws and ii) those trained on experimental or field measurements that directly link observed data to the model’s response. Multiple surrogate modelling approaches developed over the years range from data-driven machine learning models to purely mathematical approaches, such as:

• response surface method (RSE) (Box and Draper, 1987);

• kriging (Krige, 1951);

• radial basis function (RBF) (Broomhead and Lowe, 1988);

• polynomial chaos expansion (PCE);

• support vector regression (SVR) (Cortes and Vapnik, 1995); and

• artificial neural networks (ANNs).

This study focuses solely on the application of various types of ANNs for developing regression surrogate models. The training process of a neural network-based surrogate model required examination of the loss functions defined for training and validation data sets. Figure 3 shows the training and validation loss curves, with the regions of underfitting and overfitting indicated with their respective conditions. Typically, during neural network training, checkpoints are saved at each epoch and the model with the lowest validation error is retained. Figure 4 shows the steps required for developing a neural network-based surrogate model. ANNs are data-driven computational models composed of interconnected neurons arranged in layers, inspired by biological neural systems (McCulloch and Pitts, 1943). The simplest example of a single-layer ANN with a linear activation function is linear regression, which is commonly used for analysing scattered data because it provides meaningful average relationships between input and output with low computational costs. Utilisation of non-linear activation functions in multiple, mutually interconnected neurons provides ANNs with universal approximation capability, making them highly effective in capturing non-linear behaviours commonly encountered in civil and geotechnical engineering. This enables ANNs to approximate complex mappings from inputs (e.g. material properties or loading conditions) to outputs (e.g. displacements or safety factors). Despite their effectiveness, ANNs are criticised as “black-box” models due to limited interpretability, often requiring extensive data and careful hyperparameter tuning. Nevertheless, their strong predictive capabilities, versatility and availability of robust computational frameworks such as TensorFlow (Abadi et al., 2015), JAX (Bradbury et al., 2018) and PyTorch (Bradbury et al., 2018), that provide the numerical tools and infrastructure needed to efficiently build, train and apply ANNs, continue to drive their widespread adoption in engineering practice.

The following sections introduce a range of neural network-based surrogate models commonly applied in geotechnical engineering. However, the methods discussed are not exhaustive, omitting, for instance, convolutional neural networks (CNNs). The selection of an appropriate neural network architecture is highly case-dependent, as the characteristics of the available data influence which network type and activation function are most suitable for a given problem. Therefore, it is essential to define and evaluate the potential inputs and desired outputs in advance, prior to commencing surrogate modelling.

The most common variant of ANNs is the multi-layer feedforward neural network (FNN), which consists of fully connected input, hidden and output layers with backpropagation (Rumelhart et al., 1986). The neurons use weighted connections and activation functions to transform input signals into outputs. The absence of cycles or feedback loops distinguishes FNNs from other architectures, ensuring that information only moves forward, preserving the independence of the input-out-pairs, which contribute to the simplicity and stability of the model. This structure makes FNNs highly suitable for supervised learning tasks like regression and classification, where the objective is to map static, independent input parameters to a known output value. However, they are limited in handling sequential, time-dependent or spatially correlated data and may require extensive training data to achieve generalisable results.

FNNs have been successfully applied in the surrogate engineering of geotechnical phenomena. Some applications include tunnel stability in karst formations (Kovačević et al., 2021), where an FNN-based surrogate predicts the distribution of tunnel displacement, rock bolt forces and stresses in the tunnel lining, based on three parameters describing rock mass. This model has been shown to be reliable for defining design parameters for tunnelling in karstic rock masses, as validated by on-site measurements. Nguyen-Minh et al. (2025) applied an FNN using Levenberg–Marquardt backpropagation for surrogate modelling of strip anchors in frictional–cohesive soils with surcharge. Their results demonstrated that FNNs can accurately predict the seismic uplift capacity of anchors based on six input factors describing soil conditions, loading and geometry, outperforming non-neural models such as extreme gradient boosting and multi-variate adaptive regression splines. The application of FNNs for predicting soil parameters, subsequently used in a numerical model for long-term embankment settlement prediction, was investigated by Wang et al. (2023), who implemented a Monte Carlo-based dropout of connections between dense neurons to account for model uncertainty.

The physics-informed neural network (PINN) represents a computational approach that integrates the governing physical laws of a system, for example, conservation laws, that can be expressed in differential form. Typically, they are implemented as partial differential equations (PDEs), directly within the neural network training process (Lagaris et al., 1998). This integration is achieved by summing the residuals of the governing equations $LPDE$⁠, together with the residuals regarding the boundary $Lbc$ and initial conditions $Lic$⁠, into a total loss function $Ltotal$ that accounts for all residuals and other types of errors. In this manner, PINNs ensure coherence between the data-driven loss component $Ldata$ and the underlying physical principles across the entire solution domain, as presented in equation (1). Such a formulation enables the network to generate physically consistent predictions, even in scenarios where data are sparse or contaminated by noise, thereby enhancing model robustness and generalisation capability. The physical fidelity of the surrogate model can be controlled by appropriately weighting the physical components (⁠$LPDE$⁠, $Lbc$ and $Lic$⁠) of the loss function. The approach exploits automatic differentiation to efficiently evaluate the derivatives required for enforcing PDE constraints, thereby eliminating the need for conventional numerical-based discretisation methods (Raissi et al., 2019). Consequently, the selection of an appropriate activation function is essential, as it must support higher-order derivatives to ensure stable and effective PINN training (Dashtbayaz et al., 2024). PINNs effectively integrate limited data sets with governing physical laws to produce realistic and generalisable predictions beyond the training domain. However, their high computational cost, sensitivity to loss function weighting and convergence challenges in highly non-linear geotechnical problems remain significant limitations:

PINNs have emerged as a promising approach for surrogate modelling in geotechnical engineering, owing to their capacity to address challenges associated with noisy or incomplete data, which are frequently encountered in geotechnical applications. Ito et al. (2024) applied PINNs employing fully connected dense layers to simulate soil consolidation and unsaturated seepage processes, while Zhou et al. (2025) modelled large-strain consolidation using a self-adaptive variant of PINNs based on an FNN architecture, which dynamically adjusted the weighting of loss function components during training. Li et al. (2025) utilised FNN-based PINNs to predict temperature fields during artificial ground-freezing operations. Zhang et al. (2024) developed a PINN-based surrogate with fully connected layers capable of performing three-dimensional stability analyses of soil slopes subjected to seismic excitation, integrating an upper-bound limit analysis describing the kinematic failure mechanism into the loss function. Vahab et al. (2023) proposed a PINN-based surrogate model using dense, fully connected layers for piles under vertical loading, which was applied both for forward prediction of displacement fields within the model domain and for inverse analysis to identify soil material parameters from axial strain profiles along the pile. Al-Atroush (2025) incorporated Meyerhof’s bearing capacity equations into a PINN framework comprising dense layers to estimate the axial capacity of drilled shafts, emphasising that appropriate selection of loss function weights is critical to prevent the model from becoming overly constrained by the physical equations and thereby losing its inherent data-driven capability. Liu et al. (2025) enhanced a PINN model for seismic site response analysis by introducing Fourier feature embedding into a fully connected architecture for improved learning of high-frequency components. Li et al. (2026) developed a hybrid surrogate modelling framework combining a PINN based on an FNN architecture for predicting pore water pressure in slopes, incorporating governing equations, together with a CNN-based model for probabilistic slope stability assessment.

Long short-term memory (LSTM) networks are a specialised form of recurrent neural networks (RNNs) developed to efficiently handle sequential or time-dependent data (Hochreiter and Schmidhuber, 1997). Traditional RNNs struggle to process long data sequences, limiting their applicability to long-term problems (Hochreiter, 1991). The ability of LSTMs to process long, ordered data sequences makes them suitable not only for analysing time-dependent problems, but also for processing spatially ordered data, such as discretised displacements along structures. Unlike FNNs, where inputs are treated independently, and information flows only in one direction, LSTMs incorporate feedback connections and internal memory cells that retain information across time steps. Their core feature lies in a gating mechanism, comprising input, forget and output gates, that regulates the flow of information within each LSTM unit, allowing the network to preserve or discard data as required. In LSTMs, an activation function controls signal propagation to subsequent layers, while a recurrent activation function operates within the gating mechanisms, regulating information flow across time steps. This architecture more effectively addresses the vanishing and exploding gradient problems (the extreme decrease and increase of weighting updating values), by employing constant error carousel, thereby facilitating more stable training and improved learning of long-term temporal dependencies (Hochreiter, 1991). The bidirectional LSTM extends this concept by incorporating two layers that analyse data in both forward and backward directions simultaneously (Schuster and Paliwal, 1997). LSTMs are particularly advantageous for modelling long-term dependencies and capturing complex temporal or sequential relationships, but they are computationally demanding, prone to overfitting with limited data and require extensive hyperparameter tuning to achieve optimal performance. In addition, despite their design to retain information over long sequences, LSTMs can gradually forget earlier inputs due to the iterative updates of their hidden and cell states, which can reduce their effectiveness for problems requiring the retention of information from distant time steps.

LSTMs have been successfully applied to predict the real-time penetration rate of tunnel boring machines by coupling measured time-dependent data series with dense neurons representing rock mass parameters (Gao et al., 2021). The proposed model demonstrated superior performance compared with conventional RNN models and the autoregressive integrated moving average with exogenous variables. Tao et al. (2024) employed a hybrid approach combining LSTM and CNNs for time-dependent surrogate modelling of braced excavations. Their study highlighted the concept of transfer learning in surrogate modelling, where a pretrained model is fine-tuned using limited monitoring data collected during the initial stages of excavation. Pretraining was performed using both numerically generated synthetic data sets and data from an adjacent excavation constructed in advance. Cerek et al. (2026a) applied a hybrid framework for quay wall stability assessment, combining static input features related to geometry and structural stiffness, processed by FNN layers, with sequential discretised horizontal displacements processed by LSTM layers to predict the maximum bending moment of a sheet pile wall. Qi and Fourie (2018) utilised LSTMs in rock engineering to perform real-time back-analysis of rheological parameters of rock masses based on tunnel response, while Cerek et al. (2025a) applied LSTMs to predict the Young’s modulus of soil from stress–strain data. Xu and Niu (2018) used them to forecast landslides by predicting the cumulative displacement of slopes from historical records of six triggering factors. The application of LSTM-based constitutive models was investigated by Guan et al. (2023) in finite element geotechnical analysis and by Zhang et al. (2024) for modelling the stress–strain behaviour of soil. Building on these advancements, Cerek et al. (2024, 2025b) proposed an approach for predicting the stress–strain curve of soils based on initial results from constant rate of strain tests.

Transformer models represent a class of deep learning architectures originally introduced by Vaswani et al. (2017), characterised by their use of self-attention mechanisms to capture complex dependencies within data sequences. Unlike recurrent or convolutional neural networks, which process information sequentially or within local neighbourhoods, transformers operate on entire sequences simultaneously, allowing each input element to attend to all others through learnable attention weights. This parallel processing capability not only enhances computational efficiency but also enables the model to learn long-range relationships more effectively than architectures such as RNNs or LSTMs. The self-attention mechanism dynamically assigns importance to different parts of the input sequence, thereby improving learning and contextual understanding. Transformers have demonstrated outstanding performance in natural language processing, time series forecasting and spatial–temporal modelling, and are increasingly being explored in engineering domains. Transformers offer superior scalability, the ability to handle complex dependencies and efficient parallel computation. However, these benefits come at the cost of high memory requirements and computational complexity, which may challenge applications with small data sets or physics-constrained problems.

Transformer models present considerable potential for geotechnical engineering, though further research is required to realise and validate their broader applicability. For instance, Zhang et al. (2024) utilised transformers to predict the load envelope of monopile foundations for offshore wind turbines. Their study focused on developing a transformer-based surrogate model that integrates static input features with time-dependent data to predict the horizontal force and bending moment acting on a monopile, ensuring that these do not exceed a maximum allowable rotation of 0.25°. Zhang et al. (2024) based their study on synthetic data generated using a numerical model which was validated against centrifuge tests. Chang et al. (2024) proposed a dual-framework approach for predicting the posture of tunnelling shield machines, combining an FNN-based surrogate model conveying physical information, also trained on numerically generated data for the studied project, with a Transformer model trained on data from previous comparable projects. Their method leverages transfer learning by integrating historical and numerical data sets to improve predictive performance.

The DoE is crucial for developing high-quality surrogate models that accurately capture input–output relationships using a limited number of data points (Blondet et al., 2018). For simple, low-dimensional problems, exhaustive exploration of the parameter space is feasible. However, as dimensionality k increases, the required observations $nk$ grow exponentially (Bowman et al., 1993), a phenomenon known as the curse of dimensionality (Bellman, 1957). In high-dimensional spaces, machine learning models require exponentially more data to achieve reliable performance, making methodological sampling strategies essential.

Preliminary sensitivity analysis (screening) can reduce the number of design variables by identifying parameters that minimally influence model outputs. Methods such as Morris’ elementary effects (Morris, 1991) provide a computationally inexpensive way to qualitatively screen for influential variables, while Sobol indices (Sobol, 1993) quantify the contribution of each input and their interactions to output variance, allowing low-impact variables to be fixed or removed. More recently, Shapley additive explanations (Lundberg and Lee, 2017) can be applied after generating a data set to quantify feature importance for a trained surrogate model, guiding further simplification by highlighting which variables most strongly affect predictions.

Once relevant variables are identified, dimensionality reduction techniques can be applied to further simplify the model. Linear methods like principal component analysis (Jolliffe, 2002) or proper orthogonal decomposition (Berkooz et al., 1993) project data onto a lower-dimensional space, retaining the most significant modes of variation. Non-linear approaches, such as autoencoders, use neural networks to learn low-dimensional latent representations of complex data sets. These reduced-order representations improve computational efficiency, enhance convergence and reduce overfitting, especially in high-dimensional or highly correlated geotechnical data.

Proper data preprocessing is essential to ensure effective processing by neural networks. Common techniques include data normalisation, which scales inputs to a fixed range (e.g. [0, 1] or [−1, 1]), and standardisation, which centres data to zero mean and unit variance. Normalisation is particularly useful when input features have different units or magnitudes, as it improves training stability, prevents features with larger scales from dominating gradient updates and facilitates faster convergence. Standardisation is advantageous for models that assume normally distributed inputs and helps maintain consistent gradient propagation. However, both methods have limitations: normalisation is sensitive to outliers and can become problematic if input or output values exceed the predefined scaling range, while standardisation can be distorted by extreme values. Logarithmic or Box–Cox transformations can be applied to handle skewed distributions or highly non-linear scales. Careful preprocessing ensures that surrogate models are computationally efficient, train reliably and can accurately capture the underlying physics of complex geotechnical systems, while mitigating common training challenges associated with heterogeneous or high-magnitude data sets.

Careful selection of a sampling strategy is crucial for surrogate model development based on synthetic data, as it affects both predictive accuracy and computational efficiency. In high-dimensional geotechnical problems, exhaustive exploration is impractical, so well-distributed samples are needed to capture key trends and interactions while minimising redundancy. Efficient sampling reduces computational cost, improves neural network convergence and enhances model reliability across the input space. The following sections discuss some of the commonly used sampling methods in geotechnical engineering.

Grid search offers a systematic approach by exhaustively evaluating all possible combinations of hyperparameters, which are the external configuration settings defined prior to training, within a predefined grid (Bergstra and Bengio, 2012). However, the number of combinations grows exponentially with the number of hyperparameters, rapidly exceeding practical computational limits.

Random search provides a straightforward alternative to manual sampling and grid search by randomly sampling variable combinations and evaluating their performance. This method has been shown to be more efficient in high-dimensional spaces, which are common in engineering problems, as grid search often allocates excessive resources to less relevant parameters (Bergstra and Bengio, 2012). In addition, random search is inherently parallelisable, requiring no modifications to the sampling algorithm and can incorporate prior knowledge by specifying probability distributions from which to sample.

Latin hypercube sampling (LHS) is a stratified random sampling method designed to ensure space-filling coverage of the input domain (McKay et al., 1979). The range of each input variable is divided into evenly spaced strata and samples are selected such that each stratum is represented exactly once for each variable. This prevents clustering of points and ensures that projections of the sample set onto any individual input axis are well distributed. LHS is widely used in surrogate modelling because it provides better uniform coverage of high-dimensional spaces than purely random or grid-based sampling (Kamath, 2022; McKay et al., 1979). LHS maintains a more balanced representation of the input space while preserving statistical properties such as variance. In addition, LHS can be enhanced through optimisation techniques such as maximin LHS or centred LHS (Stein, 1987), which further improve space-filling properties. Despite these advantages, LHS can face challenges in very high-dimensional problems, where maintaining uniform coverage across all dimensions becomes larger and sample correlations may inadvertently arise (Helton and Davis, 2003). General options to address these limitations include increasing the sample size, employing optimised LHS variants such as orthogonal or scrambled LHS and combining LHS with sensitivity analysis to prioritise the most influential variables (Iooss and Lemaître, 2014). By doing so, high-dimensional input spaces are adequately explored while mitigating computational cost.

LHS has been widely applied in geotechnical engineering. For instance, Yang et al. (2026) utilised LHS for sampling generation for their two-step surrogate model for high-dimensional geotechnical outputs. Ruiz López et al. (2024) used LHS for generating samples for their FNN-based surrogate model of vertical shafts in London clay. Similarly, Ferrero et al. (2023) applied LHS for generating samples for an FNN-based surrogate model of braced excavations for real-time back-analysis. Yang et al. (2024) used LHS for generating a surrogate model of a monopile based on studies from the PISA Dunkirk project (Taborda et al., 2020). LHS has also been applied for conducting numerical studies of geotechnical problems, such as thermal integrity profiling of concrete piles (Sánchez Fernández et al., 2025).

Quasi-random low-discrepancy sequences, such as Sobol (1967) and Halton (1960) sequences, provide an alternative approach to space-filling sampling. Unlike purely random sampling, they deterministically generate points that uniformly cover the input space while minimising clustering. Quasi-random designs serve as an alternative to LHS and are particularly useful for global sensitivity analysis and uncertainty propagation, where uniform coverage of the input domain is critical. However, their performance can deteriorate in very high-dimensional problems due to correlations between dimensions, and generating large numbers of points while maintaining low discrepancy can be computationally challenging. The manual initialisation of the sampling process might lead to bias or non-uniform sampling within complex domains. Wang et al. (2023) utilised a weighted low-discrepancy sequence and a Sobol sequence to generate data for the probabilistic reliability analysis of soil slopes. Similarly, Hu et al. (2024) applied both Halton and Sobol sequences to construct the input parameter space for a reduced-order model of slope stability.

Adaptive sampling methods, introduced by Sacks et al. (1989), iteratively add training points in regions where the surrogate model needs improvement. Unlike one-shot designs, they dynamically select points based on current evaluation metrics, focusing on areas with high uncertainty, prediction error or critical thresholds. This approach is especially beneficial for non-uniform response surfaces with sharp gradients, discontinuities or local extrema. By refining the model in complex regions, adaptive sampling achieves higher accuracy with fewer expensive simulations compared to fixed design strategies. The main challenge of adaptive sampling lies in the implementation. The coupling between the neural network surrogate model and the sampling algorithm, which must iteratively analyse the current predictions and select regions for additional sampling.

Lyu et al. (2025) reported that an adaptive sampling strategy improved the accuracy of a random forest (RF)-based surrogate model for long-term slope deformation prediction. Their approach incorporated transfer learning, starting with an initial model trained on 500 synthetic samples and iteratively adding 10 additional synthetic samples per step, with parameters calibrated using field measurements. Chen et al. (2023) employed adaptive sequential sampling to develop a surrogate modelling framework combining an RF-based model for sensitivity analysis of input parameters with an ANN-based model incorporating Monte Carlo dropout as the primary surrogate. This framework was applied to system reliability analysis of slopes, back-calculation of soil parameters for a deep cement-mixed column-supported embankment and estimation of penetration rates in diamond drilling for assessing rock mass properties. Liang et al. (2023) conducted a comparative study of various surrogate modelling techniques utilising adaptive sampling to predict results of numerical models of a dam.

Given that the training data provide reliable observations of the primary system response, selecting the surrogate model’s architecture and training hyperparameters carefully is crucial for achieving both high accuracy and computational efficiency. The automated tuning methods discussed in the following sections aim to improve the aforementioned aspects by optimising the model accuracy, generalisation and computational performance. They do so by systematically evaluating various configurations and selecting the most suitable one based on evaluation metrics and predefined objectives. Tuning may involve selecting hyperparameters or the number of neurons and layers for ANNs. As hyperparameters are typically discrete and non-differentiable, gradient descent is not considered a suitable optimisation method in this context.

Bayesian optimisation (Mockus et al., 2014) frames hyperparameter tuning as a learning problem by constructing a surrogate model of the hyperparameter space, where the output represents the objective function evaluated on a validation set. The algorithm selects the next set of hyperparameters by balancing exploration and exploitation through an acquisition function, progressively converging towards configurations that maximise accuracy and generalisation (Snoek et al., 2015).

Bayesian optimisation has been increasingly adopted in geotechnical engineering for efficient calibration and tuning of neural network-based surrogate models. Cao et al. (2025) employed Bayesian optimisation to calibrate a CNN-based surrogate model for slope stability analysis, while Liu et al. (2024) applied it for tuning the hyperparameters of FNN- and LSTM-based models developed for virtual calibration chamber CPTs. Yang et al. (2026) utilised Bayesian sequences to optimise their two-step surrogate model for high-dimensional geotechnical problems. Asem and Gardoni (2021) adopted a Bayesian framework for predicting rock mass parameters used in the assessment of deformation and stability of underground structures. Buckley et al. (2023) applied Bayesian optimisation to dynamically update uncertain variables during pile driving, and Dai et al. (2025) used a Bayesian approach for slope back-analysis based on monitoring data.

Hyperband is a tuning algorithm that accelerates hyperparameter optimisation by dynamically allocating training resources and implementing early-stopping strategies for poorly performing trials (Li et al., 2018). It systematically explores numerous configurations by initially training multiple models for a limited number of epochs, then selectively continues only the most promising candidates for additional training. This bandit-based approach efficiently addresses hyperparameter tuning as a resource allocation problem, substantially outperforming traditional methods such as Bayesian optimisation in terms of speed. Particularly beneficial for neural network surrogates, Hyperband quickly identifies optimal network architectures without extensively training suboptimal configurations, significantly reducing computational efforts and accelerating the discovery of high-performing models.

The Hyperband tuning algorithm has found growing application in geotechnical engineering for the optimisation of neural network hyperparameters. Fan et al. (2025) applied Hyperband to tune a neural network model developed for rockburst prediction, resulting in enhanced accuracy and training efficiency. Ruiz López et al. (2024) employed the algorithm for hyperparameter tuning of an FNN-based surrogate model predicting ground displacements caused by shaft excavation. Cerek et al. (2026a) optimised both a hybrid LSTM–FNN model using Hyperband to predict the maximum bending moment of quay walls. Furthermore, Thapa and Ghani (2025) utilised Hyperband to optimise the architecture and learning parameters of neural networks for predicting the unconfined compressive strength of stabilised soils.

Neural architecture search (NAS) automates the design of neural network architectures (Hutter et al., 2019). NAS algorithms utilise methods such as reinforcement learning, evolutionary algorithms or gradient-based techniques to systematically update network structures, including layers, neuron counts and connectivity, to optimise performance for specific objectives. In surrogate modelling, NAS can identify innovative architectures tailored explicitly to geotechnical engineering challenges, potentially surpassing human-designed networks in accuracy and computational efficiency. While NAS can yield high-performing neural architectures, its computational intensity often restricts its application to research-oriented scenarios rather than routine engineering practice. Recent developments suggest combining NAS with physics-informed neural networks, highlighting its potential integration into advanced engineering solutions (Eslamlou and Huang, 2024). Guo et al. (2023) applied NAS for optimising the architecture of their PINN used in heat transfer analysis, while Hoseini Karani et al. (2025) employed NAS to tune a neural network designed for the efficient optimisation of wave energy converters. Deka et al. (2025) utilised the NSGA-II optimisation framework to refine an RNN architecture developed for predicting residential heating load.

Activation functions are a crucial element in neural networks by introducing non-linearity into the model, enabling it to capture complex relationships between input and output. They allow the network to decompose a complex non-linear problem into multiple interconnected neurons, each contributing to the computation of accurate output values based on the given inputs. Figure 5 presents various activation functions discussed in this section.

Activation functions can be grouped by their properties. Bounded functions (sigmoid, tanh and softmax) are suited for probability-like outputs, while the rectified linear unit (ReLU) family (ReLU, Leaky ReLU, parametric ReLU, exponential linear units [ELU] and Softplus) provides unbounded outputs for efficient gradient flow. Smooth functions like Swish improve learning in complex networks, and periodic functions (sine/sinusoidal representation network [SIREN]) are effective for oscillatory or highly non-linear patterns.

Vanishing and exploding gradients are key challenges in training deep networks. Exploding gradients, often caused by deep architectures, large initial weights or unbounded activations, can lead to numerical instability and model divergence. They can be mitigated via gradient clipping, appropriate weight initialisation or bounded activations such as tanh or Softplus. Conversely, vanishing gradients, typically caused by saturating functions such as sigmoid or tanh, can be lessened using non-saturating activations like ReLU, Leaky ReLU or Swish. Combining suitable activations with techniques such as batch normalisation further stabilises gradient flow, ensuring effective learning in both shallow and deep networks (Geron, 2019).

The logistic sigmoid function squashes outputs to the range between 0 and 1 (equation (2)), providing smooth transitions and interpretable probabilistic outputs (Dubey et al., 2022). It was common in early neural networks, particularly suited for binary classification tasks or any scenario requiring output values to represent probabilities. In surrogate modelling, sigmoid functions are specifically advantageous when outputs must be constrained between 0 and 1, such as in normalised quantities or failure probability predictions. Moreover, in LSTM networks, the sigmoid activation is usually used within the gating mechanism to regulate information flow and act as soft-switches for including, retaining or discarding data. However, in deep networks, sigmoid activation can suffer from saturation, leading to vanishing gradients for large positive or negative inputs, consequently slowing training in deeper layers. Therefore, the sigmoid function remains effective in shallow networks or specialised output layers requiring probability-like interpretations. The sigmoid function is infinitely differentiable, as it is based on the exponential function, and provides smooth gradients, making it suitable for use in PINNs for computing PDE residuals.

Khanlari et al. (2012) conducted a comparative study of FNN-based models using sigmoid and tanh activation functions for predicting the shear strength parameters of soils. Their results indicated that, for the investigated cases, the sigmoid function outperformed tanh for neural shallow network architectures 3-3-1 and 4-4-1. Kovačević et al. (2021) employed the sigmoid function in a neural network with 34 distinct weightings to predict tunnel displacements and loads. In the context of LSTM networks, Gao et al. (2021) used the sigmoid activation within their LSTM-based surrogate model for predicting tunnel boring machine penetration rates, while Tao et al. (2024) applied the sigmoid in an LSTM-based surrogate model of braced excavations:

The hyperbolic tangent (tanh) activation function transforms input values into an output range between −1 and 1, producing a zero-centred output beneficial for training efficiency and convergence by effectively centring data. Despite its advantages, tanh is susceptible to the vanishing gradient problem. For significantly positive or negative inputs, the function saturates, leading to near-zero gradients. This saturation can hinder training in deeper neural networks by slowing or halting weight updates during back propagation.

In geotechnical engineering, early neural network applications often employed the tanh function in hidden layers due to its ability to effectively model complex, non-linear behaviours like soil-stress relationships and structural damage patterns. The zero-centred characteristic of tanh facilitated more efficient convergence, making it a preferred choice in these contexts (Eslamlou and Huang, 2024). Nonetheless, its vulnerability to vanishing gradients remains an important consideration during network design. Nguyen-Minh et al. (2025) applied the tanh activation function to predict the seismic uplift capacity of strip anchors. In the context of PINNs, tanh provides an infinite number of derivatives and smooth gradients, making it suitable for computing PDE residuals. In LSTMs, tanh is typically used for controlling cell non-linearity, as reported by Tao et al. (2024) and Gao et al. (2021).

The ReLU activation function is characterised by outputting zero for negative inputs and linearly increasing for positive inputs, as presented in equation (3) (Parhi and Nowak, 2020). Its simplicity and effectiveness have made it a prevalent choice in modern deep networks for regression tasks. ReLU mitigates the vanishing gradient problem by maintaining a constant gradient for positive inputs, enabling faster training and improved initialisation in deeper networks. It promotes sparse activation patterns, as neurons receiving negative inputs output zero, indirectly contributing to network regularisation. Controlled gradient magnitudes in ReLU, particularly when combined with appropriate weight initialisation or complementary techniques such as batch normalisation, help prevent excessively large weight updates, thereby reducing the risk of overfitting, for example, overly network fitting noise in the data, and encouraging better generalisation. ReLU is particularly suitable for surrogate models with deep or multi-layered architectures. However, neurons can become “dead” if consistently exposed to negative inputs, outputting zero continually. To address this issue, leaky ReLU (Maas et al., 2013) is another type of activation function where the negative section allows a small gradient instead of being completely zero. This adaptation helps to reduce the risk of overfitting in neural networks. There are also parametric ReLU (He et al., 2015), ELU (Clevert et al., 2015) and others that are less commonly specifically mentioned in surrogate modelling literature, however, can be used for minor improvements to mitigate possible convergence issues. Because the first derivative of the ReLU function is piecewise constant and its second derivative is zero, its applicability in PINNs is severely limited.

In surrogate modelling for geotechnical engineering, ReLU has been applied by Gao et al. (2024) for predicting slope stability, while Wang et al. (2023) utilised ReLU for a surrogate model of embankment settlement. ReLU has also been employed in surrogate modelling of rectangular tunnels (Nguyen et al., 2025), geotechnical analysis of underground construction (Salam et al., 2025) and optimisation of in situ site investigations (Zuada Coelho et al., 2025). Xu et al. (2026) applied parametric ReLU for surrogate modelling of strip footing bearing capacity and slope stability, improving training efficiency and mitigating the dying neuron problem. Wang et al. (2025) used leaky ReLU to address dead neurons in geotechnical parameter prediction based on sparse data sets:

The softmax function transforms numerical inputs into a probability distribution where all values are between 0 and 1 and sum to 1. It calculates this by taking the exponential of each input value and dividing by the sum of all these exponentials, as shown in equation (4). In surrogate modelling, softmax is valuable when creating probabilistic approximations of complex systems and classification tasks. Softmax helps quantify uncertainty in predictions by converting model outputs into probability distributions over possible outcomes. The softmax function enables these models to express relative confidence in different predictions, allowing to understand not just what the model predicts, but how certain it is about those predictions (Wolfe et al., 2017). Softmax has been applied in multi-class CNN-based classification of unstructured geotechnical reports to generate confidence scores (Park et al., 2025). Qian and Shi (2025) employed softmax within the attention mechanism of their large language model for automated geotechnical site planning and geological characterisation. Chen et al. (2024) incorporated the softmax function into their model predicting spatially varying soil properties. The application of softmax in PINNs is limited, as its outputs sum to 1, whereas PDE solutions typically involve continuous numerical values:

The softplus activation function transforms inputs into positive outputs through a smooth, continuous logarithmic–exponential relationship (equation (5)). As illustrated in Figure 6, it acts as a differentiable counterpart to ReLU, maintaining non-zero gradients across the entire input range and thereby lessening the vanishing gradient issues often encountered with sigmoid activations. A major advantage of softplus, particularly within PINNs, is its analytical smoothness and the existence of higher-order derivatives, which facilitate accurate enforcement of differential equation constraints during training. In contrast to the sigmoid function, softplus produces unbounded positive outputs, making it appropriate for representing quantities that cannot take negative values. Nevertheless, the absence of an upper bound may limit its usefulness in applications requiring bounded or normalised outputs, such as probability-based predictions (Dashtbayaz et al., 2024). The application of the Softplus activation function in surrogate modelling of geotechnical structures has so far been limited to a few studies. Zhang et al. (2024) employed the Softplus function in an FNN-based model for predicting the true triaxial strength of sandstone, where it outperformed other activation functions. Their model used compressive strength, intermediate stress and minimum principal stress as inputs to predict the maximum principal stress at failure, based on 417 data sets from true triaxial tests on sandstone. Cerek et al. (2026b) conducted a comparative study between three FNN-based surrogate models and a PINN-based model for predicting cofferdam deformation. In their work, Softplus was incorporated as a differentiable activation function in both the FNN and PINN models, demonstrating the benefits of integrating physical equations into surrogate modelling frameworks:

The Swish activation function provides a smooth, non-monotonic mapping that retains non-zero outputs across the entire input domain, including negative values, which can contribute useful information during learning (Ramachandran et al., 2017). As presented in equation (6), swish is formulated as a sigmoid-based function with a trainable parameter $β$ that controls the curvature, ranging from nearly linear behaviour to ReLU-like characteristics at higher $β$ values. Its continuous differentiability ensures stable gradient flow and mitigates vanishing gradient problems common in deep architectures, while the ability to output negative values enhances the representation of physical processes fluctuating around zero. These properties make Swish particularly promising for PINNs, where accurately resolving fine-scale variations in PDE solutions is essential. Nevertheless, Swish is computationally more demanding than simpler functions such as sigmoid or tanh, and small gradients for negative inputs may lead to slower convergence. In geotechnical applications, Swish has been used in FNN-based surrogate modelling of true triaxial tests on sandstone (Zhang et al., 2024) and in neural network models for assessing underground excavation stability in mining engineering, trained on 175 historical case studies (Szmigiel et al., 2024). Butrat and Thongpong (2025) conducted a comparative study of several activation functions for generating normal lines for compaction curves leveraging an FNN-based model, where swish emerged as more effective than other functions like ReLU, sigmoid and tanh:

Periodic activation functions, such as sine and cosine, provide smooth, differentiable and periodic mappings capable of representing highly non-linear or oscillatory behaviours (Parascandolo et al., 2017). The SIREN utilises a sine-based activation with a frequency scaler $ω0$ (equation (7)) to capture fine-scale, wiggly patterns with high fidelity (Sitzmann et al., 2020). By producing both positive and negative outputs and maintaining continuous gradients, these functions are particularly effective for modelling repeating or wave-like phenomena, such as soil–structure interactions under cyclic loading or seismic waves in geotechnical problems. The derivatives of sine functions are cosine, making them attractive for PINNs, when solving boundary value problems. However, their non-monotonicity and inherent periodicity can create local minima and slow convergence. In addition, careful input scaling and frequency selection are essential, as large input magnitudes can lead to highly oscillatory outputs that impede learning stability, and weight initialisation schemes must be chosen carefully. To date, applications of sine-based activations in surrogate modelling in geotechnical engineering have been limited, with Ma et al. (2025) comparing FNN models for predicting gas permeability of sandstone, granite and bentonite:

Evaluation metrics provide an objective assessment of model performance for the user. Different metrics vary in their sensitivity to specific types of errors and are better suited for particular evaluation scenarios. Because these metrics often determine whether a model has achieved the desired level of accuracy, careful selection of appropriate evaluation criteria is essential for reliable model assessment and in-depth understanding of model behaviour.

The coefficient of determination R² is an essential scale-free metric for assessing the holistic performance of regression models, including surrogate models. As shown in equation (8), it measures the proportion of variance in the dependent variable explained by the model, where $yi$ is the actual value, $yi$ is the predicted value and $y-$ is the mean value. R² values range from 0 to 1, where a value of 1 represents a perfect prediction, while 0 indicates that the model is not better than simply predicting the mean value. In the context of neural networks, R² is widely used to evaluate how effectively the model captures underlying data patterns (Chicco et al., 2021):

Root mean square error (RMSE) is a widely used performance metric that quantifies the aggregate error between actual values $yi$ and predicted values $y^$ obtained from a model (Willmott, 1982). By applying the square root to the mean of squared errors (MSE), RMSE provides an interpretable measure expressed in the same units as the target variable, as shown in equation (9). However, RMSE is particularly sensitive to large errors, which can be problematic when values are not normalised:

Mean absolute percentage error (MAPE) is a widely used metric for evaluating the performance of regression models in a scale-free and interpretable manner. As shown in equation (10), the error is measured relatively, making MAPE particularly useful for comparing model performance across data sets or scales where mean absolute error may be less informative due to its dependence on the magnitude of values (de Myttenaere et al., 2016). However, MAPE is highly sensitive to very small actual values $yi$⁠, where minor prediction errors can result in disproportionately large percentage errors. Analysing the relative distribution of percentage errors across the predicted test data provides additional insight into the model’s performance. Ideally, a power-law distribution would be observed, indicating accurate model predictions, as shown in Figure 7:

This study presents a state-of-the-art review of neural network-based surrogate modelling in geotechnical engineering, providing an efficient alternative to conventional modelling techniques by approximating input–output relationships of a primary system using a secondary model trained on available data. Each step in developing a surrogate model that replicates soil behaviour and soil–structure interactions offers a range of possibilities for geotechnical engineers. Careful examination of the data sets intended for training strongly influences the selection of an appropriate neural network architecture. FNN-based models have proven efficient for processing and predicting complex non-linear data sets. However, they are limited to static input–output mappings. When the governing physical laws of a geotechnical problem are known in a differentiable form and can be incorporated into the loss function, PINNs can improve the accuracy and robustness of surrogate models by constraining solutions within a physically feasible domain. Applications of PINNs based on FNN architectures exist in the geotechnical context and have been reported in this review. Sequential data exhibiting temporal or spatial dependencies can be effectively handled by LSTMs, although their gating mechanisms may gradually forget earlier information in very long-term problems. In contrast, transformer models employ attention mechanisms and parallel processing of the entire input sequence, making them a promising architecture for complex time- and spatially-dependent surrogate modelling tasks.

Data sampling, alongside the exploration of the variable space and potential dimensionality reduction, is essential for the efficient development of surrogate models. Conventional sampling techniques, such as grid search, are generally unsuitable for complex, high-dimensional systems. LHS has therefore emerged as a more practical and easily implementable alternative, stratifying the input domain and sampling each variable once within its respective stratum. Quasi-random deterministic methods provide more uniform coverage of the variable space, but their implementation can be challenging. Adaptive sampling offers an effective means of covering the input domain and refining samples in regions of high uncertainty. However, its iterative coupling with the surrogate model makes it more complex to implement, restricting its use to cases where sampling is not performed entirely in advance.

The architecture of a surrogate model requires careful tuning to achieve optimal accuracy for a given data set. Hyperparameter optimisers, such as Hyperband and Bayesian optimisation, adjust the hyperparameters of a predefined architecture to identify the configuration that minimises the loss function. NAS represents a more advanced approach, generating multiple candidate architectures and optimising their hyperparameters to achieve the highest possible predictive performance. All of these methods have been successfully applied within the field of geotechnical engineering.

The choice of activation function is a critical factor in designing neural networks, as it directly affects model accuracy, training stability and convergence. Activation functions should reflect the characteristics of the output domain: for instance, classification tasks often require bounded outputs representing probabilities, making functions such as sigmoid and softmax appropriate. In contrast, regression tasks benefit from activations with unbounded or approximately linear behaviour for positive values, such as ReLU, Swish and Softplus, which facilitate gradient flow and efficient learning. Ultimately, selecting an activation function that aligns with the data characteristics and modelling objectives is essential for achieving robust and reliable surrogate models in geotechnical engineering.

The final evaluation of model performance relies on a range of evaluation metrics, each with its own sensitivity to specific types of errors, underscoring the importance of selecting appropriate metrics to avoid exaggerating the impact of very small or large errors for comprehensive assessment.

With widely used techniques now gaining broad adoption in geotechnical engineering as well, surrogate modelling enhances geotechnical workflows, supporting sensitivity analysis, reliability assessment, design and stability evaluation. However, rigorous quality assurance is essential at each stage to ensure model reliability and robustness.

Despite notable advances in neural network-based surrogate modelling for geotechnical engineering, several areas remain underexplored. Systematic benchmarking of various activation functions, including for instance Swish, Softplus, but also periodic functions like sine, across different architectures and problem types is needed to guide their optimal selection for complex, high-dimensional or time- and spatially-dependent geotechnical data sets. Establishing standardised benchmarks and more reliable databases of field and experimental data (ISSMGE TC304, 2025; Plante St-Cyr et al., 2025; Tang and Phoon, 2024) will provide ground-truth solutions, improve model training and enable meaningful comparisons between surrogate approaches.

Integration of hybrid architectures and physics-informed frameworks also warrants further investigation. While PINNs have shown promise by embedding governing equations into the loss function, current implementations mainly rely on fully connected FNNs. Exploring advanced architectures, such as LSTMs or transformers could improve modelling of sequential or spatial–temporal phenomena, including cyclic loading, slope deformation and long-term consolidation. Transfer learning strategies that combine numerical simulations and field measurements could further enhance predictive accuracy and reduce the data requirements for new applications (Hu et al., 2024).

Additional research inspiration may be drawn from related fields. For instance, a physics-informed attention-augmented LSTM has been developed to solve multi-degree-of-freedom systems, demonstrating its effectiveness for SHM of steel structures (Mahar et al., 2025). Other architectures, such as temporal convolutional networks and graph neural networks, can effectively capture spatial dependencies and enhance the surrogate modelling process in structural dynamics compared with conventional ANNs and LSTMs (Wen et al., 2023).

Finally, moving towards fully data-driven approaches, where neural networks directly map field observations to structural behaviour represents a promising direction. Coupled with adaptive sampling, dimensionality reduction and uncertainty quantification, such approaches could accelerate model development and increase generalisability. Overall, future research should focus on combining robust data sets, benchmarked activations and advanced architectures to create more reliable, efficient and physically consistent surrogate models for geotechnical engineering applications.

The authors gratefully acknowledge Kishore Yogaraj for their assistance with the initial literature research that supported the early development of this work.