Predictive Models to Evaluate the Interaction Effect of Soil-Tunnel Interaction Parameters on Surface and Subsurface Settlement

Nowadays, the need for subway tunnels has increased considerably with urbanization and population growth in order to facilitate movements. In urban areas, subway tunnels are excavated in shallow depths under densely populated areas and soft ground. Its associated hazards include poor ground conditions and surface settlement induced by tunneling. Various sophisticated variables influence the settlement of the ground surface caused by tunneling. The shield machine's operational parameters are critical due to the complexity of shield-soil interactions, tunnel geometry, and local geological parameters. Since all elements appear to have some effect on tunneling-induced settlement, none stand out as particularly significant; it might be challenging to identify the most important ones. This paper presents a new model of an artificial neural network (ANN) based on the partial dependency approach (PDA) to optimize the lack of explainability of ANN models and evaluate the sensitivity of the model response to tunneling parameters for the prediction of ground surface and subsurface settlement. For this purpose, 239 and 104 points for monitoring surface and subsurface settlement, respectively, were obtained from line Y, the west bond of Crossrail tunnels in London. The parameters of the ground surface, the trough, and the tunnel boring machine (TBM) were used to categorize the 12 potential input parameters that could impact the maximum settlement induced by tunneling. An ANN model and a standard statistical model of multiple linear regression (MLR) were also used to show the capabilities of the ANN model based on PDA in displaying the parameter's interaction impact. Performance indicators such as the correlation coefficient (R 2 ), root mean square error (RMSE), and t-test were generated to measure the prediction performance of the described models. According to the results, geotechnical engineers in general practice should attend closely to index properties to reduce the geotechnical risks related to tunneling-induced ground settlement. The results revealed that the interaction of two parameters that have different effects on the target parameter could change the overall impact of the entire model. Remarkably, the interaction between tunneling parameters was observed more precisely in the subsurface zone than in the surface zone. The comparison results also indicated that the proposed PDA-ANN model is more reliable than the ANN and MLR models in presenting the parameter interaction impact. It can be further applied to establish multivariate models that consider multiple parameters in a single model, better capturing the correlation among different parameters, leading to more realistic demand and reliable ground settlement assessments. This study will benefit underground excavation projects; the experts could make recommendations on the criteria for settlement control and controlling the tunneling parameters based on predicted results.


Introduction
Rapid urbanization has become a significant issue, particularly in densely populated areas with limited land. Many tunnel excavation projects are taking place concurrently around the world. Tunnel construction is becoming more common in urban areas, which poses a risk of damage to existing structures due to tunnel-induced ground deformation [1][2][3]. Recent advances in tunneling technology have allowed it to take advantage of the underground space to a greater extent. Various techniques and machines are now available to efficiently excavate under almost any soil condition. As the zone of influence caused by tunneling-induced ground deformations affects nearby structures, foundations, and services, some tunnel-soil structure interaction is unavoidable. Although such conditions are inherently undesirable, tunnel construction in such areas may be compelled by geographical and economic constraints.
Tunneling is affected by various unknowns that influence the cost and time of an underground project, such as geology, geometry, and the construction process. For example, the geological conditions at the tunnel level were primarily unknown before construction [1]. Therefore, selecting soil parameters is crucial to analyzing tunnelinginduced ground deformation. Therefore, evaluating the settlement hazards induced by tunneling is recommended since they have high significance for the long-and short-term performance of the tunnel. The effects of tunneling on the ground surface have been studied by several researchers [2][3][4]. However, their conclusions have been marred by uncertainty since it is unclear how the various elements affect the tunneling process. The optimum use of safety and economy in the design of underground excavations is possible by selecting the appropriate system and correctly modeling the chosen system.
The back-analysis procedure showed promising results in solving several geotechnical problems, particularly in tunnel engineering. Sakurai et al. [5] investigated a series of back-analysis procedures. They concluded that backanalysis would continue to be a key to achieving scientific tunnel construction in years to come. Recently, in further efforts to make the best use of advancing technologies, researchers enhanced the accuracy of the back-analysis procedure by integrating Machine Learning (ML) technology. ML technology is a subset of Artificial Intelligence (AI) that uses historical data as input to predict new output values. It has become a powerful tool in geotechnical engineering and has been used by many researchers [6][7][8].
Soil-tunnel interaction is affected by a wide range of factors that have interacted together in a manner that cannot present only by empirical, analytical, or numerical models. Although various empirical relationships in the form of regression equations have been derived from experimental results and are widely used, they do not provide accuracy or predictability when the interactions among many predictor variables are unknown, complex, or nonlinear. The inputs' limitations and the outputs' uncertainty are among the shortcomings of conventional methods. Furthermore, it does not clearly show the influential parameters and their interaction, which influences the decision-making process regarding which parameter has physical meaning in the prediction model.
Most studies that discussed ground deformation in the form of surface and subsurface settlement focused on the engineering properties representing most constitutive models, such as the elastic and strength properties. The soil index properties, specifically, Atterberg's limits and moisture content, have a significant effect on the engineering properties and, subsequently, the soil behavior, which is not wise to ignore. Das et al. [9] evaluated soil swelling pressure, considering moisture content, liquid limit, and plasticity index. Their sensitivity analysis found that the plasticity index is a more critical parameter that affects the swelling pressure, followed by the moisture content. Therefore, studying the interaction effect between various soil properties is crucial to estimating soil behavior accurately.
This article aims to develop a predictive model of the real-life surface and subsurface ground settlement induced by tunneling based on particular tunneling parameters to see how these parameters and their interactions affect the surface of the ground and the settlement of the subsurface. The current study used Multiple-Linear Regression (MLR) and Artificial Neural Networks (ANNWs) as a more effective alternative to regression models to demonstrate how much better some models can overcome these difficulties. Model development is combined with the partial dependency approach (PDA). The PD approach is not yet well established in geotechnical engineering. However, it has been utilized in other scientific domains to improve the sensitivity of the model response to parameter interaction effects. It is a robust visualization approach to capture the interaction between the predictor variables and provides a helpful guide to engineers. This enables the user to identify the inconsistencies and irregularities of the model's response.

