The Hoek–Brown (HB) failure criterion is a well-recognized rock failure model accounting for the non-linearity of the minor principal (compressive) stress. The form of a power–law relationship between the major and minor principal stresses (i.e., σ ₁ and σ ₃ ) is the mathematical expression of the HB failure criterion. Taking tensile normal stresses as positive, the HB failure criterion can be expressed as (Hoek et al., 2002) σ 3 = σ 1 − ( − m b σ 1 ( − σ c i ) ( 1 − a ) / a + s ( − σ c i ) 1 / a ) a , (1) where σ _ci is the uniaxial compressive strength of intact rock mass, and the parameters m _b , s, and a are expressed in Equations 2–4. a = 1 2   + 1 6 { exp ( − G S I 15 ) − exp ( − 20 3 ) } , (2) m b =   m i ⁡ exp ( G S I − 100 28 − 14 D F ) , (3) s =   exp ( G S I − 100 9 − 3 D F ) . (4)

In the previous equations, the geological strength index (GSI) has typical values from 10 to 100 (extremely poor rock mass to a perfectly intact rock mass). Also, DF represents the degree of disturbance, and it has typical values from 0 (undisturbed in-situ rock masses) to 1 (extremely disturbed in situ rock masses). Parameter m _i is a material constant that is related to the frictional strength of an intact rock mass and has typical values from 5 to 30.

Figure 1 shows the problem definition of 2D unlined rectangular tunnels in a plane strain condition. Due to the assumption of the plane strain condition, the problem represents a very long unlined rectangular tunnel. The tunnels have a width (D), a length (B), and a cover depth (C) above their crown. The parameters of rock masses based on the Hoek–Brown failure criterion consist of GSI, DF, m _i , σ _ci , and unit weight, γ. A uniform surcharge pressure (σ _s ) is applied over the rock surface as shown in Figure 1.

According to the aforementioned parameters, there are eight design parameters in this study (i.e., C, D, B, σ _ci , GSI, DF, m _i , and γ). The dimensionless output parameter of this problem is the stability factor denoted by σ _s /σ _ci . Eq. 5 represents the stability factor as a function of six dimensionless parameters as follows: σ s σ c i = f ( C D , B D , γ D σ c i , m i , G S I , D F ) , (5) where B/D represents the width ratio, C/D represents the cover depth ratio, DF represents the degree of disturbance, GSI represents the geological strength index, mi represents a material constant related to the frictional strength, and γD/σci represents the normalized uniaxial compressive strength.

In this study, a non-linear input–output mapping of the system of tunnels in rock masses is constructed using a neural network trained by an extreme learning algorithm. The training data are the FELA solutions of the tunnel stability factor.

The finite element limit analysis (FELA), which is the computational method based on a perfectly plastic material with an associated flow rule, employs the plastic bound theorems, finite element discretization, and mathematical optimization (Sloan, 2013; Keawsawasvong and Ukritchon, 2017; Ukritchon and Keawsawasvong, 2017; Krishnan et al., 2019; Ukritchon et al., 2019; Ukritchon and Keawsawasvong, 2020a; Ukritchon and Keawsawasvong, 2020b; Ukritchon et al., 2020; Keawsawasvong and Ukritchon, 2021). This FELA technique is carried out to derive the bracket of the true limit load from the targeted upper Bound (UB) and lower Bound (LB) solutions. A new computer software, namely, OptumG2 (OptumCE, 2020), is employed to compute the active collapse pressure (σ _s ) of unlined rectangular tunnels in rock masses. The results will be used as datasets for constructing a machine learning model.

Figures 2A–C show three numerical models of unlined rectangular tunnels in rock masses for the cases of B/D = 0.5, 1, and 2, respectively, where the others are C/D = 2, γD/σ _ci = 0, GSI = 80, DF = 0, and m _i = 20. Due to the symmetry of the problem, only half of the domain is used in the modeling. The boundary conditions at the left plane of symmetry and the right plane are set to move only in the vertical direction. The boundary condition at the bottom plane is not allowed to move in both vertical and horizontal directions. The sizes of the domain are chosen to be sufficiently large in order to avoid any error from the problem size on the computed bound solutions. The tunnel is unlined, and there is no pressure applied on the periphery of the tunnel. At the rock surface, the uniform surcharge σ _s is applied overall the area. This surcharge at the active collapse state will be optimized and used as the output from the LB and UB FELA using OptumG2.

