Surrogate models, also known as metamodels or emulators, are empirically generated models (or mathematical functions) that link the inputs and outputs of a computer model or a complicated system. They originated in the idea of design of experiments (DOE), which employed polynomial functions as a response surface methodology (RSM) [122]. Other metamodels have since been created and used, including kriging [18], artificial neural networks, radial basis functions (RBF) by [51]. RBF later has been extended by [76], and been called extended radial basis functions (ERBF). This method employs more than one basis function for each input which leads to a linear combination of radial- and non-radial basis functions to conduct the approximation. Unlike classic RSM, which simply calculates a least squares distance between data points, these approaches give an interpolation surface that spans all training data points. Furthermore, unlike RSM, which assumes a certain shape for the approximation as polynomials, these approaches often include a series of functions, each associated with individual points in the feature space. [14] presented a surrogate modeling process that combines the proper orthogonal decomposition (POD) method (also known as Karhunen Loe've decomposition [41]) with radial basis functions abbreviated as POD-RBF. This technique showed promising results for approximating finite element simulation responses in [10].

After reviewing the state of the art in metamodeling approaches and comparing their characteristics and considering tunneling simulation as an application [45] presented a hybrid technique that integrated POD and ERBF approaches. Khaledi et al. enhanced the POD-RBF method by replacing the radial basis functions with their extended version [76] and abbreviated it as POD-ERBF. Miro et al. [70] suggested a framework for evaluating and selecting the appropriate surrogate model for the approximation purpose.

Among many researches which used metamodeling techniques to replace complex FE models for system identification purposes as [35, 46], Zhao et al. [126] employed metamodels to develop a hybrid estimation concept. In which a small scale submodel is taken out of a larger complex model and continuous simulations are run in this submodel to reduce computing effort. Metamodels are used in three levels to i) correlate soil parameters and nodal displacements on submodel boundaries, ii) correlate soil parameters and tunneling-induced ground movements according to measurements, and iii) correlate tunneling process parameters and tunneling-induced building settlements.

The purpose of a surrogate model is to approximate the response u of a computational model f .**x**/ D u with an approximation uQ. As a result, the model output u may be given as

$$
u = \tilde{u} + \epsilon,\tag{2.29}$$

where  is the approximation error. The common approach to surrogate modeling is to perform computer runs for n sets of the input parameters f**x1**; **x2**;:::; **xn**g, and then using the generated pair points .**xi**; u.**xi**//, as *training data* to construct an approximation to the original model.

Three steps are necessary to do this, namely:

Testing data are often created by running the computational model with input parameters different from those used to create the training data. The third step is discussed further below.

#### **<IP_ADDRESS> Accuracy Measures of Surrogate Models**

The most significant aspect of a metamodel, as previously stated, is its agreement with the original data, or prediction goodness. One may become aware of the links between model complexity, sample size, and approximation method based on experience obtained from multiple applications of various techniques, but it is difficult to know the metamodel's estimation accuracy in advance. As a result, every metamodel must be evaluated before being used to guarantee that it is trustworthy.

Mahmoudi et al. [71], applied many measures to identify the extent of perfection in a surrogate model, as Root Mean Square Error (RMSE), normalized Root Mean Square Error (NRMSE), Mean Absolute Error (MAE), Mean Bias Error (MBE), Bias Factor, Nash-Sutcliffe Efficiency (NS), Variance Account Factor (VAF), Performance Index (PI), RMSE To Standard Deviation Ratio (RSR), Normalized Mean Bias Factor (NMBE), Mean Absolute Percentage Error (MAPE), Wilmott's Index (WI), Legate and McCabe's Index (LMI), Expanded Uncertainty (U95), t-Statistic, Global Performance Indicator (GPI), and Weighted Mean Absolute Percentage Error (WMAPE).

Additional techniques to investigating the correctness of the metamodel may be found in the literature. The coefficient of prognosis, which is a model-independent statistic to assess model quality, is suggested in [74]. It is used to generate a so-called metamodel of optimum prognosis, which is a rough estimate that excludes irrelevant factors. Atamturktur et al. [4] proposes a coverage metric, in which the coverage of the design space is examined in terms of appropriateness for metamodel development, and where extra samples in the parameter space should be placed to increase the metamodel's accuracy.

#### **2.3.3 Sensitivity Analysis**

In the concept of system analysis, the large number of input factors involved in a sophisticated geotechnical computational model is a challenge. The essential elements that control system reaction and must be calibrated using experimental or field data can be identified using sensitivity analysis.

The approaches used in sensitivity analysis are either statistical or deterministic. Two main categories are the so-called local and the global sensitivity analyses (GSA). To address complicated models including non linear phenomena like a 3D numerical simulation of TBM model advancement, GSA methods are recommended [68]. The performance of different global sensitivity analysis techniques was reviewed in [62]. This study examines three distinct global sensitivity analysis approaches, namely Sobol'/Saltelli, Random Balance Design, and Elementary Effect method. A decision graph is presented to aid in the selection of a proper technique, taking into account several model properties such as complexity, processing costs, and the number of input parameters. Based on the model features, the user may use it as a reference to choose the most appropriate and cost-effective strategy. Sophisticated SA approaches such as Sobol'/Saltelli give additional information as detecting non-linearity and interaction effects for a better understanding of the system. Therefore the designer may be ready to spend for their high computational costs in order to get such knowledge.

One major drawback of global sensitivity analysis methodologies is the requirement to construct a relatively large number of input-output data sets. As a result, the need of advanced numerical models, which are computationally expensive, makes such analysis impractical. The metamodeling notion presented in 2.3.2 is used as a solution.

Among the others, the Sobol'/Saltelli technique, has shown to be particularly useful for a variety of geotechnical systems [69, 127] since it prioritizes parameter relevance and highlights any interactions among them. Its main goal is to determine how the variance of input parameters affects the variance of the model output. The total effect sensitivity index S**T**i is a comprehensive index that considers the interaction of variables [91]. In order to employ Sobol'/Saltelli technique, two randomly chosen K 2n matrices **C**1 and **C**2 are constructed for a generic model with n parameters, where K is the number of samples taken from the model space.

Afterwards, a new matrix **R**i is defined by copying all columns from **C**2 except its ith column, which is taken from **C**1. Then, model outputs f .-/ for **C**1 and **C**2 are evaluated as

$$\mathbf{y\_{C\_1}} = f(\mathbf{C\_1}), \quad \mathbf{y\_{C\_2}} = f(\mathbf{C\_2}) \text{ and } \mathbf{y\_{R\_i}} = f(\mathbf{R\_i}). \tag{2.30}$$

Finally, the variance-based indices for model inputs are evaluated as

$$S\_{\mathbf{T}\_i} = 1 - \frac{\mathbf{y}\mathbf{c}\_2 \cdot \mathbf{y}\mathbf{a}\_i - f\_0^2}{\mathbf{y}\mathbf{c}\_1 \cdot \mathbf{y}\mathbf{c}\_1 - f\_0^2},\tag{2.31}$$

where **yC**1 , **yC**2 and **yR**i are vectors containing model evaluations for matrices **C**1, **C**2 and **R**i , respectively. The mean value f0 is further defined as

$$f\_0 = \left(\frac{1}{K} \sum\_{j=1}^{K} \mathbf{y}\_{\mathbf{C}\_l}^{(j)}\right)^2. \tag{2.32}$$

Variance-based sensitivity method is not merely used to identify the effectiveness of input factors, but also later employed in Sect.2.3.4 to optimize sensor localisation.

#### **2.3.4 Optimal Experimental Design**

The following explanation from [123] can describe the key motivation to pay attention to monitoring design in general: Experiments are expensive and time-consuming to conduct in order to gather a large enough set of experimental data. The goal of optimal experimental design (OED) is to create the necessary dynamic experiments in such a way that the parameters are estimated with the best possible statistical quality from the resulting experimental data, which is usually a measure of the estimated parameters' accuracy. In other words, the goal is to create the best feasible tests based on model candidates in order to make system identification straightforward. Fisher's work [27] is widely regarded as the first methodical and scientific investigation of the subject. Fisher implies that the variance of a distribution obtained as an outcome of an experimental design may be used to explain the discovered information content in order to link data to information. As a result, if some components of the monitoring or experimental design are ignored, the resulting distribution will be different and the information content will be reduced. In this way, the designer may compare which experimental design comes closest to a fully monitored scenario while requiring the least amount of experimental work.

To describe the information content I of the experimental data **y**Q with respect to the parameter mi , the formulation

$$I = \frac{\partial^2 \tilde{\mathbf{y}}}{\partial^2 m\_i}, \ i = 1, \ldots, s \tag{2.33}$$

can be employed, where s is the dimension of the considered parameter space. The matrix that accrues when considering all parameters and their correlations is often called the Fisher-information matrix FIM, whereby several similar formulations exist. In [6], the state of the art up to the 1970s is presented and the concepts of [27] are transferred to more systematic and fundamental approaches. This FIM is still employed in recent publications such as [120] that investigated where to place sensors to survey air pollution. In [52], a more statistical investigation is performed on sensor placements on dams, wherein the mean square error of parameter identification results is employed to define the quality of an experimental design. This approach was further developed by the employment of the bootstrap method in [94], where again repeated parameter identification is performed to obtain a covariance matrix **Cm** of all considered parameters **m**. Another aspect in the framework of OED is that one should differentiate between problems where the parameters of interest need to be back-calculated as it is the case in the present thesis and those where they can directly be measured. An example for the latter case is given in [75] where it is described how to find optimal sensor placements to detect contaminants in a water network. As the water demand in such a network is uncertain, finding these optimal placements is still a complex optimization problem, but it does not require the additional step of parameter identification. With different deterministic and stochastic optimization approaches, including different uncertainty scenarios, the optimal sensor locations are identified exhibiting quite different results and showing the need to consider effects of uncertainty in OED applications. The experimental settings can be modified to make the calibration of an unknown (or uncertain) input component as simple as possible. Furthermore, collecting data from tunnel construction measurements enables for the validation of numerical models and the reduction of parameter uncertainties. The utility UQ (as inverse equivalent to the cost function) defines how a certain design ı reduces the uncertainty and is to be maximized. Random output samples y are generated according to the defined initial distribution of the input parameter m and the considered "experimental design," y D f .m; ı/.

Hölter et al. [34] used Bootstrap technique to evaluate and compare the specific design scenarios and the aspect of including measurement error in a hydromechanical structure. Among all the various approaches to conduct OED, in this study exclusively employing sensitivity measures for performing OED, and Bayesian updating concept were applied. However, GSA is still performed as first essential step to identify the most relevant model parameters on which the specific OED approach is applied afterwards.

#### **<IP_ADDRESS> Spatial GSA**

The benefit of sensitivity analysis in relation to the issue of OED is initially determining which parameters are most relevant. The experimental design should concentrate on these parameters since discovering a parameter that has no impact would be wasteful of the generally limited resources. The conceptual relationship between GSA and OED is that changing the parameters to which the model response is most sensitive results in a wider bandwidth of outcomes. The greater the output variation, the easier it is to identify and distinguish the effect of particular parameters from uncertainty-related errors. As a result, using a back analysis, as explained before in Sect. 2.3.1 back, it is feasible to determine the most sensitive parameters with the greatest accuracy. First ideas to the concept of sensitivity analysis in a spatially distributed manner for geotechnical purposes can be found in [34]. Modern FE programs enable model answers to be provided for any point of the FE model geometry. A GSA may be done for each of these spots to utilize the position with the highest sensitivity to a given parameter to insert a sensor to identify this specific parameter. However, it should be noted that the sensitivity index STi calculated using the variance-based technique yields normalized values in each spot for all the considered parameters. Therefore, they may not be compared in various locations. This issue can be overcome by including the importance of an output via its variance and using the updated sensitivity index S Ti ,

$$S\_{T\_{i,j,k}}^{\*} = \frac{S\_i \ \sigma\_{j,k}}{\max\_k \ \sigma\_{j,k}},\tag{2.34}$$

where, S Ti;j;k denotes the well-known sensitivity index as it was introduced in Sect. 2.3.3 for a certain parameter mi , but specifically for the j th model response obtained at the kth position. j;k denotes the standard deviation of the j th model response obtained at the kth position and maxk j;k corresponds to the maximum of the standard deviation of the j th model response obtained from all of the considered positions. Contour plots may be constructed via interpolation after getting S Ti;j;k in each of the grid points, allowing for improved visualization and explanation of the data. In terms of planning an experimental setup, the contour plot areas with the greatest values of S Ti;j;k are the optimum places to install sensors of the matching response type j . The various values of S Ti;j;k should be employed as objective function weighting factors wi in the cost functions defined by Eq. 2.28.

#### **<IP_ADDRESS> Bayesian OED**