Review the Applications of Artificial Neural Networks in Geotechnical Engineering
ANNWs are a subset of ML that have recently improved prediction performance by eliminating defects in other numerical, empirical, and analytical methods. A gradual improvement was achieved in the use of ANNWs in geotechnical problems to improve the accuracy and reliability of the results. Many studies have used ANNWs to predict tunneling-induced ground settlement [10,11]. Suwansawat & Einstein [12] have studied the impact of twin tunnels on the surrounding area focusing on studying the effect of the operational tunnel boring machine (TBM) parameters such as face pressure, penetration rate, pitching angle, and grouting quality on the ground movement using neural network model. However, the study did not address the geological conditions of soil behavior, although soil properties control the magnitude of settlement.
Boubou et al. [13] used ANNWs to predict ground surface movements induced by a tunnel in Toulouse in France. Their study focused on the TBM operation parameters as input in ANNW, and they concluded that ANNWs could be a helpful tool to predict the displacements induced by the TBM. However, their study did not consider the geological variability along the tunnel drive. A later study, Boubou et al. [14] proposed a methodology for predicting ground surface movements based on TBM operation parameters and geological data. They concluded that the ANNWs adapt well to ground surface movement description and correlation with TBM parameters.
Torabi et al. [15] evaluated the effect of geotechnical parameters on TBM performance in the Tehran-Shomal highway project using ANNW and Statistical Package for the Social Sciences (SPSS) software. By comparing the results for correlation coefficients (R 2 ) obtained by ANNWs and SPSS software, it was concluded that ANNWs give more reliable results than SPSS. Khatami et al. [16] applied a computational intelligence method to predict the ground surface settlement during the construction of twin-tunnel, including the surface buildings using supervised feed-forward backpropagation neural networks. Their study depended on the geometrical parameters and structure-tunnel interaction parameters such as the building's width, bending stiffness, building weight, and distance to the tunnel centreline to feed the network. Their simulation results indicated that the proposed ANNW models could accurately predict surface settlement; however, the input parameters also ignored the geological parameters. As a result, considering the soil properties in predicting the impact of tunneling on ground deformation received later significant attention from many researchers [10,17,18]. Among the very recent investigations are the work by Chen et al. [18], Khalili et al. [19], and Moghaddasi & Noorian-Bidgoli [20], who attempted to integrate the shield's operational parameters, the tunnel's geometry, and the geological conditions in the NNW models to improve the prediction results. Gradually some efforts [2,3,8] have been made to improve the accuracy of the prediction models by developing a hybrid approach of numerical modeling and NNWs.
However, each type of soil behaves differently under specific circumstances, particularly clay soil. Therefore it is crucial to fully understand the geological condition before tunneling to ensure the safety and efficiency of tunnel construction. Predicting the short-term impact of tunneling in clay is the most challenging task for geotechnical engineers and decision-makers. Most previous studies that predicted the behavior of clay soil or soft ground during tunneling ignored the index parameters, which are usually influenced by variations in moisture content due to climatic changes or fluctuations in groundwater levels [21]. According to Moghaddasi & Noorian-Bidgoli [20], and Shreyas & Dey [22] the tunneling on soft ground is associated with some hazards, including poor ground condition, the water table above the tunnel, shallow overburden, surface, and subsurface settlement. Therefore, predicting this kind of settlement induced by tunneling is challenging. However, their investigation to predict maximum surface settlement ignored index properties such as Atterberg limits and moisture content, which seriously impact the mechanical behavior of the soil, focusing only on cohesion and Young's modulus to represent soil properties with other parameters. In a similar study, Neaupane & Adhikari [23] used ANNW to predict the ground movement around tunnels with the help of input variables such as depth to the tunnel axis, normalized volume loss, soil strength, groundwater characteristics, and construction methods. However, the characteristics of groundwater may not provide an accurate indicator of soil behaviour.
Although various attempts have been made to understand soil behavior during tunneling, their prediction model s could still not capture many aspects. In particular, the clay soil behavior, which has been consistently tied to its Atterberg limit, strongly depends on its moisture content more than other properties, as mentioned by Attah & Etim [24]. On the other hand, the previous models were insufficient to describe the sensitivity of ground deformation to the complex interactions expected between tunneling parameters. Therefore, the current study makes an effort to incorporate soil index properties, such as natural moisture content (MC), plastic limit (PL), and plasticity index (PI). In addition to mechanical properties, such as shear strength and Young modulus to represents ground parameters. The data set used to develop the prediction models in this paper was obtained from a high-quality field investigation in Hyde Park within the formation of London clay (LC) for the Crossrail project [25]. The reason behind the choice of LC is that a large amount of construction has been done on the LC outcrop. So, there are a relatively large number of site investigation reports that give details of its properties. Furthermore, the availability of high-quality field data for both short-and long-term behavior presented in detail made the LC an ideal case study for tunneling-induced ground settlement.

