Modelling Geotechnical Heterogeneities Using Geostatistical Simulation and Finite Differences Analysis

Abstract

Modelling a rock mass in an accurate and realistic way allows researchers to reduce the uncertainty associated with its characterisation and reproduce the intrinsic spatial variability and heterogeneities present in the rock mass. However, there is often a lack of a structured methodology to characterise heterogeneous rock masses using geotechnical information available from the prospection phase. This paper presents a characterization methodology based on the geostatistical simulation of geotechnical variables and the application of a scenario reduction technique aimed at selecting a reduced number of realisations able to statistically represent a large set of realisations obtained by the geostatistical approach. This type of information is useful for a further rock mass behaviour analysis. The methodology is applied to a gold deposit with the goal of understanding its main differences to traditional approaches based on a deterministic modelling of the rock mass. The obtained results show the suitability of the methodology to characterise heterogeneous rock masses, since there were considerable differences between the results of the proposed methodology, mainly concerning the theoretical tunnel displacements, and the ones obtained with a traditional approach.

1. Introduction

Existing tools and methodologies used to characterise rock masses have shown limitations in identifying the natural heterogeneities and in reducing uncertainties related to high spatial variability. Also, their impact on the geomechanical parameters and, consequently, on the overall behaviour of underground works is difficult to predict. In the past years, the use of probabilistic techniques, such as Taylor series and Monte Carlo simulation [1,2], have been able to mitigate these limitations; however, the impact of high levels of heterogeneity in the behaviour of rock masses is still a complex task and not duly solved. More recently, several authors have proposed the use of geostatistical techniques that, in contrast to traditional probabilistic techniques, take account of the spatial dependence of the data and allow for the accurate prediction of the unknown values of geomechanical parameters, such as Rock Quality Designation (RQD), Rock Mass Rating (RMR), and discontinuities parameters [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]. Furthermore, with geostatistical simulation techniques, several scenarios or realisations used to represent the variable(s) of interest are obtained, all of which reproduce the spatial correlation structure of the sample data [18].

This paper presents a methodology developed to model the geomechanical parameters of heterogeneous rock masses. The methodology starts with a probabilistic technique and a geostatistical simulation coupled with a clustering technique in order to obtain a reduced number of representative realisations, following works applied in the mining and petroleum industries [19,20,21,22,23]. With the application of the aforementioned methodology, it is possible to map, simultaneously, the uncertainty associated with the geomechanical parameters’ quantification, their spatial variability, and the heterogeneities existing in the rock masses at different spatial scales.

To validate the methodology, a numerical analysis of the rock mass mechanical behaviour during the excavation of a tunnel will be examined. For this purpose, data from a gold deposit located in the Atacama Region, northern Chile, will be used. The data include the Rock Mass Rating (RMR) empirical system and the Uniaxial Compressive Strength (UCS) of the intact rock, all extracted from mechanical boreholes and laboratory tests. Using this information, a theoretical tunnel will be modelled, using the proposed methodology to simulate the rock mass properties, and the results will be compared with a similar analysis performed using a deterministic approach that assumes the rock is a homogeneous mass.

The proposed methodology intends to take a step ahead by implementing geostatistical techniques, which are often assumed to be complex and impractical, mainly in relation to the geomechanical characterisation of rock masses.

2. Proposed Methodology

2.1. General

To apply this methodology it is necessary, in a first step, to collate all the available geotechnical information originating from different sources, such as boreholes and in situ and laboratory tests (Figure 1: step 1). Then, a statistical analysis of the available data should be performed in order to outline the behaviour of the variable(s) of interest through the calculation of the statistical parameters required for simulation.

A further preparatory step concerns the development of the grid model selected to represent the rock mass block and the underground structures. Once all of the shapes and elements of the model under study are defined, the corresponding two-dimensional (2D) or three-dimensional (3D) spatial information from all zone centroids is exported so that it can be used in the following rock mass characterisation methodology (Figure 1: step 2).

Subsequently, the chosen variables are simulated at the target grid points (zone centroids of the grid). Each variable represents either a geomechanical parameter or an empirical system value used to characterise the rock mass. If the variable is an empirical system value, the use of empirical formulas to obtain the values of the geomechanical parameters is an additional required step (Figure 1: step 3).