For several decades, applications of Bayesian updating have been applied in various fields. Miro [68] and Nguyen [78] employed Bayesian interference to update 3D simulation models of tunnel excavations. The uncertainty of the prior estimation of the soil parameters is therein reduced stepwise using measurement data that is obtained after each new excavation step of the tunnel, but each time affected with artificial random noise. As the measurement data used for back calculation is at the same time an extension of the current data knowledge, one might be interested to "update" the previous knowledge. This concept was first formulated in [7] as

$$p(m|\mathbf{y}) = \frac{p(\mathbf{y}|m)\ p(m)}{p(\mathbf{y})}.\tag{2.35}$$

Herein, p.y/ describes the probability distribution function (PDF) of model output that is used to normalise the result. Here, this would typically be noisy measurement data **y**Q that is somehow distributed around its mean value yN. p.m/ accordingly is the PDF of the model parameters. These would usually be soil parameters that are often very uncertain in geotechnical engineering. Within the Bayesian concept, this distribution p.m/ is called as "prior knowledge." Accordingly, p.mjy/ is called the posterior knowledge as it includes the new information y. The link between prior and posterior distribution is the likelihood p.yjm/ that describes how probable the new data y is, given the known probability distribution of m.

#### **2.3.5 Case Study: The Metro Line 5 in Milan**

The Metro Line No. 5 in Milan is selected as an exemplary realistic case study to show the use of previously mentioned approaches. This route, which runs from the north to the west of the city, is composed of 2 12.6-kilometer tunnels with a center-to-center distance of 16.7 meters between them. The first tunnel was completed in May 2012, and the second tunnel was completed a month later in June. The outside diameter of the tunnel tubes is D D 6:7 m. Tunnel excavation was carried out utilizing EPB machines (tunnel boring machine with earth pressure balance shield). In the segment under consideration, the tunnel axis is 15 meters below ground level. [24] developed an FE investigation to analyse the soil-structure interaction behavior that is induced by the construction of the twin-tunnel. This simulation of three soil layers using the HSsmall model [8], was used in our study as well. The subsoil consists of a 20-meter-deep gravelly sand layer, lies beneath a 5-meterthick sandy silt soil layer, according to the geotechnical soil research. Finally, a gravelly sand with a thickness of 5 m (i.e. 25–30 m below ground level) was discovered. For the finite element simulation, drained conditions are considered because of the high permeability of the soil (gravelly sand). This model was created with the FE-code Plaxis 3D and is presented in Fig. 2.24 without the soil clusters. Table 2.2 lists the soil parameters derived from the soil study and subsequent laboratory testing, as reported by [23, 24]. The considered twin tunnel is excavated in an urban area and partially underpasses a nine-story building. The building that is specifically considered in this model is selected because an extensive monitoring program has been applied to it. Figure 2.25 depicts the building's location and relationship to the two tunnel tubes. The shallow foundation of the residential structure comprises five strip footings and intermediate concrete slabs. With a height of 0.65 m, the strip footings are placed around 4 m below ground level. Figure 2.25 also shows the monitoring arrangement, comprising 13 sensors placed throughout the building. The sensors capture local vertical displacements. At the same time, the location of the TBM, applied face pressure, and grouting pressure are recorded. This allows the numerical model to be validated, as well as, to some extent, the outcomes of an OED procedure.

The ultimate goal of design, monitoring, and back calculation is to prevent compromising the integrity and serviceability of subsurface and above-ground structures. If an appropriate model that can very accurately anticipate building behavior is produced, countermeasures that lead to improved building safety may be more effectively assessed. As

**Fig. 2.24** FE model of considered interaction of twin tunnel bypassing the nine-storey building, [103]

**Fig. 2.25** Top view of the considered building and location of applied measurement points [103]

a result, the goal becomes to produce a trustworthy model by identifying the important soil features using the approach of back analysis outlined in Sect. 2.3.1, which will benefit from incorporating OED principles.

In the previous studies by the authors, e.g., in [103] the in-situ measurement data were employed to back-calculate the constitutive parameters of the surrounding soil. In that study, the grouting pressure pv and face pressure ps, the site-related volume loss factor VL, the friction angle '1, the small-strain stiffness G0;1, and the secant stiffness Eref 50;1 of the first and the second soil layer, Eref 50;2 were evaluated to be the most important parameters. It should be noted that the secant stiffness Eref 50 is supposed to be associated with the tangent stiffness Eref oed and the unloading-reloading stiffness Eref ur , by a factor of three. When the term "stiffness" is used in the following, it refers to those three factors. Only the measurement data received until the first tunnel tube is completed is used for model validation, so that a forecast for the excavation of the second tunnel tube may be made and compared to the measurement data. As a consequence, the findings shown in Fig. 2.26a, are achieved. The measured settlements obtained after excavation of the first tube in the thirteen measurement points shown in Fig. 2.25 are displayed here, along with the settlements that correspond to the identified parameters using the actual FE-model and the corresponding metamodel. The determined set of parameters is utilized for subsequent experiments and for the prediction of the settlements produced by the excavation of the second tube, as illustrated in Fig. 2.26b. The findings of the FE-model are compared to the data when the second tube is completed. With the exception of R4, T1, T2, and T3 points, there is still a strong agreement. One should keep in mind that, in the considered urban region, there are other buildings just adjacent to the one being studied, that produce pre-stressing of the subsoil and consequently an increase in stiffness, resulting in the lower settlements observed in reality.

Several GSAs are done in [103] to gain a better knowledge of the system behavior, linking the aforementioned soil and system characteristics to average settlements and tilting along different axes of the building. The results of the GSAs performed following excavation of the first and second tunnel tubes are presented in Figs. 2.27 and 2.28, linking the seven examined input parameters to the vertical settlements of each of the thirteen measurement locations. When looking at the shown findings, it is clear that they change with time and space, but that the parameters Eref 50;1, G0;1, and VL have the most impact on the analyzed outcomes. As a result, all subsequent research, have focused only on these three parameters.

In a further step, a second GSA was run for the previously discovered key parameters to examine the evolution of the sensitivity of the model response in time and space domains with respect to each of the key parameters defined from GSA. The main idea of spatial GSA as a tool for OED concept is introduced in [37]. This GSA is used in the same way as previously to find geometrical sections that are particularly promising for sensor positions, but with the output values of different places in the relevant area. Schoen et al. [102] analysed the spatial and temporal variations in the sensitivity of different parameters in various construction stages, namely during the passage of first tunnel, at the end of **Fig. 2.26** Comparison of measured and **a** back-calculated data of first tunnel excavation and **b** measured and predicted data of second tunnel excavation [33]

ST*i*

ST*i*

excavation of the first tube, and finally when both tunnels are bored. They concluded that the influence of the secant stiffness Eref r on the settlements increases with progress of the tunnel excavation. Moreover, VL has a high influence on the settlements around the TBM as surface settlements are directly related to the volume loss. Zhao et al. [125] also showed that the sensitivity of the volume loss VL depends on both the location of the TBM, and the magnitude of the surface loads above the tunnel. The spatial GSA findings reveal that the impact of each input parameter changes geographically and temporally depending on the excavation time step and location. As a result, the usefulness of monitored data for parameter identification changes depending on the position of the sensor and the period of monitoring. Therefore, in parameter identification process, it is far more necessary to utilize data from those sensor locations that are particularly impacted by the respective input parameter at the most. Afterwards, the GSA findings are utilized as weighting factors for each monitoring spot and each input parameter in the cost function for back analysis. These weighting variables are chosen based on spatial GSA results to provide more weight to those measuring places where the effect of the corresponding input parameter is considerable, and significant tunneling-induced settlement can be accurately recorded. Therefore, their values vary in different construction phases. Including these weighing factors in the back analysis improves the simulation model's accuracy, bringing it closer to the in-situ records.

Continuation of this research in relation to monitoring optimization, Hölter et al. [36] employed the introduced Bayesian OED using the Approximate Coordinate Exchange (ACE) in Sect. 2.3.4 to use obtained prior knowledge for improving the further searches. After excavating the first tunnel tube, the vertical settlements are acquired at a grid of locations separated by 5 m within the constraints of the employed FE-model. Latin hypercube sampling technique is used to build a sample set of 120 combinations of the three parameters of interest, to generate a metamodel for the FE-model. The metamodel is generated based on quadratic polynomial regression, allowing the settlement result to be received in any of the 259 parameter combinations. To calculate settlement data between these points, the polynomial interpolation approach described in [2] is used.

It is specified that six sensors should be located in this region. Because each sensor location may be swapped in both dimensions of the surface, determining the ideal design ı becomes a twelve-dimensional optimization issue. To begin, a random 62 matrix inside the provided geometric limits should be simply constructed. It should be noted that the initial design ı0, which determines where the ACE algorithm begins its search, has a significant impact on the method's performance. For example, if certain places of the model's boundaries which are not affected by the building or tunnel excavation at all are chosen, the calculations result in high utilities and the model reaction will always be the same. Accordingly, none of the investigated current designs is accepted as new best design ı and the initial design ı0 is always kept. As a result, it is reasonable to position the starting coordinates of the six sensors anywhere within the tunnel's area of impact, such that the different back analysis runs result in varying responses of the utility function u.

The obtained results are shown in Fig. 2.29. The settlement at the ground surface in the whole model domain is illustrated as a contour plot between red and green for one typical parameter combination, with the position of the building noted as a black block. The locations of the sensors in the first design ı0 are represented by white squares. They are clearly located in locations where big settlements occur for the reasons stated above. The black crossmarks represent the ultimate design to which the algorithm converges, i.e. the design with the highest utility. The black markings may be seen moving away from the initial design and towards the border of the settlement basin. It should be noted here that the settlement distribution refers to a single arbitrary realization of the input parameters **m**. The settlement distribution may be broader or smaller for various combinations. The resulting sensor configuration does not correlate to a specific set of parameters, but it is expected to be the most reliable throughout the studied range of values.

#### **<IP_ADDRESS> Reliability Measures Updating**

As mentioned before, engineers use inverse analysis techniques to increase the dependability of proposed infrastructure projects such as mechanized tunnels by evaluating the uncertainty propagation of soil parameters using recorded data. During different excavation phases, the excavation-induced settlement of the surface structures impacted by tunneling projects is commonly monitored. However, the strength characteristics of geomaterials, on the other hand, shows spatial variability owing to their natural origin. Model adaptation, which is primarily concerned with using previously acquired dataset in further prediction, may be used to precisely analyze the impact of uncertainties on dependability measurements. In this study, a recently developed approach for Bayesian updating of reliability is used while applying random field techniques to account for the geographical variability of the subsoil. The revised reliability measures may be calculated. They can be used to identify the right steering settings or to select adequate countermeasures in order to minimize long-term damage.

In this procedure, the spatial random variables are adjusted depending on the measured data from the mid-term excavation phase. After that, the updated dependability metrics for subsequent phases are being estimated. The introduced method is applied on the finite element simulation of the twin tunnel project introduced in 2.3.5. The collected findings demonstrate how the reliability update concept may be used to determine the long-term viability of an initial geotechnical design for a mechanized tunneling project.

#### **<IP_ADDRESS> Random Field**

In general, mechanized tunneling is often regarded of as a massive geotechnical project over a broad region. Almost all attributes of subsoil must be handled as heterogeneous materials owing to variation in lithological composition, stress history and sedimentation, weathering and erosion, probable faults and fractures. As a result, for every stochastic analysis, peculiarities of the spatial distribution of geological and geotechnical should be characterized and taken into consideration.

Regression analysis, geostatistics, and random field theory are three basic kinds of statistical approaches for estimating spatial variability that have been established. Regression techniques presume that all sample values within a medium have an equal and independent probability, i.e., that they are not autocorrelated, which contradicts field observations show that samples obtained closer together have a larger association than ones taken further apart. As a result, regression approaches are better suitable for preliminary analysis when just a small number of samples with wide distances between them are collected from the site [87]. However, the autocorrelation approach in geomaterials is taken into account in the other two categories.

Here, the Karhunen-Loève expansion method, which is extensively utilized in random field theory, is applied [29, 109]. The spectral decomposition of the covariance function is used in the Karhunen-Loève expansion, which may be written as

$$\hat{H}(\mathbf{x},\boldsymbol{\theta}) = \mu\_H(\mathbf{x}) + \sigma\_H(\mathbf{x}) \sum\_{i=1}^{M} \sqrt{\lambda\_i \,\,\Phi\_i(\mathbf{x}) \,\, \xi\_i(m)}. \tag{2.36}$$