Project Overview
The Crossrail (CR), or the Elizabeth line as it is now known, is a railway line in southeast England. It runs from Essex in the east to Berkshire in the west, tunneling underground through central London, a s shown in Figure 1. The current study considered the extension part of this project which runs from Maidenhead in the west, through and under London, to Sheffield in the northeast and Abbey Wood in the southeast and a spur line connecting Heathrow, see Figure 2. The project involves the construction of 21 km of twin-bore tunnels and nine underground stations. The Earth Pressure Balance (EPB) method has been used to construct the CR with a total diameter of 7.1m. Both pairs of tunnels are located at 34·5 m within the London Clay Formation (LCF). The instrumented site decided to be in Hyde Park lies on the alignment of the first drive, which runs from the Royal Oak portal west of Paddington Station to Farringdon Station, to investigate the ground response in the greenfield condition. The first tunnel boring machine (TBM1) was used to construct the Westbound Tunnel (WB) and passed beneath the Hyde Park Instrumented Site in late November 2012. The other boring machine (TBM2 was used to construct the Eastbound (EB) tunnel that passed underneath the site in early February 2013. The average driving speed was about 110 m per week.
Two lines of surface monitoring points (SMPs), the X-line and the Y-line, are shown in Figure 3 for the longitudinal section of the instrumentation site to measure surface displacements by precise leveling. Each line accommodates one or more instruments to monitor the ground responses induced by TBMs. Vertical settlement of the subsurface is measured using 14 rod extensometers, each having eight anchors that extend from ground level to the top of the tunnel. More details about the instrument's site, installation, and monitoring accuracy process are described in Standing and Selemetas [26] and Wan & Standing [27]. The scope of this paper will be limited to ground surface and subsurface settlement caused by the tunneling of TBM1-WB. Only the SMPs within line Y will be considered for surface settlement ( Figure 4). The subsurface settlement data was provided as recorded from rod extensometer anchors Figure 5).

Site Geology
London Clays One of the clays has received the most attention in geotechnical engineering. Several researchers De Freitas & Mannion [28], Hight et al. [29] and Gasparre [30] have studied the mechanical behavior of LC. These studies emphasized the similarity of the silty clays that comprise most of the LCF. They also concluded that LC is heavily over consolidated, has a high undrained shear strength, is highly impermeable, has high plasticity marine clay with an important silt content, and has a long stand-up time. They were characterized by anisotropic behavior, making their stiffness exhibit significant variability vertically and horizontally. Five major transgressive-regressive cycles are recognized within LCF, used to define five layers from A to E, described in detail by De Freitas & Mannion [28]. Most of the 7.1 m diameter of the Crossrail tunnel was excavated within the lower level of the LCF unit B (B2), see Figure 6.

Preparation of Database
The real-time data was collected from the CRT, London, during the TBM1-WB tunnel excavation in late November 2012. The framework adopted to develop the prediction model is illustrated in Figure 7. The proposed framework in Figure 7 is based on the proposed categories of Suwansawat and Einstein [12]. According to their classification, the features used to estimate the maximum settlement of the ground are separated into geological properties, tunnel geometry, and operation factors. The operation factors in the current study are represented by the TBM parameters, the trough parameters incorporated here instead of the geometry parameters, which can be significant factors when the database contains varying diameters [31]. Recently Ding et al. [7] investigated the ground settlement induced by tunneling using the categories of Suwansawat and Einstein [12]. Although they emphasized that the geological index can be used as one influencing factor, they neglected the geological features. The index properties (the moisture content and the plasticity of soil) were incorporated into the current model to improve the models of Suwansawat and Einstein [12] and Ding et al. [7] for predicting ground settlement in clay soil during tunnelling.

Modeling Approach
The data set for the current study contained 239 points for surface settlement and 104 points for subsurface settlement based on the systematic monitoring data for the Crossrail tunnels. The input data related to LC properties used in this paper were collected based on an extensive review of laboratory studies [32][33][34][35][36][37][38] that have been done for LC. The methodology includes the design of two models to predict the response of the ground to the tunneling parameters: one for surface settlement and the other for subsurface settlement. Multiple linear regression (MLR) and artificial neural networks (ANNWs) based on PDA were strategies used to create unique interactions between the selected variable and ground settlement.

Multiple Linear Regression (MLR) Models
As previously stated, settlement depends on more than one parameter, so the MLR technique is used to learn more about the relationship between several independent or predictor variables summarised in Table 1 and a dependent or criterion variable. The input parameters are organized under three categories, as shown in Figure 7, including 12 parameters indicated by (X1 to X12) to facilitate interpretation of the results, as shown in Table 1. Table 2 shows the maximum and minimum values of these parameters. The MLR procedures estimate a linear equation, as shown in Equation 1. SPSS, a powerful statistical software package, was used to calculate a linear regression between the independent and dependent variables where is the error rate of the model, subscript represents the regression coefficients or the independent contributions of each independent variable to the prediction of the dependent variable. The regression line expresses the best prediction of the dependent variable ( ), given the independent variables ( ).

Artificial Neural Network Models
The raw data in the NNW model typically cannot be used directly. Therefore, it must be pre-processed before being used to increase the processing and convergence rate during the training process and minimize the prediction error. Thus, in the first step, all input candidates designed were standardized from -1 to 1, following Equation (2). The scaling is carried out to improve the accuracy of subsequent numeric computation and obtain better output as the predictor variables have various units and ranges, as shown in Table 2, which could impact the final results. The advantage of scaling the data is preserving exactly all relationships in the data, and it does not introduce bias. The standardisation scaling technique involves centering the values on the mean with a unit standard deviation. The attribute's mean becomes zero, and the resulting distribution has one standard deviation.
where is the mean of the feature values and is the standard deviation of the feature value.
The NNWs create models that connect input and output. The smallest entity from which a NN is constructed is known as a perceptron. Each perceptron has the same number of weight factors as inputs, denoted as in Figure 8. Finally, each perceptron has an activation function that the user must select. Therefore, each perceptron has an activation function that the user must choose and is thus considered known, enhancing NN's performance and efficiency.

Figure 8. Single perceptron with and without bias and variables
In the multi-layer perceptron (MLP), each perceptron in a feed-forward network can be connected to every perceptron in the next layer, see Figure 9. To meet the demands of system response, NNW is supplemented by an Inputand Output-Layer. The training function should identify unknown weight factors during the training process [2]. Training an ANNW model involves modifying weights and biases until the model satisfies the stop criteria or the error converges to the initially determined minimum value [38]. After creating an ANN model, the number of hidden layers, hidden nodes, and the type of transfer function are optimized to create the optimal model presented in the following section. The sigmoid function (SIG) is highly recommended as a training function [16,38,39]. Accordingly, this study adopted SIG as a transfer function, including TANSIG and PURELIN for the hidden layers and the output, respectively, to develop a sensible model using MLP.

Develop an Optimum Architecture of the MLP Model
To create an optimal network that can be used to predict the ground surface and subsurface settlement, start by training and testing the NNWs on a subset of all data sets. The data set from the current study was divided into 70% used as a training set, and the other 30% remained for testing. The same NNW architectures were used for the two models with a change in the dependent variable (Y) in the subsurface model, representing the monitoring value of the settlement. The ANNW code and the training and testing process have been applied using the MATLAB-R2019b routine. The input layer is transmitted through one hidden layer that multiplies weight by a hyperbolic (Tangent) sigmoid function. The output presented by(Y) is transmitted through a second connection multiplied by weight by the Purlin function.
There is no generally accepted theory to determine how many hidden neurons are needed to approximate any given function in a single hidden layer [3,38,39]. Therefore, trial and error are the only way to arrive at a suitable learning rate and the optimal number of hidden nodes. Thus, 18 neural network models with different scenarios for the NNW parameters were created according to Table 3 to find the optimal network.

Table 3. Parameters of the proposed NNWs models
The parameter The proposed scenarios The Training Algorithm Levenberg Marquardt algorithm (LMA), Bayesian regularization (BR.) The number of nodes 5, 10, 15 The epochs 1500, 2000, 3000 The performance criteria, the Root Mean Square Error (RMSE) and Coefficient of Determination (R2) used as objective functions to assess the accuracy of prediction models through the calibration (training) and verification (testing) stages, as shown in Equations 3 and 4. RMSE is useful when significant errors are undesirable and R 2 measures the predictive capacity of models.
where: ′ and Represent the K th predicted and observed values of the target, respectively; N is the number of observations computed for the error where N represents the total number of observations; ̅ is the average of observations.

Prediction of Ground Settlement Using the MLR Model
To investigate the impact of the interaction of the trough, the TBM, and ground parameters on tunneling-induced surface and subsurface settlement, an MLR model, was performed using these parameters as predictor variables and ground settlement monitoring points as a response variable. Table 4 shows how well the regression model can 'fit' the data set. The coefficient of determination is written as R 2 , a multiple of R-Square. It denotes the percentage of variance in the response variable that the predictor variable can contribute. In model (1), the R-square is 0.8193, indicating that the tunneling parameters can explain the total variance of 82% of the surface settlement values, whereas the same tunneling parameters can explain 90% of the variance in the subsurface settlement values in the model (2). The adjusted R-squared is a modified version of the R-squared adjusted for the number of predictors in the model. It is always lower than the R-squared; however, it showed a slight difference from Multiple R-Square in the two models in Table 4. The adjusted R-squared can help compare the fit of different regression models with one another.
P-values and coefficients work together to indicate which relationships in the model are statistically significant and the nature of those relationships. The coefficients express the mathematical relationship between independent and dependent variables and which variables contribute to the model. The predictor variables X1, X10, X4, X7, X8, and X9 have physical meaning in the prediction model; Volume loss (X3) showed a minor degree of significance, as evidenced by the P-value in Table 5 for model (1) -surface settlement. They're statistically significant because they're smaller than 0.05, and there's enough evidence to reject the null hypothesis and favour the alternative hypothesis. This means the independent factors have a strong influence on the dependent variable. Predictor variables X1, X3, and X8 negatively influence surface settlement, while moisture content (X10) and shear strength (X7) favourably affect surface settlement. The ground properties significantly impacted the ground surface settlement, among other variables. However, the subsurface settlement -model (2) was strongly linked to the ground properties and showed slightly higher performance in Table 4. Most of the predictor variables were statistically significant, with a P-value less than the usual significance level; only X6 and X8 were statistically insignificant, as shown in Table 6.

Prediction of Ground Settlement from the ANN Model
The proposed models were trained using different training algorithms, and the minimum error controlled the stopping criteria. As shown in Table 7, the BR algorithm was more efficient than LMA. Although the LMA is often the fastest than the BP algorithm in MATLAB-R2019b, it is considered one of the Early Stopping (ES) techniques and interpreting, which is why it gives a low value of epochs when training. The model (A-1) based on 1500 epochs and five neurons for the hidden layer utilizing the LMA approach was the best. However, the BR algorithm, especially in a model (B-3), produced better results with lower RMSE. On the other hand, the network reached more epochs, making it the optimal model for training the dataset of this study. It should be noted that the BR algorithm does not require a separate validation data set from the training data set; instead, it uses all data, leading to better generalization than utilizing the ES validation set, which increases the accuracy [40][41][42]. This algorithm has been recently used by Ding et al. [7], and they concluded that the BR algorithm has advantages in both model selection and uncertainty measurement. These advantages enable the choice of appropriate input parameters from a list of designed candidates.
According to Table 7, the optimal scenario for building the architecture of the network for surface and subsurface settlement was (B-3), which gives the lowest RMSE. The model's features are summarized in Table 8 and presented in Figure 10.  The NNW model was trained using MATLAB-R2019b software and the Scikit-learn library in Python3.7.0. The NNW model showed good performance with 90% and 88% accuracy for models (1) and (2), respectively, as shown in Figures 11 and 12. Furthermore, there is no difference in the performance of the two prediction methods. However, the ANN models outperformed the MLR model in representing the interaction between factors, making the ANN models more valuable as a prediction tool, which corresponded well with the findings of Kaczmarek [43], Mohammadi et al. [17], and Sheil et al. [4]. In comparison, Chen et al. [18] predict the maximum surface settlement caused by EPB shield tunneling using three methods of artificial neural networks, the backpropagation neural network, the radial basis function neural network, and the general regression neural network (GRNN). According to their results, the model of GRNN showed the best performance with RMSE = 1.35. However, the current model showed less RMSE (0.133), indicating that using the BR algorithm improved the network performance. By improving the machine learning speed, Zhang et al. [44] overcame the shortcomings of single machine learning and other statistical techniques. To predict surface settlement caused by the The model's evaluation based on the RMSE showed that the XGBoost algorithm performed very well, with a value of 0.11. However, they concluded that the successful application of soft computing methods requires high-quality datasets and well-extracted features closely associated with the dependent responses. Applying the Bayesian Regularization technique in the current model optimized the network compared to Zhang et al. [44]. It showed performance that was roughly comparable to the XGBoost algorithm and had a higher coefficient of determination (R = 0.997) and less RMSE than in the ANN model of Zhang et al. [44], which was 0.23. Huang et al. [6] created one of the most recent reliable models to evaluate the risk assessment of tunneling in soft soils. They applied NNW, and the results were compared to those obtained using more traditional linear regression models. They concluded that using less computational effort, the suggested NNW model can generate a reliable model with capabilities comparable to conventional techniques. It is believed to improve risk management and decision-making.
However, black-box supervised learning models (e.g., neural networks, random forests, closest neighbour approaches, and support vector regression) also have some limitations regarding interpretability or visibility [45]. As a result, using the sci-kit-learn library's Partial Dependence approach can assist in resolving this issue and enhancing the variables selection. PDA is a robust tool and the most common approach for visualizing machine learning models' data. It is a model-agostic method usually used for global explainability between inputs and outputs and always shows a linear relationship with the following Equation: where , are the parameters for which the partial dependence function should be plotted, and , are the other parameters that are used in the machine learning model ̂. The parameter(s) in S are those for which we want to know the effect on the prediction. By marginalizing the other parameters, we obtain a function that depends only on the S parameters, including interactions with other parameters. The partial function ̂ it can then be estimated by taking averages from the training data as follow: where ( ) is actual parameter values from the dataset, and is the number of instances in the dataset. The partial function tells us about the average marginal effect on the prediction for a given value (s) of parameters S. More information on the computational theory of partial dependence plots (PDPs) can be found in [45].
The interaction of plasticity index (X12) and transverse trough width factor (X5) had no significant influence on surface settlement, according to the MLR model in Table 5. However, the PDP showed in Figures 13-a to 13-b that these parameters have a tangible impact on surface settlement. The dependent variable Y reaches its maximum value for the minimal values of both parameters; as X12 and X5 increase simultaneously, the Y and interaction impact decrease. As a result, the plasticity index (X12) is becoming increasingly crucial for the dependent variable; see Figure 13-c.

Figure 13. Partial dependence curves for model (1): (a) PDP for X12 (b) PDP for X5 (c) 2D interaction between X12 &X5
On the contrary, compared to the interaction effect of non-significant parameters in Figure 13, the interaction effect of different pairs of independent parameters that significantly impact the predictive model, such as X7 and X10, appears to be negligible, as shown in Figure 14-c. Furthermore, X7 and X10 positively impact the predictive model, which explains their minimal interaction effect because they have similar correlations with the dependent variable, allowing one of them to be ignored. The interaction of X12 and X1 was investigated in Figure 15-c and revealed that the lower values (≤ 0) had the most significant effect on surface settlement. The interaction effect decreases as X1 increases, and surface settlement becomes more dependent on X1. As a result, although the MLR model considered the degrading plasticity index (X12) as a nonsignificant parameter, it can significantly increase surface settlement, proving that traditional statistics and black-box methods have many limitations in interpreting and understanding the effects of predictors. This limitation was shown by Neaupane & Adhikari [10], Suwansawat & Einstein [12], and Chen et al. [18]. They applied the ANNWs to predict the surface settlements, and according to their findings, there is an interaction between the tunneling parameters, which the NNW cannot capture. Thus, the integration shown by PDPs improved the performance of NNW, enabling the current model to outperform the previous models and accurately demonstrate the level of parameter interaction. According to the subsurface settlement model, vertical settlement increases with depth as the monitoring point approaches the tunneling zone, as revealed in the Crossrail project investigation by Wan [25]. As a result, the interaction between tunneling parameters will appear more precisely in the subsurface zone than in the surface zone, explaining the high significance of most of the parameters observed in Table 6.  Although all of these parameters have the same significance level in the prediction of the subsurface settlement, the interaction impact decreases as the rate of face pressure increases, making the subsurface settlement strongly dependent on the soil's index properties. In particular, the moisture content X10 is attributed to the fluctuation of the groundwater table during tunneling processes, which affects its level. The dependent variable Y reaches its maximum value for the minimal value of X12. The model's final 2D partial dependence curve, shown in Figures 16-a to 16  Regarding the trough parameters, the volume loss (X4) and the transverse trough width (X5) showed the same degree of significance in the MLR model in Table 6. However, when studying their real effect using the PDPs in Figures 20-a to 20-c, the volume loss seems less effective in the prediction model subsurface settlement. It is noteworthy that the impact of the interaction of X9 and X7 with the trough parameters was similar; therefore, the presentation of one of them was enough. The interaction of these parameters with the mechanical properties of soil X9 and X7 indicated that the width of the transverse trough greatly affects the prediction model and eliminates the effect of the mechanical properties more than the volume loss, as presented in Figures 21 and Figures 22-a to 22-c. That can be explained as the width of the transverse trough produced by a mechanical process revealing the distribution of stresses caused by tunneling due to the relaxation of the soil due to loss of volume, as represented by the volume loss parameter (X4). As a result, the mechanical properties of the soil are less affected by the width of the transverse trough than by the loss of volume that degrades the properties of the soil. The subsurface reaches its maximum in minimum values of X7 and the maximum values of X4, see Figure 21. These observations are consistent with C. Zhao et al. [46]. They found that volume loss generates soil deformations and plays a crucial role in determining the ground settlement. Their investigation shows that the friction angle is the most prominent soil parameter affecting soil behavior. It is worth noting that the undrained shear strength ( ) is one of the crucial parameters in assessing soil behavior in any problem involving clayey soils, as noted by Persson [47]. However, the significance of one parameter can be influenced by another nonsignificant parameter. That cannot be captured numerically or by black-box methods, making the earlier models untrustworthy for decision-makers to consider. On the contrary, the index properties represented by the moisture content and the plastic limit showed a high interaction with the transverse trough width X5 and high independence from volume loss X4; however, the impact of the interaction on the prediction model was higher than the mechanical properties, as shown in Figures 23 to 25 and Figures 26-a to 26-c. The decrease in the plasticity index X12 and the increases in the moisture content X10 indicated a large transverse trough width and the anticipated higher subsurface settlement, as presented in Figures 23 and 25-c. Therefore, geotechnical engineers in general practice should attend closely to index properties to reduce the geotechnical risks related to tunneling induce ground settlement. In comparison, Singh et al. [48] used numerical modeling and regression analysis to examine the impact of geological and geotechnical factors on subsurface and surface settlement. According to their findings, mechanical shear properties may accurately forecast tunnel settlement by more than 85%, which agrees well with the results of the global sensitivity analysis by Zhao et al. [46]. However, their study neglected the index properties that showed a significant effect on the mechanical properties and the whole model. Wan (2014) [25] investigated the relationships between TBM operation variables and observed ground response during the CR project; no clear relationships were observed. However, the current NNW model based on PDA revealed that the face pressure as a TBM parameter and trough parameters (the volume loss and transverse trough width) were the most critical parameters affecting the mechanical properties of the soil (X9 and X7), and subsequently the settlement of the subsurface. In presenting the parameter's interaction influence, the comparative findings showed that the suggested PDA-ANN model is more trustworthy than the earlier models of Neaupane & Adhikari [10], Suwansawat & Einstein [12], and Chen et al. [18]. However, it is not yet well established in the field of geotechnical engineering. It can be further applied to develop multivariate models that consider multiple parameters in a single model, better capturing the correlation among different parameters, leading to more realistic demand and reliable ground settlement assessments. Thus, this study provided a promising model to estimate the sensitivity of ground deformation to the tunneling parameters.

Conclusions
The soil-tunneling interaction is a 3D interaction problem that can be affected by more than one parameter simultaneously. Therefore, this paper aims to develop a predictive model based on particular tunneling parameters to evaluate the relative importance of system parameters in accordance with their influence on model responses. The current study used Multiple-Linear Regression (MLR) and Artificial Neural Networks (ANNWs). Model development is combined with the partial dependency approach (PDA). The study adopted real-time data from London Crossrail tunnels under Hyde Park within the London Clay Formation. The data set from the current study contained 239 points for surface settlement and 104 points for subsurface settlement. The ANNW code and the training and testing process have been applied using the MATLAB-R2019b routine, and the MLR analysis was created using the SPSS and R software. The tunneling parameters are grouped as inherent soil properties (e.g., mechanical and index), operation parameters (e.g., face pressure), and trough parameters (e.g., tunnel volume loss) on the settlement at the ground surface and subsurface during the tunneling process. The most influential parameters proposed for the above groups included 12 parameters. The following remarks can be made based on the main findings of this study:  The ANN models outperformed the MLR model in representing the interaction between factors, making ANN models more valuable as prediction tools;  The NNW model performed well, with 90% and 88% accuracy for surface and subsurface settlement models, respectively. In contrast, the MLR indicates that the tunneling parameters can explain the total variance of 82% of the surface settlement values and 90% of the variance in the subsurface settlement values;  The Bayesian regularization (BR) showed high performance and enhanced the network's performance with a lower RMSE;  The Bayesian regularization (BR) showed high performance and enhanced the network's performance with lower RMSE;  The ground properties significantly impacted the ground surface settlement, among other variables. However, the settlement of the subsurface was strongly related to the properties of the ground, with slightly higher performance, and most of the predictor variables were statistically significant;  The interaction of two features that have a different effect on the target feature could change the overall impact of the entire model;  The partial dependence plots showed an excellent way to understand the models. Furthermore, they are simple to implement and apply to any model, allowing the user to understand the tendencies of their model better and expose inconsistencies and irregular model behaviour;  The interaction between tunneling parameters was observed more precisely in the subsurface zone than in the surface zone, explaining the high significance of most of the parameters in the prediction of the subsurface settlement;  The effect of soil properties appears more in subsurface settlement than in surface settlement;  Face pressure as a TBM parameter and trough parameters (the volume loss and transverse trough width) were the most critical parameters affecting the mechanical properties of the soil (X9 and X7), subsequently the settlement of the subsurface;  The results showed that the index properties (MC, PL, and PI) were significant parameters in predicting the ground surface and subsurface settlement and have a higher interaction impact on the prediction model than the mechanical properties ( , ℎ );  Mechanized Tunneling by PDA is an excellent criterion to estimate the sensitivity of ground deformation to tunneling parameters, reduce risk, and tackle unforeseen scenes that cannot be predicted with geotechnical investigations and back box techniques.

Data Availability Statement
The data presented in this study are available on request from the corresponding author.