Based on the simulation results, an important and required step is the application of a dimension reduction methodology, whose goal is to avoid the computational burden of the following step, the finite differences analysis of the geostatistical realisations. The goal is to guarantee that the subset of selected realisations is statistically representative of the full set. As a final step, a numerical analysis of the underground structure behaviour using finite difference elements software [24] is performed based on the selected realisations. (Figure 1: step 4).

2.2. Geostatistical Simulation

Geostatistics provides a set of techniques that allow for the prediction of rock mass properties and, at the same time, the identification of different scales of heterogeneities and an assessment of the spatial uncertainty of the properties of interest. Therefore, in this paper, geostatistical simulation that results in multiple realisations of regionalised variables was applied. Prior to simulation, it is important to model the spatial correlation structure of the data through a covariance or a variogram. Some recent works have been using geostatistics to simulate the geotechnical parameters measured in continuous or discrete scales using the following general steps [16,17,18]:

• Primarily, the data are analysed and basic statistics are calculated;
• For each variable of interest, the data are transformed into data with a standard Gaussian distribution (normal scores);
• The normal scores’ experimental variograms are computed and the spatial behaviour, including the anisotropy, of the transformed parameters is determined. This stage is followed by model fitting;
• Using the zone centroids as a grid, Gaussian random fields, with covariance functions given by the models fitted in step 3, are simulated at the target points using the turning bands algorithm [25]. This algorithm is chosen because of its ability to accurately reproduce the distribution and spatial correlation structure of the Gaussian random field to be simulated and because of its low central processing unit (CPU) time and memory storage requirements, which outperform other simulation algorithms. In addition, a post-processing stage based on kriging is used to condition the simulation to the normal scores;
• The Gaussian values of each realisation are back-transformed to their original scale.

2.3. Scenario Reduction

As an output of the simulation, several realisations of the variables of interest are obtained; however, the time required to analyse each realisation is significant, especially in geotechnics where a subsequent numerical analysis is performed, which renders this option unviable from a practical point of view. An option would be to compute and use the realisations’ average values for each point of the target grid, but this would result in smoothing the values and excluding the extreme values (in particular, the lower values that are quite important in geotechnics). As an alternative, a subset of selected realisations can be used to characterise the variables of interest, although the selection should not be made randomly and without any type of criteria.

This has led to the necessity of developing a scenario reduction methodology to allow researchers to find a reasonable number of realisations for practical applications and, at the same time, to realistically characterise the spatial uncertainty in the values of the variables of interest. The methodology proposed here to reduce the number of realisations relies on clustering techniques, more precisely on a partitioning type of clustering, the k-medoid algorithm [19]. The main goal is to find an optimum number of realisations that are capable of statistically representing the whole set of realisations without needing a high computational time to enable a field application in underground works. The general steps outlined for the methodology are as follows:

• Calculate the Euclidean distance of any pair of realisations.
• Define a dissimilarity matrix between the n realisations using the Euclidean distances resulting from step 1.
• Apply Multi-Dimensional Scaling (MDS) to represent the dissimilarity matrix in a 2D space, and herewith a spatial representation of the n realisations is accomplished. The scaling chosen was classical multidimensional scaling [26].
• Apply a kernel transform to obtain a featured space. In detail, a scale factor is chosen so that the nonlinearity of the Euclidean distance values in the 2D space is reduced.
• Study and find the optimal number of clusters by performing for a first time the kernel k-medoid clustering technique for different numbers of clusters (from 1 until a predefined maximum). To choose the optimal number of clusters, an evaluation criterion is used, such as the average silhouette method.
• Perform the kernel k-medoid clustering for a previously defined optimum number of clusters.
• Plot (in 2D) the optimum number of clusters and their respective medoids (centres).
• Perform the post-processing analysis to understand the validity of each cluster.

As a note, it is important to underline that each point represented in the 2D space corresponds to one of the n geostatistical realisations. Also, in what concerns steps 6 and 7 above, the optimum number of clusters is not assumed as an input variable but is defined using optimisation methods [27], resulting in the computation and plotting of each cluster configuration while adopting as an evaluation criterion the average silhouette width [26]. After plotting the average silhouette width for all of the clusters, the optimum number corresponds to the configuration that yields the highest value (>0.5). Figure 2 displays a summarised scheme of all of the steps required to apply the proposed scenario reduction methodology.

3. Case Study