The recent adaptivity meshing technique is also employed in this study to improve computational efficiency (e.g., Ciria et al., 2008). Using this adaptivity meshing technique, a number of elements are automatically added in the zones that contain large plastic shear strain. As a result, the differences between UB and LB solutions become smaller after a few iteration steps. The setting of the adaptivity meshing technique in this study is set as an initial mesh of 5,000 elements at the first step and then will be increased to 10,000 elements after three iterations of mesh adaptivity. Examples of typical adaptive meshes can be seen in Figures 3A–C for the cases of B/D = 0.5, 1, and 2, respectively. Examples of absolute velocity contours of an unlined rectangular tunnel in the rock mass are also shown in Figures 4A–C for the cases of B/D = 0.5, 1, and 2, respectively, which depict the associated failure mechanisms of this tunnel stability problem. All results of the tunnel stability obtained from the LB and UB FELA (about 2,160 solutions) will be averaged and then used as the training datasets in the machine learning approach as described later in the next section. The results of the stability index obtained by the FELA will be used in the output layer datasets in order to construct a machine learning model.

Table 1 concludes the ranges of dimensionless parameters considered in the FELA. It is important to note that these ranges cover realistic geometrical and geological properties of the tunnel and rock masses. The input includes major parameters: γD/σ _ci , GSI, m _i , DF, C/D, and B/D affecting the stability index that is analyzed using the FELA and applied as an output variable in machine learning models.

A multilinear regression model is first conducted, and the regression coefficients or weights are obtained. The coefficients are optimized by minimizing the error in WEKA software. The multiple linear regression equation is shown as follows: y = − 8.9111 x 1 + 0.2889 x 2 + 0.4117 x 3 − 1.7415 x 4 + 2.4081 x 5 − 3.4669 x 6 − 25.4051 , (10) where y represents the stability factor σ _s /σ _ci , whereas x 1 , x 2 , x 3 , x 4 , x 5 ,   a n d   x 6 are the dimensionless input parameters, namely, γD/σ _ci , GSI, m _i , DF, C/D, and B/D, respectively.

Figure 6 presents the comparison between the actual stability factor (average values from UB and LB solutions) obtained from the FELA and the predicted value obtained from the multiple linear regression. The performance of the developed equation can be accessed via statistical tests, R ², MAE, and RMSE which are found to be 0.8466, 3.2983, and 4.5924, respectively (see Table 2).

In order to maximize the accuracy of the ANN models, the number of hidden layers and neurons should be optimized. In this study, one hidden layer is considered, while the number of hidden neurons is varied. Figure 7 presents the performance of ANN models for the stability factor of rectangular tunnels in rock masses. It is clear that the performance of ANN models is increased if the number of hidden neurons is increased. It should be noted that R ², MAE, and RMSE must be calculated as the R ² value alone cannot really show the performance of the models as the variation of MAE and RMSE can still be seen clearly. Figure 8 presents the performance of the models against the number of hidden neurons. It is found that when the number of hidden neurons is over 11, the performance of ANN models is likely to be stabilized. In this case, ANN with the architecture of 6-11-1 is chosen to be the optimal MLP model as it shows the lowest MAE and RMSE values among the other models, while R ² is the highest among the models. Table 2 also compares the performance between the MLR and MLP models. It is clear that the MLP model performs much better than the MLR model. This optimal MLP model with the architecture of 6-11-1 is used in the next section for the sensitivity analysis.

After obtaining the optimal ANN architecture, the approximate general functions can be employed considering the weighted inputs and the transfer function to create the outputs. The layer number defines the superscript on the weight matrix in multiple-layer networks, as seen in Figure 9. The proper notation is used in the two-layer tansig/purelin network. This network is useful for approximating general functions. Given a sufficient number of neurons in the hidden layer, it can arbitrarily approximate well any function with a finite number of discontinuities. In this section, the final weights of each parameter have been calculated in order to study the effects of each parameter on the stability index. Figure 9 shows an example of the dimension of weight matrix and bias of the optimal ANN model for the heading tunnel. Predictive Equation 11 can be developed based on the tansig function, weight, and bias from the ANN model. P r e d i c t e d   v a l u e =   [ ∑ i = 1 N I W 2 i t a n s i g ( ∑ j = 1 J I W 1 i j x j + b 1 i ) + b 2 ] , (11) where N is the number of hidden neurons, X is the number of input variables, and J is the number of input variables. The weight matrix (IW1 and IW2) and bias (b_1i and b₂) in the hidden and output layers corresponding to the optimal ANN models are obtained. Hidden weight (IW1) is obtained based on the number of input parameters (J) and hidden neurons (N). There is one weight for every input to neuron connection between the layers. Each neuron in the hidden layer has its own bias constant (b_1i). As for the output weight matrix (IW2), the number of rows matches the number of hidden layer neurons (N) and the number of columns matches the number of output layer neurons (k). There is one column for every neuron in the output layer. In this case, the output layer contains only one column. Table 3 presents the neural network constants of the optimal ANN model including the weight matrix and bias for the stability index calculation of rectangular tunnels. These values obtained from the optimal ANN networks can be used to develop predictive equation functions and test on new datasets with different variations of parameters within required ranges.