The eigenvalues and eigenfunctions of a Fredholm integral equation, where the autocovariance function (i.e., autocorrelation function multiplied by the variance of the random field) is specified as a kernel, are i and ˚i.x/, respectively. The mean and standard deviation are represented by H .x/ and H .x/, respectively. M is the number of series terms in the equation, and i. / is a vector of uncorrelated random variables that represents the unknown parameters' randomness. This vector will be subjected to modifications in a Bayesian way later on in the model adaptation context. The further details regarding random discretisation can be found in [61, 63]. A random field discretisation example of Elastic modulus is shown in Fig. 2.30.

**Fig. 2.30** Random field exemplar

#### **<IP_ADDRESS> Bayesian Back Analysis**

With a given performance function Gx, the failure event is defined as

$$P\_F = P[G\_x \le 0] = \int f\_{(x)} \mathrm{d}x,\tag{2.37}$$

where f.x/ is the joint probability density function of x, x being the vector of uncertain input variables.

Straub [106] presented the Bayesian back analysis approach, which is used here.

Employing the aforementioned Milan Metro tunnel in Sect. 2.3.5, Mahmoudi et al. [59], performed a reliability analysis at the end of the second tunnel excavation, merely considering the spatial uncertainty of material parameters. Afterwards, the Bayesian reliability updating concept is utilized to take into account the settlement recording at the end of first tube excavation. The ground measurement is assumed to be recorded equal to 12, 9 and 2 mm at the measurement point from left to right after the excavation of the first one is completed. These values are less than average amounts observed for the homogeneous cases. Figure 2.30 depicted a typical random field discretization generated by Karhunen-Loe've expansion. The spatial arrangement of random field will be changed after applying the Bayesian updating technique. The results show a particular change in the reliability measures prediction through considering the mid-term measurements. In addition, the effect of the arrangement and accuracy of the sensors on the reliability measures were investigated. Table 2.3 shows the change in the failure probability indices before and after Bayesian updating. Here the Mid-term measurement represent the monitored ground settlement at the middle of two tunnels. The considered error for the measurement is noted by m, the updated reliability index and correlated probability of failure were shown by ˇ and Pf , respectively. The reader is referred to [59] for detailed information about the applied techniques to obtain the reliability measures and updating them.

#### **2.3.6 Pattern Recognition**

Geotechnical construction has historically employed sensing technology for a variety of objectives, including identifying geological irregularities, health structure monitoring, and ensuring integrity. Using common instruments in geotechnics as e.g., strain gauges, load cells or piezometers, a vast amount of data in the form of time series will be provided, specifically, in large projects like tunneling. The high dimensionality and feature correlation of large time series databases, along with noises, is an evident challenge when studying them. This challenge can be overcome by using automated pattern recognition algorithms to time series for translating a dataset representation of an item or relationship to an output category. The data logs of pore water pressure and ground settlements recorded by a TBM during tunnel excavation are assessed in this work. Due to the different stress distribution in the soil domain as a result of different geological formation, the resultant vertical displacements and water pressure logs may vary substantially during tunneling. As a result, the changes in sensor data trend or motif may be identified and labeled across different geological conditions through supervised learning methods. Pattern recognition is referred to as a supervised learning process when a given pattern is assigned to one of the pre-defined classes using pre-labeled data to construct a model. A supervised machine learning (SML) algorithm can identify correlations and patterns in the fluctuation of observed data and correlate it to potential geological formation changes.

The Support Vector Method (SVM) [17], is used in this work as a well-known approach for conducting SML. Several scholars have used SVM in geotechnical engineering problems, tunnel excavation design, and reliability analysis [38]. Yao et al., [124] used SVM to forecast the displacement of the surrounding rock of a railway tunnel for tunnel safety estimation. Mahdevari et al. [58] used this technique as a nonlinear regression approach to forecast tunnel convergence after excavation. The results of performing SVM and K-Nearest Neighbor (KNN) analysis were presented by the authors in [60]. The learning process in human brains is the motivation for neural nets. A neural network is made up of a linked graph of perceptrons, each of which receives input from previous levels and passes output to subsequent ones. The weight supplied to each individual link, which is determined by the cost function and the optimizer, determines how each layer output becomes the input for the following layer. The ultimate findings of a classification task are really probabilities. The input dataset will be assigned with the label that received the highest probability. A neural network, on the other hand, has a lot of hyperparameters to adjust. To successfully comprehend diverse datasets, multiple sets of hyperparameters are required. However, the enormous number of hyperparameters makes it difficult to identify their values to achieve an accurate prediction/classification results. Therefore, finding the appropriate settings of hyperparameters to create the model from a certain dataset requires hyperparameter tuning. For instance, the number of neurons, activation function, optimizer, learning rate, batch size, and epochs are the hyperparameters to be modified. Although, there are several techniques concentrating on optimizing the hyperparameter evaluation, the problem of determining the right combination of hyperparameters remains a challenge. In this study, a network for optimization using the stochastic gradient descent learning technique established in [89] is used. It includes five hyperparameters, namely, the number of neurons in hidden layers, epochs, learning rate, batch size, and lambda value. In order to tune these parameters a grid search is conducted.

In the following, the neural network concept as another supervised learning method, is applied to predict any geological distortion ahead of a TBM excavating a 8.5 m tunnel in the depth of 8.5 m. In order to train and validate the network, the gathered data from 10 ground settlement and ten pore water pressure sensors are employed. The details regarding the system configuration is explained by the authors in [88].

In order to measure the accuracy level of trained neural network, various indications can be used, namely accuracy, recall, precision and F1-score. They can be obtained as

$$\text{Accuracy} = \frac{T\_P}{T\_S},\tag{2.38}$$

$$\text{Precision} = \frac{T\_P}{T\_P + F\_P},\tag{2.39}$$

$$\text{Recall} = \frac{T\_P}{T\_P + F\_N},\tag{2.40}$$

$$\text{Fl-score} = 2 \frac{1}{\frac{1}{\text{Precision}} + \frac{1}{\text{Recall}}},\tag{2.41}$$

where TP indicates the number of cases in which machine predicted the right class, and TS is the number of all the validation test data. FP and FN denotes the false positive and false negative cases. Regarding multi-class cases as the present study, the precision and recall values are firstly calculated for all the classes, and finally use the mean value as a general metric.

The calculated metrics in various distance between the TBM face and the geological anomaly is represented in Table 2.4. In this neural network, 5% white error is assumed for the pore water and settlement gauges. The obtained results show acceptable accuracy level for all the cases, however, the prediction accuracy increased when the distance decreased.

The highest accuracy in classification, according to the data in [59], belongs to the combination of kriging and KNN, which improves recognition accuracy even when the excavation is taking place at a distance from the anomaly. However, the SVM and neural

**Fig. 2.31** Overview of the three-stage exploration approach. (Reproduced with permission from Elsevier from C. Riedel et al. "A hybrid exploration approach for the prediction of geological changes ahead of mechanized tunnel excavation". In: *Journal of Applied Geophysics*, 203, p. 104684 (2022))

network accuracy is sufficient for scenarios when the anomaly is predicted in the near field, despite the fact that it requires less computing work.

#### **2.4 Three-Stage Concept**

The previously explained methods of the current chapter can be combined into one hybrid method to improve the prediction of the subsoil properties as well as to reduce computational demand. In [88], a three-stage concept is presented which combines the pattern recognition approach (Sect. 2.3.6), UHSA (Sect. <IP_ADDRESS>), and the adjoint approaches in time (Sect. 2.2.4) and frequency (Sect. 2.2.6) domain. The concept is shortly explained in the following, where an overview of the three-stage exploration approach is given in Fig. 2.31.

In the first stage, pore water pressures and ground settlements are continuously recorded while the tunnel is excavated. The acquired data is processed by the pattern recognition approach in order to identify a potential anomaly. To do so, various machine learning algorithms are trained for different scenarios, e.g. a fully homogeneous model, an interlayer, a layer change or a boulder. The applied supervised methods, namely SVM and KNN shown high accuracy to classify the geological condition ahead of TBM, merely based on the monitored settlement and pore water pressure variation. If the SML approach successfully identifies an anomaly, the TBM is stopped in order to perform a seismic survey. Subsequently, the second stage is initiated, which is FWI with UHSA. Here, all available knowledge derived from the SML approach is processed, where the predicted scenario of the SML approach is used to implement a parametrization of the anomaly. UHSA might already bring a sufficient prediction; however, it is probable that some parts of the anomaly are not reconstructed in its entirety because the chosen parametrization might not be able to describe all features. Therefore, the third stage is executed in the form of adjoint FWI, where the final model of UHSA is used as the initial model. Both the time domain and the frequency domain approach can be used in order to finalize the prediction of the subsoil properties and to image further details. In [88], the three-stage concept is successfully tested with synthetic data, where the reference model contains a boulder. Each of the three stages can improve the image of the anomaly and on the basis of the final prediction, the location, shape and material properties of the boulder can be precisely estimated.

#### **2.5 Soil-Shield Machine – Interactions**

This chapter summarizes some result from a research focused on identification of ground conditions by analyzing machine data during excavation of shield driven tunnels in soils. The goal of this research is to utilize TBM machine data for estimating soil properties. Particularly data about movement of the cutting wheel and the shield with respect to subsoil conditions have been analyzed. Towards this goal, the geological and geotechnical properties of subsoil have been continuously evaluated and the raw machine data collected at a hydro-shield project. The raw data have been separated between active and passive parameters: Active parameters are set by the operator and include the cutting wheel rpm and pressure of thrust cylinders. Passive parameters (thrust force on cutting wheel, torque on cutting wheel, advance speed and rate of penetration) are the result of active parameters. The latter are termed excavation-specific data components which may serve as an indicator for the influence from the ground's resistance to excavation (ground conditions) and for the influence of the state of the cutting wheel (design of the wheel, state of tool wear and clogging), respectively.

#### **2.5.1 Excavation-Specifc Data Components**

Table 2.5 distinguishes between excavation-specific and excavation-independent raw TBM data. Those data must be revised for arriving at meaningful interpretations of interactions between a TBM and the ground [22, 32]. For example, the total thrust force FTh of the shield machine (raw data) must be stripped from excavation-independent data as shown in Fig. 2.50 to arrive at data of tool forces FT to be compared with the soil condition. The detailed procedure is given in [22].

To connect machine data with soil conditions the following should be considered. The penetration rate Penraw is a target value in shield tunneling and requires a certain contact force Fcon of the cutting wheel for soil with a given resistance to be excavated. Additionally, the torque at cutting wheel Mraw is a function of the specified penetration rate Penraw and the resistance of the soil to be excavated. It is defined accordingly

There are some general connections of the resistance of the soil with the specific torque moment and the specific penetration as shown in Table 2.6.

**Table 2.5** Excavation-specific and excavation-independent components of machine raw-data

**Table 2.6** Tendency of excavation resistance for the main soil types (SPT-values [67, 112])

**Fig. 2.32** Schematic presentation of forces acting on cutting wheel of a hydro-shield machine: FT : tool forces; FS : resulting force of face support pressure acting on area of cutting wheel drive; FB : internal friction forces of main bearing; FC : internal friction force of hydraulic cylinders; FTh: thrust force on cutting wheel raw-data [21]

#### **2.5.2 Case Study: Application to a Hydro Shield Project**

A reference project, involved a 2.27 km long tunneling project of diameter 9.5 m driven through terraced rubble of the quaternary near the river Rhine is considered. Details about the project, the TBM and the geology and can be found in [21]. The soil consists of sand and gravel as shown in Fig. 2.33.

The raw TBM data were processed as described above and plotted in Fig. 2.34. The wide spread of the specific penetration can be attributed to different geological conditions at the face as shown in Fig. 2.35.

To compare different projects with different shield machines the specific contact force (SCF), which is contact force Fcon divided by the area of the cutting wheel ACW, is computed. The Force Index FI is introduced as

$$\text{FI} = \text{Pen}\_{\text{raw}} / \text{SCF},\tag{2.42}$$

and the Torque Index TI,

$$\text{TI} = M\_{\text{spec}} / A\_{\text{CW}},\tag{2.43}$$

and show the plots in Fig. 2.35.

**Fig. 2.33** Typical soil (left) and screenlines (right) at the hydro-shield drive [21]

**Fig. 2.34** Plot of specific torque vs. specific penetration at the hydro-shield project [21]

**Fig. 2.35** Correlation of data with site conditions [21]

The different lines of the TI-FI plots reflect the different designs of the cutting wheel and TBM specifics. Analysis of the soil conditions at the various projects leads to Fig. 2.36. The torque of a shield TBM is mainly attributed to the shear strength of the soil and contact force is attributed to the packing density or consistency of soil.

For future studies it is possible to employ soil characteristics from site investigation to estimate the magnitudes of FI and TI, respectively. For a given penetration rate (the target parameter) the torque and the contact force of the cutting wheel may be tentatively estimated.

**Fig. 2.36** Force-Index vs. Torque-Index plots for a wide range of case studies [21]