The methodology is applied to a theoretical tunnel modelled using geotechnical data gathered in a gold deposit, namely information regarding the parameters of the Rock Mass Rating (RMR) empirical system [28] and Uniaxial Compressive Strength (UCS) of the intact rock. The system is used to classify the quality of the rock mass into five classes (very poor, poor, fair, good, and very good). Five of the six parameters that compose the system are summed, with the exception of the sixth parameter that is subtracted as it takes into account the relationship between the orientation of discontinuities and the structure, resulting in a continuous scale that varies from 0 to 100. The mentioned RMR parameters are: (a) Uniaxial compressive strength of rock material (P1); (b) RQD (P2); (c) Discontinuity spacing (P3); (d) Condition of discontinuities (P4); (e) Groundwater conditions (P5); and (f) Orientation of discontinuities (P6).

3.1. Available Data

The data set used in this case study is composed of information recovered from exploration boreholes in a gold deposit located in the Cordillera de Los Andes, Atacama Region, Chile. The regional geology of the area is characterized by a group of intrusive, volcanic, and sedimentary rocks that are affected by fault zones that control the mineralization, which allow for the identification of four main lithological units of sedimentary rocks.

Approximately 4000 samples are available with information about the RMR’s five initial parameters. The samples are located on a regular grid of 40 m × 40 m × 20 m. At this stage, it is important to mention that the RMR under consideration for this study is the basic RMR that considers the sum of parameters P1 to P5. The uniaxial compressive strength (P1) was measured in the laboratory through uniaxial compressive tests using cylindrical rock samples (according to the standard ASTM D4543–08) and then normalized to a diameter of 50 mm (UCS_50mm). The remaining parameters (P2 to P5) were estimated directly from borehole logging. According to the results of the rock mechanics laboratory tests and the interpreted RMR values from the borehole samples, the rock mass is classified as fair to good (mostly in the range of 50 to 60) in terms of geotechnical quality.

3.2. Three-Dimensional Numerical Model

A 3D numerical model was developed using the Flac3d software [24] in order to investigate the differences in the tunnel behaviour when a heterogeneous approach is used to obtain the geomechanical parameters of the rock mass. The modelled tunnel presents a length of 90 m and is composed of a 6 m radius arch with vertical sidewalls of the same height. The model grid is 108 m wide along the $\mathrm{x}$ direction, 90 m along the $\mathrm{y}$ direction, 200 m along the $\mathrm{z}$ direction, and is composed, mostly, of brick elements, whose size increases as one moves further away from the tunnel, resulting in a finer mesh near the excavation. A total of 100,800 zones and 105,742 grid points compose the model. The tunnel’s central axis is located 96 m below the ground surface (Figure 3). In what concerns the model boundary conditions, the horizontal displacements in the model’s vertical boundaries were fixed, along with all of the displacements in the lower boundary. Note that this theoretical tunnel has the proposed methodology’s validation as its main purpose.

The adopted tunnel support system is composed of 0.20 m thick shotcrete simulated with shell elements with a linear elastic isotropic behaviour, with an elastic modulus of 20 GPa and a Poisson ratio of 0.25. The construction process begins with the excavation of the tunnel arch in a 3 m length followed by the application of shotcrete in the arch. Then, the remaining part of the tunnel is excavated and finally the shotcrete is applied in the walls of the tunnel. A gravitational initial stress field with a horizontal-to-vertical stress ratio (${\mathrm{K}}_0$) of 0.5 is adopted. A unit weight of 25 kN/m³ and a Poisson ratio of 0.20 were assumed for the rock mass.

3.3. Exploratory Data Analysis

The basic statistics of the RMR and UCS variables are presented in

Table 1. Rock Mass Rating (RMR) and UCS (in MPa) statistics for the samples’ values.

Statistical Parameter	RMR	UCS (in MPa)
Number of samples	3969	3969
Minimum	48	138
Maximum	78	208
Mean	66.7	176.9
Variance	14.2	1054.7

. The RMR shows a minimum of 48 and a maximum of 78 and its geomechanical quality varies from fair to good. Regarding the UCS, although it has a higher variance, within the RMR rating scale this parameter only presents a short range with only two different scales (rating 10 and rating 12). An additional analysis, prior to the geostatistical simulation, was made in order to analyse the cross-correlation between both variables. Hence, the Pearson product-moment correlation coefficient was computed, resulting in a weak correlation with a value of 0.09 (less than 0.3). Therefore, both variables will be simulated separately.

3.4. Model Spatial Continuity