The sensitivity analysis of the stability of rectangular tunnels in rock masses using the ANN approach is presented next to portray the influences of all considered input dimensionless parameters (e.g., B/D, C/D, DF, GSI, m _i , and γD/σ _ci ) on the stability factor (σ _s /σ _ci ). Figures 10A,B, respectively, show the effect of the width ratio B/D on the stability factor σ _s /σ _ci for the cases of (γD/σ _ci = 0, DF = 0, GSI = 80, and m _i = 5 and 30). The non-linear relationship between σ _s /σ _ci and B/D can be seen in Figure 10. It is found that an increase in B/D can decrease the geometrical arching effect resulting in the reduction of tunnel stability. As a result, a larger B/D ratio causes a smaller value of σ _s /σ _ci . In Figure 10, the stability factor of the lowest case (B/D = 0.5) is higher than that of the highest case (B/D = 2), which is about 120–300%. The impact of cover depth ratio C/D on the stability factor σ _s /σ _ci for the cases of γD/σ _ci = 0, DF = 0, m _i = 20, and GSI = 60 and 100 are shown in Figures 11A,B, respectively. A non-linear relationship between σ _s /σ _ci and C/D is observed in Figure 11. A larger C/D value also yields an incensement of the geometrical arching effect which is positive in improving the tunnel stability factor (σ _s /σ _ci ). Generally, when B/D = 0.5 to 1, the case of the deepest cover depth (C/D = 5) has a larger stability factor about three to five times of the case of C/D = 1, as shown in Figure 11. However, when B/D = 2, the difference between the stability for the cases of C/D = 1 and 5 becomes very large about 40 times.

The influences of Hoek–Brown parameters for rock masses including DF, GSI, m _i , and γD/σ _ci are presented next. From Figures 12A,B, the impact of geological strength index GSI on the stability factor σ _s /σ _ci is illustrated for the cases of γD/σ _ci = 0, C/D = 3, B/D = 0.75, and DF = 0.25 and 0.5, respectively. Numerical results in Figure 12 have shown that there is an exponential relationship between GSI and σ _s /σ _ci , where an increase of GSI results in a non-linear increase of σ _s /σ _ci for all m _i values. Such results can be referred to the exponential function used in the model of the Hoek–Brown failure criterion as expressed in Eqs 2–4. It is due to the fact that a large GSI value is a highly undisturbed rock mass which yields an increase in the stability of rock tunnels. From Figure 12, it is also found that when the value of m _i is large, the non-linearity of the σ _s /σ _ci and GSI relationship is also high. The effect of the m _i parameter on the stability factor σ _s /σ _ci is presented in Figures 13A,B for the cases of DF = 0.25 and 0.5, respectively. The plots are for γD/σ _ci = 0, C/D = 3, and B/D = 0.75. The results in Figure 13 show that the relationship between σ _s /σ _ci and m _i is linearly increasing. An increase of m _i results in an increase of σ _s /σ _ci for all B/D values. In a physical meaning, m _i depends upon the mineralogy, composition, and grain size of the intact rock. From Figure 13, it can be observed that the slope of the σ _s /σ _ci and m _i line becomes higher when the value of GSI becomes larger, meaning that the impact of m _i is more prominent when the GSI value is high. The influence of the degree of disturbance DF on the stability factor σ _s /σ _ci is demonstrated in Figures 14A,B for the cases of GSI = 60 and 100, respectively, where the others are γD/σ _ci = 0, C/D = 3, and m _i = 20. A high value of DF represents a larger degree of disturbance meaning that the stability factor σ _s /σ _ci is reduced. The slopes of the lines of DF effect decrease as the values of GSI increase (see Figures 14A,B). For the case of GSI = 60 (see Figure 14A), the structure of rock masses for this case is very blocky and has partially distributed masses so that the impact of DF is very high due to the low quality of rocks. On the other hand, when GSI = 100 (see Figure 14B), the structure of rock masses is perfectly intact with few very widely spaced discontinuities. As a result, the effect of DF on this perfectly intact rock on the stability factor is then smaller than that of blocky rocks. Finally, the influence of the normalized unit weight and the uniaxial compressive strength ratio γD/σ _ci is presented in Figures 15A,B for the cases of DF = 0, GSI = 60, C/D = 3, and m _i = 5 and 30, respectively. Obviously, all data are plotted horizontally, meaning that the increase of σ _ci /γD does not significantly affect the results of σ _s /σ _ci . This is because the stress induced by the unit weight γ of rock masses is low in comparison to that of σ _ci for rock masses. Thus, the impact of γD/σ _ci is negligible for the tunnel stability in rock masses.