**Fig. 2.37** Dependency of the torque index to the force index with respect to the shear strength and packing density [21]

#### **2.6 Monitoring of Building Damages by Radar Interferometry**

One important aspect of the settlement damage research was the monitoring of building damages using both traditional terrestrial- and advanced satellite-based monitoring techniques. Also, a combination of terrestrial and extra-terrestrial proved to be of valuable service to TBM operators. Of course, this presupposes proper damage assessment methods with validated mathematical and engineering theories are readily available.

#### **2.6.1 On the Importance of Settlement Monitoring for the Management of Damage for Inner-City Tunnel Construction**

The continuous increase in traffic volumes in metropolitan areas gives rise to further construction of underground transport facilities. Typical examples are road or subway tunnels usually driven by mechanized shield tunneling. This trend densifies the traffic facilities in the underground and consequently reduces the coverage of existent building structures [82].

With inner-urban tunneling projects reports on actual or potential damages due to settlements of existing structures often predominate in the press. Such reports unmistakably reflect the concerns of citizens and residents and let problems of acceptance arise that must be expected to extend the construction period already in advance. Thus, in particular the risk of damage to existent structures above-ground comes to the attention of developers and construction companies at an early stage [66].

A look into the past shows that this topic was already taken up in the 1950s when some rudimentary recommendations for the prediction of structural damage were first developed. However, significant influences from the settlement event itself, as well as the building materials and key properties of the structures, were not taken into account. It was not until the 1970s that Burland & Wroth [16] published their approach that is still current today, which made a semi-empirical prediction of potential damages to buildings due to tunneling-induced settlements possible.

Up to now published semi-empirical approaches to compute vertical settlements can be differentiated with reference to Fig. 2.38 into

In addition to its stiffness, the effective displacement of the structure is of central importance in the evaluation of structural damage caused by tunneling-induced settlements. Due to (too) far-reaching simplifications, e.g. with regard to the often highly heterogeneous building stock, the approach according to Burland & Wroth shows deficits that

may lead to conservative but also to uncertain evaluation results. In tunneling practice, in particular (too) conservative evaluation results, intended to protect the structures, lead to expensive and time-consuming auxiliary measures (e.g. compensation injections).

#### **2.6.2 Recent Methods of Damage Assessment Due to Tunneling Induced Settlements**

Skempton & MacDonald [105] developed the first purely empirical approach to evaluate damages to above-ground structures due to settlements. The evaluation criterion for damages therein are the so-called angular distortions =L or ˇ, i.e., purely shear-based structural deformations. Limiting values of these distortions classify the extent of damage according to Table 2.7.

Burland & Wroth (1974) [16] and Burland et al. (1977) [15] recognized the deficits of a purely shear-based approach and published an enhanced semi-empirical approach to damage assessment, the Limiting Tensile Strain Method (LTSM). The basis of their approach is the transformation of the structure into an equivalent mass-less beam with a shear extension according to Timoshenko (1955) [113]. They recommend the evaluation

**Fig. 2.39** Definition of structural rotation and translation: deflection ratio =L, tilt !, rotation , angular strain ˛ and angular distortion ˇ [95]

of damage by means of the limiting tensile strain "lim, which indicates the onset of cracking and is defined as the ratio of the relative deflection of the beam to its length L, the so called deflection ratio (cf. Fig. 2.39). Assuming the neutral axis at the center of gravity of the beam, they derived the analytical equations Eqs. 2.44 and 2.45. Thereby "b;max is the maximum bending strain due to direct tensile loading and "d;max is the maximum shear strain caused by diagonal tensile loading (cf. Fig. 2.40, middle).

$$\frac{\Delta}{L} = \left(0.167 \frac{L}{H} + 0.25 \frac{H}{L} \frac{E}{G}\right) \varepsilon\_{b,\text{max}},\tag{2.44}$$

$$\frac{\Delta}{L} = \left(0.65 \frac{L^2}{H^2} \frac{G}{E} + 1\right) \varepsilon\_{d,\text{max}}.\tag{2.45}$$

In Eq. 2.44 and Eq. 2.45, the stiffness of the structure is estimated by the moment of inertia of the equivalent substitute beam with length B, height H, shear modulus G, and Young's modulus E. The dimensions of the structure and the substitute beam are identical. Roof superstructures are not considered in the determination of H (Fig. 2.39). The influence of the estimated bending stiffness EI on the shape of the settlement trough is neglected in their approach, which makes the evaluation of the structural damage again conservative. The green-field settlement trough itself – as the damage-inducing event – can be determined using the approaches from Sec. 2.6.1 [97].

**Fig. 2.40** Damage assessment scheme acc. to Burland and Boscardin [95]

The damage categorization proposed by Boscardin & Cording [11] finally relates Burland's limiting tensile strain "lim to the expected damage at a structure. Moreover, damaging contributions of horizontal strains could be integrated to the LTSM approach which practically comprises four steps (cf. Fig. 2.40):

To limit the settlement trough's region of influence, Mair et al. [64] proposed the 1 mm settlement line. For long structures that span several regions of the green-field settlement trough, only those parts within that boundary are to be considered.

#### **2.6.3 Overview and Options for the Monitoring of Settlements**

As mentioned above, the two most important approaches to settlement monitoring include traditional terrestrial monitoring and advanced satellite-based monitoring. A strategic combination of terrestrial and extra-terrestrial monitoring proved to be even more advantageous when the synergetic aspects of these methods were exploited.

#### **<IP_ADDRESS> Terrestrial Monitoring**

Above all, the risk of damage can be more accurately assessed, controlled, and also predicted by monitoring or surveillance of above-ground structures. State of the art are terrestrial methods that record mainly vertical displacements of structures due to shield driving at the foundation or at the street level (e.g., barometric by tube water leveling). However, three-dimensional monitoring, i.e. the measurement of displacements in all coordinate axes including rotations, among others, with fully automatic total stations (Fig. 2.41), has already been proven suited several times. A recent example is the new metro crossing Amsterdam from north to south, where fully automatic measurement systems were used to monitor the facades of settlement-sensitive historic buildings. Regardless the devices used, all methods must

The accuracy of terrestrial methods is usually ˙0.2 mm per 100 m for total stations; leveling (e.g., barometric by tube water leveling) achieves less than 0.2 mm (per 100 m)

**Fig. 2.41** Left: reflector applicable on tripods or an a facade, right: Measurement with tachymeter (total station) acc. to [95]

or up to 0.1 mm. Terrestrial settlement monitoring does not consist of individual measurements, but of the interaction of all (settlement) measurements before, during and after a tunnel construction measure. The boundary conditions and the scope of measurements, including the methods used, are defined in a so-called measurement program. This takes into account the expected settlement sensitivity (vulnerability) of individual structures [66, 82].

#### **<IP_ADDRESS> Satellite-Based Monitoring by Radar Interferometry (Synthetic Aperture Radar, SAR)**

Since the launch of the German radar satellite TERRASAR-X, remote sensing systems for settlement monitoring in large urban areas have become competitive due to the increased geometric resolution of the radar system [5]. On-site installations, e.g. in buildings, are no longer required. The method records complex raw data signals (amplitude and phase information) of radar waves (electromagnetic waves) reflected at the Earth surface (Fig. 2.42). Suitable natural radar reflectors are metallic objects or parts of buildings, e.g. a window sill. These are visible as bright spots in the radar image (SAR image, Fig. 2.43). Water surfaces and vegetation in forests or open spaces reflect the radar waves poorly or even not at all and thus stay dark in radar images. Therefore, radar-interferometric settlement monitoring is particularly suited in urban areas. For open spaces, additional measures as installation of artificial reflectors becomes necessary.

The signals recorded by the radar sensor must be spatially and temporally stable to allow for superposition (interference) of two wave signals (coherent radar signals). By com-

**Fig. 2.42** Principle of radar interferometry (simplified) according to [19]

**Fig. 2.43** Radar scan by SAR with reflections of different strength [19]

paring two phase positions i recorded from slightly different satellite positions (SAR1 and SAR2) at the same position of an object, the interferometric phase ,

$$
\phi = \phi\_{\text{flat\\_Earth}} + \phi\_{\text{topo}} + \phi\_{\text{disp}} + \phi\_{\text{path}} + \phi\_{\text{noise}} + n2\pi,\tag{2.46}
$$

can be calculated (phase shift). If the other components of the interferometric phase remain sufficiently stable, the phase fraction disp can be used to infer the displacement, e.g. of a natural radar reflector like a window sill. The displacement is then to be converted from the satellite's "line of sight" into the sought components of the vertical settlement or uplift. With the SAR method for acquisition and evaluation, the physically limited resolution in azimuth direction is increased by combining multiple images from different satellite positions.

However, in general, the interferometric phase is composed of several components that can have a significant impact on the accuracy of the method. Beside the target deformation disp the further components are assigned to the impact of the Earth curvature flat\_Earth and to the height model of the Earth surface topo. path and noise capture perturbations from the atmosphere and noise.  If the expected settlements are larger than the wavelength used, the evaluation becomes complex and erroneous. Further details on the technique used and the procedures (D-InSAR, PS-InSAR) for the evaluation of radar images can be found in the literature [3, 40].

Comprehensive analyses in [100] prove that an accuracy of less than 1.5 mm (standard deviation) can be reached with a sufficiently large sample. For the TERRASAR-X radar satellite the maximum repetition rate today, i.e., the time interval between two subsequent overflights, is 11 days. Thus, live monitoring of settlement data accompanying the excavation process is not possible yet. But, future satellite missions are already scheduled that will reduce the repetition rate down to 4 days by combining several satellites. The remarkable density of settlement information at so-called persistent scatterers (PS), i.e. at long-term stable natural or artificial radar reflectors, amounts to approximately 25,000 PS/km2 in inner-city regions.

**Fig. 2.44** conceptual combination of terrestrial and radar-interferometric techniques in tunneling and its spatial coverage acc. to [19]

#### **<IP_ADDRESS> Combination of Terrestrial and Extra-Terrestrial Techniques**

If the advantages of both terrestrial and radar-interferometric methods are strategically combined, both performance and cost-effectiveness can be significantly increased. This also comprises the advantages of the satellite-based method in terms of gathering evidence. Without on-site installations, depending on the spatial resolution of the radar image on the ground, the monitoring region can be expanded almost at free will and without great additional expense. A potential combination of both methods is shown in Fig. 2.44.

The fundamental principle is the reduction of terrestrial measurements in favor of radarinterferometric settlement monitoring. In the close range of the tunnel axis, terrestrial measurements are still indispensable at the time of the shield passage, since short measurement intervals are urgently needed for checking the alarm or threshold values. This is different in the long-distance range in space as well as in time, areal satellite-based methods are preferred which require some time for stochastic evaluation of several overflights. The best division of regions for terrestrial and radar-interferometric methods finally depends on the topology, the geology, the damage risk of existent buildings and potential compensation measures [98].

#### **2.6.4 Data Handling and Visualization Techniques**

Besides the actual settlement monitoring, visualization plays a decisive role [30]. Visualized in space and time, critical conditions such as torsional deformations when underpassing structures in parallel are recognized easier. Then, countermeasures can be taken more quickly and in a targeted manner. Numerous commercial data management systems offer closed solutions for visualization. This holistic and visual settings allow settlements to be interpreted significantly easier. Visualization supports different perspectives, provides fading of additional information and various zooming options for detailed viewing.

**Fig. 2.45** VR-Lab at the Ruhr-Universität Bochum, Germany

Time-related data such as settlement data can be animated which offers a high degree of interoperability between users. At best, a three-dimensional presentation of the results in a virtual reality (VR) environment supports the aforementioned analysis options (Fig. 2.45).

Virtual Reality offers the possibility of creating a realistic computer-generated threedimensional and dynamic image of the reality. The generation takes place in real time and can be converted to the user's viewing angle. To perceive the 3D effect, tools such as shutter glasses are required. A so-called flystick enables navigation. Both components can be tracked to account for the current position and viewing direction of the user in the computation of the image and perspective (Fig. 2.45).

If data is not only to be visualized but also made available for further analyses, e.g. for FE computations, integral product models on a BIM (Building Information Modeling) and IFC (Industry Foundation Classes) basis are needed [96]. Figure 2.46 exemplifies the

visualization of structures, tunneling, settlement data and the soil during construction of the Wehrhahn line in Dusseldorf, Germany [30].

#### **2.6.5 Case Study: Wehrhahn line in Dusseldorf**

The Earth-observation TERRASAR-X satellite of the German Aerospace Center (DLR) was launched in June 2007 to provide high-quality radar images from an altitude of 514 km. In the best available mode (Spotlight) the spatial resolution is about 11m[28]. However, in the region of interest (Fig. 2.47), the tunnel construction site in Dusseldorf, Germany [30], another mode (Stripmap) with a resolution up to 33 m was used.