Taking into account the available data of the case study, a geostatistical simulation is performed for the RMR and the UCS on the previously constructed grid mesh (i.e., considering the zone centroids). The variograms of both variables (previously transformed into normally distributed variables) are built along two main directions of anisotropy (regarding the spatial variation of the data), namely the horizontal ($\mathrm{xy}$) plane and the vertical ($\mathrm{z}$) direction, considering a lag distance multiple of 40 m for the horizontal plane and 20 m for the vertical direction (Figure 4).

It is possible to observe that the RMR shows an almost isotropic behaviour, since the variogram is similar for both analysed directions. However, this does not happen for the UCS, with the vertical direction showing a smaller correlation range than the horizontal plane, i.e., confirming the distinctive spatial behaviour in both directions (higher variability along the vertical direction). Nonetheless, neither of the variograms shows a nugget effect, proving the spatial continuity of the variables at short scale. The variograms are modelled with exponential and Gaussian basic structures as presented below (the distances between brackets represent the practical correlation range for the horizontal and vertical directions, respectively, and the number preceding the basic structure indicates the sill of the structure):

▪ RMR variogram model:

$$\gamma =0.61Exponential\left(250m;250m\right)+0.05Gaussian\left(450m;250m\right)$$

▪ UCS variogram model:

$$\gamma =0.44Exponential\left(450m;150m\right)+0.66Gaussian\left(500m;105m\right)$$

3.5. Conditional Simulation Results

Both variables were simulated using the turning bands algorithm (with a total of 1500 turning lines) and conditioned to the sample data. The number of realisations was set to 100, resulting in 100 RMR and UCS realisations.

The first realisation of both variables is presented in

Table 2. RMR and UCS (in MPa) statistics for realisation number 1.

Statistical Parameter	Realisation n. 1 of the RMR	Realisation n. 1 of the UCS
Number of grid points	100,800	100,800
Minimum	64	137
Maximum	77	207
Mean	71.1	159.5
Variance	2.5	391.5

and Figure 5 in order to demonstrate the existing spatial variability. The simulated values’ basic statistics are significantly different from those of the initial data, namely regarding the minimum and mean values that are higher for the simulation results. This discrepancy between the minimum and mean simulated values and, as a consequence, the variance value, can be justified with the size and location of the used grid since the grid area is significantly smaller than the data area.

To demonstrate the spatial variability existing in the rock mass, the RMR’s first realisation and the average of 100 realisations are mapped in Figure 5. Overall, the scale of variation for the RMR is low and the rock mass can be assumed to be nearly homogeneous.

3.6. From Geotechnical Data to Geomechanical Parameters

To perform the mechanical behaviour analysis of the theoretical tunnel, it was necessary to import into the finite differences software the geomechanical parameters of the rock mass instead of the geotechnical parameters obtained so far. Therefore, this crucial step consists in obtaining the geomechanical parameters, namely the rock mass deformation modulus, since a linear elastic behaviour is adopted to represent the rock mass. To perform this task, the application of empirical formulas that use as input the geotechnical variables was carried out. The number of formulas that could be used with the RMR systems and the UCS value is large, so that not all can be applied; however, in order to reduce the uncertainty associated with the use of a single formula, a set of four formulas are selected and a statistical methodology applied [29]. It is important to note that the used formulas have been developed using data from similar rock types to the one used in this case. Since the information regarding the Uniaxial Compressive Strength (UCS) exists, some of the selected formulas also consider this value (

Table 3. Empirical expressions used to obtain ${\mathrm{E}}_{\mathrm{m}}$ and the corresponding authors.

Number	Authors	Required Parameters	Limitations	Equation (${\mathit{E}}_{\mathit{m}}$ in GPa)
1	[30]	RMR	RMR > 50	${\mathrm{E}}_{\mathrm{m}}=2\times \mathrm{RMR}-100$
2	[31]	RMR	-	${\mathrm{E}}_{\mathrm{m}}=0.0003\times {\mathrm{RMR}}^3-0.0193\times {\mathrm{RMR}}^2+0.315\times \mathrm{RMR}+3.4064$
3	[32]	UCS, GSI	UCS > 100 MPa	${\mathrm{E}}_{\mathrm{m}}=\left(1-\frac{\mathrm{D}}{2}\right)\times {10}^{\left(\mathrm{GSI}-10\right)/40}$
4	[33]	RMR	-	${\mathrm{E}}_{\mathrm{m}}=0.1\times {\left(\frac{\mathrm{RMR}}{10}\right)}^3$

GSI: Geological Strength Index.

). Likewise, the same happens for the Geological Strength Index (GSI) that can be obtained from the RMR using a correlation formula [30]. This is a linear relation and the GSI is obtained by subtracting five units from the RMR value. Table 3 shows the selected ${\mathrm{E}}_{\mathrm{m}}$ expressions along with their limitations of application and their corresponding authors.

