A new model for predicting the advance rate of a Tunnel Boring Machine (TBM) in hard rock conditions

The prediction of the advance rate of a Tunnel Boring Machine (TBM) in hard rock conditions is one of the most important concerns for estimating the time and costs of a tunnel project. In this paper, in the first step, a model based on Rock Engineering Systems (RES) is proposed to predict geotechnical risks (representing media characteristics) in rock TBM tunnelling. Fifteen main parameters that influence the geotechnical hazards were used in the modelling. In establishing an interaction matrix and also a parameter rating, the views of five experts were taken into account. The Vulnerability Index (VI) (geotechnical risk levels) for 2058 datasets out of 2168 sets of data from 53 geological zones in 11 km of the Zagros long tunnel was obtained. In the second step, based on the machine operating parameters such as torque, cutter head rotation per minute, cutter normal force and media characteristics (represented by VIs), which were used as input parameters and advance rate was used as an output parameter, while using 2058 datasets, linear and non-linear multiple regression analyses were carried out. 110 datasets (out of 2168 datasets), which were not used in the modelling, were applied to evaluate the performance of regression models and other models in literature and the results were compared. The obtained results showed that the new linear model proposed with R2=0.83 and RMSE=0.12 has a better performance than the other models.


Introduction
In tunnel construction by TBM, the estimation of time and cost of a project are one of the most important parameters, which has always been a challenge between contractors and owners. Penetration and advance rates are key parameters in the performance estimation of a TBM, which results in estimating the time and cost of a project. In recent years, many researchers devoted their work to predict the performance of a TBM. There are two categories of research works while considering modelling of predicting performances of hard rock TBMs. One category is theoretical models and the other empirical models (Rostami, 2015). In principle, theoretical models are developed, by use of tests of indentation or full-scale cutting tests in a laboratory, which provide an estimation of cutting forces based on the type of cutter, geometry of cutting, spacing, and the cut penetration (Crow, 1975;Roxborough and Phillips, 1975;Ozdemir, 1978;Snowdon et al., 1982;Sanio, 1985;Sato et al., 1991;Rostami and Ozdemir, 1993; Rostami, 1997; Balci et al. 2009;Bilgin et al. 2005Bilgin et al. , 2007Bilgin et al. , 2012; Tumac et al., 2012). The main weakness of these models is that the tests used do not entirely describe the facts regarding rock mass characteristics which the disc cutters of TBM face in the field. This type of equipment for tests may not be easily accessible in all centers for research around the world (Khademi Hamidi et al. 2010). The first models were developed, using intact rock parameters such as tensile strength and Uniaxial Compressive Strength (UCS) (Graham, 1977;Farmer and Glossop, 1980 (Farrokh, 2018). Different laboratory tests such as Taber abrasion, Shore hardness, Schmidt hammer, point load index, and drilling rate index tests were applied for the prediction of PR (Howarth et al. 1986;McFeat-Smith, 1977;Ozdemir, 1977;Pang et al., 1989;Yagiz, 2002;Yagiz and Ozdemir, 2001;Fattahi and Moradi, 2017). Other prediction models have been improved in the last three decades to consider a wide range of rock properties and rock mass conditions. Palmstrom (1995) introduced a forecast model based on the Rock Mass index (RMi) and the Q TBM model was developed by Barton (1999Barton ( , 2000 The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2020, pp. 57-74, DOI: 10.17794/rgn.2020.2.6 based on the developed Q system. The NTNU model was introduced by the Norwegian University of Science and Technology (Bruland, 2012;Blindheim, 1979 Hamidi et al. (2010). In 2008, Yagiz developed a relationship between rock mass properties and PR. Hassanpour et al. (2009,2010,2011,2016) investigated TBM performance based on the Field Penetration Index (FPI). Zhao (2007, 2009)    The objective of the present piece of work is to introduce a new model to improve the accuracy of advance rate prediction in different geotechnical conditions in rock TBM tunnelling, considering both operational parameters of TBM (machine parameters) as well as media characteristics. Media characteristics are represented by the geotechnical risk levels in the modelling. The geotechnical risk levels are obtained, using an RES based model. In the new model, machine parameters and geotechnical risk levels are input parameters and the advance rate is an output parameter. The field data used for the modelling was collected from section two of the Zagros long tunnel in Iran.

Description of the site and data collection
The Zagros long tunnel with a circular cross-section, an excavation diameter of 6.73 m, a final diameter of 6 m, and a length of 52 km is a water transfer tunnel which consists of two sections. This tunnel is located in the northwest of Iran. It was planned to transfer about 70 m 3 /s of water through this tunnel. The data was collected from the second section of the tunnel, which has a length of 26 km. This section of the tunnel was excavated by a double-shield TBM. The machine specifics are shown in Table 1.
The tunnel passes through three formations. The Pabdeh (PEPd) formation includes a mix of greenish-grey limestone argillaceous and dark grey limy shale. The Gurpi (KGu) formation combines argillaceous limestone and limy shale. Furthermore, the Ilam (Ki) formation consists of brownish-grey limestone as illustrated in The database for this research work (2168 sets of data) was obtained by combining the data in the preconstruction and construction phases. During the construction phase, the geological and geomechanical data, which were obtained before construction, were exam- During construction, the operational data of TBM such as Revolutions per Minute (RPM), torque, thrust, and advance rate were recorded. The basic descriptive statistics of the sets of data are revealed in Table 2. Hudson (1992) introduced the Rock Engineering Systems (RES) approach in which interactions between rock mechanics parameters in a rock system can be obtained. This approach was applied to different rock engineering fields such as an investigation of underground excavation stability (Lu and Hudson, 1993), with rockfall risk assessment (Cancelli and Crosta, 1993), and characteristics of a rock mass to indicate natural slope stability (Mazzoccola and Hudson, 1996). (Latham and Lu, 1999) introduced an assessment system to determine the rock mass's blastability,  rock fragmentation by blasting and predicted the advance rate of face, while determination of the operation efficiency in retreat longwall mining panel was studied by (Aghababaei et al., 2019). In the approach of RES, establishing the matrix of interaction plays an important role. By generating the interaction matrix, the weighting of the parameters of the rock mass system can be evaluated. In this matrix, the principal parameters, which affect the system, are placed on the leading diagonal of the matrix and the interactions (effects of each parameter on any other parameter) are located on the off-diagonal cells. The assignment of values to off-diagonal cells is called coding of the interaction matrix. The simplest of the interaction matrix, a two parameter system, is shown in Figure 2(a). Beside that, a view of interaction matrix coding is shown in Figure 2(b). The row passing through Pi denotes the impact of Pi on all the other parameters in the system, while the column through Pi shows the impacts of other parameters, or the remaining of the system, on the Pi. In principle, there is no limit to the number of parameters that may be included in an interaction matrix.

Rock engineering system
For numerically coding the interaction matrix, different procedures such as the 0-1 binary, Expert Semi-Quantitative (ESQ) (Hudson, 1992), and the Continuous Quantitative Coding (CQC) (Lu and Latham, 1994) were applied. Among the suggested coding procedures, the most commonly used is the ESQ coding. Based on this coding procedure, the interaction intensity is presented by the values from 0 for no interaction to 4 for critical interaction as revealed in Table 3. (1) Where: C i -Cause of i th parameter, E i -Efect of i th parameter, a i -Weighting of the i th parameter.

Prediction of the level of geotechnical risks, using an RES based model
The level of geotechnical risks was used to reflect the effects of media characteristics in the modelling of TBM performance prediction. In this research, the principles of the work carried out by Benardos and Kaliampakos, (2004) were adopted to define a model for predicting the level of geotechnical risks in rock TBM tunnelling. The merit of this modelling is that all of the important and easily obtainable parameters affecting geotechnical risks in rock TBM were taken into account.
There are two main phases for defining the RES based model. In the first phase, the geotechnical parameters (hazards), which cause the occurrence of risk in the case of rock TBM tunnelling, are defined. Also, their behaviors are analyzed and the impact (weight) that each one has in the whole risk conditions is evaluated. Further, the RES principles can be applied to evaluate the weighting of the parameters that influence the risk system.
In phase two, the Index of Vulnerability (VI), representing the level of risk, can be determined, using Equation 2 (Hudson, 1992): (2) Where: VI -Index of Vulnerability, a i -Weighting of the i th parameter, P i -Value (rating) of the i th parameter, P max -Maximum value assigned the i th parameter (normalization factor). In this equation, the level of total risk can be calculated based on the weight and individual influence of each parameter on the total risk level (the parameter's value rating). The parameter value ratings can be carried out based on the effects of these parameters on the vulnerability conditions.
Based on the VI obtained, (using Equation 2) and the classification of VI, the level of total risk in each zone of the tunnel can be identified. In the classification for VI, it is divided into three main classes with a normalized scale of 0 to 100, (see Table 4). Class I is including, small-scale problems. These problems cannot significantly affect tunnel driving. In class II, a challenging region might occur, which must be taken into consideration. In the matrix of interaction, the summation of a row is called the "cause" value ( ) and the summation of a column is the "effect" value ( ), which are the coordinates (C, E) for a parameter in the system. In the cause and effect space, the coordinate values for each parameter in the system can be plotted, which is called the C-E plot. The interactive intensity value of each parameter in the system is indicated as the summation of the C and E values ( ) and it can be applied as an indication of the significance of a parameter in the system. Also, a parameter's weighting factor ( ) is defined using the percentage value of ( ) as shown in Equation 1.
In class III, it is expected that in certain individual regions, several difficulties during tunnel excavation, are to be observed.

Parameters influencing geotechnical risks of TBM
In reviewing the TBM geotechnical risk assessment in hard ground conditions and the pieces of literature pub-  Table 5 are the most important ones from available parameters in the database and can be obtained easily. These parameters cover the most important hazards which might be encountered during rock TBM tunnelling such as instability, fault zone, squeezing, sticking due to muddy condition, and water ingress. In the selected parameters, the components of RMR such as UCS, RQD, joint spacing, joint orientation, joint condition, and water inflow are preferred instead of RMR as they describe other issues of risk raised as well as instability as a hazard.
UCS shows the compressive strength of intact rocks and also has significant effects on RMR and the associated stability (Goodman., 1989). Moreover, rocks with very high UCS cause wear and tear of cutters as well as a considerable reduction in the advance rate. Also, in rocks with low UCS, gripping problems might occur. The condition of the joints significantly influences an excavation's stability. It also affects the gripping of TBM. The total performance of TBM in joints with very rough surfaces of a limited extent and a hard wall is much better than open joints filled with more than 5 mm of gouge, or open more than 5 mm, or when joints extend more than several meters.
The orientation of the joints comparative to the tunnel axis has an impact on the chipping behavior of rocks during the excavation by TBM. Based on Bieniawski (1984Bieniawski ( , 1989, the orientation of the joints perpendicular to the tunnel axis is the most favorable for TBM tunnelling. The spacing of joints affects the rating in RMR and therefore the stability of the tunnel. On the other hand, based on the recommendations by DAUB (German Tunnelling Committee) (DAUB, 1997), in the rocks with joint spacing of 0.2 m-0.6 m, gripping problems might be encountered, for which shield-TBM suggests. In rocks with joint spacing of 0.06 m-0.2 m, the possible falling rocks might trap the shield-TBM causing serious delays.
Groundwater can strongly influence rock mass behavior (Goodman., 1989). Also, the pore pressure in rocks depends on the level of groundwater. Higher groundwater level means more pore pressure, resulting in decreasing stability, and in the case of high permeability and the presence of joints, an increase in water inflow. In return, the high inflow of water causes delays and possible damage.
Gas emission such as Methane (CH 4 ), Hydrogen Sulfide (H 2 S), Carbone Monoxide (CO) and Hydrogen Cyanide (HCN) into a tunnel might cause safety concerns such as human loss and damage due to an explosion, respiratory problems as well as corrosion damage, resulting in delays in tunnelling. In most cases, the main concern is CH 4 and H 2 S. CH 4 is a flammable gas, lighter than air and odorless. The Lower Explosive Limit (LEL) and Upper Explosive Limit (UEL) for CH 4 are 5% and 15% respectively (Industrial Training Branch of the National Coal Board, 1981).
H 2 S is a colorless and toxic gas with LEL and UEL of 4% and 44% respectively. It is slightly heavier than air with a strong smell of rotten eggs. H 2 S dissolves in water, making a weakly acidic solution. The main dangers are related to its toxic effect and corrosive property to metals. It can be recognized by smelling below 1 ppm. The 8 hour occupational exposure limit is 10 ppm, for short times (15 min), 15 ppm is acceptable. Above 15 ppm, full face masks with filters are used. This restricts the ability to work to a great extent. Above 100 ppm, H 2 S deadens a person's sensitivity to smell it. Also,  Class  I  II  III  0-33 33-66 66-100 RQD describes the fracturing degree of a rock mass. It is one of the most important components of RMR, which has a significant role in the excavation's stand-up time, indicated by RMR (Goodman., 1989). Also, rocks with low RQD cause difficulties through TBM gripping, putting excavation at risk and causing a reduction in the advance rate.

High-Very high
The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2020, pp. 57-74, DOI: 10.17794/rgn.2020.2.6 above 100 ppm it is immediately dangerous to life or health. Therefore, filters should only be used up to 100 ppm. At higher concentrations, masks with positive pressure self-supplying respirators shall be used.
Concentrations between 50 and 100 ppm have to be negotiated for longer periods. H 2 S also causes corrosion to metals, in particular to the electrical installations on the TBM. To reduce the effects, some electrical cabinets, the operator cabin, and a rest cabin have been connected to a fresh air supply (Industrial Training Branch of the National Coal Board, 1981).
Squeezing ground may lead to inadmissible deformations of the tunnel, damage to the support or in the case of mechanized excavation, to the immobilization of TBM due to sticking cutter heads or jamming of the shield. Depending on the number and the length of the critical stretches, squeezing conditions may even put the feasibility of the TBM drive into question ( The occurrence of karstic zones through tunnel alignment results in decreasing rock strength, which in return reduces the load bearing of rock, while gripping is required. On the other hand, in the case of having gases such as CH 4 and H 2 S, karstic zones facilitate the leakage of such gases into the tunnel, which raises the TBM tun-nelling risks (Marinos, 2001). The requirements, which cause karstic zones include carbonate rocks, joints, and flow of water. In the Zagros long tunnel, in limited zones, the needs for the occurrence of karstic zones are met.
Encountering the alternatives of hard and soft rocks in the tunnel face, the so-called mixed face condition might cause serious damage to the cutters, which results in a delay in tunnelling. Also, in the case of facing abrasive rocks, which are identified by their quartz contents, cutter wear might occur, which decreases the penetration rate, causing significant delays to tunnelling and economic loss.
The presence of water and rocks containing clay minerals such as shale and marl might cause swelling, which in return increase the stress on the support system used. Also, in the presence of water and such rocks, a sticky condition might be encountered, causing delays and sometimes trapping of the cutterhead.
The width of crushed zones, relevant water ingress and RQD affect the risk associated with the fault zones. With a wider crushed zone, lower RQD and in the presence of high water inflow, a higher degree of risk is expected, which in return causes jamming of the TBM, delays and economic loss.

Interaction matrix
The 15 principal parameters affecting the geotechnical risk in TBM tunnelling are located along the leading diagonal of the matrix and the effects of each parameter on any other parameter (interactions) are placed on the off-diagonal cells. The assigning values to off-diagonal cells, coding the matrix, were carried out by five experts, using the ESQ coding as proposed by Hudson (1992).  The interaction matrix for geotechnical hazards in TBM tunnelling is presented in Table 6. Table 7 gives the cause (C), effect (E), interactive intensity (C+E), dominance (C-E), and weight of each parameter (a i ) calculated by Equation 1 for the parameters affecting the hazards in TBM tunnelling. Also, the weight for each parameter is shown in Figure 3. As can be seen in Table 7 and Figure 3, the most interactive parameters are permeability, RQD, joints condition karstic zones, whereas the least interactive parameter is the mixed face condition.

Rating of parameters
The rating of the parameters' values was carried out based on their effect on the vulnerability conditions. Totally 4 classes (except UCS, which has five classes), from 0 to 3 were considered, where 0 denotes the worst case (most unfavorable) and 3 the best (most favorable condition). In the case of TBM tunnelling, the rating of each parameter is presented in Table 8. The ranges of parameters in Table 8 were proposed based on empirical results, practical limits and the experiences of different experts (five experts in this research work).

Geotechnical risk estimation of the Zagros long tunnel
Geotechnical risk estimation of 2058 sets of data from 11 km of section two of the Zagros long tunnel was carried out using the vulnerability index (Equation 2). To have a clear understanding, the parameter values and the corresponding VI obtained for one set of data are shown in Table 10. In this table, Q i and Q Max for each parameter can be obtained based on its value or description through Table 8.

Multiple regression modelling for the prediction of the advance rate of TBM in hard rock condition
In this paper, based on 2058 datasets collected from section two of the Zagros long tunnel (out of 2168 datasets), linear and non-linear regression analyses were carried out for developing a new model to predict the advance rate of TBM in hard rock conditions. Torque (T q ), cutter head rotation per minute (RPM), cutter normal force (F n ) as machine parameters, and vulnerability in- : Very coarse surfaces of limited extent, hard wall rock; A 2 : Slightly coarse surfaces, opening less than 1 mm, hard wall rock; A 3 : Slightly coarse surfaces, opening less than 1 mm, soft wall rock; A 4 : Open joints with smooth surfaces, 1 to more than 5mm aperture, filled with gouge filling or 1 to more than 5mm aperture, joints extended more than several meters dex (VI) as the effect of media characteristics were used as input parameters and Advance Rate (AR) was the output parameter. Data distributions for these parameters   are shown in Figure 4. Also, by using the database, a simple regression analysis was performed to identify the relationship between each parameter with AR as shown to carry out a multiple linear regression analysis. According to the statistical analysis, the predictive model is as shown in Equation 3: ( Where: AR -Advance rate (m/h), T q -Torque (nominal) of the cutterhead (kN.m), RPM -Cutter head rotation per minute (rotation/minute), F n -Cutter normal force (kN), VI -Index of Vulnerability.

Multiple linear regression analysis
Tq, RPM, F n , and VI as independent variables and AR as a dependent variable, using SPSS software, were used

Model Relation
Graham (1977) Farmer and Glossop (1980) Cassinelli et al. (1982) Innaurato et al. (1991) Barton, (1999) Hassanpour et al. (2009Hassanpour et al. ( , 2010 Khademi et al. (2010) PR: penetration rate, PRev: penetration per revolution, FPI: field penetration index, F n : cutter normal force, UCS: uniaxial compressive strength of intact rock, TS: tensile strength, RQD: rock quality designation, RMCI: rock mass cuttability index, σ c = uniaxial compressive strength of rock, α: the angle between the tunnel axis and the planes of weakness, Jc: RMR joint condition partial rating, RSR: rock structure rating, RMR: rock mass rating, Q TBM : Barton rock mass quality rating for TBM driven tunnels The model summary for Equation 3 is shown in Table  11. Also, to check the degree of correlation between input variables, a multicollinearity analysis was carried out. The Variance Inflation Factor (VIF), as one of the most useful tools for finding the level of multicollinearity, was used in this research. Its range varies from one to infinity. In general, VI greater than 10 shows a possibility of multicollinearity (Montgomery and Peck, 1992).
The VIF values for independent variables in Equation 3 were determined and can be observed in Table 12. As can be seen from Table 12, the VIF for each independent variable in Equation 3 is less than 10. It shows that there is no high correlation between input independent variables. Furthermore, analysis of variance (ANOVA) for Equation 3 is illustrated in Table 13. The model statistic value F and significance (Sig.) are used to provide enough proof to discard the theory of ''no effect''. From Table 13, F of 0.004 and Sig. of 0.000 (less than 0.05) were obtained, which reveals that the null theory can be discarded. It means that at least one of the input parameters significantly influences the AR.

Multiple non-linear regression modelling
For multiple non-linear regressions modelling, polynomial, power, exponential, and logarithmic models with the same datasets were used. For the polynomial model with R 2 =0.86 obtained, the mathematical equation is as shown in Equation 4: (4) Finally, for the exponential model with R 2 =0.87, the relation is as shown in Equation 7: (7)

Evaluation performance of the models
The performance of regression models were compared with the models from literature such as Khademi  Table 14), using 110 randomly selected datasets (out of 2168 datasets) from the Zagros long tunnel, which were not used in the modelling. Two signs, the determination coefficient (R 2 ) and the Root Mean Square Error (RMSE) (Equations 8 and 9) were used to do the models performance evaluation and the results obtained are illustrated in Table 15. The predicted AR from all models for 110 datasets was compared with the recorded AR as shown in Figure 6. The predicted AR from these models for five datasets was compared with the recorded AR as illustrated in Figure  7. As it can be seen from Figure 7 and Table 15 the linear regression model with R 2 =83 and RMSE=0.12 shows the best performance among the models.
Where: R 2 -Coefficient of determination, x iactual -The i th measured element, x ipredicted -The i th predicted element, n -The number of datasets.
Where: RMSE -Root Mean Square Error, x iactual -The i th measured element, x ipredicted -The i th predicted element, n -The number of datasets.

Conclusions
In this research work, using parameters of machine and media characteristics, different models were developed based on linear and non-linear multiple regression analysis to determine the advance rate in rock TBM tunnelling. In the modelling, torque, cutter head rotation per minute, and cutter normal force are machine parameters and the level of geotechnical risks (VI), representing media characteristics, which were used as input parameters and AR was the output parameter. For predicting VI, an RES based model was proposed. The merit of this modelling is that all of the important and easily obtainable parameters affecting geotechnical risks in rock TBM were taken into account.
Validation was carried out, using 110 randomly selected datasets (out of 2168 datasets) from the Zagros long tunnel, which were not used in the modelling. For 110 datasets, recorded AR were compared with the AR estimated through regression models and also with the AR obtained from the models such as Khademi et al., Hassanpour et al., QTBM, NTNU, Farmer and Glossop, Cassinelli et al., Innaurato et al., and Graham. The obtained results showed that the multiple linear regression model with R 2 of 0.83 and RMSE of 0.12 has a better performance among the different models.
As it is evident, the linear regression model was constructed based on site-specific parameters and limited datasets from the Zagros long tunnel. This model can be used to estimate the advance rate in the Zagros long tunnel or tunnels with almost the same ground conditions. However, it cannot be generalized for other tunnel pro- jects. For tunnelling projects with other geological environments, a similar approach with other parameters can be developed.