**Results of displacements monitoring** To monitor displacements of existent infrastructures induced by tunneling a set of 16 TERRASAR-X images are evaluated, gaining already high density of persistent scatterers. A block of buildings near the starting shaft (Fig. 2.48 and 2.49) illustrates the settlement effects in a very clear way. It is especially appropriate, as the tunnel boring machine undergoes the buildings with only a small overburden between tunnel roof and existent foundations. Moreover, extended initial ground injections are arranged here to "prestress" the remaining ground and compensate subsequent settlements.

**Results of uplift monitoring** To monitor the process of intended uplift by injections and the following settlements a conventional rubber tube leveling system was established at the foundations of the buildings close to the tunnel. The terrestrial data is used to verify and calibrate the satellite data in space and time. For this purpose terrestrial and PS measuring points are geo-referenced in identical coordinate systems and superimposed on the 3D terrestrial laser scan model of the intra-urban area (Fig. 2.49) including positions on roofs, facades and the ground surface.

**Fig. 2.47** Radar-Image and tunneling axis (Dusseldorf, Germany) [99]

**Fig. 2.48** Tunnel lining and region of compensating injections [99]

**Fig. 2.49** Distribution of persistent scatterer along with displacements after eight months [99]

Figure 2.50 shows a representative time-series of displacements obtained from PS measurements compared to terrestrial values close to the PS points. All values are related to initial values prior to the tunneling workout to ensure an overall identical basis. The vertical ground movements derived from the PS-InSAR technique correspond very well to the terrestrial values. This especially holds true, if the very small extents of movements in

**Fig. 2.50** Comparison of time-series of terrestrial determined displacements (solid lines) and satellite determined displacements (dashed lines) during compensating injections and shield tunneling crossing date [99]

a millimeter range are taken into account as well as certain inevitable local differences in the single points. Figure 2.50 clearly illustrates the intended initial uplift of the injections and its compensating effect on the subsequent settlements. Thus, an almost unaffected state of the buildings is fulfilled. It should be noted that the displacements measured by the satellite are given in line of sight (LOS). They have to be converted with the angle of incidence to end up with purely vertical displacements that are comparable to the terrestrial measurements on ground. However, an exact evaluation of error rates is still in progress.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **3 Excavation Simulations and Cutting Tool Wear**

Lukas Brackmann , Arne Röttger, Hoang-Giang Bui , Sahir Butt , Golnaz Hoormazdi , Abdiel Ramon Leon Bal, Sebastian Priebe , Dennis Wingender, Hongwei Yang, Daniel Balzani , Klaus Hackl , Günther Meschke , Inka Mueller , and Jörg Renner

L. Brackmann (-) Materials Technology, Ruhr-Universität Bochum, Bochum, Germany e-mail<EMAIL_ADDRESS>
A. Röttger New Production Technologies and Materials, Bergische Universität Wuppertal, Wuppertal, Germany e-mail<EMAIL_ADDRESS>
H.-G. Bui - G. Hoormazdi - K. Hackl Mechanics of Materials, Ruhr-Universität Bochum, Bochum, Germany e-mail<EMAIL_ADDRESS>
G. Hoormazdi e-mail<EMAIL_ADDRESS>
K. Hackl e-mail<EMAIL_ADDRESS>
S. Butt - A. R. Leon Bal - S. Priebe - G. Meschke Structural Mechanics, Ruhr-Universität Bochum, Bochum, Germany e-mail<EMAIL_ADDRESS>
A. R. Leon Bal e-mail<EMAIL_ADDRESS>
S. Priebe e-mail<EMAIL_ADDRESS>
G. Meschke e-mail<EMAIL_ADDRESS>
© The Author(s) 2023 93 G. Meschke et al. (Eds.), *Interaction Modeling in Mechanized Tunneling*, https://doi.org/10.1007/978-3-031-24066-9\_3

#### **Abstract**

The mechanized tunnel construction is carried out by tunnel boring machines, in which the soil in front of the working face is removed, and the tunnel lining is carried out with shotcrete or the setting of segments and their back injection. Advancements in this field aim towards increase of the excavation efficiency and increase of the tool lifetime, especially in rock-dominated grounds. The latter is achieved by understanding the wear mechanisms abrasion and surface-fatigue, and by knowledge of the microstructureproperty relation of the utilized materials. Improvements for tool concepts are derived, based on experiments and simulations. A key parameter towards efficient rock excavation is the shape of the cutting edge of the utilized disc cutters. Sharp cutting edges have proven to generate higher rock excavation rates compared to blunt ones. The compressive strength of the utilized steel has to be high, to inhibit plastic deformation and thereby to maintain sharp cutting edges. This requirement competes with the demand for toughness, which is necessary to avoid crack-growth in the case of cyclic loading. Solutions for this contradiction lie in specially designed multiphase microstructures, containing both hard particles and ductile microstructural constituents. Besides adapting the alloying concept, these required microstructures and the associated properties can be adjusted by specific heat-treatments.

D. Wingender -D. Balzani

D. Balzani e-mail<EMAIL_ADDRESS>
H. Yang School of Civil Engineering, Sun Yat-sen University, Guangzhou, China e-mail<EMAIL_ADDRESS>
I. Mueller

Technical Mechanics, Bochum University of Applied Sciences, Bochum, Germany e-mail<EMAIL_ADDRESS>
J. Renner

Experimental Geophysics, Ruhr-Universität Bochum, Bochum, Germany e-mail<EMAIL_ADDRESS>
Continuum Mechanics, Ruhr-Universität Bochum, Bochum, Germany e-mail<EMAIL_ADDRESS>
#### **3.1 Introduction**

The excavation of soil and rock in mechanized tunneling is fundamentally different from that in other tunnel construction methods, such as open construction pits or blast excavation. In mechanized tunneling, tools are attached to a rotating shield, which also defines the tunnel's diameter. Thereby, the tools are selected taking into account the geological conditions. For example, chisels and scrapers are used in the case of non-cohesive soil, as they remove the geomaterial in a similar manner as a shovel. In contrast, disc cutters are employed in soft and hard rock, where material degradation occurs by crack initiation, crack propagation, and spalling of rock fragments. The profitability of tunnel construction is largely defined by the efficiency and speed of soil excavation. The excavation efficiency decreases due to the wear of the tools, so they have to be replaced upon reaching their wear limit. Precise knowledge of the tool wear and the associated excavation efficiency is essential for planning optimal tool change intervals and estimating construction costs. This chapter deals with excavation of soil and medium-strength rock, and the wear of the tools used for their excavation. It is important to understand both the soil-tool interaction to describe the degradation mechanisms and the tool wear on the microscale, thus deriving optimal material concepts for tunneling tools.

At the beginning of this chapter, the fundamental principles of soil and rock excavation are explained, and the results of recent advancements in the understanding of the rock excavation process are presented. The results of rock indentation tests and numerical simulations form the base for future design optimizations of tunneling tools, aiming towards improved excavation efficiency. Subsequently, the results of numerical simulations of the interaction between tunneling tools, such as disc cutters and chisel tools, and various soils and rocks are discussed with a focus on the determination of evolving cutting forces, especially at material interfaces, as is the case if tunneling occurs in heterogeneous ground conditions.

In the second part, the material concepts and wear mechanisms of tunneling tools are examined. Tool wear is discussed on macro- and micro-scale by explaining the interaction of the different tool microstructures and the abrasive particles in the ground to be excavated. In addition to abrasive wear, special emphasis is placed on the influence of cyclic mechanical loading and the connected wear mechanism of surface spalling. Several material concepts, such as hardfacing alloys and tool steels, are considered concerning their abrasive wear resistance and fatigue resistance. Furthermore, existing test methods for abrasiveness of soil and rock, as well as fatigue tests are described and critically evaluated. The numerical simulations investigate the previously introduced wear mechanisms on meso- and microscale. Mesoscale simulations address abrasive wear and focus on the interaction of abrasive particles with the tool surface during a scratching process. A bridge is built toward the scratch behavior of multiple particles, as well as particle mixtures from single-scratch simulations. The microscale simulations explore the material response of hard phase containing materials under static and cyclic loading conditions. The underlying computational methods are explained in detail, and the results of crack propagation simulations in different hard-phase containing materials are presented. Subsequently, novel techniques of in-situ monitoring of tool wear and damages using vibration analysis are presented. An experimentally based proof of concept for detecting damaged disc cutters is provided, after a short digress into the fundamental principles of wear detection and vibration analysis. A transfer of the gained results concludes the chapter with practical recommendations on increasing tunneling efficiency and improving tool lifetime.

#### **3.2 Excavation of Geomaterials in Mechanized Tunneling: Experiments and Simulations on Failure Mechanisms**

Mechanized tunneling relies on tunnel boring machines (TBMs) that have to be equipped with appropriate tools for the excavation of the material along the anticipated path. The terminology regarding the material "below our feet", for which we use the term "geomaterials" in an overarching sense, is unfortunately convoluted. Conventionally, the material below our feet is classified as either soil or rock. Appealingly simple, with only two groups, the classification and its use in different scientific communities bears complications. In geo-engineering, for example, the term "soil" is at times used synonymous with "ground," the latter meaning "solid material below us" in an objective way. Furthermore, either group is diverse in composition and mechanical behavior. Soils represent loose and unconsolidated sediments produced by the deposition of particles after their transport through air or water over vastly variable distances. A first subgroup of soil, the organic soils, with the adjective often dropped in agricultural context, is characterized by a substantial fraction of organic matter and forms the pedosphere. The introduction of the two further subgroups, granular and cohesive soils, rests on a mixture of structural and mechanical characteristics. For soils, cohesiveness, actually at conflict with the defining "loose," results from the presence of water that mediates and amplifies electrical forces between soil particles to the extent that they form aggregates with some shape stability, in cases associated with significant plasticity. While the mechanical property, cohesiveness, cannot be used alone to distinguish soils and rocks, the cause for the cohesion and its mechanical ramifications are strikingly different for them. The cohesiveness of rocks results from welding of grains, either by a "substance", the cement, acting as a glue leading to lithification of sediments or by the formation of "intact" grain boundaries (held together by interatomic forces as in metals or ceramics) during their genesis (crystallization from melts or recrystallization by solid-state reactions).

The overall success and the predictability of costs of excavation projects critically hinge on an (in the best of cases) a-priori knowledge or an experience-based principal understanding of the deformation behavior of the subsurface material to be penetrated. The selection of an appropriate excavation technique and associated tools requires at least a material classification, but the more is known about strength parameters the better. The excavation techniques and the activated failure processes shift from soil to rock. Cutting

**Fig. 3.1** Two basic processes in excavation: indentation (**a**) and cutting (**b**)

or scraping suffices to disintegrate soils, while indentation and fragmentation is required for rock (Fig. 3.1).

As a rule of thumb, the contents of organic matter decreases while the cohesiveness increases with depth. On an absolute scale, organic soils and cohesive soils exhibit compressive strengths as low as tens of kPa. The compressive strength of granular soils and highly porous sedimentary rocks overlap (1 to 10 MPa); crystalline rocks may exhibit strength up to a few hundreds of MPa. The significant changes in porosity during deformation and their effect on strength was first noted for soils; the related critical-state concept [38, 104] also has some applicability to porous rocks [22]. Likewise, the concept of effective stresses, introduced by Terzaghi for soils [93, 94], applies to rocks as well [71].

#### **3.2.1 Excavation of Soft Soils**

The mechanized excavation of soils by means of tunnel boring machines is characterized by the engagement of the scraper tools mounted on the cutting wheel with the ground. The TBM advancement actually leads to a strong and transient tool-soil interaction, during which the excavation tools penetrate the tunnel face and push the soil mass away, leading to destructuration and failure of the material that ultimately enters the excavation chamber. The excavation operation of soils involves the development of large displacements and deformations, as well as a strong coupling between skeleton deformation and pore pressure variations causing significant fluid flow in partially and fully saturated soils [63, 64]. The intrinsic complexity of mechanized excavation problems, associated to a great extent with the interaction between the excavation instrument and the penetrated soil, has stimulated experimental and numerical investigations [6–8, 52, 53] of tool-soil interaction processes. Here, we report on the development of an experimental device that allows for measuring the topology of the soil during the excavation and the evolution of the reaction forces on the tool when penetrating dry and partially saturated soils, and on the validation and application of single-phase [7, 8] and two-phase [53] numerical models based on the Particle Finite Element Method (PFEM) [68] coupled with a standard hypoplastic formulation for the modeling of granular materials [103].

**Fig. 3.2** Setup and main components of the developed excavation device (re-drawn after: [6])

#### **<IP_ADDRESS> RUB Excavation Device**