In the third formula, an additional parameter represented by the D letter is required. This parameter is dependent on the disturbance degree of the rock mass due to the excavation process and, in this case, an intermediate value was adopted (D = 0.5).

The admissible interval (AI) methodology [29] uses as an input the data obtained when different formulas are applied and combines them in order to get an interval, which can be calculated using the following equation:

$$\mathrm{AI}=\mathrm{mean}\pm \mathsf{\sigma }$$

(1)

where the mean term corresponds to the average values of ${\mathrm{E}}_{\mathrm{m}}$ using all of the chosen formulas and $\sigma$ represents the standard deviation of the average values. The range of the interval is given by the ${\mathrm{E}}_{\mathrm{m}}$ mean plus and minus one standard deviation.

Subsequently, a filter on the ${\mathrm{E}}_{\mathrm{m}}$ values for each one of the four formulas is applied, which is aimed at identifying the values located outside the range of the AI. These outside values are then eliminated and a new average of ${\mathrm{E}}_{\mathrm{m}}$ is computed using the remaining accepted values. As a consequence, a single ${\mathrm{E}}_{\mathrm{m}}$ value is calculated for each grid point. In the case of all of the values being rejected, the ${\mathrm{E}}_{\mathrm{m}}$ value is obtained with the average of the four formulas. Note that the methodology is applied to each point of the grid mesh and then to each one of the 100 realisations in order to obtain ${\mathrm{E}}_{\mathrm{m}}$.

The basic statistics of ${\mathrm{E}}_{\mathrm{m}}$ are presented in

Table 4. Statistical analysis of the ${\mathrm{E}}_{\mathrm{m}}$ (in GPa) values obtained through the application of empirical formulae using the geotechnical information.

Statistical Parameter	${\mathit{E}}_{\mathit{m}}$ (in GPa)
Total of points with information	100,800
Mean	39.50
Variance	7.30
Standard deviation	2.70
Minimum	19.74
Maximum	71.53

and the histogram of all of the 100 realisations is displayed in Figure 6. Most points show a ${\mathrm{E}}_{\mathrm{m}}$ value between 30 GPa and 45 GPa, corresponding to a rigid rock mass. Regarding the parameter dispersion, the range of values varies from 19 to 70 GPa; however, the number of points with these extreme values is low.

3.7. Scenario Reduction

As a first step, the 100 realisations of the rock mass deformation modulus, obtained using geostatistical simulation, are represented in a 2D space using the Euclidean distance as a metric (Figure 7). In order to apply the kernel algorithm, the Euclidean space should be transformed to a more linear space, the Feature Space. This intermediate stage can be achieved through the application of the kernel algorithm with a bandwidth of approximately 20% of the maximum Euclidean distance. To choose the optimal number of clusters, a silhouette average width test is performed [26]. This test allows for an assessment of the number of clusters that best represents the overall group of realisations and is performed from 2 until 10 clusters (i.e., the maximum number of clusters is predefined as 10). In this case, the optimal number of clusters is 4, which, from a geotechnical point of view, should be enough to cover the more conservative, the intermediate, and the optimistic scenarios. The final division into four clusters is obtained after applying kernel k-medoid clustering, which was run 500 times. Figure 7b displays the resulting configuration, where each cluster medoid is identified with a square and the corresponding elements with points in the same colour.

In order to confirm that the selected realisations statistically represent the full set, a validation is performed based on the comparison between the results for percentiles 10, 50, and 90 for all of the realisations. These percentiles are computed for each point of the grid point and then are averaged over the grid points. In this way, it is possible to obtain a single (average) value to represent each percentile (10, 50, and 90). Secondly, the percentiles 10, 50, and 90 computed for each grid point are summed point-by-point in order to understand the percentile trend for each realisation set. From Figure 8, it is possible to verify that, for a four-cluster configuration, the correspondence between the selected realisations and the full realisation set is slightly higher in the case of the selected realisations. Indeed, the difference in the ${\mathrm{E}}_{\mathrm{m}}$ average value is less than 1 GPa, which can be considered a low value taking into account the variation scale of the parameter (between 19 GPa and 70 GPa). Following the same trend, through the analysis of the percentile values for each point of the grid (Figure 8b), it is possible to notice that the clusters are able to achieve almost the same values as the full realisation set with quasi-coincident percentiles values.

4. Numerical Modelling Results

In order to understand the advantages of using the proposed methodology, the numerical models resulting from the individual realisations are compared with a deterministic approach, which considers the rock mass as a homogeneous medium (deterministic values for the geomechanical parameters). Therefore, for the deterministic model, here referred to as model 1, the rock mass deformation modulus is obtained firstly by averaging the RMR values in all of the 100,800 zone centroids and, subsequently, by applying the four aforementioned empirical formulas and averaging the ${\mathrm{E}}_{\mathrm{m}}$ results. Likewise, a model referred to as model 2 contemplates the four medoids (individual realisations) selected after applying the scenario reduction explained above. Finally, and since the computational time of each model is modest due to the simplicity of the theoretical tunnel, the results from the previous models are compared with the ones obtained for each one of the 100 realisations. It is important to point out that this analysis will also serve as a validation for the scenario reduction methodology and, subsequently, for the proposed approach.

For the sake of clarification for all of the analysed models, Figure 9 presents a scheme illustrating the path followed on the mentioned process, starting with the variable selection and its geostatistical simulation, then the application of the scenario reduction methodology, and finally the ${\mathrm{E}}_{\mathrm{m}}$ output values. Regarding model 1, the deterministic value adopted to represent ${\mathrm{E}}_{\mathrm{m}}$ is 32.73 GPa.

In this section, the numerical results obtained from the finite difference calculations are presented. Emphasis will be given to the results in terms of displacements, since the deformation modulus affects them in a direct way; however, some reference will be made to the maximum and minimum principal stress values calculated in the rock mass. This analysis is performed for all of the mentioned models. In detail, regarding the displacements on the rock mass, the maximum values obtained on the top of the tunnel arch, at the mid-point of the invert and at mid-height on the left and right sidewalls, are registered.

Primarily,

Table 5. Displacements and stresses values obtained for the 100 individual realisations.

Statistical Parameter	Vertical Displacement (mm)		Horizontal Displacement (mm)		Max. Stress (MPa)	Min. Stress (MPa)
	Arch	Invert	Left Wall	Right Wall
Minimum	5.20	6.82	1.50	1.42	6.67	1.77
Maximum	8.34	11.31	2.75	2.63	7.91	1.84
Mean	6.60	8.62	2.00	2.00	6.97	1.80
P10	5.92	7.76	1.73	1.71	6.81	1.77
P50	6.60	8.62	2.00	2.00	6.97	1.79
P90	7.49	9.80	2.28	2.33	7.21	1.82

presents a summary of the maximum, minimum, and percentile values obtained from the calculations for the 100 individual realisations. Since model 2 concerns individual realisations, they are also included in the 100 realisations values. Note that these values should serve as a reference for further comparisons.

Comparing the 100 realisations results, some differences are worthy of reference, namely the differences between the maximum and minimum values for the vertical and horizontal displacements, where a significant range of values is obtained. In terms of principal stresses, these differences are subtler with variations smaller than 1 MPa. Regarding the displacements’ mean values, the differences between the left and right sidewalls are zero, confirming that the heterogeneities present in this rock mass are small.

Looking in detail at models 1 and 2,

Table 6. Summary of the displacements and stresses results obtained in models 1 and 2.

Model	Corresponding Realisation	Vertical Displacement (mm)		Horizontal Displacement (mm)		Max. Stress (MPa)
		Arch	Invert	Left Wall	Right Wall	Arch	Invert
Model 1	-	5.86	9.36	1.79	1.79	6.91	1.77
Model 2-Cluster 1	Realisation 28	6.31	9.25	2.00	2.14	6.85	1.80
Model 2-Cluster 2	Realisation 46	6.07	8.24	1.82	1.87	7.05	1.80
Model 2-Cluster 3	Realisation 53	6.95	10.1	2.11	2.22	6.99	1.80
Model 2-Cluster 4	Realisation 74	6.02	7.03	1.56	1.84	7.91	1.80