The setup of the excavation device consists of a rectangular container filled in with sand and enclosed by four Plexiglass walls. The relevant dimensions are indicated in Fig. 3.2. The excavation tool is a rectangular Plexiglass panel with a height of 25 cm, displaced along the excavation container by a stepper motor and a trolley mounted on a dual linear guide rail system with sliding bearings. The excavation device permits adjusting the penetration depth of the cutting blade. Four force transducers, attached to the cutting blade, register reaction forces and torque during the experiment.

The excavation tests were carried out in dense Haltener Silbersand [76], a siliceous sand with rounded grains. The following conditions were considered in the experiments: initial height of sand in the container H0 D 30 cm, penetration depth of the tool dp D 10 cm, and horizontal velocity of the tool vOx D 1.2 cm/s. After each test run, the sand was leveled to the same initial height (30 cm). Water-saturation is reached via a hose entering the container; water is added until a level is attained that slightly exceeds the soil surface.

#### **<IP_ADDRESS> Computational Analysis of Tool-Soil Interactions**

The computational plane-strain and 3D PFEM simulations of tool-soil interactions for dry and water-saturated conditions, involving single and multiple cutting tool, used the configurations displayed in Fig. 3.3. Excavation analyses of dry soils were performed utilizing dimensions extracted from [20]: L D 2:2 m, width w = 20 cm, initial height of granular material H0 = 30 cm, penetration depth of the tool dp = 20 cm, height of the tool Htool = 25 cm and its thickness ttool = 2.5 cm. The tool moves with a prescribed horizontal velocity of vOx = 1.0 cm/s. In simulations with multiple cutting tools, the horizontal velocity of the tools was vOx = 15 cm/s. For simulation of the excavation of water-saturated soil, the model dimensions were: L = 1.4 m, w = 23 cm, H0 = 30 cm, dp = 10 cm, and ttool = 2.5 cm.

**Fig. 3.3** Initial geometry, model setup and boundary conditions of the 2D (top, re-drawn after: [53]) and 3D (bottom, re-drawn after: [7]) PFEM excavation models