shows a summary of the displacements and stresses obtained from the numerical analysis. Likewise, in order to understand the full data set representation using the clusters, an individual analysis is completed for each displacement and stress value. The aim is to represent the maximum and minimum values of all of the analysed models. Therefore, in Figure 10 and Figure 11, all of the model values for the displacements and principal stresses obtained after the tunnel excavation are represented.

In what concerns the maximum displacements results, in a more general analysis, the following comments can be made: (1) there are significant differences between the displacements’ maximum values from cluster to cluster, mainly in the tunnel invert; (2) differences are registered in the horizontal displacements of the tunnel sidewalls; and (3) model 1 (homogeneous approach) shows, in general, lower displacement values in comparison with each individual cluster.

In a detailed analysis, the maximum vertical displacements of model 1 (homogeneous approach) are observed in the centre of the tunnel invert with a magnitude of 9.36 mm, while the central point of the tunnel arch shows a smaller value of 5.86 mm (Figure 10).

Concerning model 2, cluster 3 is the one that registers the highest displacement values in all of the analysed points, while cluster 4 registers the lowest ones. The maximum vertical displacement occurs in the centre of the tunnel invert with, approximately, 10.1 mm (cluster 3), whereas the maximum displacement on the tunnel arch admits a magnitude of 6.95 mm (cluster 3). For the sidewalls, the values vary between 1.56 mm (left) and 1.84 mm (right) for cluster number 4 and 2.11 mm (left) and 2.22 mm (right) for cluster number 3. The observed disparity between the horizontal displacements can be explained by the assumed variability in the ${\mathrm{E}}_{\mathrm{m}}$ values (Figure 11). Looking at Figure 10 and Figure 11, it is possible to conclude that there are higher values and more dispersion in the displacements obtained in the tunnel invert in comparison with the tunnel arch. In addition, model 2 tends to slightly underestimate the displacement values in comparison with the 100 realisations results (Figure 10 and Figure 11), namely in the tunnel sidewalls where the dispersion of values is higher. Moreover, the homogeneous approach tends to underestimate the tunnel arch displacements and slightly overestimate the displacements in the tunnel invert. Regarding the tunnel sidewalls, the resulting displacement values show a lower dispersion and the homogeneous approach shows a tendency to be not on the safe side, that is, with underestimated values.

Note that all of the previously presented values were obtained after the tunnel’s total excavation and in the same section of the tunnel, with a $\mathrm{y}$ value equal to 0. Figure 12 shows the displacements contour in model 1 (Figure 12a) and all of the clusters of model 2 (Figure 12b–e) for comparison. The overall displacements are small due to the good geomechanical quality of the rock mass, resulting in high ${\mathrm{E}}_{\mathrm{m}}$ values.

In terms of stress distribution, and based on the maximum and minimum principal stresses, it is possible to conclude that the rock mass surrounding the tunnel is predominantly in compression for both models. In model 2-cluster 4, the tunnel walls are the most compressed zones, with values of 7.9 MPa and 1.8 MPa for the maximum and minimum principal stresses, respectively. Moreover, the tunnel arch and invert show compression stresses ranging between 4.0 MPa and 0.4 MPa. Similarly, the stresses of model 1 present the same range of values as model 2. Indeed, this similarity in the stress values can be justified by the lower influence of the ${\mathrm{E}}_{\mathrm{m}}$ variation on the rock mass principal stresses. As expected, the stresses are smaller on the tunnel invert and arch where the maximum displacements are observed. Regarding the tensile stresses, they only appear in the centre of the walls and invert with values under 0.2 MPa.

5. Discussion

In this section, an analysis will be made comparing the obtained values for the displacements and principal stresses of this tunnel case, with the aim to address three main points: (1) the heterogeneity representation of the individual realisations when compared with a traditional approach that assumes the rock mass is a homogeneous medium; (2) comparing the results, in terms of values and computational time, of using the scenario reduction optimal cluster configuration with respect to all of the 100 realisations, i.e., giving an answer to the question if the optimal clusters configuration is able to integrate the most conservative and optimistic values and embody a good solution to integrate into the numerical analysis; and (3) the computational time needed to analyse the 100 geostatistical realisations and the ones obtained from the proposed methodology.

Regarding the first point, it is possible to observe from Figure 12 that the heterogeneity of the rock mass is manifested by the asymmetrical pattern of the displacements contour in model 2. Indeed, the disparity of the values of the horizontal displacements confirms this statement. In relation to model 1, model 2 shows percentage differences in displacements values ranging from 1% to 25%, being either positive or negative. In terms of principal stresses, the maximum difference between the heterogeneous (model 2) and the homogeneous approach was 14% for the maximum principal stresses and 2% in the case of the minimum principal stresses.

It seems important to mention that, even if they are small, the differences between a heterogeneous and a homogeneous model can be relevant in the geotechnical analysis of underground works; in fact, these differences should be higher as the rock mass spatial variability increases and as the rock mass shows worse characteristics.

To perform the analysis regarding point number 2, histograms containing all the models’ results have been computed (Figure 13 and Figure 14). Similar to the graphical representation of the displacements and principal stress values, using the histograms it is possible to confirm the same tendency observed for model 2; that is, the ability of the model, introducing the concept of clustering, to show more optimistic and more conservative results. Also, it is possible to state that many of the 100 realisations show values in-between the referred clusters, confirming the capability of the scenario reduction methodology to represent, statistically, the full set of realisations. To complete this analysis, the first three moments (mean, standard deviation, and skewness) obtained from the 100 realisations distribution fitting can be consulted in

Table 7. Distribution curve fitting details (first three moments) regarding the 100 realisations’ obtained values for displacements.

Displacement	Mean	Standard Deviation	Skewness
Maximum horizontal displacement (mm): left sidewall	2.01	0.67	0.091
Maximum horizontal displacement (mm): right sidewall	2.00	0.24	0.100
Maximum vertical displacement (mm): centre of the arch	6.63	0.67	0.006
Maximum vertical displacement (mm): centre of the invert	8.71	0.83	0.003

. For all of the analysed zones, the displacements show a symmetrical distribution with some tendency to the right side (positive skewness); that is, higher displacements. About the third point concerning computational time, it was approximately 2 h for each individual realisation. Hence, to analyse the 100 realisations a total of 200 h is needed. If a reduction methodology is adopted, the time can go down to a maximum of 21 h (2 h per realisation for a maximum of 10 realisations, plus 1 h to apply the scenario reduction methodology). However, we strongly believe that this time difference can even be bigger in the case of more complex numerical models that require more computational time in the finite difference software.

6. Conclusions

In this work, a methodology for the numerical modelling of heterogeneous rock masses that considers uncertainty quantification, variability of the geomechanical parameters, and the presence of heterogeneities was presented. The methodology is outlined to combine geostatistical simulation and finite difference analysis. Therefore, using geotechnical information commonly obtained from field surveys, it is possible to simulate the variables of interest while honouring the values observed or measured at the data locations. A scenario reduction methodology allows researchers to obtain a smaller number of realisations without losing relevant information so that a numerical analysis of the underground structure can be made using a smaller number of possible scenarios. A prior step requires the conversion of the geotechnical information into the deformation modulus (E_m) through the use of empirical formulas. If detailed in situ data based on large-scale tests existed, this step could be avoided.

In order to validate the proposed methodology, a real data set containing geotechnical data about a gold deposit was used. The geotechnical set contains information of the empirical system RMR and values of the UCS measured by laboratory tests. The main goal of this work was not only to validate the proposed methodology, but also to compare the main differences in the mechanical behaviour of the structure when modelled using the heterogeneous approach and a deterministic, i.e., an isotropic and homogeneous, model. The heterogeneous approach proved to be encouraging for heterogeneous rock mass modelling due to the lower computational time required and the capability to statistically represent a full set of realisations (resulting in different scenarios, from conservative to optimistic ones).

Therefore, the main findings of this work are: (1) the developed methodology allows for the simulation of a rock mass’s geomechanical properties and the export of that information into numerical analysis software to perform a rock mass mechanical behaviour analysis; (2) the use of geostatistical information allows for the identification of a rock mass’s spatial variability and heterogeneities, which results in more accurate and realistic rock mass models; (3) although the rock mass used here had lower variability, some interesting differences were detected between the proposed approach and the homogeneous one, with the homogeneous one having a tendency to underestimate the displacement value (which is not on the safe side); and (4) the scenario reduction methodology allows researchers to obtain a reduced and practical number of realisations that produces results that are representative of the full realisations set.

In future developments, the proposed methodology should be validated with a different type of rock where the geomechanical variability is higher.