In this numerical model, the height of the tool is 5 cm larger than in the experiments (i.e. Htool = 30 cm), to avoid soil over-passing the tool at the end. A horizontal tool velocity of 1.2 cm/s was applied. A frictional tool-soil interface discretized with triangular contact elements (shown as a red layer in Fig. 3.3, top, was considered in plane-strain excavation analyses, whereas a no-slip tool-soil interface was assumed for 3D simulations.

Excavation simulations for dry soil were performed assuming corn kernels [20]. The geotechnical and hypoplastic parameters adopted for corn can be found in [7]. The initial void ratio was set to e0 D 0:82 (relative density of Id D 0:35). Simulations of watersaturated soil, were carried out for Silbersand, assuming an initial void ratio of e0 D 0:66 (Id D 0:70). The hypoplastic parameters adopted for this sand are contained in [6, 53]. The constitutive model was initialized using at-rest later pressure conditions.

**Fig. 3.4** Spatial distribution of the total velocities of the particles (in [m/s]) in the excavation according to the **a** 2D and **b** 3D excavation models, for a tool horizontal displacement of Sx D 20 cm

In both 2D and 3D numerical analyses assuming partially or fully saturated soil, the sides of the excavation models were impervious (Fig. 3.3), while the ground free surface was allowed to drain freely (i.e., pOw D 0). In plane-strain simulations of partially saturated sand, the saturated hydraulic conductivity Ksat w was estimated from the Kozeni-Carman permeability model [18]; the Soil-Water Characteristic Curve (SWRC) was defined via van Genuchten's model [97].

The spatial distributions of the total velocities of the particles in the deformed ground, according to the 2D and 3D excavation models (Fig. 3.4), show higher velocities near the cutting tool and within the heap of excavated material, than in the rest of the domain. The simulation results also show a shear slip plane in the velocity field, propagating from the bottom of the excavation tool towards the ground free surface.

**Fig. 3.5 a** Profiles of the ground free surface for a tool horizontal displacement of Sx D 40 cm. **b** Evolution of the tool reaction forces, according to experiments and DEM simulations [20] as well as PFEM analyses [7] (re-drawn after: [7])

The profile of the ground free surface computed with the PFEM is compared with experimental and DEM results reported in [20] (Fig. 3.5a). In general, a very good correlation between laboratory results and numerical predictions from the hypoplastic PFEM model is observed for the assessed tool displacement. The evolution of the horizontal reaction force generated on the cutting tool vs. its horizontal displacement, displays an initial increase with a steep slope, where predictions obtained from PFEM analyses exceed DEM and laboratory results (Fig. 3.5b). As the excavation process continues, the reaction forces steadily increase due to the accumulation of granular material in front of the tool. For the evaluated horizontal displacement range of the tool, PFEM predictions [7] agree well with the results obtained from the DEM and the experiments presented in [20].

Tool-soil interaction analyses considering three cutting tools simultaneously excavating in soil are analyzed. The setup consists of one leading tool (tool 1) located ahead of the other two trailing tools (tool 2 and 3), which are positioned at the same level (Fig. 3.3, top). Similar to the PFEM results concerning excavations with a single tool (Fig. 3.4), distinctive shear slip lines emerge from the lower part of each cutting tool (Fig. 3.6a).  Computed reaction forces on tools 1 (solid line) and 2 (dashed line) are very close, due to the similarity in the topologies of their associated excavation fronts, while for tool 3 (line with marks), higher reaction forces are calculated (Fig. 3.6b).

Excavation analyses in initially fully saturated Silbersand, are now considered. A staggered topology of the ground free surface occurs, i.e., bumps of soil develop ahead of the scrapper (Fig. 3.7a). In general, the computed and measured profiles of the free surface are in good agreement (Fig. 3.7b). For the assessed horizontal displacement of 35 cm, the maximum height attained by the heap of material in the test (denoted by HOmax) was nearly HOmax D 20:26 cm, at a horizontal distance of 11.03 cm from the tool. According to the PFEM simulations, the maximum height in the ground from the bottom of the tool (Hmax)

**Fig. 3.6 a** Spatial distribution of the total velocities of the particles (in [m/s]) for a horizontal displacement of the tools of Sx D 20 cm. **b** spatio-temporal evolution of the reaction forces in tools 1, 2 and 3 (re-drawn after: [7])

**Fig. 3.7 a** Measured and **b** computed profiles of the ground free surface for a tool horizontal displacement of Sx D 35 cm, in water-saturated excavations

is circa Hmax D 22:63 cm computed at 10.29 cm ahead of the tool. It can also be observed a denser nodal distribution in shear deformation zones of the soil (i.e. shear bands) as compared to the rest of the excavation domain. This is achieved by means of an adaptive re-meshing procedure based on a soil dilation criterion, incorporated into the hypoplastic PFEM formulation for the improved capture of strain localization zones in the ground.

The capabilities of the proposed two-phase PFEM formulation are further assessed by means of the the computed and measured profiles of the ground free surface at selected tool displacements (Fig 3.8a). Experimental data (solid gray lines) show maximum heights of the excavation heap of 18.6 cm and 22.75 cm, for tool displacements of 25 cm and 45 cm, respectively. The maximum heights of the ground computed with the PFEM model (solid and dashed black lines), for corresponding tool displacements, are 18.9 cm and 25.5 cm. The computed and measured topology of the ground for the assessed tool displacements, are in good agreement.

**Fig. 3.8** Computed and measured **a** excavation profiles for tool horizontal displacements of Sx D 25 and Sx D 45 cm and **b** reaction force-tool displacement curves obtained from excavation analyses in dry and water-saturated sand (re-drawn after: [53])

We evaluate the tool reaction force-displacement curves generated during excavations in Silbersand (Fig. 3.8b). To complement the experiments performed in water-saturated sand, dry excavation tests were also carried out. Test results in initially fully saturated sand (red solid lines), show an initial increase of the reaction forces, followed by strong oscillations in the reaction forces, where a maximum reaction force of around 1100 N, at around 27.5 cm of tool displacement, is registered by the force transducers. Excavation experiments performed in dry conditions (aqua solid line), on the contrary, show a maximum reaction force of 490 N at the final displacement of the tool, corresponding to 45 cm. In this cases, no strong oscillations in the force plot are detected. Although the proposed model (results in dark red and dark aqua solid lines) is not able to fully reproduce the large oscillations observed in the reaction forces obtained from excavation tests in initially saturated sand, the numerical results lie within the experimental range.

Finally, 2D and 3D excavation analyses in fully saturated Silbersand, are presented. For these simulations, no-slip conditions at the tool-soil interface and a constant saturated hydraulic conductivity of Ksat w D 1 104 m/s, are considered. Furthermore, the tool is horizontally displaced at constant velocity of 10 cm/s. The spatial distributions of pore water pressures in the deformed configuration of ground at a tool horizontal displacement of 36 cm, assuming initially loose i.e. e0 D 0:79 (Id D 0:32) and dense i.e. e0 D 0:70 (Id D 0:58) sands, are investigated (Figs. 3.9 and 3.10). For the evaluated tool displacement, pore pressures in the ground are mostly positive for initially loose sand, whereas for the denser specimen, negative pore pressures develop ahead of the excavation tool while the rest of the soil domain (e.g. near the left and bottom boundaries) undergoes positive pressures. Negative pore pressures are normally associated to the excavation [63] and strain localization phenomena [83] of dense, dilatant materials. Comparable distributions of pore pressures are computed with the 2D PFEM excavation model, for similar initial soil densities.

**Fig. 3.9 a** Three- and **b** two-dimensional spatial distributions of pore pressures (in [Pa]) in excavations performed on fully saturated loose sand (e0 D 0:79) at a tool displacement of Sx D 36 cm

The topology of the ground free surface for a horizontal displacement of the tool of 36 cm, computed with both the 2D and 3D PFEM models, is analyzed (Fig. 3.11a). In general, slightly higher elevations of the free surface are observed in simulation results pertaining to dense sand (black solid line, blue dots), as compared to results involving loose sand. For loose sand, the 2D excavation model (red line) predicts a higher elevation of the ground free surface in comparison to its 3D counterpart (yellow dots). In the case of dense sand, the predicted curves remain close for the most part. Predictions obtained from the 2D and 3D PFEM models agree well. The computed evolution of tool reaction forces with respect to the horizontal displacement traveled by the tool, is assessed (Fig. 3.11b). Larger reaction forces are obtained for the tool-soil interaction carried out in the denser sand. For this set of simulations, the difference in the force levels between dense and loose sand is nearly sixfold at the end of the excavation process. For the selected soil parameters and tool displacement range, predictions from both versions of the excavation model, are in close agreement.

**Fig. 3.10 a** Three- and **b** two-dimensional spatial distributions of pore pressures (in [Pa]) in excavations performed fully saturated dense sand (e0 D 0:70) at a tool displacement of Sx D 36 cm

**Fig. 3.11 a** Computed free surface profiles of the ground for a tool horizontal displacement of Sx D 36 cm. **b** Reaction force-displacement curves computed during the tool-soil interactions performed with the 2D and 3D PFEM models, assuming dense and loose sands

#### **3.2.2 Experimental and Simulation based Investigation of Rock Fragmentation**

After briefly discussing the deformation characteristics of rocks, we present a suite of laboratory indentation test on a variety of intermediate-strength rock types. Finally, a peridynamic simulation model [13, 17] is presented that was used to simulate the indentation processes.

#### **<IP_ADDRESS> Deformation Characteristics of Rocks**

Here, we briefly review the state of knowledge of the deformation characteristics of rocks focusing on aspects relevant for modeling purposes and intermediate-strength rocks. For modeling purposes, two questions are immanent: "How many parameter are needed to describe the deformation behavior? Where can the values of the relevant parameters be found?". We address the first one below and regarding the second we refer to data collections in [3, 86]. For further reading, we suggest the compact overview [59] and the extensive treatments [41, 71].

Rocks are aggregates of minerals, naturally occurring, in most cases crystalline compounds of elements. The chemical composition of the most common rock-forming minerals is actually restricted to a limited number of elements, i.e., Si, O, Ca, Mg, Fe, Al, Na, K, and H. The abrasiveness, the extent to which rock fragments scratch a (metallic) tool, is controlled by the hardness of the minerals, well known for the common rock-forming minerals and classified by Mohs (relative) scale. In contrast, fracturing of rocks, the elementary step of excavation, is controlled by the structure formed by the minerals. Planning as well as substantial modeling of an excavation project obviously requires a quantitative description of deformation behavior.

Elastic deformation of rocks is probably of subordinate relevance for their excavation. Elastic strains resulting from removal of material, the creation of openings, will remain well below 1% in most cases, because elastic moduli of rocks range from a few to a few tens of Gigapascal [34]. Elastic in-situ parameters are constrained by surveys using elastic waves. However, such dynamic parameters often significantly exceed the relevant static parameters [29]. This discrepancy originates from rock-mass heterogeneity associated with fractures and faults that also causes a scale dependence of elastic parameters. The arrival of elastic waves reflects the fastest paths associated with the least damaged rock sections, while static parameters represent bulk behavior dominated by the weakest sections. Non-linearity and inelasticity even for modest stress perturbations associated with excavation (typically orders of magnitude smaller than Young's moduli) are further consequences of the damage inventory.

As for any polycrystalline solid, the strength of rocks, i.e., a measure of the maximum stress they can bear in a specific loading configuration, depends on state variables, such as stress tensor components, temperature, and chemical milieu, and internal variables, such as grain size and porosity, including crack density. Under compression, the prevailing condition in the subsurface, geomaterials tend to fail on localized planes exhibiting obtuse angles to the maximum principal (compressive) stress. This morphology and orientation relation led to address the failure planes as shear faults.

The presence of cracks and pores leads to the prominent dependence of compressive strength on mean stress. Failure of low porosity rocks in compression is accompanied by dilation, a relative increase in their volume, because failure results from nucleation, growth, and interaction of microfractures; therefore their strength increases with "confinement." The higher the porosity is the more shear-enhanced compaction may occur, counterbalancing the dilation. Porous rocks may develop localized compaction bands normal to the largest principal stress or deform in a ductile manner by non-localized cataclastic flow [96] at sufficiently high confinement. Pore collapse can even be induced by isostatic loading owing to the stress concentrations on grain contacts.

The faster the loading the stronger the rock appears; however, the effect of loading rate on strength is irrelevant for most technical applications since a change of one order of magnitude in strength requires more than 10 orders of magnitude change in strain rate. An increase in temperature tends to weaken rocks. With the potential exception of some carbonate rocks, evaporates, and claystones, the reduction stays well below an order of magnitude as long as temperatures stay below about 300 °C, the current limit for engineering subsurface projects. The limited number of data on samples of significantly different size indicates a reduction on strength with increasing size that likely reflects a scaling between the size of pre-existing flaws and sample size. The presence of water tends to have a weakening effect; the identification of the underlying physical and chemical processes has proven difficult.

Before representing descriptions of the failure of rocks, it seems mandated to emphasize that the user of such relations (and associated empirical parameters) has to answer the question what is to be modeled. The immediate interaction of a cutting tool and a rock is probably dominated by the strength as determined on rock samples that are intact before the experiment, in the sense that pre-existing interfaces do not completely dissect them. In-situ, the pre-existing inventory (size, density) of discontinuities (joints, faults) affect failure progression and the fragment size. Studies on rock properties distinguish (intact) rock and rock mass (including mesoscopic structure, joints, faults).

Frictional strength poses a lower limit for the compressive strength of rock masses. Frictional strength of interfaces in rocks varies significantly with their roughness and the acting normal stress. Frictional sliding is associated with local wear of asperities and thus strongly depends on deformation history. When described by a conventional linear relation between shear stress and normal stress, Amonton's law, the intercept, often addressed as cohesion, tends to increase with roughness. The slope, the coefficient of friction, converges to values between about 0.6 and 0.8 for many rocks, an empirical observation today known as Byerlee's law or, probably more appropriately, rule [11].

In contrast to metals, rocks -not unlike concrete- exhibit a uniaxial compressive strength C that is about an order of magnitude larger than their tensile strength T [59]. The determination of tensile strength is technically cumbersome and thus such experiments are seldom performed. In contrast, measurements of the mode I fracture toughness KIc , the resistance of a material against propagation of a single tensile fracture, are simple and thus often performed. From a fracture mechanical perspective, tensile strength obeys

$$T = K\_{Ic} / \sqrt{\pi c\_{\rm ini}} \tag{3.1}$$

and thus constitutes a measure of the size of the pre-existing crack cini that leads to the macroscopic failure.

Fracture criteria of rocks, surfaces in the three-dimensional space of principal stresses, have been extensively investigated for more than half a century [71]. The commonly used linear Mohr-Coulomb criterion postulates that failure occurs on a plane, for which a critical shear stress

$$|\tau| = S + \mu\_{\text{int}} \sigma\_n \tag{3.2}$$

is reached, where S denotes an intrinsic shear strength and int the coefficient of internal friction. The shear strength is related to the uniaxial compressive strength as S D C = 2 q 1 C 2 int C int. Such a linear description is often acceptable for the limited range in normal stress relevant for technical applications. However, this formulation implies that the slope of the criterion int and the orientation of the failure plane with respect to the least principal stress, ˇ, are related by

$$
\beta = \frac{\pi}{4} \pm \frac{\arctan \mu\_{\text{int}}}{2}. \tag{3.3}
$$

Experimental evidence does not support this assertion but documents an increase in the failure angle and a decrease in the coefficient of internal friction with increasing mean stress [59]. The non-linearity is reflected by the empirical Hoek-Brown criterion, in its most general form requiring the determination of three parameters [27]. Quite some effort has been spent on relating the parameters to rock-mass indices to allow for modeling of failure on the meter to decameter scale. Costamagna et al. [21] related failure to the Cardano condition for the existence of three real-valued eigenvalues of the characteristic equation of the stress tensor to arrive at linked friction and fracture criteria involving three parameters.

Murrell [67] extended the micromechanical Griffith concept of failure [37] to overall compressive stress states to arrive at a non-linear criterion with a single parameter,

$$
\sigma\_{\rm oct}^2 = 8T\sigma\_{\rm oct} = \frac{2}{3}C\sigma\_{\rm oct},\tag{3.4}
$$

where oct D q .
1 2/ 2 C .
1 3/ 2 C .
2 3/ 2 /=3 and oct D .
1 C 2 C 3/=3 denote octahedral shear and normal stress, the latter being identical to mean stress.

The second equality in Eq. 3.4 reflects the relation between tensile and uniaxial compressive strength C D 12T , inherent in the criterion, in good agreement with experimental evidence. The fundamental criticism of fracture-mechanics based criteria as Eq. 3.4 addresses the notion that compressive failure is not caused by the propagation of a single critical defect. In contrast, compressive stresses eventually penalize growth [26] and the formation of the shear fault results from complex interaction and coalescence of multiple micro-cracks [58]. Subsequent fracture-mechanical treatments of compressive failure related the parameters of Mohr-Coulomb type linear criteria to the central parameters fracture toughness and initial crack length amended by a microscopic friction coefficient addressing the sliding of closed microfractures [5, 62, 89].

#### **<IP_ADDRESS> Laboratory Indentation Tests**

Rock fragmentation in mechanized tunneling involves in general two basic processes: indentation and cutting (Fig. 3.1). In both processes, the fragmentation depends on the penetration depth of the used tool, though they differ in whether the tool travels normal (indentation) or parallel (cutting) to the rock surface. When penetration increases up to a critical depth, the rock behavior transits from non-localized deformation to localized deformation, i.e., fracturing, the latter of which is favorable for rock fragmentation. The simple, well-reproducible and accurate indentation test has long been utilized to measure different material properties, such as hardness [43], yield stress [91], and fracture toughness [85], though mostly for glasses and ceramics, and to assess rock properties, such as drillability and cuttability [90, 92]. The maximum indentation pressure pmax is a pivotal parameter in interpretation of indentation tests; it represents the specific energy for removing a unit rock volume in cutting experiments with a Non-Truncated Tip Indenter (NTTI) [92].

The indentation process can be continuously monitored using advanced experimental techniques, such as electron scanning microscopy [54], digital image correlation [107], acoustic emission [19, 106], infrared thermography [57], and electronic speckle interferometry [19]. Numerical simulations, based, for example, on the discrete element method [40] and the finite element method [56], have also been performed. These experimental and numerical studies revealed that rock indentation involves several operating processes including volumetric compaction, plastic deformation, and macro fracturing.

Many simplified models have been proposed to explain the stresses and the deformation under indenters in brittle solids, perhaps most notably in glasses, among which is the cavity expansion model (CEM) [42, 51]. Ever since its transfer to rock, a frictional geomaterial [19, 39], the model has gained popularity, but formulations were limited to the Non-Truncated Tip Indenter (NTTI). Such ideally sharp indenters suffer from severe wear due to the high stress concentration at their tips. The tool wear lowers the fragmentation efficiency as energy is invested in the deformation of the tools instead of in rock breakage, which motivates usage of cutting tools with truncated tips. The use of Truncated Tip Indenters (TTI) changes the failure mechanism; a phase of compressive failure below the flat indenter precedes the tensile splitting. Often, indentation tests conducted using a TTI were erroneously interpreted using relations developed for NTTI. The incorporation of tip truncation by [1] does not account for the ultimate fracturing. Recently, Yang et al. [105] proposed a model based on the CEM for TTI that captures plastic deformation in compression and brittle tensile fracturing, and finite sample size. Here, we use the term "plastic" for macroscopically non-localized deformation irrespective of the deformation mechanisms on grain scale.

**Mechanical and acoustic emission responses during indentation** The investigations performed focused on intermediate-strength rocks, positioned between the two endmembers, soil and rock, either due to "weak" minerals (e.g., calcite, clays, halite) or "weak" structures (e.g., porous), because they pose particular problems regarding the "right" selection of excavation tools. We describe the general response of such intermediate-strength rocks to indentation using exemplary data from tests on a variety of sandstones, limestones, and tuffs with compressive strengths (UCS) from 11 to 140 MPa (Table 3.1). The damage progression and failure process was monitored during the indentation test using an acoustic-emission (AE) system (ASC Milne, Applied Seismology Consulting, UK). The uncertainty of locating AE hypocenters is around 8 ˙ 2 mm that is comparable to the size of the used AE sensors. Testing apparatus, procedures and specimen preparation are detailed in [105].

Typically, the force increases almost linearly to a distinct peak during the indentation process (Fig. 3.12); AE hypocenters gradually form a cluster beneath the indenter during the loading (Fig. 3.13). At the end of a test, the AE hypocenters trace the macroscopic fracture closely, indicating the validity of the AE technique for monitoring of the failure process. The indentation pressure typically exhibits a plateau-like maximum preceding the peak force, indicating that growth of the damage zone, where compressive failure conditions are reached, eventually proceeds at almost constant energy input until the tensile

**Fig. 3.12** Force, indentation pressure, number and total number of AE counts over an interval of 20 s as a function of time for a Gildehaus sandstone specimen of 84 mm diameter and 100 mm height. Tests performed in displacement control with a constant piston velocity of 0.05 mm/min

3 Excavation Simulations and Cutting Tool Wear 111 specimens

three-point bending tests on specimens of 10 mm diameter and 33 mm length; UCS: uniaxial compressive strength; pmax : peak indentation pressure determined from tests on specimens of 30 mm

diameter and 50 mm height. Subscripts "dry" and "sat" indicate testing of dry or saturated specimens, respectively

**Fig. 3.13** Cumulative temporal-spatial distribution of AE hypocenters (circles) at different percentages of peak force (**a**–**d**) indicated in Fig. 3.12 by the vertical dashed lines, end of test (**e**) and photograph after indentation test (**f**). The marker size indicates the relative magnitude of AE energy. Crosses indicate sensor locations. The color bar indicates time relative to the one at peak force (with regard to [105])

tangential stress at the rim of this zone suffice to induce growth of a macroscopic tensile fracture that ultimately splits the specimen. When peak indentation pressure is surpassed, the cluster region representing the damage zone below the indenter does barely grow further (Fig. 3.13c, d), until growth of the macroscopic fracture becomes apparent.

The force-displacement curves for the "hard," Anröchter Sandstone (AS in Table 3.1), with the highest uniaxial compressive strength (140 MPa) among the tested rocks, differ from that for the "weak" to "intermediate-strength" (uniaxial compressive stress less than about 80 MPa) rocks regarding the post-peak stage of indentation. For the intermediatestrength rocks, indentation beyond peak force is associated with a drop to a finite load that is then maintained for a displacement of several tenths of millimeters, during which the median crack propagates towards the bottom of the specimen. For the hard rock, brittle fracturing after peak force is more violent and typically results in coeval splitting of the sample and total loss of load bearing capacity.

**Analysis of the results of indentation tests** Analysis of the mechanical and acousticemission response indicates that the failure associated with rock indentation involves successively initial elastic deformation, punch-in of indenter into the specimen associated

**Fig. 3.14** Correlation between indentation characteristics and the two conventional strength parameters, unixial compressive strength (**a**) and fracture toughness (**b**). Uniaxial compressive strength (UCS) was only determined for dry samples. Indentation tests were not performed on saturated GBS samples. Error bars indicate standard deviations; a minimum of 3 tests were performed for each rock variety and saturation condition

with gradual formation of a crushed zone, formation of a damaged zone, and nucleation and growth of a macroscopic fracture that eventually splits the sample in half. We simplified sample failure due to indentation with a conceptual model as detailed in [105]. The model emphasizes the role of the peak-indentation pressure pmax that indicates the resistance of the tested rock to compressive failure in the damage zone, e.g., described by the Mohr-Coulomb failure law (Eq. 3.2), that precedes the initiation of macroscopic tensile fracturing.

The experimentally observed maximum indentation pressures correlate with fracture toughness for the eight investigated rocks (Table 3.1). While saturation of samples with water systematically reduced toughness, the sample-to-sample variability is too large to identify an effect of saturation on maximum indentation pressure (Fig. 3.14a). The correlation between pmax and uniaxial compressive strength (UCS) is also positive but the highly porous Gildehaus sandstone (GBS) seems to deviate from the general trend (Fig. 3.14b). The correlation of the maximum indentation pressure with the two strength parameters is consistent with our interpretation that it marks the transition from non-localized compressive failure in the damage zone to localized tensile fracturing at its boundary. The maximum indentation pressure increases with sample size, as evidenced by results for GBS specimens with a diameter from 30 mm to 84 mm.

For a specific rock variety, the sample-to-sample variability in maximum indentation pressure is not random, but we find an inverse correlation between pmax and the associated penetration depth dpmax (Fig. 3.15a). The pairs (dpmax, pmax) do however show no systematic relation for the suite of rocks.

**Fig. 3.15** Correlation between maximum indentation pressure and associated penetration depth in dimensional (**a**) and normalized form (**b**). Markers in blue, black and cyan indicate sandstone, tuff, and limestone, respectively. All but the encircled data points represent tests on samples with 30 mm diameter. The legend holds for (**a**) and (**b**)

To address material differences, the data is presented in normalized form, according to the theoretical analysis of [105]. Uniaxial compressive strength serves as a zero-order estimate of the stress corresponding to the yield conditions for the material in the damage zone, and we thus normalize maximum penetration pressure by it. The damage zone grows until its boundary reaches a pre-existing flaw, for which the tangential stress on the boundary corresponds to the critical stress concentration, the fracture toughness. For a given material, critical penetration depth and damage zone size correlate; damage zone size and the stresses in it–determined by compressive failure conditions–control the stress concentration on pre-existing flaws. The size of initial flaws scales approximately with .KIC=UCS/ 2 [5]; thus, we use this ratio to normalize critical penetration depth. As a result, there appears to be a general decreasing trend between the ratio of pmax/UCS and dpmax=.KIC=UCS/ 2 (Fig. 3.15b); yet, some differences between the various tested rock types with compositional differences remain. A potential reason could be the materialdependence of the initial flaw-size relation beyond that accounted for by the two strength parameters. Furthermore, Zhu et al. [108] found that the initial-flaw size scaling involves a dependence on porosity beyond its effect on fracture toughness and uniaxial compressive strength.

#### **<IP_ADDRESS> Numerical Investigation of Rock Fragmentation**

Simulation of rock excavation requires tools that can handle the discontinuities associated with fracturing leading to fragmentation. We developed a peridynamic model [87, 88] providing a natural way of incorporating discontinuities in the simulation domain. All peridynamic simulations were carried out using an extended version of the open source software *Peridigm* [55, 70]. The interested reader is referred to [14, 16] for further details on the simulation model.

We performed two sets of simulations with respect to the kinematics of the disc. Simulations of indentation tests (see Sect. <IP_ADDRESS> and [105]) were performed to provide insight into the formation of the crushed zone at the tool-rock interface and of tensile fractures emanating from that zone. Linear Cutting Machine (LCM) tests [79, 80] were treated to estimate the cutting forces required in the excavation process as a function of penetration depth. Besides homogeneous materials, we also considered mixed or heterogeneous ground conditions, situations where two or more types of rock or soil with significantly different mechanical properties are exposed at the tunnel face, to assess the variability of loads on the cutting discs and associated wear (see also Sect. 3.3.5).

These analyses are motivated by the characterization of the impulse load exerted on the cutting disc when it traverses a material interface, as this can cause excess vibration in the cutter head and increased fluctuations in torque and thrust of TBM. These peak loads also complicate the extraction of damage sensitive features based on the force oscillations.

**Simulation of indentation tests** Our simulations of indentation tests cover a total of six specimen sizes, with a combination of three diameters 30, 50 and 84 mm and two heights 50 and 100 mm. The truncated-tip indenter used is the same for all simulations (Fig. 3.16). For validating the current model, we use the experimental data from indentation tests on Gildehaus Sandstone reported in [105]. A qualitative validation relates to the formation of the crushed zone, from which tensile cracks initiate, while a quantitative validation rests on comparing the force-penetration data obtained from the experiments and the simulations.

**Fig. 3.17** Temporal evolution (left to right columns) of the fracture process during the indentation test (top row), the associated damaged and cracked regions are filtered out for visualization (bottom row)

The temporal evolution of fractures occurring due to the indentation load is presented in Fig. 3.17 for a specimen with 30 mm diameter and 100 mm height. The simulations capture the experimental observation of the successive formation of the crushed zone and the initiation of a central macroscopic tensile fracture splitting the sample in half. The interaction of the growing fracture with the sample's boundary leads to crack tip bifurcation [13, 15] and the crack branches. Finally, the specimen splits into two main and several small fragments. The loading stiffness as well as the peak load predicted by the simulation are in a good agreement with the experimental records (Fig. 3.18). However, the post-peak behavior differs slightly, which can be explained by the absence of damping in the simulations. As the main fracture starts propagating, the two large fragments move opposite to each other and separate faster than in the experiments, causing the indenter to

loose contact with the rock as soon as the specimen splits and thus the force at the indenter goes to zero. The progressing fracture opening can be can be observed qualitatively by comparing the last two columns of Fig. 3.17.

**Simulation of rock excavation with cutting discs** At the level of a single cutting disc, interactions between the rock and the cutting disc can be characterized by the reaction force, which is decomposed into normal, rolling and side forces. These forces need to be estimated for an individual disc to predict the performance of a TBM, i.e., the global torque and thrust requirements. In the LCM test [35, 78] developed at the Colorado School of Mines (CSM) to predict the performance of a single cutting tool, the cutting disc moves along the rock specimen at a known penetration level and the reaction forces are measured at the cutting disc. We use this test as a benchmark for our analysis model.

We performed peridynamic simulations of the LCM test with a cutting disc of constant cross-section and a diameter of 432 mm for four levels of fixed tool penetration. The dimensions of the specimen are 0.4 m 0.4 m 0.15 m in the simulations and we ascribed properties of an elastic-brittle constitutive relation corresponding to Colorado red granite (Young's modulus 41.0 GPa, Poisson's ratio 0.234, uniaxial compressive strength 158 MPa, fracture energy 85.7 J/m2). The simulation results obtained from the rock cutting with a tool penetration depth of 7 mm (Fig. 3.19a) show an evolving damaged zone beneath the cutting disc; the temporal evolution of rolling, normal and side forces (Fig. 3.19b) are examined from 0.1 to 0.3 m to avoid the influence of the specimen boundaries. The average normal forces obtained from the peridynamic simulations [14] agree well with the experimental data from LCM tests [35] (Fig. 3.19c).

**Simulations of excavation with cutting discs in mixed ground** The excavation process in mixed ground conditions is analyzed using a scaled-down cutting disc. The analyses

**Fig. 3.19** Simulation of the LCM test showing the evolving fragmentation and crushing zone in the rock (**a**), evolution of cutting forces for a tool penetration of 7 mm (**b**) and average normal cutting forces for different tool penetrations compared against the results from the LCM experiment [35] (**c**)

consider a cutting disc moving from soft soil to hard rock medium. The soil material is modelled using an elastic-plastic constitutive relation and the rock material is modelled using an elastic-brittle constitutive relation. The soil's yield strength was 100 KPa and the rock's fracture energy Gf = 23.7 J/m2.

The simulated reaction forces at the cutting disc vary significantly with lateral displacement for the two penetration levels (Fig. 3.20b, d). The normal cutting forces are negligible (for the chosen yield strength) when the disc cutter is in the soil domain, but forces increase as it approaches the soil-rock interface (cutting length of 0.05 m) and when it thrusts through the rock domain. The peak in the cutting forces increases with an increase in tool penetration. Plastic deformations in the soil as well as the damage level in the rock also increase with increasing penetration levels (Fig. 3.20a, c).

**Fig. 3.20** Peridynamic simulation of a cutting disc working in mixed ground conditions (soft soil to hard rock domain) with a penetration of p = 2 mm and p = 3 mm. Plastic deformations and damage/fragmentation in hard rock is shown in **a** and **c** and the associated cutting forces acting on the disc are shown in **b** and **d**, respectively. (The vertical grey line represents the soil-rock interface)

#### **3.2.3 Implications for Mechanized Tunneling and Outlook**

Tool-soil interactions are a central component of mechanized tunneling excavations. The complexity of these interactions underlies the rapid evolution of the ground free surface and the development of large deformations, accompanied by the flow of the interstitial fluid in water-saturated soils. Laboratory tests performed on dry and water-saturated sand using a custom-made excavation device provided first-hand information regarding the topology of the material and the reaction forces generated on the tool during the excavation process. The presence of interstitial fluid and the development of shear bands strongly influences the reaction force-displacement curve. Presumably, this influence is associated with localized pore-pressure variations in the soil during the excavation. Numerical simulations of tool-soil interactions demonstrated that the proposed hypoplastic PFEM formulations were able to reproduce the main features of the mechanized excavation processes, including the soil deformations and the reaction forces on the tool. Furthermore, the presented models provided valuable insight into the spatial distribution of pore pressures in the ground and the effects of the initial density of the soil on the excavation process. Ongoing and future research concerning soil excavation include numerical simulations of mechanized excavations using an EPB cutting wheel model. Relevant boundary conditions, such as the required support pressure at the tunnel face, rotational speed and penetration rate of the wheel, will be considered. We will investigate the evolution of the torque generated at the wheel, the flow of excavated material at the face and the spatial distribution of pore pressures in the soil.

A crushed zone develops in rocks when the cutting discs installed on a rotating cutter head are hydraulically pressed against the tunnel face; as the disc continues to penetrate further, stresses around the crushed zone increase to the point that eventually macroscopic fractures are initiated. Fractures associated with neighboring discs coalesce leading to fragmentation and disintegration of the rock. Experimental investigation and analytical modeling of indentation tests identified the governing parameters for the involved mechanisms, non-localized damaging and localized brittle fracturing. The peak indentation pressure appears to indicate the transition from non-localized damaging below a truncated indenter to localized brittle fracturing, pre-requisite for the rock fragmentation that is aimed for in excavation procedures. Therefore, in practice, the tool has to exceed the penetration depth associated with the peak indentation pressure for efficient rock excavation. The corresponding thrust force may then be estimated through the established correlation between normalized peak indentation pressure and penetration. Thus, rock indentation tests performed in the laboratory provide relevant information on determination of operational parameters, such as thrust force and optimum penetration. In-situ, the rock fragmentation in mechanized tunneling is obviously complicated by many other factors, such as lithostatic stresses and the water-weakening effect, whose effects are poorly understood, and are therefore subject of ongoing research.

Numerical analyses performed using the framework of peridynamics theory have been proposed to predict the performance of a TBM in various scenarios. Results obtained from the indentation tests provided insight into the formation of the crushed zone from where the localized tensile fractures initiate leading to rock disintegration. Furthermore, the cutting forces obtained for a full-scale linear cutting test agree with the experimental recordings. Future work should account for different disc geometries, tool penetrations, tool spacings and cutting speeds. The cutting forces obtained from simulations of cutting discs moving from a soft ground domain to a strong rock layer yield a peak load at the soil-rock interface. The model was expanded to account for localized abrasive wear caused to the disc due to the cutting forces. Mechanized tunneling in mixed grounds results in highly variable loads on the cutting discs that may cause abnormal cutter wear leading to unexpected TBM stoppage.

#### **3.3 Tool Wear in Mechanized Tunneling – Appearances, Mechanisms, and Countermeasures**

Wear of tunneling tools causes downtimes of the TBM, both in a planned manner but more critically at unforeseen times if the wear rate of the tools exceeds preliminary wear predictions. It is crucial to understand the effective wear mechanisms that arise in the various geological conditions, to predict the rate of tool wear and its influence on the efficiency of the tunnel boring machine. Based on this assessment, the material concept for the tunneling tools can be adapted to minimize tool wear and reduce the risk of failure. This chapter illustrates the macroscopic appearance of the most important wear mechanisms, explains the underlying tribological micro-mechanisms, gives examples on how the wear behavior of materials can measured, and points out possible improvements for material concepts.